Book contents
Contents
raw Math

Eigenvalues and eigenvectors only tell the full geometric story for a square matrix that has enough independent eigenvectors to span the space. A defective matrix has too few of them, and a rectangular matrix does not even map a space back into itself, so "eigenvector" is not even defined. The singular value decomposition (SVD) answers a more modest but universally valid question instead: no matter what a matrix does, which orthogonal input directions does it map to orthogonal output directions, and by how much does it stretch each of them? Every matrix, square or not, diagonalisable or not, has an answer.

What Every Matrix Does to the Unit Circle

Take any \(2\times 2\) matrix \(\mathbf{A}\) and feed it every unit vector \(\mathbf{x}\) with \(\|\mathbf{x}\|_2=1\), i.e. every point on the unit circle. The image \(\mathbf{Ax}\) always traces out an ellipse. This is not obvious from the entries of \(\mathbf{A}\) alone, but it becomes obvious once the right pair of perpendicular directions is found: two orthogonal input directions \(\mathbf{v}_1,\mathbf{v}_2\) whose images \(\mathbf{Av}_1,\mathbf{Av}_2\) are again orthogonal, pointing exactly along the major and minor axes of the ellipse.

The unit circle with the two orthogonal directions \(\mathbf{v}_1,\mathbf{v}_2\) that will map to the ellipse's axes.
The image under \(\mathbf{A}=\left(\begin{smallmatrix}1&1\\0&1\end{smallmatrix}\right)\): an ellipse whose axes are \(\sigma_1\mathbf{u}_1\) and \(\sigma_2\mathbf{u}_2\).

The lengths \(\sigma_1\ge\sigma_2\ge 0\) of the semi-axes are the singular values of \(\mathbf{A}\), the input directions \(\mathbf{v}_i\) are its right singular vectors, and the output directions \(\mathbf{u}_i\) are its left singular vectors. Collecting the \(\mathbf{v}_i\) as columns of an orthogonal matrix \(\mathbf{V}\), the \(\mathbf{u}_i\) as columns of an orthogonal matrix \(\mathbf{U}\), and the \(\sigma_i\) on the diagonal of \(\mathbf{\Sigma}\), the relation \(\mathbf{Av}_i=\sigma_i \mathbf{u}_i\) for every \(i\) is exactly the matrix equation \(\mathbf{AV}=\mathbf{U}\mathbf{\Sigma}\), i.e.

\[\mathbf{A}=\mathbf{U}\mathbf{\Sigma}\mathbf{V}^T.\]

Reading the product right to left as it acts on a vector \(\mathbf{x}\) gives a three-step geometric recipe that works for any matrix, not just the ones with a nice eigenbasis: \(\mathbf{V}^T\mathbf{x}\) first rotates (or reflects) \(\mathbf{x}\) into the coordinate system where \(\mathbf{v}_1,\ldots,\mathbf{v}_n\) become the standard axes; \(\mathbf{\Sigma}\) then stretches each of those axes independently by \(\sigma_i\); and finally \(\mathbf{U}\) rotates (or reflects) the result into its true output orientation. Every matrix, however complicated, is nothing more than rotate, stretch along perpendicular axes, rotate again — the only question the SVD answers is exactly which axes and which stretch factors. If \(\mathbf{A}\) is singular, one or more \(\sigma_i\) are zero and the picture degenerates accordingly: a singular \(2\times 2\) matrix flattens the unit circle onto a line segment instead of an ellipse, because the corresponding axis is stretched by a factor of \(0\).

Definition

This decomposition is not merely a convenient parametrisation, it exists for every real matrix \(\mathbf{A}\in\mathbb{R}^{m\times n}\), and it can be built directly from facts already established for eigenvalues and eigenvectors. Consider the square matrix

\[\mathbf{S}:=\mathbf{A}^T\mathbf{A}\in\mathbb{R}^{n\times n}.\]

\(\mathbf{S}\) is symmetric (\(\mathbf{S}^T=\mathbf{A}^T\mathbf{A}^{TT}=\mathbf{A}^T\mathbf{A}=\mathbf{S}\)) and positive semidefinite, since for any \(\mathbf{x}\),

\[\mathbf{x}^T\mathbf{S}\mathbf{x}=\mathbf{x}^T\mathbf{A}^T\mathbf{A}\mathbf{x}=(\mathbf{Ax})^T(\mathbf{Ax}) =\|\mathbf{Ax}\|_2^2\ge 0.\]

