10.3 Least Squares, Eigen Decomposition, and SVD: From Exact Solutions to Best Approximations

Real observations rarely satisfy one exact linear system. Thousands of request-load and latency records may need to be approximated with only dozens of features, so the equations greatly outnumber the unknowns. Linear algebra turns “they cannot all be exact” into a measurable best-approximation problem.

Learning Objectives

• Understand least squares from a projection perspective;
• Explain invariant directions represented by eigenvalues and eigenvectors;
• Grasp that singular value decomposition applies to any matrix;
• Clarify the centering, dimensionality reduction, and information loss boundaries in PCA.

1. Least Squares: Minimizing Residuals

When $A\mathbf{x}=\mathbf{b}$ has no exact solution, least squares selects:

$$\hat{\mathbf{x}}=\arg \underset{\mathbf{x}}{min}\Vert A\mathbf{x}-\mathbf{b}{\Vert }_2^2.$$

The residual $\mathbf{r}=\mathbf{b}-A\hat{\mathbf{x}}$ is orthogonal to each column of $A$, so:

$$A^{\mathsf{T}}(\mathbf{b}-A\hat{\mathbf{x}})=0,$$

leading to the normal equations:

$$A^{\mathsf{T}}A\hat{\mathbf{x}}=A^{\mathsf{T}}\mathbf{b}.$$

Geometrically, $A\hat{\mathbf{x}}$ is the orthogonal projection of $\mathbf{b}$ onto the column space of $A$.

The fitted projection $A\hat{\mathbf{x}}$ is unique, but the coefficient vector is unique only when $A$ has full column rank. A nonzero null-space vector leaves predictions unchanged: $A(\hat{\mathbf{x}}+\mathbf{z})=A\hat{\mathbf{x}}$. Among multiple minimizers, lstsq returns one with minimum Euclidean norm. Least squares also applies when an exact solution exists.

For full-column-rank $A$, ${\kappa }_2(A^{\mathsf{T}}A)={\kappa }_2(A{)}^2$. The normal equations are useful for derivation, but explicitly forming them can amplify rounding problems. Common computational methods include QR and SVD. The example calls NumPy 2.x lstsq; it does not assume automatic selection of a method best suited to the application.

import numpy as np

# y ≈ slope * x + intercept
x = np.array([0.0, 1.0, 2.0, 3.0])
y = np.array([1.1, 2.9, 5.2, 6.8])
A = np.column_stack([x, np.ones_like(x)])

solution, residuals, rank, singular_values = np.linalg.lstsq(
    A, y, rcond=None
)
slope, intercept = solution
print(f"{slope:.2f} {intercept:.2f} {rank}")  # 1.94 1.09 2
residual = y - A @ solution
print(f"{residual @ residual:.3f}")            # 0.082
print(np.allclose(A.T @ residual, 0.0))       # True

Least squares implicitly assumes an error metric: large residuals are penalized more heavily due to squaring, making the method sensitive to outliers. When data violate assumptions such as independence or homoscedasticity, parameter estimates may still be computed, but statistical interpretation requires additional validation.

residuals contains sums of squared residuals, not per-sample errors. It is empty when rank is below the column count or rows do not outnumber columns; empty does not mean zero fitting error. Compute y - A @ solution directly when needed. Explicit rcond=None uses machine precision times the larger matrix dimension as the relative rank threshold. NumPy’s documentation records the default change in 2.0.

2. Eigenvectors: Directions That Remain Unchanged Under Transformation

For a square matrix $A$, if there exists a non-zero vector $\mathbf{v}$ such that

$$A\mathbf{v}=\lambda \mathbf{v},$$

then $\mathbf{v}$ is an eigenvector and $\lambda$ is the corresponding eigenvalue. Along this direction, the transformation acts only by scaling; negative eigenvalues also reverse the direction of the vector.

Eigenvalues can help us understand:

• The growth or decay of a linear dynamical system along different directions;
• The structural properties encoded in adjacency and Laplacian matrices of graphs;
• The principal directions of variation in symmetric covariance matrices.

Not all matrices have a complete set of eigenvectors in the real number domain. Matrices may possess complex eigenvalues, or may lack sufficient linearly independent eigenvectors to be diagonalized. Therefore, "computing eigen decomposition" is not a universally applicable or unconditional step for arbitrary matrices.

A 90° planar rotation has no real eigenvectors, while [[1,1],[0,1]] has only a one-dimensional eigenspace and is not diagonalizable. A real symmetric matrix does admit an orthonormal eigenbasis, for which libraries provide eigh. A repeated eigenvalue can have more than one choice of basis.

3. SVD: Directions and Magnitudes of Any Matrix

Any real matrix $A\in {\mathbb{R}}^{m\times n}$ admits a singular value decomposition:

$$A=U\mathrm{\Sigma }{V}^{\mathsf{T}}.$$

• The columns of $V$ represent orthogonal directions in the input space;
• The non-negative diagonal entries of $\mathrm{\Sigma }$ are the singular values, indicating the scaling magnitude along each direction;
• The columns of $U$ represent orthogonal directions in the output space.

The transformation can be viewed as a three-step process: first rotate or reflect the input via $V^{\mathsf{T}}$, then scale along coordinate axes via $\mathrm{\Sigma }$, and finally rotate or reflect the result into the output space via $U$.

In the full SVD, $U$ is $m\times m$, $\mathrm{\Sigma }$ is $m\times n$, and $V$ is $n\times n$; the middle step also removes or pads coordinates. The reduced form uses $r=min(m,n)$ with factor shapes $(m,r)$, $(r,r)$, and $(r,n)$. NumPy returns a one-dimensional singular-value array and Vh, which equals V.T for real matrices. Multiplying as though that result were V reverses the intended orientation.

