QR Decompositions¶
QR decompositions are a powerful tool in linear algebra and data science for several reasons. They provide a way to decompose a matrix into an orthogonal matrix $Q$ aand an upper triangular matrix $R$, which can simplify many computations and analyses.
Theorem: Let $A$ is an $m \times n$ matrix with linearly independent columns ($m \geq n$ in this case), then $A$ can be decomposed as $A = QR$ where $Q$ is an $m \times n$ matrix whose columns form an orthonormal basis for Col($A$) and $R$ is an $n \times n$ upper-triangular invertible matrix with positive entries on the diagonal.
In the literature, sometimes the QR decomposition is phrased as follows: any $m \times n$ matrix $A$ can also be written as $A = QR$ where $Q$ is an $m \times m$ orthogonal matrix ($Q^T = Q^{-1}$), and $R$ is an $m \times n$ upper-triangular matrix. One follows from the other by playing around with some matrix equations. Indeed, suppose that $A = Q_1R_1$ is a decomposition as above (that is, $Q_1$ is $m \times n$ and $R_1$ is $n \times n$). Use can use the Gram-Schmidt procedure to extend the columns of $Q_1$ to an orthonormal basis for all of $\mathbb{R}^m$, and put the remaining vectors in a $(m - n) \times n$ matrix $Q_2$. Then
$$ A = Q_1R_1 = \begin{bmatrix} Q_1 & Q_2 \end{bmatrix}\begin{bmatrix} R_1 \\ 0 \end{bmatrix}. $$
The left matrix is an $m \times m$ orthogonal matrix and the right matrix is $m \times n$ upper triangular. Moreover, the decomposition provides orthonormal bases for both the column space of $A$ and the perp of the column space of $A$; $Q_1$ will consist of an orthonormal basis for the column space of $A$ and $Q_2$ will consist of an orthonormal basis for the perp of the column space of $A$.
However, we will often want to use the decomposition when $Q$ is $m \times n$, $R$ is $n \times n$, and the columns of $Q$ form an orthonormal basis for the column space of $A$. For example, the python function numpy.linalg.qr give QR decompositions this way (again, assuming that the columns of $A$ are linearly independent, so $m \geq n$).
Key take-away. The QR decomposition provides an orthonormal basis for the column space of $A$. If $A$ has rank $k$, then the first $k$ columns of $Q$ will form a basis for the column space of $A$.
For small matrices, one can find $Q$ and $R$ by hand, assuming that $A = [ a_1\ \cdots\ a_n ]$ has full column rank. Let $e_1,\dots,e_n$ be the unnormalized vectors we get when we apply Gram-Schmidt to $c_1,\dots,c_n$, and let $u_1,\dots,u_n$ be their normalizations. Let $$ r_j = \begin{bmatrix} \langle e_1,c_j \rangle \\ \vdots \\ \langle e_n, c_j \rangle \end{bmatrix}, $$ and note that $\langle e_i,c_j \rangle = 0$ whenever $i > j$. Thus $$ Q = \begin{bmatrix} u_1 & \cdots & u_n \end{bmatrix} \text{ and } R = \begin{bmatrix} r_1 & \cdots & r_n \end{bmatrix} $$ give rise to a $A = QR$, where the columns of $Q$ form an orthonormal basis for $\text{Col}(A)$ and $R$ is upper-triangular. We can also compute $R$ directly from $Q$ and $Q$. Indeed, note that $Q^TQ = I$, so $$ Q^TA = Q^T(QR) = IR = R. $$
Example. Find a QR decomposition for the matrix $$ A = \begin{bmatrix} 1 & 1 & 1 \\ 0 & 1 & 1 \\ 0 & 0 & 1 \\ 0 & 0 & 0 \end{bmatrix}. $$ Note that one trivially see (or by applying the Gram-Schmidt procedure) that $$ \begin{bmatrix} 1 \\ 0 \\ 0 \\ 0 \end{bmatrix}, \begin{bmatrix} 0 \\ 1 \\ 0 \\ 0 \end{bmatrix}, \begin{bmatrix} 0 \\ 0 \\ 1 \\ 0 \end{bmatrix} $$ forms an orthonormal basis for the column space of $A$. So with $$ Q = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ 0 & 0 & 0 \end{bmatrix} \text{ and }R = \begin{bmatrix} 1 & 1 & 1\\ 0 & 1 & 1 \\ 0 & 0 & 1 \end{bmatrix}, $$ we have $A = QR$.
Let's do a more involved example.
Example. Consider the matrix $$ A = \begin{bmatrix} 1 & 0 & 0 \\ 1 & 1 & 0 \\ 1 & 1 & 1 \\ 1 & 1 & 1 \end{bmatrix}. $$ One can apply the Gram-Schmidt procedure to the columns of $A$ to find that $$ \begin{bmatrix} 1 \\ 1 \\ 1 \\ 1 \end{bmatrix}, \begin{bmatrix} -3 \\ 1 \\ 1 \\ 1 \end{bmatrix}, \begin{bmatrix} 0 \\ -\frac{2}{3} \\ \frac{1}{3} \\ \frac{1}{3}\end{bmatrix} $$ forms an orthogonal basis for the column space of $A$. Normalizing, we get that $$ Q = \begin{bmatrix} \frac{1}{2} & -\frac{3}{\sqrt{12}} & 0 \\ \frac{1}{2} & \frac{1}{\sqrt{12}} & -\frac{2}{\sqrt{6}} \\ \frac{1}{2} & \frac{1}{\sqrt{12}} & \frac{1}{\sqrt{6}} \\ \frac{1}{2} & \frac{1}{\sqrt{12}} & \frac{1}{\sqrt{6}} \end{bmatrix} $$ is an appropriate $Q$. Thus $$ \begin{split} R = Q^TA &= \begin{bmatrix} \frac{1}{2} & \frac{1}{2} & \frac{1}{2} & \frac{1}{2} \\ -\frac{3}{\sqrt{12}} & \frac{1}{\sqrt{12}} & \frac{1}{\sqrt{12}} & \frac{1}{\sqrt{12}} \\ 0 & -\frac{2}{\sqrt{6}} & \frac{1}{\sqrt{6}} & \frac{1}{\sqrt{6}} \end{bmatrix}\begin{bmatrix} 1 & 0 & 0 \\ 1 & 1 & 0 \\ 1 & 1 & 1 \\ 1 & 1 & 1 \end{bmatrix} \\ &= \begin{bmatrix} 2 & \frac{3}{2} & 1 \\ 0 & \frac{3}{\sqrt{12}} & \frac{2}{\sqrt{12}} \\ 0 & 0 & \frac{2}{\sqrt{6}} \end{bmatrix}. \end{split} $$ So all together, $$A = \begin{bmatrix} \frac{1}{2} & -\frac{3}{\sqrt{12}} & 0 \\ \frac{1}{2} & \frac{1}{\sqrt{12}} & -\frac{2}{\sqrt{6}} \\ \frac{1}{2} & \frac{1}{\sqrt{12}} & \frac{1}{\sqrt{6}} \\ \frac{1}{2} & \frac{1}{\sqrt{12}} & \frac{1}{\sqrt{6}} \end{bmatrix}\begin{bmatrix} 2 & \frac{3}{2} & 1 \\ 0 & \frac{3}{\sqrt{12}} & \frac{2}{\sqrt{12}} \\ 0 & 0 & \frac{2}{\sqrt{6}} \end{bmatrix}. $$
To do this numerically, we can use numpy.linalg.qr.
import numpy as np
# Define our matrices
A = np.array([[1,1,1],[0,1,1],[0,0,1],[0,0,0]])
B = np.array([[1,0,0],[1,1,0],[1,1,1],[1,1,1]])
# Take QR decompositions
QA, RA = np.linalg.qr(A)
QB, RB = np.linalg.qr(B)
Our resulting matrices are:
print(f"QA = {QA}\n")
print(f"RA = {RA}\n")
print(f"QB = {QB}\n")
print(f"RB = {RB}")
QA = [[ 1. 0. 0.] [-0. 1. 0.] [-0. -0. 1.] [-0. -0. -0.]] RA = [[1. 1. 1.] [0. 1. 1.] [0. 0. 1.]] QB = [[-0.5 0.8660254 0. ] [-0.5 -0.28867513 0.81649658] [-0.5 -0.28867513 -0.40824829] [-0.5 -0.28867513 -0.40824829]] RB = [[-2. -1.5 -1. ] [ 0. -0.8660254 -0.57735027] [ 0. 0. -0.81649658]]
How to use QR decompositions¶
One of the primary uses of QR decompositions is to solve least squares problems, as introduced above. Assuming that $A$ has full column rank, we can write $A = QR$ as a QR decomposition, and then we can find a least-squares solution to $Ax = b$ by solving the upper-triangular system.
Theorem. Let $A$ be an $m \times n$ matrix with full column rank, and let $A = QR$ be a QR factorization of $A$. Then, for each $b \in \mathbb{R}^m$, the equation $Ax = b$ has a unique least-squares solution, arising from the system $$ Rx = Q^Tb. $$
Normal equations can be ill-conditioned, i.e., small errors in calculating $A^TA$ give large errors when trying to solve the least-squares problem. When $A$ has full column rank, a QR factorization will allow one to compute a solution to the least-squares problem more reliably.
Example. Let $$ A = \begin{bmatrix} 1 & 0 & 0 \\ 1 & 1 & 0 \\ 1 & 1 & 1 \\ 1 & 1 & 1 \end{bmatrix} \text{ and } b = \begin{bmatrix} 1 \\ 1 \\ 1 \\ 0 \end{bmatrix}. $$ We can find the least-squares solution $Ax = b$ by using the QR decomposition. Let us use the QR decomposition from above, and solve the system $$ Rx = Q^Tb. $$ As $$ \begin{bmatrix} \frac{1}{2} & -\frac{3}{\sqrt{12}} & 0 \\ \frac{1}{2} & \frac{1}{\sqrt{12}} & -\frac{2}{\sqrt{6}} \\ \frac{1}{2} & \frac{1}{\sqrt{12}} & \frac{1}{\sqrt{6}} \\ \frac{1}{2} & \frac{1}{\sqrt{12}} & \frac{1}{\sqrt{6}} \end{bmatrix}^T\begin{bmatrix} 1 \\ 1 \\ 1 \\ 0 \end{bmatrix} = \begin{bmatrix} \frac{3}{2} \\ -\frac{1}{2\sqrt{3}} \\ -\frac{1}{\sqrt{6}}, \end{bmatrix} $$ we are looking at the system $$ \begin{bmatrix} 2 & \frac{3}{2} & 1 \\ 0 & \frac{3}{\sqrt{12}} & \frac{2}{\sqrt{12}} \\ 0 & 0 & \frac{2}{\sqrt{6}} \end{bmatrix}x =\begin{bmatrix} \frac{3}{2} \\ -\frac{1}{2\sqrt{3}} \\ -\frac{1}{\sqrt{6}} \end{bmatrix}. $$ Solving this system yields that $$ x_0 = \begin{bmatrix} 1 \\ 0 \\ -\frac{1}{2} \end{bmatrix} $$ is a least-squares solution to $Ax = b$.
Let us set this system up in python and use numpy.linalg.solve.
import numpy as np
# Define matrix and vector
A = np.array([[1,0,0],[1,1,0],[1,1,1],[1,1,1]])
b = np.array([[1],[1],[1],[0]])
# Take the QR decomposition of A
Q, R = np.linalg.qr(A)
# Solve the linear system Rx = Q.T b
beta = np.linalg.solve(R,Q.T @ b)
This yields
beta
array([[ 1.00000000e+00],
[ 6.40987562e-17],
[-5.00000000e-01]])
which (basically) agrees with our exact least-squares solution.
Note that numpy.linalg.lstsq still gives a ever so slightly different result.
np.linalg.lstsq(A,b)[0]
array([[ 1.00000000e+00],
[ 2.22044605e-16],
[-5.00000000e-01]])
Let's go back to the house example. While we're at it, let's get used to using pandas to make a dataframe.
import numpy as np
import pandas as pd
# First let us make a dictionary incorporating our data.
# Each entry corresponds to a column (feature of our data)
data = {
'Square ft': [1600, 2100, 1550, 1600, 2000],
'Bedrooms': [3, 4, 2, 3, 4],
'Price': [500, 650, 475, 490, 620]
}
# Create a pandas DataFrame
df = pd.DataFrame(data)
# Create our matrix X and our target y
X = df[["Square ft", "Bedrooms"]].to_numpy()
y = df[["Price"]].to_numpy()
# Augment X with a column of 1's (intercept)
X_aug = np.hstack((np.ones((X.shape[0], 1)), X))
# Perform QR decomposition
Q, R = np.linalg.qr(X_aug)
# Solve the upper triangular system Rx = Q^Ty
beta = np.linalg.solve(R, Q.T @ y)
Let's look at the output.
print(f"Q = {Q} \n\nR = {R} \n\nbeta = {beta}")
Q = [[-0.4472136 0.32838365 0.40496317] [-0.4472136 -0.63745061 -0.22042299] [-0.4472136 0.42496708 -0.7689174 ] [-0.4472136 0.32838365 0.40496317] [-0.4472136 -0.44428376 0.17941406]] R = [[-2.23606798e+00 -3.95784032e+03 -7.15541753e+00] [ 0.00000000e+00 -5.17687164e+02 -1.50670145e+00] [ 0.00000000e+00 0.00000000e+00 7.27908474e-01]] beta = [[-3.05053797e-13] [ 3.00000000e-01] [ 5.00000000e+00]]
As we can see, the least-squares solution agrees with what we got by hand and by other python methods (if we agree that the tiny first component is essentially zero).
The QR decomposition of a matrix is also useful for computing orthogonal projections.
Theorem. Let $A$ be an $m \times n$ matrix with full column rank. If $A = QR$ is a QR decomposition, then $QQ^T$ is the projection onto the column space of $A$, i.e., $QQ^Tb = \text{Proj}_{\text{Col}(A)}b$ for all $b \in \mathbb{R}^m$.
Let's see what our range projections are for the matrices above. Note that the first example above will have the orthogonal projection just being $$ \begin{bmatrix} 1 \\ & 1 \\ & & 1\\ & & & 0 \end{bmatrix}. $$ Let's look at the other matrix.
Example. Working with the matrix $$ A = \begin{bmatrix} 1 & 0 & 0 \\ 1 & 1 & 0 \\ 1 & 1 & 1 \\ 1 & 1 & 1 \end{bmatrix}, $$ the projection onto the column space if given by $$ QQ^T = \begin{bmatrix} 1 \\ & 1 \\ & & \frac{1}{2} & \frac{1}{2} \\ & & \frac{1}{2} & \frac{1}{2} \end{bmatrix}. $$ This is a well-understood projection: it is the direct sum of the identity on $\mathbb{R}^2$ and the projection onto the line $y = x$ in $\mathbb{R}^2$.
Now let's use python to implement the projection.
import numpy as np
# Create our matrix A
A = np.array([[1,0,0],[1,1,0],[1,1,1],[1,1,1]])
# Take the QR decomposition
Q, R = np.linalg.qr(A)
# Create the range projection
P = Q @ Q.T
P
array([[1.00000000e+00, 2.89687929e-17, 2.89687929e-17, 2.89687929e-17],
[2.89687929e-17, 1.00000000e+00, 7.07349921e-17, 7.07349921e-17],
[2.89687929e-17, 7.07349921e-17, 5.00000000e-01, 5.00000000e-01],
[2.89687929e-17, 7.07349921e-17, 5.00000000e-01, 5.00000000e-01]])
The output gives
array([[1.00000000e+00, 2.89687929e-17, 2.89687929e-17, 2.89687929e-17],
[2.89687929e-17, 1.00000000e+00, 7.07349921e-17, 7.07349921e-17],
[2.89687929e-17, 7.07349921e-17, 5.00000000e-01, 5.00000000e-01],
[2.89687929e-17, 7.07349921e-17, 5.00000000e-01, 5.00000000e-01]])
As we can see, the two off-diagonal blocks are all tiny, hence we treat them as zero. Note that if they were not actually zero, then this wouldn't actually be a projection. This can cause some problems.
Let's write a function to implement this, assuming that columns of A are linearly independent.
import numpy as np
def proj_onto_col_space(A):
# Take the QR decomposition
Q,R = np.linalg.qr(A)
# The projection is just Q @ Q.T
P = Q @ Q.T
return P
We'll come back to this later. We should really be incorporating some sort of error tolerance so that things are super super tiny can actually just be sent to zero.
Remark. Another way to get the projection onto the column space of an $n \times p$ matrix $A$ of full column rank is to take $$ P = A(A^TA)^{-1}A^T. $$ Indeed, let $b \in \mathbb{R}^n$ and let $x_0 \in \mathbb{R}^p$ be a solution to the normal equations $$ A^TAx_0 = A^Tb. $$ Then $x_0 = (A^TA)^{-1}A^Tb$ and so $Ax_0 = A(A^TA^{-1})A^Tb$ is the (unique!) vector in the column space of $A$ which is closest to $b$, i.e., the projection of $b$ onto the column space of $A$. However, taking transposes, multiplying, and inverting is not what we would like to do numerically.
Singular Value Decomposition¶
The SVD is a very important matrix decomposition in both data science and linear algebra.
Theorem. For any matrix $n \times p$ matrix $X$, there exist an orthogonal $n \times n$ matrix $U$, an orthogonal $p \times p$ matrix $V$, and a diagonal $n \times p$ matrix $\Sigma$ with non-negative entries such that $$ X = U\Sigma V^T. $$
- The columns of $U$ are left left singular vectors.
- The columns of $V$ are the right singular vectors.
- $\Sigma$ has singular values $\sigma_1 \geq \sigma_2 \geq \cdots \geq \sigma_r > 0$ on its diagonal, where $r$ is the rank of $X$.
Remark. The SVD is clearly a generalization of matrix diagonalization, but it also generalizes the polar decomposition of a matrix. Recall that every $n \times n$ matrix $A$ can be written as $A = UP$ where $U$ is orthogonal (or unitary) and $P$ is a positive matrix. This is because if $$ A = U_0\Sigma V^T $$ is the SVD for $A$, then $\Sigma$ is an $n \times n$ diagonal matrix with non-negative entries, hence any orthogonal conjugate of it is positive, and so $$ A = (U_0V^T)(V\Sigma V^T). $$ Take $U = U_0V^T$ and $P = V\Sigma V^T$.
By hand, the algorithm for computing an SVD is as follows.
- Both $AA^T$ and $A^TA$ are symmetric (they are positive in fact), and so they can be orthogonally diagonalized; one can form an orthogonal basis of eigenvectors. Let $v_1,\dots,v_p$ be an orthonormal basis of eigenvectors for $\mathbb{R}^p$ which correspond to eigenvectors of $A^TA$ in decreasing order. Suppose that $A^TA$ has $r$ non-zero eigenvalues. Let $V$ be the matrix whose columns contain the $v_i$'s. This gives our right singular vectors and our singular values.
- Let $u_i = \frac{1}{\sigma_i}Av_i$ for $i = 1,\dots,r$, and extend this collection of vectors to an orthonormal basis for $\mathbb{R}^n$ if necessary. Let $U$ be the corresponding matrix.
- Let $\Sigma$ be the $n \times p$ matrix whose diagonal entries are $\sigma_1 \geq \sigma_2 \geq \cdots \geq \sigma_r$, and then zeroes if necessary.
Example. Let us compute the SVD of $$ A = \begin{bmatrix} 3 & 2 & 2 \\ 2 & 3 & -2 \end{bmatrix}. $$ First we note that $$ A^TA = \begin{bmatrix} 13 & 12 & 2 \\ 12 & 13 & -2 \\ 2 & -2 & 8 \end{bmatrix}, $$ which has eigenvalues $25,9,0$ with corresponding eigenvectors $$ \begin{bmatrix} 1 \\ 1 \\ 0 \end{bmatrix}, \begin{bmatrix} 1 \\ -1 \\ 4 \end{bmatrix}, \begin{bmatrix} -2 \\ 2 \\ 1 \end{bmatrix}. $$ Normalizing, we get $$ V = \begin{bmatrix} \frac{1}{\sqrt{2}} & \frac{1}{3\sqrt{2}} & -\frac{2}{3} \\ \frac{1}{\sqrt{2}} & -\frac{1}{3\sqrt{2}} & \frac{2}{3} \\ 0 & \frac{4}{3\sqrt{2}} & \frac{1}{3} \end{bmatrix}. $$ Now we set $u_1 = \frac{1}{5}Av_1$ and $u_2 = \frac{1}{3}Av_2$ to get $$ U = \begin{bmatrix} \frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} \\ \frac{1}{\sqrt{2}} & -\frac{1}{\sqrt{2}} \end{bmatrix}. $$ So $$ A = \begin{bmatrix} \frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} \\ \frac{1}{\sqrt{2}} & -\frac{1}{\sqrt{2}} \end{bmatrix} \begin{bmatrix} 5 & 0 & 0 \\ 0 & 3 & 0 \end{bmatrix} \begin{bmatrix} \frac{1}{\sqrt{2}} & \frac{1}{3\sqrt{2}} & -\frac{2}{3} \\ \frac{1}{\sqrt{2}} & -\frac{1}{3\sqrt{2}} & \frac{2}{3} \\ 0 & \frac{4}{3\sqrt{2}} & \frac{1}{3} \end{bmatrix}^T $$ is our SVD decomposition.
We note that in practice, we avoid the computation of $X^TX$ because if the entries of $X$ have errors, then these errors will be squared in $X^TX$. There are better computational tools to get singular values and singular vectors which are more accurate. This is what our python tools will use.
Let's use numpy.linalg.svd for the above matrix.
import numpy as np
#Define our matrix
A = np.array([[3,2,2],[2,3,-2]])
# Take the SVD
U, S, Vh = np.linalg.svd(A)
Our SVD matrices are
print(f"U = {U}\n\nS = {S}\n\nVh.T = {Vh.T}")
U = [[-0.70710678 -0.70710678] [-0.70710678 0.70710678]] S = [5. 3.] Vh.T = [[-7.07106781e-01 -2.35702260e-01 -6.66666667e-01] [-7.07106781e-01 2.35702260e-01 6.66666667e-01] [-6.47932334e-17 -9.42809042e-01 3.33333333e-01]]
Because the eigenvalues of the hermitian squares of $$ \begin{bmatrix} 1 & 1 & 1\\ 0 & 1 & 1 \\ 0 & 0 & 1 \\ 0 & 0 & 0 \end{bmatrix} \text{ and } \begin{bmatrix} 1 & 0 & 0 \\ 1 & 1 & 0 \\ 1 & 1 & 1 \\ 1 & 1 & 1 \end{bmatrix} $$ are quite atrocious, an exact SVD decomposition is difficult to compute by hand. However, we can of course use python.
import numpy as np
# Define our matrices
A = np.array([[1,1,1],[0,1,1],[0,0,1],[0,0,0]])
B = np.array([[1,0,0],[1,1,0],[1,1,1],[1,1,1]])
# SVD decomposition
U_A, S_A, Vh_A = np.linalg.svd(A)
U_B, S_B, Vh_B = np.linalg.svd(B)
The resulting matrices are
print(f"U_A = {U_A}\n\nS_A = {S_A}\n\nVh_A.T = {Vh_A.T}\n\nU_B = {U_B}\n\nS_B = {S_B}\n\nVh_B.T = {Vh_B.T}")
U_A = [[ 0.73697623 0.59100905 0.32798528 0. ] [ 0.59100905 -0.32798528 -0.73697623 0. ] [ 0.32798528 -0.73697623 0.59100905 0. ] [ 0. 0. 0. 1. ]] S_A = [2.2469796 0.80193774 0.55495813] Vh_A.T = [[ 0.32798528 0.73697623 0.59100905] [ 0.59100905 0.32798528 -0.73697623] [ 0.73697623 -0.59100905 0.32798528]] U_B = [[-2.41816250e-01 7.12015746e-01 -6.59210496e-01 0.00000000e+00] [-4.52990541e-01 5.17957311e-01 7.25616837e-01 6.71536163e-17] [-6.06763739e-01 -3.35226641e-01 -1.39502200e-01 -7.07106781e-01] [-6.06763739e-01 -3.35226641e-01 -1.39502200e-01 7.07106781e-01]] S_B = [2.8092118 0.88646771 0.56789441] Vh_B.T = [[-0.67931306 0.63117897 -0.37436195] [-0.59323331 -0.17202654 0.7864357 ] [-0.43198148 -0.75632002 -0.49129626]]
Another final note is that the operator norm of a matrix $A$ agrees with its largest singular value.
Pseudoinverses and using the SVD¶
The SVD can be used to determine a least-squares solution for a given system. Recall that if $v_1,\dots,v_p$ is an orthonormal basis for $\mathbb{R}^p$ consisting of eigenvectors of $A^TA$, arranged so that they correspond to eigenvalues $\sigma_1 \geq \sigma_2 \geq \cdots \geq \sigma_r$, then $\{Av_1,\dots,Av_r\}$ is an orthogonal basis for the column space of $A$. In essence, this means that when we have our left singular vectors $u_1,\dots,u_n$ (constructed based on our algorithm as above), we have that the first $r$ vectors form an orthonormal basis for the column space of $A$, and that the remaining $n - r$ vectors form an orthonormal basis for the perp of the column space of $A$ (which is also equal to the nullspace of $A^T$).
Definition. Let $A$ be an $n \times p$ matrix and suppose that the rank of $A$ is $r \leq \min\{n,p\}$. Suppose that $A = U\Sigma V^T$ is the SVD, where the singular values are decreasing. Partition $$ U = \begin{bmatrix} U_r & U_{n-r} \end{bmatrix} \text{ and } V = \begin{bmatrix} V_r & V_{p-r} \end{bmatrix} $$ into submatrices, where $U_r$ and $V_r$ are the matrices whose columns are the first $r$ columns of $U$ and $V$ respectively. So $U_r$ is $n \times r$ and $V_r$ is $p \times r$. Let $D$ be the diagonal $r \times r$ matrices whose diagonal entries are $\sigma_1,\dots, \sigma_r$, so that $$ \Sigma = \begin{bmatrix} D & 0 \\ 0 & 0 \end{bmatrix} $$ and note that $$ A = U_rDV_r^T. $$ We call this the reduced singular value decomposition of $A$. Note that $D$ is invertible, and its inverse is simply $$ D = \begin{bmatrix} \sigma_1^{-1} \\ & \sigma_2^{-1} \\ & & \ddots \\ & & & \sigma_r^{-1} \end{bmatrix}. $$ The pseudoinverse (or Moore-Penrose inverse) of $A$ is the matrix $$ A^+ = V_rD^{-1}U_r^T. $$
We note that the pseudoinverse $A^+$ is a $p \times n$ matrix.
With the pseudoinverse, we can actually find least-squares solutions quite easily. Indeed, if we are looking for the least-squares solution to the system $Ax = b$, define $$ x_0 = A^+b. $$ Then $$ \begin{split} Ax_0 &= (U_rDV_r^T)(V_rD^{-1}U_r^Tb) \\ &= U_rDD^{-1}U_r^Tb \\ &= U_rU_r^Tb \end{split} $$ As mentioned before, the columns of $U_r$ form an orthonormal basis for the column space of $A$ and so $U_rU_r^T$ is the orthogonal projection onto the range of $A$. That is, $Ax_0$ is precisely the projection of $b$ onto the column space of $A$, meaning that this yields a least-squares solution. This gives the following.
Theorem. Let $A$ be an $n \times p$ matrix and $b \in \mathbb{R}^n$. Then $$ x_0 = A^+b$$ is a least-squares solution to $Ax = b$.
Taking pseudoinverses is quite involved. We'll do one example by hand, and then use python -- and we'll see something go wrong! There is a function numpy.linalg.pinv in numpy that will take a pseudoinverse. We can also just use numpy.linalg.svd and do the process above.
Example. We have the following SVD $A = U\Sigma V^T$. $$ \begin{bmatrix} 1 & 1 & 2\\ 0 & 1 & 1 \\ 1 & 0 & 1 \\ 0 & 0 & 0 \end{bmatrix} = \begin{bmatrix} \sqrt{\frac{2}{3}} & 0 & 0 & -\frac{1}{\sqrt{3}} \\ \frac{1}{\sqrt{6}} & \frac{1}{\sqrt{2}} & 0 & \frac{1}{\sqrt{3}} \\ \frac{1}{\sqrt{6}} & -\frac{1}{\sqrt{2}} & 0 & \frac{1}{\sqrt{3}} \\ 0 & 0 & 1 & 0 \end{bmatrix} \begin{bmatrix} 3 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}\begin{bmatrix} \frac{1}{\sqrt{6}} & -\frac{1}{\sqrt{2}} & -\frac{1}{\sqrt{3}} \\ \frac{1}{\sqrt{6}} & \frac{1}{\sqrt{2}} & -\frac{1}{\sqrt{3}} \\ \sqrt{\frac{2}{3}} & 0 & \frac{1}{\sqrt{3}} \end{bmatrix}^T. $$ Can we find a least-squares solution to $Ax = b$, where $$ b = \begin{bmatrix} 1 \\ 1 \\ 1 \\ 1 \end{bmatrix}? $$ The pseudoinverse of $A$ is $$ \begin{split} A^+ &= \begin{bmatrix} \frac{1}{\sqrt{6}} & -\frac{1}{\sqrt{2}} \\ \frac{1}{\sqrt{6}} & \frac{1}{\sqrt{2}} \\ \sqrt{\frac{2}{3}} & 0 \end{bmatrix} \begin{bmatrix} 3 \\ & 1 \end{bmatrix} \begin{bmatrix} \sqrt{\frac{2}{3}} & 0 \\ \frac{1}{\sqrt{6}} & \frac{1}{\sqrt{2}} \\ \frac{1}{\sqrt{6}} & -\frac{1}{\sqrt{2}} \\ 0 & 0 \end{bmatrix}^T \\ &= \begin{bmatrix} \frac{1}{9} & -\frac{4}{9} & \frac{5}{9} & 0 \\ \frac{1}{9} & \frac{5}{9} & -\frac{4}{9} & 0 \\ \frac{2}{9} & \frac{1}{9} & \frac{1}{9} & 0\end{bmatrix}, \end{split} $$ and so a least-squares solution is given by $$ \begin{split} x_0 &= A^+b \\ &= \begin{bmatrix} \frac{1}{9} & -\frac{4}{9} & \frac{5}{9} & 0 \\ \frac{1}{9} & \frac{5}{9} & -\frac{4}{9} & 0 \\ \frac{2}{9} & \frac{1}{9} & \frac{1}{9} & 0\end{bmatrix}\begin{bmatrix} 1 \\ 1 \\ 1 \\ 1 \end{bmatrix} \\ &= \begin{bmatrix} \frac{2}{9} \\ \frac{2}{9} \\ \frac{4}{9} \end{bmatrix}. \end{split} $$
Now let's do this with python, and see an example of how things can go wrong. We'll try to take the pseudoinverse manually first.
import numpy as np
# Create our matrix A and our target b
A = np.array([[1,1,2],[0,1,1],[1,0,1],[0,0,0]])
b = np.array([[1],[1],[1],[1]])
# Take the SVD decomposition
U, S, Vh = np.linalg.svd(A)
# Prepare the pseudoinverse
# Recall that we invert the non-zero diagonal entries of the diagonal matrix.
# So we first build S_inv to be the appropriate size
S_inv = np.zeros((Vh.shape[0], U.shape[0]))
# We then fill in the appropriate values on the diagonal
S_inv[:len(S), :len(S)] = np.diag(1/S)
# Build the pseudoinverse
A_pinv = Vh.T @ S_inv @ U.T
# Compute the least-squares solution
beta = A_pinv @ b
What is the result?
beta
array([[ 2.74080345e+15],
[ 2.74080345e+15],
[-2.74080345e+15]])
This is WAY off the mark. So what happened? Well, when we look at our singular values, we have
S
array([3.00000000e+00, 1.00000000e+00, 1.21618839e-16])
As we got this matrix numerically, the last entry is actually non-zero, but tiny. This isn't exactly what's going on since we know that the rank of A is 2. So when we invert the singular values and throw them on the diagonal, have 1/1.21618839e-16 which is a very large value. This value then messes up the rest of the computation. So how do we fix this? One can set tolerances in numpy, but we'll get to that later. Let's just note that numpy.linalg.pinv will already incorporate this. Let's see what we get.
import numpy as np
# Create our matrix A and our target b
A = np.array([[1,1,2],[0,1,1],[1,0,1],[0,0,0]])
b = np.array([[1],[1],[1],[1]])
# Build the pseudoinverse
A_pinv = np.linalg.pinv(A)
# Compute the least-squares solution
beta = A_pinv @ b
print(f"A_pinv={A_pinv}\n\nbeta={beta}")
A_pinv=[[ 0.11111111 -0.44444444 0.55555556 0. ] [ 0.11111111 0.55555556 -0.44444444 0. ] [ 0.22222222 0.11111111 0.11111111 0. ]] beta=[[0.22222222] [0.22222222] [0.44444444]]
The Condition Number¶
Numerical calculations involving matrix equations are quite reliable if we use the SVD. This is because the orthogonal matrices $U$ and $V$ preserve lengths and angles, leaving the stability of the problem to be governed by the singular values of the matrix $X$. Recall that if $X = U\Sigma V^T$, then solving the least-squares problem involves dividing by the non-zero singular values $\sigma_i$ of $X$. If these values are very small, their inverses become very large, and this will amplify any numerical errors.
Definition. Let $X$ be an $n \times p$ matrix and let $\sigma_1 \geq \cdots \geq \sigma_r$ be the non-zero singular values of $X$. The condition number of $X$ is the quotient $$ \kappa(X) = \frac{\sigma_1}{\sigma_r} $$ of the largest and smallest non-zero singular values.
A condition number close to 1 indicates a well-conditioned problem, while a large condition number indicates that small perturbations in data may lead to large changes in computation. Geometrically, $\kappa(X)$ measures how much $X$ distorts space.
Example. Consider the matrices $$ A = \begin{bmatrix} 1 \\ & 1 \end{bmatrix} \text{ and } B = \begin{bmatrix} 1 \\ & \frac{1}{10^6} \end{bmatrix}. $$ The condition numbers are $$ \kappa(A) = 1 \text{ and } \kappa(B) = 10^6. $$ Inverting $X_2$ includes dividing by $\frac{1}{10^6}$, which will amplify errors by $10^6$.
Let's look our main example in python by using numpy.linalg.cond.
import numpy as np
import pandas as pd
# First let us make a dictionary incorporating our data.
# Each entry corresponds to a column (feature of our data)
data = {
'Square ft': [1600, 2100, 1550, 1600, 2000],
'Bedrooms': [3, 4, 2, 3, 4],
'Price': [500, 650, 475, 490, 620]
}
# Create a pandas DataFrame
df = pd.DataFrame(data)
# Create out matrix X
X = df[['Square ft', 'Bedrooms']].to_numpy()
# Check the condition number
cond_X = np.linalg.cond(X)
Let's see what we got.
cond_X
np.float64(4329.082589067693)
so this is quite a high condition number! This should be unsurprising, as clearly the number of bedrooms is correlated to the size of a house (especially so in our small toy example).