By the symmetric-matrix and positive-semidefinite properties established earlier, \(\mathbf{S}\) therefore has an orthonormal eigenbasis \(\mathbf{v}_1,\ldots,\mathbf{v}_n\) with real, non-negative eigenvalues \(\lambda_1\ge\cdots\ge\lambda_n\ge 0\). Define the singular values

\[\sigma_i:=\sqrt{\lambda_i},\qquad i=1,\ldots,n,\]

and let \(r\le n\) be the number of strictly positive \(\sigma_i\) (it will turn out that \(r=\mathrm{rank} (\mathbf{A})\)). For each \(i=1,\ldots,r\), define

\[\mathbf{u}_i:=\frac{1}{\sigma_i}\mathbf{Av}_i.\]

These vectors are automatically orthonormal: using \(\mathbf{Sv}_j=\lambda_j\mathbf{v}_j\) and the orthonormality of the \(\mathbf{v}_i\),

\[\mathbf{u}_i\cdot\mathbf{u}_j=\frac{1}{\sigma_i\sigma_j}\mathbf{v}_i^T\mathbf{A}^T\mathbf{A}\mathbf{v}_j =\frac{\lambda_j}{\sigma_i\sigma_j}\,\mathbf{v}_i\cdot\mathbf{v}_j =\begin{cases}1 & i=j\\0& i\neq j\end{cases}.\]

If \(r<m\), extend \(\mathbf{u}_1,\ldots,\mathbf{u}_r\) to a full orthonormal basis of \(\mathbb{R}^m\) (any orthonormal set can be completed this way). If \(r<n\), the remaining \(\mathbf{v}_{r+1},\ldots,\mathbf{v}_n\) already satisfy \(\mathbf{Av}_i=\mathbf{0}\), because \(\sigma_i=0\) means \(\|\mathbf{Av}_i\|_2^2=\mathbf{v}_i^T\mathbf{Sv}_i=\lambda_i=0\). Collecting everything into \(\mathbf{U}=(\mathbf{u}_1,\ldots,\mathbf{u}_m)\), \(\mathbf{V}=(\mathbf{v}_1,\ldots,\mathbf{v}_n)\), and the \(m\times n\) matrix \(\mathbf{\Sigma}\) with \(\sigma_i\) on the diagonal (zero elsewhere), the relations \(\mathbf{Av}_i=\sigma_i\mathbf{u}_i\) for every column give \(\mathbf{AV}=\mathbf{U\Sigma}\), and since \(\mathbf{V}\) is orthogonal (\(\mathbf{V}^{-1}=\mathbf{V}^T\)),

\[\boxed{\mathbf{A}=\mathbf{U}\mathbf{\Sigma}\mathbf{V}^T}\]

with \(\mathbf{U}\in\mathbb{R}^{m\times m}\) and \(\mathbf{V}\in\mathbb{R}^{n\times n}\) orthogonal, and \(\mathbf{\Sigma}\in\mathbb{R}^{m\times n}\) diagonal with non-negative entries \(\sigma_1\ge\cdots\ge \sigma_{\min(m,n)}\ge 0\). Unlike an eigendecomposition, this construction never required \(\mathbf{A}\) to be square, symmetric, or diagonalisable at any point.

A Concrete Example

The shear matrix from the visualization above,

\[\mathbf{A}=\begin{pmatrix}1&1\\0&1\end{pmatrix},\]

illustrates the construction step by step. First,

\[\mathbf{S}=\mathbf{A}^T\mathbf{A}=\begin{pmatrix}1&0\\1&1\end{pmatrix}\begin{pmatrix}1&1\\0&1 \end{pmatrix}=\begin{pmatrix}1&1\\1&2\end{pmatrix}.\]

Its characteristic polynomial is \(\lambda^2-3\lambda+1\), with roots \(\lambda=\frac{3\pm\sqrt5}{2}\). Writing \(\varphi=\frac{1+\sqrt5}{2}\) for the golden ratio, and using the identity \(\varphi^2=\varphi+1\), these roots are exactly \(\lambda_1=\varphi^2\) and \(\lambda_2=\varphi^{-2}\), so the singular values are the golden ratio itself and its reciprocal:

\[\sigma_1=\varphi\approx 1.618,\qquad \sigma_2=\frac1\varphi\approx 0.618.\]

Solving \((\mathbf{S}-\lambda_i\mathbf{I})\mathbf{v}_i=\mathbf{0}\) and normalizing gives the right singular vectors

