Singular Value Decomposition for Ill-Conditioned Data

Why Use the Singular Value Decomposition?

Multivariate data often contain directions that are redundant, nearly redundant, or measured on very different scales. These features can make a calculation appear to have a unique answer when the data do not support one, or can make small numerical errors produce large changes in the answer. The singular value decomposition (SVD) gives a direct way to inspect the directions in a matrix, quantify their strength, and solve least-squares problems even when the matrix is rank deficient.

Let \(A\) be a real \(m\times p\) matrix. Its singular values describe how strongly \(A\) acts along orthogonal directions. Large singular values correspond to directions with substantial signal in the matrix; zero singular values correspond to directions that \(A\) eliminates. Very small singular values indicate directions that are difficult to distinguish reliably from numerical error or measurement noise.

The Decomposition and Its Meaning

Definition: A singular value decomposition of \(A\in\mathbb{R}^{m\times p}\) is a factorization \(A=U\Sigma V^{T}\), where \(U\in\mathbb{R}^{m\times m}\) and \(V\in\mathbb{R}^{p\times p}\) are orthogonal matrices, and \(\Sigma\) is an \(m\times p\) diagonal-shaped matrix. Its diagonal entries are the singular values \(\sigma_1\geq\sigma_2\geq\cdots\geq\sigma_q\geq0\), where \(q=\min(m,p)\). The columns of \(U\) are left singular vectors, and the columns of \(V\) are right singular vectors. The number of positive singular values is the rank of \(A\).

The right singular vectors describe directions in the input space, while the left singular vectors describe the corresponding directions in the output space. In particular, \(Av_j=\sigma_j u_j\) for each positive singular value \(\sigma_j\). Thus, \(A\) maps the unit vector \(v_j\) to a vector of length \(\sigma_j\) in direction \(u_j\). A zero singular value means \(Av_j=0\), so \(v_j\) lies in the null space of \(A\).

The SVD exists for every real rectangular matrix. One way to see why is to apply the spectral theorem to the symmetric positive semidefinite matrix \(A^T A\). Its eigenvalues are nonnegative; their square roots are the singular values. For each positive eigenvalue \(\sigma_j^2\) with unit eigenvector \(v_j\), define \(u_j=Av_j/\sigma_j\). These \(u_j\) are orthonormal, since \(u_i^T u_j=v_i^T A^T A v_j/(\sigma_i\sigma_j)=0\) for \(i\ne j\), and equals 1 for \(i=j\). Completing these vectors to orthonormal bases gives \(U\) and \(V\), and the resulting factorization is \(A=U\Sigma V^T\). The zero-eigenvalue directions satisfy \(Av_j=0\). This also shows that the positive singular values count the rank of \(A\). \(\square\)

This construction explains why singular values are useful even when the matrix is rectangular: they are defined through the action of \(A\), not by requiring \(A\) itself to be square. It also gives an important computational principle. Although \(A^T A\) can explain the theory, forming it in a numerical calculation squares the ratio between the largest and smallest nonzero singular values. For ill-conditioned data, computing from \(A\) directly by an SVD avoids that loss of information.

Rank-Deficient Least Squares

For a response vector \(y\in\mathbb{R}^{m}\), least squares seeks an \(x\in\mathbb{R}^{p}\) minimizing \(\|Ax-y\|_2\). If \(A\) has linearly dependent columns, the minimizer may not be unique. The SVD identifies both the fitted component and the ambiguity.

Theorem: Suppose \(A=U\Sigma V^T\) has rank \(r\), with positive singular values \(\sigma_1,\ldots,\sigma_r\). The minimum-Euclidean-norm least-squares solution is \(x_{\min}=\sum_{j=1}^{r}(u_j^T y/\sigma_j)v_j\). Every least-squares solution has this form plus an arbitrary vector in the null space of \(A\).

To prove the result, write any \(x\) in the orthonormal basis of right singular vectors as \(x=\sum_{j=1}^{p}c_jv_j\). In the left singular-vector basis, the residual has components \(\sigma_jc_j-u_j^Ty\) for the \(r\) positive singular values, and components \(-u_j^Ty\) in the remaining left directions. The latter components cannot be changed by choosing \(x\). Therefore the residual norm is minimized by setting \(c_j=u_j^Ty/\sigma_j\) for \(1\leq j\leq r\). For \(j>r\), the corresponding right singular vectors lie in the null space, so their coefficients do not affect the residual. Setting all those coefficients to zero gives the solution with the smallest norm, because the \(v_j\) are orthonormal. Allowing arbitrary coefficients in those null directions gives every minimizer. \(\square\)

Worked Example: A Minimum-Norm Fit With Redundant Columns

Consider \(A=\begin{pmatrix}1&1\\2&2\end{pmatrix}\) and \(y=(1,2)^T\). The two columns of \(A\) are identical, so it is rank deficient. Write \(A=ab^T\), where \(a=(1,2)^T\) and \(b=(1,1)^T\). Their lengths are \(\sqrt{5}\) and \(\sqrt{2}\), giving the sole positive singular value \(\sigma_1=\sqrt{10}\), with \(u_1=(1,2)^T/\sqrt{5}\) and \(v_1=(1,1)^T/\sqrt{2}\).

The theorem gives \(u_1^Ty=5/\sqrt{5}=\sqrt{5}\), and hence \(x_{\min}=(\sqrt{5}/\sqrt{10})v_1=(1/2,1/2)^T\). Indeed, \(Ax_{\min}=(1,2)^T=y\), so the residual is zero. Every vector with \(x_1+x_2=1\) also fits exactly, because \(Ax=(x_1+x_2)(1,2)^T\). Among these vectors, \((1/2,1/2)^T\) has the smallest squared norm: for \(x_2=1-x_1\), \(x_1^2+x_2^2=2(x_1-1/2)^2+1/2\), which is minimized at \(x_1=1/2\).