In exact arithmetic, the number of nonzero singular values equals the rank. For an invertible square matrix, ${\kappa }_2(A)={\sigma }_{max}/{\sigma }_{min}$. A singular square matrix is not invertible and conventionally has infinite condition number; discarding zero singular values does not make that full problem well conditioned. A ratio restricted to nonzero directions describes a different, restricted problem.

4. Low-Rank Approximation and Compression

Arrange the singular values in descending order and retain only the first $k$:

$$A_k=U_k{\mathrm{\Sigma }}_k{V}_k^{\mathsf{T}}.$$

$A_k$ is a low-rank approximation with rank at most $k$. For both the spectral norm and the Frobenius norm, it attains the optimal approximation error, though the minimizing matrix need not be unique. Intuitively, larger singular values correspond to stronger patterns of variation in the matrix, while smaller singular values correspond to weaker directions.

For $0\le k<min(m,n)$, the spectral-norm error is ${\sigma }_{k+1}$ and the squared Frobenius error is $\sum _{i>k}{\sigma }_i^2$. These guarantees concern matrix approximation without additional constraints. Requiring nonnegative entries, sparsity, or business rules changes the optimization problem.

This technique can be applied to image compression, latent semantic analysis, and recommendation models. However, "small" does not automatically imply "noise." Rare patterns or outliers may also lie along low-energy directions. Before compression, it is essential to examine the task context to determine whether such low-energy components are meaningful or should be discarded.

5. The Relationship Between PCA and SVD

Principal Component Analysis (PCA) identifies the orthogonal directions that maximize data variance. When the data matrix $X$ is structured such that each row represents a sample, the first step is to center the data along the columns:

$$X_c=X-\mathbf{1}{\mathit{\mu }}^{\mathsf{T}}.$$

The covariance matrix is proportional to $X_c^{\mathsf{T}}{X}_c$. Performing SVD on the centered data:

$$X_c=U\mathrm{\Sigma }{V}^{\mathsf{T}},$$

the first few columns of $V$ correspond to the principal component directions, and the squared singular values indicate the variance explained along each direction.

For $N>1$ samples, the sample covariance is $X_c^{\mathsf{T}}{X}_c/(N-1)$ and each component variance is ${\sigma }_i^2/(N-1)$. Its explained-variance ratio is ${\sigma }_i^2/\sum _j{\sigma }_j^2$. If all samples are identical, the denominator is zero and there is no nonzero variance to allocate.

import numpy as np

X = np.array([
    [2.0, 1.0],
    [3.0, 2.0],
    [4.0, 2.5],
    [5.0, 4.0],
])

mean = X.mean(axis=0)
centered = X - mean
_, singular_values, vt = np.linalg.svd(centered, full_matrices=False)

principal_direction = vt[0]
one_dimensional = centered @ principal_direction
reconstructed = np.outer(one_dimensional, principal_direction) + mean
ratio = singular_values[0] ** 2 / np.sum(singular_values ** 2)
error = np.linalg.norm(X - reconstructed, "fro")
print(f"{ratio:.4%} {error:.6f}")  # 99.0588% 0.301962
print(np.isclose(error, singular_values[1]))  # True

If the features have vastly different scales, whether to standardize beforehand depends on the specific problem: PCA is sensitive to scale, and standardization alters the meaning of "maximum variance." The mean, scale, and principal components derived from the training data must be fixed and applied consistently to new data.

Here one retained direction leaves a Frobenius reconstruction error equal to the discarded singular value. Changing a singular vector’s overall sign does not change reconstruction; equal singular values also permit rotations within their subspace. Compare reconstructions and subspaces across implementations, rather than requiring identical vector signs. See NumPy’s SVD shape conventions.

6. Which Tool to Choose

Question	Common Tools	Key Reminders
Full-rank square linear systems	LU, specialized solve	Avoid explicit inversion
Overdetermined systems and fitting	QR, SVD, lstsq	Check residuals, rank, and data assumptions
Principal directions of symmetric matrices	Eigenvalue decomposition	Use the symmetric structure
Matrix rank and low-rank approximation	SVD	Pay attention to condition number and truncation error
Data dimensionality reduction	PCA/SVD after centering	Scaling strategy belongs to model definition

Common Misconceptions

• PCA automatically identifies the most important business features: It maximizes variance, but has no understanding of business objectives.
• Retaining 95% of variance preserves 95% of information: This is a heuristic tied to task relevance, not a universal equivalence.
• The direction with the largest eigenvalue is always the most stable: Stability also depends on matrix properties, spectral gaps, and perturbation sensitivity.
• SVD is just a theoretical decomposition: It is, in fact, a foundational numerical tool underpinning least squares, pseudoinverses, dimensionality reduction, and condition analysis.

Exercise

1. Explain why the least-squares residuals are orthogonal to the column space.
2. Compute the eigenvalues and eigenvectors of the diagonal matrix diag(4, 1), and describe the geometric transformation it represents.
3. Use SVD to approximate a grayscale matrix as a rank-1 matrix, and compare the reconstruction error.
4. Perform PCA on the same dataset first with only centering, and then with standardization followed by centering. Explain why the principal directions differ between the two approaches.

Summary

The next chapter moves from static transformations to continuous change: how functions behave in a small neighborhood around a point, and how many such local changes accumulate to form the overall behavior.