Linear Algebra and Numerical Computing Questions

The idea in plain words

A matrix $A$ is a machine that takes a vector and returns another vector: it can stretch, squash, rotate and shear. Most input vectors come out pointing in a new direction. An eigenvector is a special nonzero direction that comes out pointing along the same line (or exactly the opposite way), only longer or shorter. The eigenvalue is the factor by which it is scaled:

$A v = \lambda v$

Here $v$ is the eigenvector and $\lambda$ is its eigenvalue. A negative $\lambda$ flips the vector; a $\lambda$ of 0 collapses it to the zero vector. The geometric picture: feed every point of a unit circle through a symmetric matrix (one equal to its own transpose, so entry row $i$, column $j$ equals entry row $j$, column $i$) and you get an ellipse whose axes lie along the eigenvectors, with semi-axis lengths equal to the absolute eigenvalues.

Worked example

Take the symmetric matrix with rows (3, 1) and (1, 2).

By hand first. $Av = \lambda v$ can be rewritten $(A - \lambda I)v = 0$, where $I$ is the identity matrix. A nonzero $v$ can satisfy that only if $A - \lambda I$ squashes some direction to nothing, which is the case exactly when its determinant (a single number that is zero for such a matrix; for a 2 by 2 matrix with rows $(a, b)$ and $(c, d)$ it is $ad - bc$) is zero:

$\det\begin{pmatrix} 3-\lambda & 1 \\ 1 & 2-\lambda \end{pmatrix} = (3-\lambda)(2-\lambda) - 1 = \lambda^2 - 5\lambda + 5 = 0$

The 5 in the middle is the trace (sum of the diagonal, 3 + 2) and the last 5 is the determinant of $A$ itself (3 times 2 minus 1 times 1). The quadratic formula gives $\lambda = (5 \pm \sqrt{5})/2$, which is 1.382 and 3.618. For the eigenvector of each $\lambda$, the first row of $(A - \lambda I)v = 0$ reads $(3 - \lambda)v_1 + v_2 = 0$, so $v = (1, \lambda - 3)$, and dividing by its length gives a vector of length 1.

The program repeats that arithmetic, computes the eigenpairs with a library routine, checks $Av = \lambda v$ numerically, checks that the eigenvalues sum to the trace and multiply to the determinant, and then shows the ML uses. The sign of an eigenvector is arbitrary, so the code fixes it (largest component positive); the hand-built direction for 1.382 comes out with the opposite sign, which is the same line.

import numpy as np

np.set_printoptions(precision=4, suppress=True)

# 1. Geometry: a symmetric 2x2 matrix, its eigenpairs, and what it does to vectors.
A = np.array([[3.0, 1.0], [1.0, 2.0]])
vals, vecs = np.linalg.eigh(A)           # eigenvalues in ascending order
for i in range(2):                       # fix the arbitrary sign: largest component positive
    if vecs[np.argmax(np.abs(vecs[:, i])), i] < 0:
        vecs[:, i] = -vecs[:, i]
print("eigenvalues:", vals)
for i in range(2):
    v = vecs[:, i]
    print(f"eigenvector {i}:", v, " A@v:", A @ v, " lambda*v:", vals[i] * v)
print("trace:", np.trace(A), " sum of eigenvalues:", vals.sum())
print("det:", round(np.linalg.det(A), 4), " product of eigenvalues:", round(vals.prod(), 4))
x = np.array([1.0, 0.0])
print("generic x:", x, " A@x:", A @ x, " (direction changed)")

# By hand: det(A - lambda*I) = lambda^2 - trace*lambda + det = 0, solved with the quadratic formula.
tr, dt = np.trace(A), np.linalg.det(A)
roots = (tr + np.array([-1.0, 1.0]) * np.sqrt(tr**2 - 4 * dt)) / 2
print("roots of lambda^2 - 5*lambda + 5:", roots)
for lam in roots:
    # first row of (A - lam*I) v = 0 is (3 - lam)*v1 + v2 = 0, so v = (1, lam - 3); scale to length 1
    v = np.array([1.0, lam - 3.0])
    print(f"lambda={lam:.4f}: direction (1, {lam - 3:.4f}) scaled to length 1 ->", v / np.linalg.norm(v))

# 2. ML use 1: PCA directions are eigenvectors of the covariance matrix.
rng = np.random.default_rng(0)
Z = rng.normal(size=(500, 2)) * np.array([3.0, 0.5])
c, s = np.cos(0.6), np.sin(0.6)
R = np.array([[c, -s], [s, c]])
X = Z @ R.T                              # stretched cloud, rotated by 0.6 radians
C = np.cov(X, rowvar=False)              # 2x2 sample covariance (divisor n-1)
cv, cV = np.linalg.eigh(C)
print("covariance eigenvalues:", cv, " variance share of top direction:", round(cv[1] / cv.sum(), 4))
top = cV[:, 1] * np.sign(cV[0, 1])
print("top eigenvector:", top, " true long axis:", np.array([c, s]))

# 3. ML use 2/3: Hessian eigenvalues set gradient-descent speed (f = 0.5 x^T H x, step 1/lambda_max).
def steps_to_converge(H, tol=1e-6, lr=None):
    x = np.array([1.0, 1.0])
    if lr is None:
        lr = 1.0 / np.linalg.eigvalsh(H).max()
    k = 0
    while np.linalg.norm(x) > tol:
        x = x - lr * (H @ x)
        k += 1
    return k

for name, H in [("eigenvalues 1 and 1", np.diag([1.0, 1.0])),
                ("eigenvalues 1 and 10", np.diag([10.0, 1.0])),
                ("eigenvalues 1 and 100", np.diag([100.0, 1.0]))]:
    print(f"{name}: condition number {np.linalg.cond(H):.0f}, steps to reach 1e-6: {steps_to_converge(H)}")

