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

Real-world observations are often inconsistent. You've recorded thousands of requests with their load and latency, hoping to explain the results using just a few features. The number of equations far exceeds the number of unknowns, and it's nearly impossible for all equations to hold exactly. The value of linear algebra lies here: transforming the idea of "impossible to be perfectly accurate" into a measurable problem of best approximation.

Learning Objectives

• Understand least squares from a projection perspective;
• Explain stable 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\min_{\mathbf{x}} \lVert A\mathbf{x} - \mathbf{b} \rVert_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$.

While the normal equations are elegant for derivation, they are not always suitable for direct computation. Forming $A^\mathsf{T}A$ squares the condition number, potentially amplifying numerical instability. In practice, QR decomposition or SVD is often preferred; NumPy’s lstsq automatically selects the most appropriate numerical method.

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(slope, intercept, rank)

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.

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.

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 \Sigma V^\mathsf{T}. $$

• The columns of $V$ represent orthogonal directions in the input space;
• The non-negative diagonal entries of $\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 $\Sigma$, and finally rotate or reflect the result into the output space via $U$.

The number of non-zero singular values equals the rank of the matrix. The ratio between the largest and smallest non-zero singular values relates to the condition number: the more extreme this ratio, the more certain directions amplify input errors into solution errors.

4. Low-Rank Approximation and Compression

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

$$ A_k = U_k \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 provides the strictly optimal low-rank approximation. Intuitively, larger singular values correspond to stronger patterns of variation in the matrix, while smaller singular values correspond to weaker directions.

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}\boldsymbol{\mu}^\mathsf{T}. $$

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

$$ X_c = U\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.

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

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.

6. Which Tool to Choose

Question	Common Tools	Key Reminders
Linear systems in matrix form	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	Leverage 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

Linear equations strive for exact matches, while least squares handles observations that cannot be precisely matched. Eigenvalue decomposition reveals special directions in square matrices, and SVD provides input directions, scaling factors, and output directions for any matrix. PCA is an application of this structure in data dimensionality reduction, it is not a magic spell that automatically interprets data meaning.

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.