\[\mathbf{v}_1=\frac{1}{\sqrt{2+\varphi}}\begin{pmatrix}1\\ \varphi\end{pmatrix}\approx\begin{pmatrix}0.526\\ 0.851\end{pmatrix},\qquad \mathbf{v}_2=\frac{1}{\sqrt{2+\varphi}}\begin{pmatrix}\varphi\\ -1\end{pmatrix}\approx\begin{pmatrix}0.851\\ -0.526\end{pmatrix}.\]

The left singular vectors follow from \(\mathbf{u}_i=\mathbf{Av}_i/\sigma_i\):

\[\mathbf{u}_1=\frac{1}{\sigma_1}\mathbf{Av}_1\approx\begin{pmatrix}0.851\\0.526\end{pmatrix},\qquad \mathbf{u}_2=\frac{1}{\sigma_2}\mathbf{Av}_2\approx\begin{pmatrix}0.526\\-0.851\end{pmatrix}.\]

One checks directly that \(\mathbf{u}_1\cdot\mathbf{u}_2\approx 0\), and that \(\mathbf{U}\mathbf{\Sigma}\mathbf{V}^T\) reconstructs \(\mathbf{A}\) exactly. Notice how the entries of \(\mathbf{u}_1\) and \(\mathbf{v}_2\) match up to a sign, and likewise for \(\mathbf{u}_2\) and \(\mathbf{v}_1\): this near-symmetry is a coincidence of this particular shear matrix, not a general rule.

Right and Left Singular Vectors as Eigenvectors

The construction above shows \(\mathbf{v}_i\) are eigenvectors of \(\mathbf{A}^T\mathbf{A}\), but the same relation can be read off directly from the boxed decomposition. Since \(\mathbf{U}\) is orthogonal (\(\mathbf{U}^T\mathbf{U}=\mathbf{I}\)) and \(\mathbf{\Sigma}^T\mathbf{\Sigma}\) is diagonal with entries \(\sigma_i^2\),

\[\mathbf{A}^T\mathbf{A}=(\mathbf{U\Sigma V}^T)^T(\mathbf{U\Sigma V}^T) =\mathbf{V}\mathbf{\Sigma}^T\mathbf{U}^T\mathbf{U}\mathbf{\Sigma}\mathbf{V}^T =\mathbf{V}(\mathbf{\Sigma}^T\mathbf{\Sigma})\mathbf{V}^T,\]

which is precisely the eigenvalue decomposition of \(\mathbf{A}^T\mathbf{A}\): its columns \(\mathbf{v}_i\) are eigenvectors with eigenvalues \(\sigma_i^2\). By the symmetric argument (using \(\mathbf{V}^T\mathbf{V}=\mathbf{I}\)),

\[\mathbf{A}\mathbf{A}^T=\mathbf{U}(\mathbf{\Sigma}\mathbf{\Sigma}^T)\mathbf{U}^T,\]

so the \(\mathbf{u}_i\) are eigenvectors of \(\mathbf{AA}^T\), with the very same eigenvalues \(\sigma_i^2\). This gives a second, equivalent route to the SVD: eigendecompose either \(\mathbf{A}^T\mathbf{A}\) or \(\mathbf{AA}^T\) — whichever is smaller — and derive the other side from it (see "Robust Computation" below for why deriving one from the other, rather than eigendecomposing both independently, matters in practice).

Full versus Compact SVD, and the Rectangular Case

The construction keeps all \(n\) columns of \(\mathbf{V}\) and all \(m\) columns of \(\mathbf{U}\), including those with \(\sigma_i=0\); this is the full SVD. Since a zero singular value contributes nothing to \(\mathbf{U\Sigma V}^T\), it is common to drop the last \(n-r\) columns of \(\mathbf{V}\), the last \(m-r\) columns of \(\mathbf{U}\), and the zero rows/columns of \(\mathbf{\Sigma}\), keeping only

\[\mathbf{A}=\mathbf{U}_r\mathbf{\Sigma}_r\mathbf{V}_r^T,\qquad \mathbf{U}_r\in\mathbb{R}^{m\times r},\quad\mathbf{\Sigma}_r\in\mathbb{R}^{r\times r},\quad\mathbf{V}_r\in\mathbb{R}^{n\times r},\]

the compact SVD, where \(r=\mathrm{rank}(\mathbf{A})\). For a square, full-rank matrix (\(r=m=n\)) the two coincide. For a genuinely rectangular or rank-deficient \(\mathbf{A}\), the discarded columns are not wasted information: the trailing columns of \(\mathbf{V}\) span the null space of \(\mathbf{A}\) (solutions of \(\mathbf{Ax}=\mathbf{0}\)), and the trailing columns of \(\mathbf{U}\) span the orthogonal complement of the column space of \(\mathbf{A}\) — exactly the four fundamental subspaces of \(\mathbf{A}\), all produced by one decomposition.

