Eigenvalues, Eigenvectors, and Quadratic Forms

(a) Diagonal, so the eigenvalues are the diagonal entries: \(3,5\), with eigenvectors \(\ve_1,\ve_2\). Both positive → positive definite.
(For any diagonal matrix, \(\mD\ve_i = d_{ii}\ve_i\) directly from the definition — no characteristic polynomial needed.)

(b) \(\tr=5\), \(\det = 4-4=0\). So \(\lambda^2-5\lambda+0=0\), giving \(\lambda(\lambda-5)=0\): \(\lambda=5,0\).
For \(\lambda=5\): \(\begin{bmatrix}-1&2\\2&-4\end{bmatrix}\vv=\vzero \Rightarrow v_1=2v_2\), so \(\vv=(2,1)\).
For \(\lambda=0\): \(\begin{bmatrix}4&2\\2&1\end{bmatrix}\vv=\vzero \Rightarrow 2v_1+v_2=0\), so \(\vv=(1,-2)\).
Orthogonal ✓ (\(2-2=0\)). One eigenvalue is zero → positive semi-definite but not definite. The zero eigenvalue means \(\vv=(1,-2)\) is in the null space: the matrix is singular, and \(q(\vv)=0\) there even though \(\vv\neq\vzero\).

(c) \(\tr=0\), \(\det=-1\). \(\lambda^2-1=0\), so \(\lambda=\pm1\).
\(\lambda=1\): \(v_1=v_2\), giving \((1,1)\). \(\lambda=-1\): \(v_1=-v_2\), giving \((1,-1)\).
Mixed signs → indefinite. This is the swap matrix, and \(q(\vx)=2x_1x_2\), which is clearly positive in the first quadrant and negative in the second. A saddle.

Proof of the bound. Diagonalise \(\mA=\mQ\Lambda\mQ^\top\) and set \(\vz=\mQ^\top\vx\). Since \(\mQ\) is orthogonal, \(\|\vz\|=\|\vx\|\) (orthogonal matrices preserve length: \(\|\mQ^\top\vx\|^2 = \vx^\top\mQ\mQ^\top\vx = \vx^\top\vx\)). Then

Each \(z_i^2\ge0\), so replacing every \(\lambda_i\) by the largest can only increase the sum:

Replacing every \(\lambda_i\) by the smallest gives the lower bound identically. ∎

Equality. Take \(\vx=\vq_1\), so \(\vz=\ve_1\) and the sum is exactly \(\lambda_1\), while \(\|\vx\|^2=1\). Equality attained ✓. Similarly \(\vq_d\) attains the lower bound.

Equivalently: \(\lambda_1 = \max_{\vx\neq\vzero}\frac{\vx^\top\mA\vx}{\vx^\top\vx}\), the Rayleigh quotient. This variational characterisation is what makes PCA's "maximise the variance" and "take the top eigenvector" the same instruction.

The ridge claim. If \(\mA\vq=\lambda\vq\) then

so \(\vq\) is still an eigenvector, with eigenvalue shifted by \(c\). Since the eigenvectors are unchanged and every eigenvalue shifts by the same \(c\), the ordering is preserved and \(\lambda_{\max}(\mA+c\mI)=\lambda_{\max}(\mA)+c\), \(\lambda_{\min}(\mA+c\mI)=\lambda_{\min}(\mA)+c\). ∎

Hence for \(\mA=\mX^\top\mX\succeq0\) and \(c=\lambda>0\), every eigenvalue is at least \(\lambda>0\), so none is zero, so the matrix is invertible — the promise from Lecture 10. And the condition number becomes \(\frac{\lambda_{\max}+\lambda}{\lambda_{\min}+\lambda}\), which is strictly smaller than \(\frac{\lambda_{\max}}{\lambda_{\min}}\) whenever \(\lambda_{\max}>\lambda_{\min}\). Regularisation and preconditioning are literally the same operation.

1. Check by multiplying: \(\mQ\Lambda^{-1}\mQ^\top\,\mQ\Lambda\mQ^\top = \mQ\Lambda^{-1}(\mQ^\top\mQ)\Lambda\mQ^\top = \mQ\Lambda^{-1}\Lambda\mQ^\top = \mQ\mQ^\top = \mI\) ✓. And \(\Lambda^{-1}\) is diagonal with entries \(1/\lambda_i\), which exist because positive definiteness gives \(\lambda_i>0\).

2. Induction. \(\mA^2 = \mQ\Lambda\mQ^\top\mQ\Lambda\mQ^\top = \mQ\Lambda^2\mQ^\top\), the middle \(\mQ^\top\mQ\) collapsing to \(\mI\). Repeating gives \(\mA^k=\mQ\Lambda^k\mQ^\top\). ∎
This is why powers of matrices are easy in the eigenbasis and horrible outside it, and why the eigen-decomposition is the right tool for anything involving repeated application — Markov chains, discrete dynamical systems, and the gradient-descent recursion above.