# Why those step counts: each coordinate is multiplied by (1 - lambda_i/lambda_max) per step.
for kappa in (10, 100):
    shrink = 1 - 1 / kappa
    print(f"condition number {kappa}: flat coordinate shrinks by {shrink} per step; ceil(ln(1e-6)/ln({shrink})) = {int(np.ceil(np.log(1e-6) / np.log(shrink)))}")

# A larger step helps only a little: compare 1/lambda_max with the best fixed step 2/(lambda_max + lambda_min).
for kappa in (10, 100):
    H = np.diag([float(kappa), 1.0])
    best = 2.0 / (kappa + 1.0)
    print(f"condition number {kappa}: step 1/lambda_max takes {steps_to_converge(H)}, step 2/(lambda_max+lambda_min) takes {steps_to_converge(H, lr=best)}")

Output:

eigenvalues: [1.382 3.618]
eigenvector 0: [-0.5257  0.8507]  A@v: [-0.7265  1.1756]  lambda*v: [-0.7265  1.1756]
eigenvector 1: [0.8507 0.5257]  A@v: [3.0777 1.9021]  lambda*v: [3.0777 1.9021]
trace: 5.0  sum of eigenvalues: 5.0
det: 5.0  product of eigenvalues: 5.0
generic x: [1. 0.]  A@x: [3. 1.]  (direction changed)
roots of lambda^2 - 5*lambda + 5: [1.382 3.618]
lambda=1.3820: direction (1, -1.6180) scaled to length 1 -> [ 0.5257 -0.8507]
lambda=3.6180: direction (1, 0.6180) scaled to length 1 -> [0.8507 0.5257]
covariance eigenvalues: [0.2344 8.7383]  variance share of top direction: 0.9739
top eigenvector: [0.8247 0.5656]  true long axis: [0.8253 0.5646]
eigenvalues 1 and 1: condition number 1, steps to reach 1e-6: 1
eigenvalues 1 and 10: condition number 10, steps to reach 1e-6: 132
eigenvalues 1 and 100: condition number 100, steps to reach 1e-6: 1375
condition number 10: flat coordinate shrinks by 0.9 per step; ceil(ln(1e-6)/ln(0.9)) = 132
condition number 100: flat coordinate shrinks by 0.99 per step; ceil(ln(1e-6)/ln(0.99)) = 1375
condition number 10: step 1/lambda_max takes 132, step 2/(lambda_max+lambda_min) takes 71
condition number 100: step 1/lambda_max takes 1375, step 2/(lambda_max+lambda_min) takes 709

Reading the output:

• Same direction. For both eigenvectors, A@v equals lambda*v. For the generic vector (1, 0), A@x is (3, 1), which points elsewhere.
• Hand and library agree. The roots of $\lambda^2 - 5\lambda + 5$ are the same 1.382 and 3.618, and the direction $(1, \lambda - 3)$ scaled to length 1 is (0.8507, 0.5257) for 3.618, as in the library output. For 1.382 it is (0.5257, -0.8507), the library's vector with both signs flipped.
• Sanity checks. The eigenvalues 1.382 and 3.618 sum to the trace 5 and multiply to the determinant 5, which matches the coefficients of the quadratic above.
• Symmetric matrices are the friendly case. A real symmetric matrix has real eigenvalues and mutually orthogonal eigenvectors (perpendicular: their dot product is zero; here $-0.5257 \times 0.8507 + 0.8507 \times 0.5257 = 0$). A set of vectors that are perpendicular and each of length 1 is called orthonormal. This is why eigh (the routine for symmetric matrices) is the right call for covariances and Hessians. A general matrix, such as a rotation, can have complex eigenvalues or too few independent eigenvectors.

Where eigenvalues show up in machine learning