Seen this way, the SVD is the natural completion of the picture built around the four fundamental subspaces: \(\mathbf{v}_1,\ldots,\mathbf{v}_r\) is an orthonormal basis of the row space and \(\mathbf{v}_{r+1},\ldots, \mathbf{v}_n\) one of the null space, together spanning all of \(\mathbb{R}^n\); \(\mathbf{u}_1,\ldots, \mathbf{u}_r\) is an orthonormal basis of the column space and \(\mathbf{u}_{r+1},\ldots,\mathbf{u}_m\) one of the left null space, together spanning \(\mathbb{R}^m\). The relation \(\mathbf{Av}_i=\sigma_i\mathbf{u}_i\) says precisely how \(\mathbf{A}\) carries the first basis to the second, coordinate by coordinate: any \(\mathbf{x}\in\mathbb{R}^n\) written as \(\mathbf{x}=\sum_i c_i\mathbf{v}_i\) is mapped to \(\mathbf{Ax}=\sum_i c_i\sigma_i\mathbf{u}_i\), so the entire, possibly rectangular or singular, action of \(\mathbf{A}\) reduces to independently rescaling each coordinate of one orthonormal basis to get the coordinates of another.

This correspondence between the four subspaces is often drawn as the Strang diagram (after Gilbert Strang, MIT): two axes for the row space and null space in the domain, two axes for the column space and left null space in the codomain, with the σi-scaled map \(\mathbf{v}_i\mapsto\mathbf{u}_i\) connecting the two pictures.

ℝⁿ ℝᵐ Row space Column space Null space Left null space x Ax 0 A A

Is the Eigendecomposition of a Symmetric Matrix Always Its SVD?

Sometimes. If \(\mathbf{A}\) is real symmetric, the spectral theorem already gives \(\mathbf{A}=\mathbf{Q}\mathbf{\Lambda}\mathbf{Q}^T\) with \(\mathbf{Q}\) orthogonal. When every eigenvalue is non-negative (\(\mathbf{A}\) is positive semidefinite), this is already a bona-fide SVD with \(\mathbf{U}=\mathbf{V}=\mathbf{Q}\) and \(\mathbf{\Sigma}=\mathbf{\Lambda}\), since the construction above, run on this \(\mathbf{A}\), reproduces exactly this \(\mathbf{Q}\) (because \(\mathbf{A}^T\mathbf{A}=\mathbf{A}^2=\mathbf{Q}\mathbf{\Lambda}^2\mathbf{Q}^T\) has the same eigenvectors). The symmetric example \(\mathbf{A}=\left(\begin{smallmatrix}4&1&1\\1&4&1\\1&1&4\end{smallmatrix}\right)\) from the eigenvalues chapter, with eigenvalues \(6,3,3\) and orthonormal eigenvectors \(\mathbf{q}_1,\mathbf{q}_2,\mathbf{q}_3\), is exactly such a case: its SVD and its EVD are the same decomposition, \(\mathbf{U}=\mathbf{V}=\mathbf{Q}\), \(\mathbf{\Sigma}=\mathrm{diag}(6,3,3)\).

If a symmetric matrix has a negative eigenvalue, however, the two decompositions differ. Singular values are never negative by definition, so a negative eigenvalue \(\lambda_i<0\) contributes the singular value \(\sigma_i=|\lambda_i|\), and the corresponding column of \(\mathbf{U}\) or \(\mathbf{V}\) must have its sign flipped to absorb the discrepancy. For \(\mathbf{B}=\mathrm{diag}(3,-2)\), the eigenvalues are \(3\) and \(-2\), but a valid SVD is \(\mathbf{\Sigma}=\mathrm{diag}(3,2)\) with \(\mathbf{V}=\mathbf{I}\) and \(\mathbf{U}=\mathrm{diag}(1,-1)\): indeed \(\mathbf{U}\mathbf{\Sigma}\mathbf{V}^T=\mathrm{diag}(3,-2)=\mathbf{B}\). So the eigendecomposition of a symmetric matrix is its SVD precisely when the matrix is positive semidefinite; otherwise the magnitudes agree but a sign has to move from \(\mathbf{\Lambda}\) into \(\mathbf{U}\) or \(\mathbf{V}\).