3. Define \(\mA^{1/2} = \mQ\Lambda^{1/2}\mQ^\top\) with \(\Lambda^{1/2}=\operatorname{diag}(\sqrt{\lambda_i})\), well defined since \(\lambda_i>0\). Then \((\mA^{1/2})^2 = \mQ\Lambda^{1/2}\Lambda^{1/2}\mQ^\top = \mQ\Lambda\mQ^\top=\mA\) ✓, and it is symmetric (\((\mQ\Lambda^{1/2}\mQ^\top)^\top\) is itself) with positive eigenvalues \(\sqrt{\lambda_i}\), hence positive definite. ∎
Used constantly: sampling from \(\mathcal{N}(\vzero,\mS)\) is done as \(\mS^{1/2}\vz\) with \(\vz\sim\mathcal{N}(\vzero,\mI)\), and whitening a dataset is multiplication by \(\mS^{-1/2}\).

4. Write the Newton step in the eigenbasis of the Hessian. With \(\vg=\grad f\) and \(\tilde{\vg}=\mQ^\top\vg\), the step is

i.e. component \(i\) of the gradient is divided by \(\lambda_i\). Directions of high curvature get small steps, directions of low curvature get large ones. That is exactly what is needed to fix the zig-zag in the figure above: after the rescaling, every direction has effective curvature 1, the condition number becomes 1, and the ellipses become circles. Newton converges in a single step on a quadratic for precisely this reason.

And it identifies the cost: you need \(\mH^{-1}\), which is \(O(d^3)\) to compute and \(O(d^2)\) to store. Adam's diagonal rescaling is the cheap version — it fixes the conditioning that a diagonal Hessian would cause, and ignores the rest.

import numpy as np
rng = np.random.default_rng(0)
n, d, k_true = 500, 10, 3

Z = rng.normal(size=(n, k_true))          # latent
W = rng.normal(size=(k_true, d))          # loading
X = Z @ W + 0.35*rng.normal(size=(n, d))  # observed
Xc = X - X.mean(0)

S = Xc.T @ Xc / n
vals, vecs = np.linalg.eigh(S)
o = np.argsort(vals)[::-1]; vals, vecs = vals[o], vecs[:, o]
print("eigenvalues:", np.round(vals, 4))
print("explained  :", np.round(np.cumsum(vals)/vals.sum(), 4))

def recon_err(Q):                          # Q: d x k with orthonormal columns
    P = Q @ Q.T
    return np.mean(np.sum((Xc - Xc @ P)**2, axis=1))

for k in (1, 2, 3, 4):
    Qk = vecs[:, :k]
    best = recon_err(Qk)
    rand = min(recon_err(np.linalg.qr(rng.normal(size=(d, k)))[0])
               for _ in range(1000))
    print(f"k={k}  PCA err {best:.5f}   best of 1000 random subspaces {rand:.5f}"
          f"   {'PCA wins' if best <= rand else 'PROBLEM'}")

U, sv, Vt = np.linalg.svd(Xc, full_matrices=False)
print("\\nsigma^2 / n :", np.round(sv**2/n, 6)[:5])
print("eigenvalues  :", np.round(vals, 6)[:5])
print("match?", np.allclose(sv**2/n, vals))
print("top singular vector vs top eigenvector (up to sign):",
      np.allclose(np.abs(Vt[0]), np.abs(vecs[:, 0])))

What you should see.

2. The first three eigenvalues are large and comparable; the remaining seven are small and roughly equal to the noise variance \(0.35^2\approx0.12\). That gap — the "elbow" — is the visible signature of the true rank. Cumulative explained variance jumps to about 97–98% at \(k=3\) and then creeps up slowly, which is the standard way of choosing \(k\).

3. PCA beats every random subspace at every \(k\), usually by a wide margin. That is the Eckart–Young theorem in empirical form: the top-\(k\) eigen-subspace is the optimal rank-\(k\) approximation, not merely a good one. Note that reconstruction error and retained variance are two sides of one identity — total variance is fixed, so minimising what you throw away is maximising what you keep.

4. \(\sigma_i^2/n = \lambda_i\) to machine precision, and the right singular vectors equal the eigenvectors up to sign. This is the practical route: computing the SVD of \(\mX_c\) directly is both faster and numerically better than forming \(\mX_c^\top\mX_c\), because forming that product squares the condition number — the same lesson as lstsq versus inv in Lecture 10. Every real PCA implementation, including scikit-learn's, uses the SVD.

Sign ambiguity is worth noting: if \(\vq\) is an eigenvector so is \(-\vq\), so the sign of each principal component is arbitrary and can flip between library versions or platforms. If you ever see a PCA plot mirror itself for no reason, that is why.