1. PCA (principal component analysis). The covariance matrix $C$ of centered data (data with each feature's average subtracted) has the variance of each feature on its diagonal (the average squared distance from the mean) and, off the diagonal, how each pair of features moves together; it is symmetric. Its eigenvectors are the directions of the principal components, and its eigenvalues are the variance along each direction, meaning the spread of the data after every point is projected onto that direction. In the program the cloud was built by stretching one axis and rotating by 0.6 radians; the top eigenvector (0.8247, 0.5656) recovers the long axis (0.8253, 0.5646), and that direction carries 0.9739 of the total variance.
2. The Hessian of a loss. The Hessian is the matrix of second derivatives, and it measures curvature. Near a minimum its eigenvectors are the perpendicular directions along which the loss curves independently, and its eigenvalues are the curvature along each. Large eigenvalue means a steep direction; small means a flat one. All eigenvalues positive (the matrix is called positive definite) means the loss curves upward in every direction, a bowl with a minimum at the bottom. Mixed signs mean a saddle point: a flat spot that curves up in some directions and down in others, like a mountain pass.
3. Conditioning of a linear system or a quadratic. The condition number is the ratio of the largest to the smallest eigenvalue magnitude (for a symmetric positive definite matrix, meaning one with all eigenvalues positive). A value of 1 means every direction is equally curved, a round bowl; the larger it is, the longer and narrower the valley, and the harder it is for gradient descent (and for solving linear systems accurately) to treat all directions at once. Gradient descent on $f(x) = \tfrac{1}{2} x^\top H x$ with the standard step size $1/\lambda_{\max}$ (any step below $2/\lambda_{\max}$ still converges) moves slowly along the flat directions. The output counts the steps needed to bring the starting point (1, 1) within 1e-6 of the minimum: 1 step with condition number 1, 132 with condition number 10, and 1375 with condition number 100.

Why a flat direction slows things down: for this quadratic, one step multiplies each coordinate (measured along an eigenvector) by $1 - \lambda_i/\lambda_{\max}$. The steepest direction is multiplied by 0 and vanishes at once. The flattest direction is multiplied by $1 - 1/\kappa$, where $\kappa$ is the condition number: 0.9 per step at $\kappa = 10$ and 0.99 at $\kappa = 100$. Shrinking by a factor of $10^{-6}$ takes $\ln(10^{-6})/\ln(0.9) = 131.1$, so 132 steps, and $\ln(10^{-6})/\ln(0.99) = 1374.6$, so 1375 steps, exactly the counts printed. A larger step helps only a little. Above $2/\lambda_{\max}$ the steepest coordinate blows up, and the best fixed step, $2/(\lambda_{\max}+\lambda_{\min})$, cuts the count from 132 to 71 at $\kappa = 10$ and from 1375 to 709 at $\kappa = 100$, as the last two output lines show. That is about half the steps, not a cure: the flat direction still shrinks by roughly $1 - 2/\kappa$ per step, so the count still grows in proportion to $\kappa$.

Trade-offs and common mistakes

• Saying "an eigenvector is a direction that does not change". It can flip (negative $\lambda$) or shrink to zero; what does not change is the line it lies on.
• Treating eigenvectors as unique. Any nonzero multiple is also an eigenvector, so only the direction counts; when two eigenvalues are equal, any direction in their shared plane works.
• Using a general-purpose eigen solver on a covariance matrix. Use the symmetric routine: it is faster and guarantees real, sorted eigenvalues with orthonormal vectors.
• Forgetting that eigenvalues exist only for square matrices. For a rectangular data matrix, the analogue is the singular value decomposition (SVD), which factors any matrix into a rotation, a stretch along the axes and another rotation.

The idea in plain words

A sample covariance matrix $S$ for $d$ features is estimated from $n$ samples. After centering, the data matrix has rank at most $n-1$, so when $n \le d$ the $d \times d$ matrix $S$ has at least $d-n+1$ eigenvalues that are exactly zero. A matrix with a zero eigenvalue is singular, and with a near-zero one it is ill-conditioned: the condition number (ratio of largest to smallest eigenvalue for a symmetric positive semidefinite matrix, one whose eigenvalues are all zero or positive, as every covariance matrix is) is huge, and tiny errors in the data are amplified enormously when you invert it or solve with it.

Symptoms, with the numbers

The program draws 20 samples in 50 dimensions from a Gaussian whose true covariance is diagonal with eigenvalues from 4.0 down to 0.25 (true condition number 16). It prints the rank, the conditioning, and a held-out test of several repairs. Held-out data is fresh data from the same distribution that was not used to build the estimate, so it shows how well the estimate generalizes. The score used is the mean Gaussian log-likelihood: for each held-out row it asks how plausible that row is under a Gaussian with the fitted covariance, then averages. Higher means the fitted covariance describes new data better (the constant dropped from the score shifts every value by the same amount, so only differences between scores matter, and the values are not required to be negative); the true covariance gives the best score that is attainable on average, so the gap to it is the damage done by estimation error.

import numpy as np
from sklearn.covariance import LedoitWolf

rng = np.random.default_rng(0)
n, d = 20, 50                                   # fewer samples than features
true_eigs = np.linspace(4.0, 0.25, d)           # true covariance: diagonal, condition number 16
scale = np.sqrt(true_eigs)
X = rng.normal(size=(n, d)) * scale
Xnew = rng.normal(size=(2000, d)) * scale       # held-out data from the same distribution
S = np.cov(X, rowvar=False)                     # d x d sample covariance (divisor n-1)

ev = np.linalg.eigvalsh(S)[::-1]                # descending
print("shape of S:", S.shape, " rank:", np.linalg.matrix_rank(S), " (at most n-1 =", n - 1, ")")
print("nonzero eigenvalues:", int((ev > 1e-10).sum()), " of", d)
print("true condition number:", true_eigs.max() / true_eigs.min(), " sample condition number:", f"{np.linalg.cond(S):.2e}")
print("largest sample eigenvalue:", round(ev[0], 3), "(true largest 4.0)")

def shrink(S, a):
    # (1 - a) * S + a * mu * I, with mu = trace(S) / d
    mu = np.trace(S) / S.shape[0]
    return (1 - a) * S + a * mu * np.eye(S.shape[0])

def gaussian_loglik(Sigma, Z):
    # mean Gaussian log-likelihood of the rows of Z, dropping the constant -d/2 * log(2 pi)
    sign, logdet = np.linalg.slogdet(Sigma)
    P = np.linalg.inv(Sigma)
    return float(-0.5 * np.mean(np.einsum("ij,jk,ik->i", Z, P, Z)) - 0.5 * logdet)

print("\ntrue covariance, held-out log-likelihood:", round(gaussian_loglik(np.diag(true_eigs), Xnew), 3))
print("\nshrinkage a | cond(S_a) | held-out mean log-likelihood")
for a in (0.01, 0.1, 0.3, 0.6, 0.9, 1.0):
    Sa = shrink(S, a)
    print(f"{a:11.2f} | {np.linalg.cond(Sa):9.2e} | {gaussian_loglik(Sa, Xnew):10.3f}")

lw = LedoitWolf().fit(X)                        # scikit-learn centers the data and divides by n
print("\nLedoit-Wolf shrinkage chosen from the data:", round(lw.shrinkage_, 3),
      " held-out:", round(gaussian_loglik(lw.covariance_, Xnew), 3))

Output:

shape of S: (50, 50)  rank: 19  (at most n-1 = 19 )
nonzero eigenvalues: 19  of 50
true condition number: 16.0  sample condition number: 1.35e+18
largest sample eigenvalue: 15.735 (true largest 4.0)

true covariance, held-out log-likelihood: -38.983

shrinkage a | cond(S_a) | held-out mean log-likelihood
       0.01 |  7.83e+02 |  -1437.922
       0.10 |  7.21e+01 |   -143.218
       0.30 |  1.94e+01 |    -60.995
       0.60 |  6.27e+00 |    -46.183
       0.90 |  1.88e+00 |    -43.686
       1.00 |  1.00e+00 |    -43.914

Ledoit-Wolf shrinkage chosen from the data: 0.764  held-out: -44.455

• Rank. $S$ is 50 by 50 with rank 19 = n-1, so 31 of its 50 eigenvalues are zero (to rounding). Its sample condition number is about 1e18, which is set by rounding error in the zero eigenvalues, not by the data.
• Distorted spectrum. The true largest eigenvalue is 4.0 but the sample one is 15.735. Sample eigenvalues spread out: the large ones are overestimated and the small ones underestimated. This happens even when it is full rank.

Reading the helper code: np.linalg.slogdet returns the sign and the logarithm of the determinant without overflowing, and np.einsum("ij,jk,ik->i", Z, P, Z) computes $z^\top P z$ for every row $z$ of the held-out matrix Z at once, where P is the inverse covariance (also called the precision matrix). The score is the average of $-\tfrac12 z^\top P z - \tfrac12 \log\det \Sigma$.

What goes wrong downstream

Check these in order; the first three all contain the inverse or the determinant of $S$.

1. Inversion fails or returns garbage. All of these use $S^{-1}$: linear and quadratic discriminant analysis (classifiers that model each class as a Gaussian and score a point with the inverse covariance), Mahalanobis distance (distance measured in units of the data's own spread, $(x-\mu)^\top S^{-1}(x-\mu)$), Gaussian likelihoods, generalized least squares (regression that down-weights noisy, correlated errors using the inverse of their covariance) and minimum-variance portfolio weights (asset weights proportional to $S^{-1}\mathbf{1}$, the mix with the lowest variance). With a zero eigenvalue the inverse does not exist; with a tiny one, the inverse has huge entries (a direction with eigenvalue 0.001 counts a unit move 1000 times as heavily as one with eigenvalue 1) and the result is dominated by noise directions.
2. The log-determinant diverges. $\log\det S = \sum \log \lambda_i$ goes to $-\infty$ as an eigenvalue goes to 0, which wrecks any likelihood comparison.
3. Noisy small-eigenvalue directions get all the weight. Those directions have the least data support but the largest weight in $S^{-1}$, so the model chases noise (overfits).
4. Iterative solvers slow down, since convergence depends on the condition number.

Fixes, in the order to try them

The first two rows are the ones to reach for first because they cost one line and always give an invertible matrix. PCA comes next when size is the obstacle, structure when there is domain knowledge, and the truncated pseudo-inverse is a last resort because its threshold is a guess. More data is the cleanest fix but rarely on offer.

In the output, the held-out mean log-likelihood (higher is better, constant dropped) is -1437.922 at $a = 0.01$, which is nearly unregularized and badly conditioned (condition number 7.83e2), and improves steadily to -43.686 at $a = 0.9$. The true covariance scores -38.983, the ceiling. Full shrinkage to the identity ($a = 1$) scores -43.914, slightly worse than $a = 0.9$, because the true covariance is not the identity. Ledoit-Wolf chose $a = 0.764$ and scored -44.455. That is a one-seed comparison, so what it supports is that heavy shrinkage helps enormously when $n < d$, not a ranking of the last three settings. scikit-learn's LedoitWolf centers the data and divides by $n$ rather than $n-1$, a minor difference from np.cov.

Recommendation

Default to Ledoit-Wolf shrinkage for a Gaussian-style model, because it needs no tuning loop and its output is always invertible and well conditioned. Use PCA first when $d$ is so large that even a $d \times d$ matrix is too big to store. Choose a structured estimate when you have domain knowledge (for example independent sensors). If you only need $S^{-1}y$, solve the linear system with a Cholesky factorization of the shrunk matrix (writing it as $LL^\top$ with $L$ triangular, after which two cheap triangular solves give the answer) instead of forming the inverse.

Common mistakes

• Adding a tiny constant like 1e-10 to the diagonal and calling it fixed. It makes the code run but the smallest retained eigenvalue is then about that constant, so the condition number stays enormous and the noise directions still dominate.
• Using np.linalg.pinv silently, which hides the problem and still amplifies the smallest retained directions.
• Judging the fix on the training data, where the unregularized estimate looks perfect. Always judge on held-out data.

The idea in plain words

The singular value decomposition (SVD) writes any matrix $X$ (any shape, no square or symmetric requirement) as a rotation, a stretch along the axes, and another rotation:

$X = U \Sigma V^\top$

For a data matrix with $n$ rows (samples) and $d$ columns (features), after dropping empty directions:

• $U$ is $n \times r$ with orthonormal columns (each column has length 1 and any two columns are perpendicular, so $U^\top U = I$); these are the left singular vectors, one coordinate per sample.
• $\Sigma$ is an $r \times r$ diagonal matrix of singular values $\sigma_1 \ge \sigma_2 \ge \dots \ge 0$ (how much the matrix stretches along each axis).
• $V$ is $d \times r$ with orthonormal columns (the right singular vectors, directions in feature space).
• $r$ is the rank: the number of independent directions the matrix really uses, which equals its number of nonzero singular values and is at most $\min(n, d)$.

PCA (principal component analysis) finds the orthogonal directions in feature space that carry the most variance. The insight: if $X_c$ is the data with each column's mean subtracted, its SVD gives those directions directly. The columns of $V$ are the principal directions, and each singular value gives the variance along its direction.

Mechanism, in order

1. Center. Subtract each column's mean: $X_c = X - \mathbf{1}\mu^\top$, with $\mu$ the vector of column means. This step is not optional (see the last example below).
2. Decompose. $X_c = U \Sigma V^\top$.
3. Relate to the covariance. The sample covariance matrix is $C = \frac{1}{n-1} X_c^\top X_c$, and substituting the SVD gives

$C = \frac{1}{n-1} V \Sigma U^\top U \Sigma V^\top = V \,\frac{\Sigma^2}{n-1}\, V^\top$

since $U^\top U = I$. That is an eigendecomposition of $C$: the eigenvectors are the columns of $V$ and the eigenvalues are $\sigma_i^2 / (n-1)$.

Step 3 on four points: the program below takes four centered points in two features, computes $X_c^\top X_c/(n-1)$ directly, and checks that $V\,\mathrm{diag}(\sigma_i^2/(n-1))\,V^\top$ gives the same matrix. It then builds a matrix with singular values 1 and 1e-9 to show the condition-number effect discussed under the cost heading.

import numpy as np

# four centered points in 2 features: check C = V diag(S^2/(n-1)) V^T by hand-sized numbers
Xc = np.array([[3.0, 2.0], [1.0, 0.0], [-1.0, 0.0], [-3.0, -2.0]])
n = len(Xc)
U, S, Vt = np.linalg.svd(Xc, full_matrices=False)
C = Xc.T @ Xc / (n - 1)
print("Xc^T Xc / (n-1) =", C.tolist())
print("singular values:", np.round(S, 4), " S^2/(n-1):", np.round(S**2 / (n - 1), 4))
print("eigenvalues of C:", np.round(np.linalg.eigvalsh(C)[::-1], 4))
print("V diag(S^2/(n-1)) V^T equals C:", np.allclose(Vt.T @ np.diag(S**2 / (n - 1)) @ Vt, C))
print("columns of V have length 1 and are perpendicular (V^T V = I):", np.allclose(Vt @ Vt.T, np.eye(2)))
print("rank (number of nonzero singular values):", int((S > 1e-10).sum()))

# squaring the condition number: singular values 1 and 1e-9
rng = np.random.default_rng(0)
Q1, _ = np.linalg.qr(rng.normal(size=(50, 2)))
Q2, _ = np.linalg.qr(rng.normal(size=(2, 2)))
X = Q1 @ np.diag([1.0, 1e-9]) @ Q2.T
s = np.linalg.svd(X, compute_uv=False)
print("SVD of X, singular values:", s, " condition number:", f"{s[0] / s[1]:.2e}")
ev = np.linalg.eigvalsh(X.T @ X)[::-1]
print("eigenvalues of X^T X:", ev, " (true values 1 and 1e-18)")
print("sqrt of those eigenvalues:", np.sqrt(np.abs(ev)))

Output:

Xc^T Xc / (n-1) = [[6.666666666666667, 4.0], [4.0, 2.6666666666666665]]
singular values: [5.2361 0.7639]  S^2/(n-1): [9.1388 0.1945]
eigenvalues of C: [9.1388 0.1945]
V diag(S^2/(n-1)) V^T equals C: True
columns of V have length 1 and are perpendicular (V^T V = I): True
rank (number of nonzero singular values): 2
SVD of X, singular values: [1.00000000e+00 1.00000001e-09]  condition number: 1.00e+09
eigenvalues of X^T X: [1.00000000e+00 6.37266093e-17]  (true values 1 and 1e-18)
sqrt of those eigenvalues: [1.0000000e+00 7.9828948e-09]

Reading it: the covariance matrix is $\begin{pmatrix} 6.667 & 4 \\ 4 & 2.667 \end{pmatrix}$, its eigenvalues 9.1388 and 0.1945 equal $\sigma_i^2/(n-1)$ for the singular values 5.2361 and 0.7639 ($5.2361^2/3 = 9.1388$), the columns of $V$ have length 1 and are perpendicular, and the rank is 2. The last block of the output is used in the cost section below.

4. Variance explained. Variance along component $i$ is $\lambda_i = \sigma_i^2/(n-1)$, and the explained variance ratio is $\lambda_i / \sum_j \lambda_j$. The divisor $n-1$ is the sample-variance convention (the one np.cov and scikit-learn's PCA use); with divisor $n$ every variance scales by the same constant, so the ratios do not change.
5. Project. The component scores (coordinates of each sample in the new basis) are $X_c V_k = U_k \Sigma_k$, using the first $k$ columns.

Worked example

The program builds 200 samples in 5 features from 2 hidden factors plus small noise and a nonzero mean. It computes PCA three ways (SVD, eigendecomposition of the covariance, scikit-learn) and checks that they agree. It also prints a wrong-divisor check that is expected to fail, and shows what skipping the centering does.

import numpy as np
from sklearn.decomposition import PCA

rng = np.random.default_rng(0)
n, d = 200, 5
latent = rng.normal(size=(n, 2)) * np.array([4.0, 2.0])
W = rng.normal(size=(2, d))
X = latent @ W + 0.5 * rng.normal(size=(n, d)) + np.array([10.0, -3.0, 0.0, 5.0, 1.0])

# PCA through SVD of the CENTERED data matrix (n x d)
Xc = X - X.mean(axis=0)
U, S, Vt = np.linalg.svd(Xc, full_matrices=False)    # U: n x d, S: d, Vt: d x d
print("shapes U, S, Vt:", U.shape, S.shape, Vt.shape)

# variance along component i = S_i^2 / (n - 1)   (sample variance, divisor n-1)
var_svd = S**2 / (n - 1)
ratio = var_svd / var_svd.sum()
print("variance from SVD:       ", np.round(var_svd, 4))
print("explained variance ratio:", np.round(ratio, 4))

# same thing via eigendecomposition of the covariance matrix (also divisor n-1)
evals, evecs = np.linalg.eigh(np.cov(X, rowvar=False))
evals = evals[::-1]
print("eigenvalues of covariance:", np.round(evals, 4))
print("S^2/(n-1) equals eigenvalues:", np.allclose(var_svd, evals))
print("S^2/n equals eigenvalues:", np.allclose(S**2 / n, evals))

# compare with scikit-learn (it centers internally; components match up to sign)
pca = PCA(n_components=d).fit(X)
print("sklearn explained_variance_:", np.round(pca.explained_variance_, 4))
print("matches SVD variances:", np.allclose(pca.explained_variance_, var_svd))
same_axes = np.allclose(np.abs(pca.components_), np.abs(Vt))
print("components match up to sign:", same_axes)

# scores: projecting onto the top k components, two equivalent ways
k = 2
scores_a = Xc @ Vt[:k].T
scores_b = U[:, :k] * S[:k]
print("scores via Xc @ V equal U*S:", np.allclose(scores_a, scores_b))

# skipping the centering step changes the answer
U2, S_unc, Vt2 = np.linalg.svd(X, full_matrices=False)
print("top singular value, centered vs uncentered:", round(S[0], 3), round(S_unc[0], 3))
mean_dir = X.mean(0) / np.linalg.norm(X.mean(0))
print("|cosine| of first component with the mean direction, centered:  ", round(abs(Vt[0] @ mean_dir), 3))
print("|cosine| of first component with the mean direction, uncentered:", round(abs(Vt2[0] @ mean_dir), 3))

# truncated / randomized SVD: only the top 2 components, as used for large data
rp = PCA(n_components=2, svd_solver="randomized", random_state=0).fit(X)
print("randomized top-2 variances:", np.round(rp.explained_variance_, 4))
print("close to the full SVD's top 2:", np.allclose(rp.explained_variance_, var_svd[:2], rtol=1e-3))

Output:

shapes U, S, Vt: (200, 5) (5,) (5, 5)
variance from SVD:        [139.5319   8.4273   0.2688   0.2384   0.2119]
explained variance ratio: [0.9385 0.0567 0.0018 0.0016 0.0014]
eigenvalues of covariance: [139.5319   8.4273   0.2688   0.2384   0.2119]
S^2/(n-1) equals eigenvalues: True
S^2/n equals eigenvalues: False
sklearn explained_variance_: [139.5319   8.4273   0.2688   0.2384   0.2119]
matches SVD variances: True
components match up to sign: True
scores via Xc @ V equal U*S: True
top singular value, centered vs uncentered: 166.634 169.18
|cosine| of first component with the mean direction, centered:   0.032
|cosine| of first component with the mean direction, uncentered: 0.718
randomized top-2 variances: [139.5319   8.4273]
close to the full SVD's top 2: True

What the output shows:

• Shapes: $U$ is 200 by 5, $\Sigma$ is stored as a length-5 vector, and $V^\top$ is 5 by 5.
• $\sigma_i^2/(n-1)$ equals the covariance eigenvalues (True), and equals scikit-learn's explained_variance_ (True). Dividing by $n$ instead is False, which is the check that the divisor matters.
• The first component carries 0.9385 of the variance and the second 0.0567, so two components keep about 0.9952 of it, as expected from 2 hidden factors.
• Components match scikit-learn up to sign. The sign of a singular vector is arbitrary, so comparisons use absolute values.
• Without centering, the top singular value changes from 166.634 to 169.18 and the first direction tilts toward the mean. Its absolute cosine with the mean direction is 0.718 uncentered versus 0.032 centered; how far it tilts depends on how large the mean is compared with the spread of the data. Uncentered, the first "component" is partly just the direction from the origin to the data cloud, not the direction of greatest spread.

Cost for large datasets

• Full SVD of an $n \times d$ matrix costs on the order of $n d \min(n,d)$ operations. Forming the covariance first costs about $n d^2$ and then a $d \times d$ eigendecomposition costs about $d^3$. Squaring the data this way also squares the condition number (the ratio of the largest to the smallest singular value, a measure of how lopsided the matrix is), which loses precision for small singular values, so the SVD of the data matrix is the more accurate route. In the last block of the four-point program's output, a matrix with singular values 1 and 1e-9 (condition number 1e9) is decomposed correctly by the SVD (1 and 1.00000001e-09), but the eigenvalues of $X^\top X$ come out as 1 and 6.37e-17 instead of 1e-18. The second value is rounding noise: 64-bit floats resolve about 16 digits, so against an eigenvalue of 1 anything near 1e-16 cannot be told from zero, and the squared condition number 1e18 is beyond that reach. Its size and even its sign change from matrix to matrix (200 other random matrices with the same singular values gave values between about -3e-16 and +3e-16, negative roughly half the time, which is why the code takes the absolute value before the square root). Here the square root, 7.98e-09, is nearly 8 times the true 1e-9, and on another matrix it could be larger, or zero.
• For ordinary data (up to thousands of columns) the full SVD of the centered matrix is the right default. When you need only the top $k$ components (the number you want to keep) with $k \ll \min(n,d)$, use a truncated SVD (iterative Lanczos-type methods, which find the few leading singular vectors using only products of the matrix with vectors) or a randomized SVD. Randomized SVD multiplies $X_c$ by a random $d \times (k+p)$ matrix, where $p$ is a small number of extra random columns (the oversampling, for example 5 to 10) that make a good approximation more likely; the product is a sketch, a much smaller matrix $n \times (k+p)$ that captures the dominant directions of $X_c$. An exact SVD of the small sketch then gives the answer. Cost drops to roughly $n d (k+p)$. In the program, scikit-learn's randomized solver with random_state=0 returns top-2 variances that match the full SVD to the printed precision.
• Randomized methods are approximate. They are most accurate when singular values decay quickly; with a flat spectrum (singular values that decay slowly), add oversampling $p$ or power iterations (multiplying by $X_c^\top X_c$ a few times before sketching, which widens the gap between large and small singular values).

Common mistakes

• Forgetting to center (shown above), or centering but not scaling when features have different units.
• Confusing singular values with variances: the variance is the singular value squared, divided by $n-1$.
• Assuming eigenvectors and singular vectors come with a fixed sign.
• Choosing $k$ by a fixed variance target alone: look at the spectrum for a gap, and judge the downstream task.

The idea in plain words

$f(x) = \|Ax - b\|^2$ is the sum of squared residuals of a linear model: $A$ is an $m \times n$ matrix (m equations or samples, n unknowns or features), $x$ has n entries, and $b$ has m entries. It is a quadratic bowl in $x$. Setting its gradient to zero gives the normal equations, and the bowl has a single lowest point exactly when the columns of $A$ are linearly independent (no column can be built as a combination of the others).

Derivation

Write the squared norm as a product, $\|v\|^2 = v^\top v$, and expand:

$f(x) = (Ax - b)^\top (Ax - b) = x^\top A^\top A x - 2\, b^\top A x + b^\top b$

Use two standard rules for the gradient with respect to $x$ (the vector of partial derivatives): the gradient of $x^\top M x$ is $(M + M^\top)x$, which is $2Mx$ when $M$ is symmetric (and $A^\top A$ is), and the gradient of $c^\top x$ is $c$. With $M = A^\top A$ and $c = A^\top b$ (so $b^\top A x = c^\top x$):

$\nabla f(x) = 2 A^\top A x - 2 A^\top b = 2 A^\top (Ax - b)$

Where the two rules come from, on the smallest cases. In one dimension ($n = m = 1$, $A = a$ a single number) the loss is $f(x) = (ax - b)^2$, and the ordinary chain rule gives $f'(x) = 2a(ax - b) = 2a^2 x - 2ab$, which is the formula above with $A^\top A = a^2$ and $A^\top b = ab$. For general $n$, write the quadratic form as a double sum, $x^\top M x = \sum_{i,j} M_{ij} x_i x_j$. The partial derivative with respect to $x_k$ picks up the terms where $x_k$ appears as $x_i$ (giving $\sum_j M_{kj} x_j = (Mx)_k$) and as $x_j$ (giving $\sum_i M_{ik} x_i = (M^\top x)_k$), so the gradient is $(M + M^\top)x$. For the linear term, $c^\top x = \sum_i c_i x_i$ has partial derivative $c_k$ with respect to $x_k$, so its gradient is $c$.

Shapes: $A^\top$ is $n \times m$ and the residual $Ax - b$ has m entries, so the gradient has n entries, like $x$. Read it as "each column of A, dotted with the residual".

Setting it to zero gives the normal equations:

$A^\top A\, x = A^\top b$

They say the residual is orthogonal to every column of $A$: the best fit leaves nothing that a column could still explain. The Hessian (matrix of second derivatives) is $2 A^\top A$.

A hand-sized example

Three points $(t, b) = (0, 1), (1, 2), (2, 4)$, fitted by a line $b \approx x_0 + x_1 t$. Each row of $A$ is $(1, t)$ and the unknowns are the intercept $x_0$ and slope $x_1$. The normal equations $A^\top A\, x = A^\top b$ become a 2 by 2 system that can be solved with pencil and paper, and the program checks the arithmetic and the orthogonality statement.

import numpy as np

# Three points (t, b): (0, 1), (1, 2), (2, 4); fit b = x0 + x1 * t
A = np.array([[1.0, 0.0], [1.0, 1.0], [1.0, 2.0]])
b = np.array([1.0, 2.0, 4.0])
AtA, Atb = A.T @ A, A.T @ b
print("A^T A =", AtA.tolist())
print("A^T b =", Atb.tolist())
x = np.linalg.solve(AtA, Atb)
print("solution of A^T A x = A^T b:", x, " 5/6 =", round(5 / 6, 4))
r = A @ x - b
print("residual:", np.round(r, 4), " A^T r:", np.round(A.T @ r, 12))

Output:

A^T A = [[3.0, 3.0], [3.0, 5.0]]
A^T b = [7.0, 10.0]
solution of A^T A x = A^T b: [0.83333333 1.5       ]  5/6 = 0.8333
residual: [-0.1667  0.3333 -0.1667]  A^T r: [0. 0.]

By hand: $A^\top A = \begin{pmatrix} 3 & 3 \\ 3 & 5 \end{pmatrix}$ (its determinant is $15 - 9 = 6$), $A^\top b = (7, 10)$, and solving gives $x_0 = 5/6$ and $x_1 = 3/2$. The residuals sum to zero and are orthogonal to the slope column, which is exactly "$A^\top$ times the residual is zero".

When the solution is well defined

• Full column rank (the n columns of $A$ are linearly independent, which needs $m \ge n$): $A^\top A$ is positive definite (all eigenvalues positive), so it is invertible, the bowl is strictly convex (it curves upward in every direction, so there is exactly one lowest point), and the unique minimizer is $x^\star = (A^\top A)^{-1} A^\top b$.
• Rank deficient (the rank of $A$ is its number of independent columns, and here it is less than $n$: some column is a combination of the others, or $m < n$): $A^\top A$ is singular (not invertible, with a zero eigenvalue). A minimizer still exists (the normal equations are always solvable), but there are infinitely many, differing by vectors in the null space of $A$ (the vectors $z$ with $Az = 0$, which change $x$ without changing any prediction). The pseudoinverse $A^+$ (computed from the SVD, the factorization $A = U\Sigma V^\top$) returns the solution with the smallest norm, and adding a ridge penalty $\lambda \|x\|^2$ (which turns $A^\top A$ into $A^\top A + \lambda I$) makes the answer unique again.

Worked example

The program fits a line to six points, checks the analytic gradient against central finite differences, solves the normal equations and compares with np.linalg.lstsq. It then duplicates a column to make $A$ rank deficient and shows two different solutions with identical loss, adds ridge, and shows the condition number effect.

import numpy as np

np.set_printoptions(precision=4, suppress=True)
rng = np.random.default_rng(0)

# Full-column-rank case: 6 equations, 2 unknowns
A = np.array([[1.0, 0.0], [1.0, 1.0], [1.0, 2.0], [1.0, 3.0], [1.0, 4.0], [1.0, 5.0]])
b = np.array([1.1, 2.9, 5.2, 7.1, 8.8, 11.2])
f = lambda x: np.sum((A @ x - b) ** 2)
grad = lambda x: 2 * A.T @ (A @ x - b)

# 1. Check the derived gradient against central finite differences at a random point
x = rng.normal(size=2)
eps = 1e-6
num = np.array([(f(x + eps * e) - f(x - eps * e)) / (2 * eps) for e in np.eye(2)])
print("analytic gradient:   ", grad(x))
print("finite differences:  ", num)
print("match:", np.allclose(grad(x), num, atol=1e-5))
print("match if the factor 2 is dropped:", np.allclose(grad(x) / 2, num, atol=1e-5))

# 2. Set the gradient to zero: normal equations A^T A x = A^T b
x_star = np.linalg.solve(A.T @ A, A.T @ b)
print("normal equations solution:", x_star)
print("largest |gradient| there below 1e-9:", np.abs(grad(x_star)).max() < 1e-9)
print("lstsq agrees:", np.allclose(x_star, np.linalg.lstsq(A, b, rcond=None)[0]))
print("Hessian 2 A^T A eigenvalues (all > 0):", np.linalg.eigvalsh(2 * A.T @ A))

# 3. Not well defined: duplicate a column, so A^T A is singular
B = np.hstack([A, A[:, [1]]])                     # 6 x 3, columns 2 and 3 identical
print("rank of B:", np.linalg.matrix_rank(B), " columns:", B.shape[1])
print("smallest eigenvalue of B^T B is ~0:", abs(np.linalg.eigvalsh(B.T @ B)[0]) < 1e-9)
x_min = np.linalg.pinv(B) @ b                      # minimum-norm minimiser
x_alt = x_min + np.array([0.0, 1.0, -1.0])         # shift along the null space of B
print("min-norm solution:", x_min, " objective:", round(np.sum((B @ x_min - b) ** 2), 6))
print("shifted solution: ", x_alt, " objective:", round(np.sum((B @ x_alt - b) ** 2), 6))

# 4. Ridge makes it well defined again; and the normal equations square the condition number
lam = 0.1
x_ridge = np.linalg.solve(B.T @ B + lam * np.eye(3), B.T @ b)
print("ridge solution:", x_ridge)
C = np.vander(np.linspace(0, 1, 8), 6, increasing=True)   # polynomial design matrix, nearly dependent columns
print("cond(C):", f"{np.linalg.cond(C):.3e}", " cond(C^T C):", f"{np.linalg.cond(C.T @ C):.3e}",
      " cond(C)^2:", f"{np.linalg.cond(C) ** 2:.3e}")

Output:

analytic gradient:    [ -75.0544 -262.3596]
finite differences:   [ -75.0544 -262.3596]
match: True
match if the factor 2 is dropped: False
normal equations solution: [1.0429 2.0029]
largest |gradient| there below 1e-9: True
lstsq agrees: True
Hessian 2 A^T A eigenvalues (all > 0): [  3.5457 118.4543]
rank of B: 2  columns: 3
smallest eigenvalue of B^T B is ~0: True
min-norm solution: [1.0429 1.0014 1.0014]  objective: 0.134857
shifted solution:  [1.0429 2.0014 0.0014]  objective: 0.134857
ridge solution: [1.0046 1.0057 1.0057]
cond(C): 3.405e+03  cond(C^T C): 1.159e+07  cond(C)^2: 1.159e+07

Reading the output:

• The analytic and numeric gradients agree, and the version with the factor 2 dropped does not match, so the check can fail and the derivation passes it.
• The solution is intercept 1.0429 and slope 2.0029; the gradient there is zero to within 1e-9; the Hessian eigenvalues 3.5457 and 118.4543 are positive, so the minimum is unique.
• With a duplicated column, the rank is 2 for 3 columns and the smallest eigenvalue of $B^\top B$ is about 0. The two printed solutions differ by the vector (0, 1, -1) in the null space, and both give objective 0.134857. Ridge with $\lambda = 0.1$ picks one answer, splitting the weight evenly between the identical columns (1.0057 each).

Trade-offs and pitfalls

• Squared condition number. The condition number (largest singular value divided by smallest, a measure of how much input errors get amplified) is what matters here. For the polynomial design matrix, $\operatorname{cond}(C) = 3.405\text{e}3$ while $\operatorname{cond}(C^\top C) = 1.159\text{e}7$, which equals the square. Forming $A^\top A$ loses roughly twice as many digits as working with $A$, so numerical libraries solve least squares through a QR factorization ($A = QR$ with $Q$ having orthonormal columns and $R$ triangular) or an SVD of $A$ rather than inverting $A^\top A$.
• Never write the inverse explicitly. Solve the linear system; the inverse is slower and less accurate.
• Large n. Forming $A^\top A$ costs about $m n^2$ and solving costs about $n^3$; for huge problems use gradient-based or iterative methods (conjugate gradient, which improves the solution step by step using only products with $A$ and $A^\top$, or SGD) that never form $A^\top A$.
• Near-dependent columns do not make the problem undefined, but make the solution unstable to noise: small changes in $b$ move $x^\star$ a lot. Diagnose with the condition number or the singular values of $A$.