Singular Values versus Eigenvalues in General

For a general matrix, singular values and eigenvalues are not the same numbers at all — they answer different questions (how much does \(\mathbf{A}\) stretch orthogonal directions, versus which directions does \(\mathbf{A}\) leave invariant) and coincide only in the positive-semidefinite-symmetric case above. Two examples already met in this book underline how far apart they can get:

Uniqueness of the Singular Value Decomposition

The singular values \(\sigma_1\ge\cdots\ge\sigma_{\min(m,n)}\ge 0\) are always uniquely determined (they are the square roots of the eigenvalues of \(\mathbf{A}^T\mathbf{A}\), which are unique). The singular vectors, however, are not.

Let \(\mathbf{S}=\mathrm{diag}(s_1,\ldots,s_r)\) with each \(s_i\in\{-1,+1\}\) be any sign matrix, so that \(\mathbf{S}^2=\mathbf{I}\). If \(\mathbf{A}=\mathbf{U}_r\mathbf{\Sigma}_r\mathbf{V}_r^T\) is a compact SVD, then \(\mathbf{U}_r'=\mathbf{U}_r\mathbf{S}\) and \(\mathbf{V}_r'=\mathbf{V}_r\mathbf{S}\) give another valid SVD of the same matrix, since \(\mathbf{\Sigma}_r\) is diagonal and therefore commutes with any diagonal matrix:

\[\mathbf{U}_r'\mathbf{\Sigma}_r\mathbf{V}_r'^T =\mathbf{U}_r\mathbf{S}\mathbf{\Sigma}_r\mathbf{S}\mathbf{V}_r^T =\mathbf{U}_r(\mathbf{S}\mathbf{\Sigma}_r\mathbf{S})\mathbf{V}_r^T =\mathbf{U}_r\mathbf{\Sigma}_r\mathbf{V}_r^T=\mathbf{A}.\]

Every pair \((\mathbf{u}_i,\mathbf{v}_i)\) may therefore be flipped simultaneously without changing \(\mathbf{A}\) — exactly the freedom already used above to write the SVD of \(\mathrm{diag}(3,-2)\) with \(\mathbf{U}=\mathrm{diag}(1,-1)\) instead of \(\mathbf{U}=\mathbf{I}\), \(\mathbf{V}=\mathrm{diag}(1,-1)\).

A second, larger source of non-uniqueness appears whenever a singular value repeats. If \(\sigma_i=\sigma_j=\sigma\) for \(i\neq j\), then for any orthogonal \(2\times 2\) matrix \(\mathbf{R}\) rotating the corresponding pair of columns jointly in \(\mathbf{U}\) and \(\mathbf{V}\), the identity \(\mathbf{R}^T(\sigma\mathbf{I})\mathbf{R}=\sigma\mathbf{I}\) (since \(\sigma\mathbf{I}\) commutes with everything) shows the rotated pair reconstructs \(\mathbf{A}\) just as well. The repeated singular value \(\sigma_2=\sigma_3=3\) of \(\mathbf{A}=\left(\begin{smallmatrix}4&1&1\\1&4&1\\1&1&4 \end{smallmatrix}\right)\) is exactly such a case: any orthonormal basis of the plane spanned by \(\mathbf{q}_2,\mathbf{q}_3\), used simultaneously for the corresponding columns of \(\mathbf{U}\) and \(\mathbf{V}\), is an equally valid choice.

So the SVD is unique up to (a) independent sign flips of paired columns for distinct singular values, and (b) arbitrary simultaneous rotations of paired columns within a repeated singular value's subspace.

Robust Computation: Deriving \(\mathbf{V}\) from \(\mathbf{U}\)

The freedom just described causes a real numerical hazard. Eigendecomposing \(\mathbf{A}^T\mathbf{A}\) and \(\mathbf{AA}^T\) independently can hand back sign- or rotation-inconsistent \(\mathbf{U}\) and \(\mathbf{V}\) that individually diagonalise the right matrix, but whose combination \(\mathbf{U\Sigma V}^T\) no longer equals \(\mathbf{A}\): if the two eigensolvers happen to pick different sign matrices \(\mathbf{S}_j\neq\mathbf{S}_k\) for \(\mathbf{U}\) and \(\mathbf{V}\) respectively, the reconstruction becomes \(\mathbf{U}_j\mathbf{\Sigma}\mathbf{V}_k^T=\mathbf{U}\mathbf{S}_j\mathbf{\Sigma}\mathbf{S}_k\mathbf{V}^T\), which differs from \(\mathbf{A}\) whenever \(\mathbf{S}_j\neq\mathbf{S}_k\).