Principal Components From the Data Matrix

Let \(Z\) be an \(n\times p\) matrix of observations after each variable has been centered by subtracting its sample mean. Its SVD, \(Z=U\Sigma V^T\), gives principal components directly. The right singular vectors are directions in variable space; the scores along those directions are the columns of \(ZV\). Since \(Zv_j=\sigma_j u_j\), the \(j\)th score vector is \(\sigma_j u_j\). With the usual sample covariance matrix \(Z^TZ/(n-1)\), the variance of that component is \(\sigma_j^2/(n-1)\).

This relationship follows from the squared length of the score vector. Because \(u_j\) is a unit vector, the sum of squared scores is \((\sigma_j u_j)^T(\sigma_j u_j)=\sigma_j^2\). Centering makes the sample mean score zero, so dividing that sum by \(n-1\) gives the sample variance. The SVD therefore gives the component directions, their scores, and their variances without first forming a covariance matrix.

Worked Example: Principal Components of Two Redundant Variables

Suppose the centered data matrix is \(Z=\begin{pmatrix}-1&-1\\0&0\\1&1\end{pmatrix}\). The two columns are identical. Set \(v_1=(1,1)^T/\sqrt{2}\). Then \(Zv_1=(-\sqrt{2},0,\sqrt{2})^T\), whose length is 2, so \(\sigma_1=2\) and \(u_1=(-1,0,1)^T/\sqrt{2}\). The orthogonal direction \(v_2=(1,-1)^T/\sqrt{2}\) is in the null space, so its singular value is zero.

The first principal component scores are \((-\sqrt{2},0,\sqrt{2})^T\). There are \(n=3\) observations, so their sample variance is \(\sigma_1^2/(n-1)=4/2=2\). The second component has zero scores and zero variance. As a check, \(Z^TZ=\begin{pmatrix}2&2\\2&2\end{pmatrix}\); its sample covariance matrix is \(Z^TZ/2\), whose trace is 2, equal to the total component variance.

Singular-Value Decay, Conditioning, and Numerical Rank

For a full-column-rank matrix, the 2-norm condition number is \(\kappa_2(A)=\sigma_1/\sigma_p\). A large value means that some input directions are compressed far more than others. Consequently, a small perturbation in the response or in the matrix can produce a much larger change in a least-squares solution. If \(A\) is exactly rank deficient, its usual condition number for solving a unique full-column-rank problem is infinite: there is no unique solution in the null directions.

The pattern of singular values conveys more than a single condition number. A sharp drop followed by values near zero signals that the data matrix has directions carrying very little variation compared with its leading directions. The exact rank is the number of strictly positive singular values in exact arithmetic. In finite-precision computation, however, a singular value may be positive but too small to distinguish reliably from rounding error. Numerical rank is therefore an operational assessment of how many singular directions are treated as meaningful at the scale and precision of the calculation. It is not an intrinsic integer independent of context: measurement noise, variable scaling, and computational precision all matter.

Worked Example: A Small Singular Value Amplifies Error

Take \(A=\begin{pmatrix}1&0\\0&10^{-8}\end{pmatrix}\). Its singular values are 1 and \(10^{-8}\), so its condition number is \(10^8\). For \(y=(1,10^{-8})^T\), the exact solution is \(x=(1,1)^T\). Now perturb only the second response entry by \(10^{-8}\), giving \(\tilde y=(1,2\cdot10^{-8})^T\). Solving \(A\tilde x=\tilde y\) gives \(\tilde x=(1,2)^T\).

The response perturbation has Euclidean norm \(10^{-8}\), while the solution changes by norm 1. The small singular direction divides the second response component by \(10^{-8}\), magnifying the perturbation. The singular values also show why blindly treating every nonzero direction as equally reliable can be misleading: whether a direction of size \(10^{-8}\) is practically informative depends on the precision and noise scale in the application.

Common Mistake and Practical Interpretation

A common mistake is to say, “The matrix is almost singular, so the least-squares solution is nearly unique and reliable.” This confuses existence with stability. A very small but nonzero singular value may technically yield a unique full-rank solution, yet the division by that value can make the solution highly sensitive. Conversely, an exactly rank-deficient matrix has many least-squares solutions, but the SVD still identifies a well-defined minimum-norm one. The correction is to inspect the singular values and the intended interpretation of the coefficients, rather than relying only on whether a numerical routine returned a vector.

In principal component analysis, a small singular value means little observed variation along its right singular vector. It does not by itself prove that the direction is scientifically irrelevant: the scale of measurement error and the variables’ units matter. Centering is essential when the goal is covariance-based PCA; otherwise the leading direction may mainly describe the overall mean rather than variation around it.

Key takeaway: The SVD separates a matrix into orthogonal directions and their scales. It reveals exact rank deficiency, provides the minimum-norm least-squares fit, and computes principal components directly. Singular-value decay and condition numbers warn when a formally available solution may be numerically sensitive.

Check Your Understanding

Answer each question by referring to the directions and scales in the SVD.

What does a zero singular value imply about the corresponding right singular vector?
Why can a rank-deficient least-squares problem have multiple minimizers but still have a unique minimum-norm minimizer?
For a centered data matrix, how is the sample variance of the \(j\)th principal component related to its singular value?
What does a large condition number indicate about sensitivity to small perturbations?
Why can numerical rank depend on measurement noise and computational precision?