This is easy to underestimate until it is seen going wrong. Return to \(\mathbf{B}=\mathrm{diag}(3,-2)\) from above, whose correct SVD is \(\mathbf{U}=\mathrm{diag}(1,-1)\), \(\mathbf{\Sigma}=\mathrm{diag}(3,2)\), \(\mathbf{V}=\mathbf{I}\). Both \(\mathbf{B}^T\mathbf{B}=\mathrm{diag}(9,4)\) and \(\mathbf{BB}^T=\mathrm{diag}(9,4)\) happen to be the same diagonal matrix here, and any eigensolver is free to return \(\pm\mathbf{e}_1,\pm\mathbf{e}_2\) as eigenvectors for either one, independently. Suppose the solver for \(\mathbf{BB}^T\) happens to return \(\mathbf{U}=\mathbf{I}\) (a perfectly valid choice, since \(\mathbf{I}^T(\mathbf{BB}^T)\mathbf{I}\) is diagonal), while the solver for \(\mathbf{B}^T\mathbf{B}\), run independently, happens to flip a sign and return \(\mathbf{V}=\mathrm{diag}(-1,1)\) (equally valid on its own). Naively combining them with the shared \(\mathbf{\Sigma}=\mathrm{diag}(3,2)\) gives

\[\mathbf{U}\mathbf{\Sigma}\mathbf{V}^T=\mathbf{I}\cdot\mathrm{diag}(3,2)\cdot\mathrm{diag}(-1,1) =\mathrm{diag}(-3,2)\neq\mathbf{B}.\]

Both \(\mathbf{U}\) and \(\mathbf{V}\) were legitimate eigenvector matrices for their respective sides, yet combining them reconstructs the wrong matrix, purely because the two independent eigensolvers broke the sign ambiguity differently. The fix below removes this hazard entirely by never eigendecomposing both sides.

The fix is to compute only one side by eigendecomposition and derive the other directly from \(\mathbf{A}=\mathbf{U}_r\mathbf{\Sigma}_r\mathbf{V}_r^T\) itself. Multiplying both sides by \(\mathbf{U}_r^T\) and using \(\mathbf{U}_r^T\mathbf{U}_r=\mathbf{I}\),

\[\mathbf{U}_r^T\mathbf{A}=\mathbf{U}_r^T\mathbf{U}_r\mathbf{\Sigma}_r\mathbf{V}_r^T =\mathbf{\Sigma}_r\mathbf{V}_r^T \quad\Longrightarrow\quad \boxed{\mathbf{V}_r^T=\mathbf{\Sigma}_r^{-1}\mathbf{U}_r^T\mathbf{A}}\]

(valid since \(\mathbf{\Sigma}_r\) contains only strictly positive singular values and is therefore invertible). This is automatically consistent: if an eigensolver returns \(\mathbf{U}_j=\mathbf{U}_r\mathbf{S}_j\) for some sign matrix \(\mathbf{S}_j\), then computing \(\mathbf{V}_j^T\) from the boxed formula gives

\[\mathbf{V}_j^T=\mathbf{\Sigma}_r^{-1}\mathbf{U}_j^T\mathbf{A} =\mathbf{\Sigma}_r^{-1}\mathbf{S}_j^T\mathbf{U}_r^T\mathbf{A} =\mathbf{S}_j\mathbf{\Sigma}_r^{-1}\mathbf{U}_r^T\mathbf{A} =\mathbf{S}_j\mathbf{V}_r^T\]

(using that \(\mathbf{S}_j\) is diagonal and therefore commutes with the diagonal \(\mathbf{\Sigma}_r^{-1}\), and \(\mathbf{S}_j^T=\mathbf{S}_j\)). So \(\mathbf{V}_j=\mathbf{V}_r\mathbf{S}_j\) picks up exactly the same sign matrix as \(\mathbf{U}_j\), and \(\mathbf{U}_j\mathbf{\Sigma}_r\mathbf{V}_j^T=\mathbf{A}\) holds regardless of which sign convention the eigensolver used for \(\mathbf{U}_r\). Applying this formula to the cautionary example above, with the same arbitrarily-chosen \(\mathbf{U}=\mathbf{I}\), gives \(\mathbf{V}^T=\mathrm{diag}(\frac13,\frac12)\cdot\mathbf{I}\cdot\mathrm{diag}(3,-2)=\mathrm{diag}(1,-1)\) instead of the mismatched \(\mathrm{diag}(-1,1)\) — and indeed \(\mathbf{I}\cdot\mathrm{diag}(3,2)\cdot\mathrm{diag}(1,-1)=\mathrm{diag}(3,-2)=\mathbf{B}\), correctly, no matter which sign the eigensolver had picked for \(\mathbf{U}\). This is why practical SVD algorithms eigendecompose only one of \(\mathbf{A}^T\mathbf{A}\)/\(\mathbf{AA}^T\) (typically the smaller one) and derive the other side algebraically, rather than eigendecomposing both independently.

Truncated SVD and the Eckart–Young Theorem

Expanding the matrix product \(\mathbf{U\Sigma V}^T\) column by column shows that any matrix decomposes into a sum of rank-one pieces, ordered by shrinking singular value:

\[\mathbf{A}=\sum_{i=1}^{r}\sigma_i\,\mathbf{u}_i\mathbf{v}_i^T.\]

Keeping only the first \(k<r\) terms gives the rank-\(k\) matrix

\[\widetilde{\mathbf{A}}_k:=\sum_{i=1}^{k}\sigma_i\,\mathbf{u}_i\mathbf{v}_i^T.\]

The Eckart–Young theorem states that \(\widetilde{\mathbf{A}}_k\) is the best possible rank-\(k\) approximation of \(\mathbf{A}\) in the Frobenius norm: no other matrix of rank \(\le k\) is closer to \(\mathbf{A}\), and the resulting approximation error is exactly the singular values that were dropped,

\[\|\mathbf{A}-\widetilde{\mathbf{A}}_k\|_F^2=\sigma_{k+1}^2+\cdots+\sigma_r^2.\]

The symmetric example \(\mathbf{A}=\left(\begin{smallmatrix}4&1&1\\1&4&1\\1&1&4\end{smallmatrix}\right)\) makes this concrete. Its SVD equals its EVD (it is positive definite, as shown above), so the best rank-\(1\) approximation is \(\widetilde{\mathbf{A}}_1=6\,\mathbf{q}_1\mathbf{q}_1^T\), and the theorem predicts an error of \(\sqrt{\sigma_2^2+\sigma_3^2}=\sqrt{3^2+3^2}=3\sqrt2\) — which matches a direct computation of \(\|\mathbf{A}-\widetilde{\mathbf{A}}_1\|_F\). Truncation is the mathematical basis of low-rank compression: an \(m\times n\) matrix needs \(mn\) numbers to store, but its rank-\(k\) approximation needs only \((m+n+1)k\) numbers (\(k\) columns of \(\mathbf{U}\), \(k\) columns of \(\mathbf{V}\), and \(k\) singular values) — a large saving whenever \(k\) is small compared to \(m\) and \(n\), as is typical for images and other structured data.

Properties

Applications

Solving linear systems and least squares. For an overdetermined or ill-conditioned system \(\mathbf{Ax}\approx\mathbf{b}\), the vector \(\mathbf{x}=\mathbf{A}^{+}\mathbf{b}\) is the minimum-norm least squares solution: it minimizes \(\|\mathbf{Ax}-\mathbf{b}\|_2\) and, among all minimizers, has the smallest \(\|\mathbf{x}\|_2\). Substituting \(\mathbf{A}=\mathbf{U\Sigma V}^T\) confirms this satisfies the normal equations \(\mathbf{A}^T\mathbf{Ax}=\mathbf{A}^T\mathbf{b}\) from linear regression, but computing it through the SVD avoids ever forming \(\mathbf{A}^T\mathbf{A}\) explicitly — a numerically safer route whenever \(\mathbf{A}\) is ill-conditioned, since forming \(\mathbf{A}^T\mathbf{A}\) squares the condition number (\(\kappa(\mathbf{A}^T\mathbf{A})=\kappa(\mathbf{A})^2\)).

Written out via the pseudo-inverse's own definition, this minimum-norm solution is the explicit sum

\[\mathbf{x}=\mathbf{A}^{+}\mathbf{b}=\sum_{i=1}^r\frac{\mathbf{u}_i^T\mathbf{b}}{\sigma_i}\,\mathbf{v}_i,\]

using only the nonzero singular values. The system \(\mathbf{Ax}=\mathbf{b}\) has an exact solution precisely when \(\mathbf{b}\) lies in the column space, i.e. when it is orthogonal to every left null-space vector \(\mathbf{u}_{r+1},\ldots,\mathbf{u}_m\); in that case every solution has the form \(\mathbf{x}=\mathbf{A}^{+}\mathbf{b}+\mathbf{n}\) for \(\mathbf{n}\) in the null space, an affine subspace of dimension \(n-r\). Among all of these, \(\mathbf{A}^{+}\mathbf{b}\) is the shortest: it is built entirely out of \(\mathbf{v}_1,\ldots,\mathbf{v}_r\) (the row space), which is orthogonal to the null space \(\mathbf{n}\) lives in, so by the Pythagorean theorem \(\|\mathbf{A}^{+}\mathbf{b}+\mathbf{n}\|_2^2=\|\mathbf{A}^{+}\mathbf{b}\|_2^2+\|\mathbf{n}\|_2^2\) is minimized exactly at \(\mathbf{n}=\mathbf{0}\), which is exactly why \(\mathbf{A}^{+}\mathbf{b}\) rather than any other point of the solution set was singled out above.

Data compression. The truncated SVD from the previous section is used directly to compress images and other tabular data: replacing an \(m\times n\) matrix by its rank-\(k\) approximation keeps the entries that matter most (largest singular values first) while shrinking the storage from \(mn\) to \((m+n+1)k\) numbers.

Principal component analysis. Given a mean-centered data matrix \(\mathbf{X}\) (one row per observation), the sample covariance matrix is \(\mathbf{\Sigma}_{\text{cov}}=\frac{1}{n-1}\mathbf{X}^T\mathbf{X}\), and its eigenvectors are exactly the right singular vectors of \(\mathbf{X}\), with eigenvalues \(\sigma_i^2/(n-1)\). Principal component analysis therefore reduces to computing the SVD of \(\mathbf{X}\) directly, which is both faster and numerically more stable than first forming \(\mathbf{X}^T\mathbf{X}\) and then eigendecomposing it — the same squared-condition-number issue as above. The resulting principal axes are exactly the ellipse axes drawn when visualizing a covariance matrix as a confidence ellipse; the full derivation from variance maximization, the equivalent reconstruction-error view, and how to read explained-variance percentages are worked out in principal component analysis.

Homogeneous least squares. Ordinary least squares assumes \(\mathbf{b}\neq\mathbf{0}\), but many fitting problems instead ask for a nontrivial \(\mathbf{x}\) making \(\mathbf{Ax}\) as close to \(\mathbf{0}\) as possible, since \(\mathbf{x}=\mathbf{0}\) always solves \(\mathbf{Ax}=\mathbf{0}\) exactly but carries no information. Forbidding that trivial answer with the constraint \(\|\mathbf{x}\|_2=1\) turns the problem into minimizing \(\|\mathbf{Ax}\|_2^2=\mathbf{x}^T(\mathbf{A}^T\mathbf{A})\mathbf{x}\) over the unit sphere, which is exactly the Rayleigh quotient of \(\mathbf{A}^T\mathbf{A}\) again. Its minimum is the smallest eigenvalue of \(\mathbf{A}^T\mathbf{A}\), attained at the corresponding eigenvector — so the best-fitting unit vector is simply \(\mathbf{v}_n\), the right singular vector belonging to the smallest singular value of \(\mathbf{A}\), and the achieved minimum is exactly \(\sigma_n\).

For \(\mathbf{A}=\left(\begin{smallmatrix}1&1\\1&-1\\1&0\end{smallmatrix}\right)\), the columns are already orthogonal, so \(\mathbf{A}^T\mathbf{A}=\left(\begin{smallmatrix}3&0\\0&2\end{smallmatrix}\right)\) is diagonal with no computation needed: its eigenvalues are simply \(3\) and \(2\), read off the diagonal, giving singular values \(\sigma_1=\sqrt3\) and \(\sigma_2=\sqrt2\) with \(\mathbf{v}_1=\mathbf{e}_1\), \(\mathbf{v}_2=\mathbf{e}_2\). The unit vector minimizing \(\|\mathbf{Ax}\|_2\) is therefore \(\mathbf{x}=\mathbf{v}_2=(0,1)\), giving \(\mathbf{Ax}=(1,-1,0)\) with norm \(\sqrt2=\sigma_2\) — and no other unit vector does better, since every other direction has a component along \(\mathbf{v}_1\), which is stretched by the larger factor \(\sqrt3\) instead. This is the same technique used, for instance, to fit a plane through noisy 3D points or a line through the origin in total least squares: collect the constraints as rows of \(\mathbf{A}\) and read off \(\mathbf{v}_n\).