Eigenvalues and SVD: Calculate, Reconstruct, and Approximate

Decomposing a matrix into its action along a few directions reveals how it stretches vectors, which information the transformation loses, and how much error remains when fewer directions are kept. The linear algebra overview introduces vectors, rank, orthogonal projection, and least squares. Start here with a 2×2 eigenvalue calculation, then use a 3×2 matrix to work through SVD, reconstruction, and low-rank approximation. All matrices are real, vectors are columns, and $T$ denotes transpose.

Eigenvectors stay on an invariant line

For a square matrix $B$, a nonzero vector $v$ satisfying $Bv=\lambda v$ is an eigenvector, and $\lambda$ is its eigenvalue. MIT 18.06SC's eigenvalue notes begin with this relationship: the line through $v$ is invariant under the transformation. A positive $\lambda$ preserves orientation, a negative $\lambda$ reverses it, and $\lambda=0$ maps the vector to zero. The zero vector is excluded because it would satisfy the equation for every $\lambda$.

Take

$$B=\begin{bmatrix}2&1\\1&2\end{bmatrix}.$$

The equation $(B-\lambda I)v=0$ has a nonzero solution only if $B-\lambda I$ is singular. First solve the characteristic equation:

$$\det(B-\lambda I)=(2-\lambda)^2-1 =\lambda^2-4\lambda+3=(\lambda-3)(\lambda-1)=0.$$

This gives $\lambda_1=3$ and $\lambda_2=1$. Find the corresponding null spaces:

• For $\lambda=3$, $B-3I=\begin{bmatrix}-1&1\\1&-1\end{bmatrix}$, so the coordinates satisfy $v_2=v_1$.
• For $\lambda=1$, $B-I=\begin{bmatrix}1&1\\1&1\end{bmatrix}$, so the coordinates satisfy $v_2=-v_1$.

Normalize the vectors to obtain

$$q_1=\frac{1}{\sqrt2}\begin{bmatrix}1\\1\end{bmatrix},\qquad q_2=\frac{1}{\sqrt2}\begin{bmatrix}1\\-1\end{bmatrix}.$$

Direct multiplication checks $Bq_1=3q_1$ and $Bq_2=q_2$. Every $x=c_1q_1+c_2q_2$ therefore satisfies $Bx=3c_1q_1+c_2q_2$: the first coordinate is stretched by three, while the second keeps its length.

A real $n\times n$ matrix admits $B=P\Lambda P^{-1}$ with real $P,\Lambda$ exactly when it has $n$ linearly independent real eigenvectors. A real matrix can have complex eigenvalues; the invariant-line description above applies to real eigenvectors. For example, $J=\begin{bmatrix}1&1\\0&1\end{bmatrix}$ has the repeated eigenvalue 1, but the null space of $J-I$ is only the line spanned by $(1,0)^T$. Its eigenvectors cannot form a two-dimensional basis.

Symmetry gives an orthogonal decomposition

The spectral theorem in MIT's symmetric-matrix notes states that a real symmetric matrix has real eigenvalues and a complete orthonormal eigenbasis. For distinct eigenvalues $\lambda_i\ne\lambda_j$, orthogonality follows directly:

$$\lambda_j q_i^Tq_j=q_i^TBq_j=(Bq_i)^Tq_j =\lambda_i q_i^Tq_j, \qquad q_i^Tq_j=0.$$

For a repeated eigenvalue, arbitrary vectors in the same eigenspace need not be orthogonal, but an orthonormal basis can be chosen within it. With $Q=[q_1\ q_2]$, we have $Q^TQ=I$ and $Q^{-1}=Q^T$. Our example becomes

$$B=Q\begin{bmatrix}3&0\\0&1\end{bmatrix}Q^T =3q_1q_1^T+q_2q_2^T.$$

Each term on the right projects onto one axis and scales by its eigenvalue. Symmetry does not require positive eigenvalues; both happen to be positive in this example.

Rectangular SVD and factor shapes

A 3×2 matrix maps two-dimensional inputs to three-dimensional outputs, so $Av=\lambda v$ cannot directly compare its input and output. MIT's SVD notes establish that every real matrix has an SVD, using two sets of orthonormal directions to connect the spaces:

$$A=U\Sigma V^T,\qquad Av_i=\sigma_i u_i.$$

The right singular vector $v_i$ lies in the input space, and the left singular vector $u_i$ lies in the output space. Singular values are ordered $\sigma_1\ge\cdots\ge\sigma_p\ge0$, where $p=\min(m,n)$. Multiplication by $V^T$ extracts input coordinates, $\Sigma$ stretches them, and $U$ maps them into output coordinates.

Form, $A\in\mathbb R^{m\times n}$	$U$	$\Sigma$	$V^T$
Full SVD	$m\times m$	$m\times n$	$n\times n$
Economy SVD, $p=\min(m,n)$	$m\times p$	$p\times p$	$p\times n$
Nonzero terms only, $r=\operatorname{rank}(A)$	$m\times r$	$r\times r$	$r\times n$

In the full form, $U,V$ are square orthogonal matrices. In the other two forms, their columns are orthonormal, but the matrices may be rectangular. Economy SVD can still contain zero singular values: $p$ is not necessarily the rank $r$.

The SciPy linear algebra tutorial names the returned objects U, s, Vh. Here s is a one-dimensional singular-value array, and Vh is $V^T$ for real matrices, rather than $V$. The scipy.linalg.svd reference specifies the economy shapes obtained with full_matrices=False.

Reconstructing a matrix from two singular components

Take

$$A=\begin{bmatrix}2&0\\0&2\\-2&-2\end{bmatrix},\qquad A^TA=\begin{bmatrix}8&4\\4&8\end{bmatrix}=4B.$$

Because $A^TA v_i=\sigma_i^2v_i$, we can reuse the eigenvectors: $v_1=q_1$ and $v_2=q_2$, with eigenvalues 12 and 4. Thus

$$\sigma_1=\sqrt{12}=2\sqrt3,\qquad \sigma_2=2.$$

For each nonzero singular value, compute $u_i=Av_i/\sigma_i$:

$$\begin{aligned} Av_1&=\begin{bmatrix}\sqrt2\\\sqrt2\\-2\sqrt2\end{bmatrix},& u_1&=\frac{1}{\sqrt6}\begin{bmatrix}1\\1\\-2\end{bmatrix};\\[6pt] Av_2&=\begin{bmatrix}\sqrt2\\-\sqrt2\\0\end{bmatrix},& u_2&=\frac{1}{\sqrt2}\begin{bmatrix}1\\-1\\0\end{bmatrix}. \end{aligned}$$

These vectors have unit length and $u_1^Tu_2=0$. Setting $U_p=[u_1\ u_2]$, $\Sigma_p=\operatorname{diag}(2\sqrt3,2)$, and $V_p=[v_1\ v_2]$ gives an economy SVD with shapes $(3\times2)(2\times2)(2\times2)$. For the full form, add $u_3=(1,1,1)^T/\sqrt3$ and a bottom row of zeros to $\Sigma$; $A^Tu_3=0$.

Expand the product into outer products. Each nonzero component has rank 1:

$$\begin{aligned} A&=\sigma_1u_1v_1^T+\sigma_2u_2v_2^T\\ &=\begin{bmatrix}1&1\\1&1\\-2&-2\end{bmatrix} +\begin{bmatrix}1&-1\\-1&1\\0&0\end{bmatrix} =\begin{bmatrix}2&0\\0&2\\-2&-2\end{bmatrix}. \end{aligned}$$

Changing the signs of both $u_i,v_i$ leaves their outer product unchanged. For a real symmetric matrix, singular values are the absolute values of eigenvalues; negative signs are carried by the relative orientations of the left and right singular vectors. Singular values equal the eigenvalues themselves only in the positive semidefinite case.

Truncation and approximation error

For $1\le k<r$, keep only the largest $k$ components:

$$A_k=\sum_{i=1}^{k}\sigma_i u_i v_i^T.$$

Section 5 of Stanford CS168's low-rank approximation notes gives the optimality result: among matrices of rank at most $k$, $A_k$ minimizes Frobenius error. The norm is defined by $\|M\|_F^2=\sum_{a,b}M_{ab}^2$. Distinct outer products are orthogonal in this inner product because

$$\langle u_iv_i^T,u_jv_j^T\rangle_F =(u_i^Tu_j)(v_i^Tv_j).$$

The squared residual norm is therefore the sum of the discarded squared singular values:

$$\|A-A_k\|_F=\sqrt{\sum_{i=k+1}^{r}\sigma_i^2},\qquad \|A-A_k\|_2=\sigma_{k+1}.$$

The matrix $2$-norm measures the greatest length amplification of a unit input vector. The residual attains the stated value along $v_{k+1}$. If $k\ge r$, reconstruction error is zero.

Keeping the first component of our example gives

$$A_1=\begin{bmatrix}1&1\\1&1\\-2&-2\end{bmatrix},\qquad A-A_1=\begin{bmatrix}1&-1\\-1&1\\0&0\end{bmatrix}.$$

Hence $\|A-A_1\|_F=2$, $\|A\|_F=\sqrt{12+4}=4$, and relative Frobenius error is $1/2$. The retained fraction of squared norm is $12/(12+4)=75\%$. Retaining 75% of squared norm does not mean 25% relative norm error: the discarded squared error is 25%, and taking its square root gives the relative norm error.

PCA turns squared singular values into variance

The scikit-learn PCA documentation describes orthogonal components that explain the most variance in centered data. Centering does not automatically scale features to the same magnitude.

Let each row of $X\in\mathbb R^{N\times d}$ be an observation, with column means subtracted and $N>1$. Using the sample covariance convention with denominator $N-1$, as in the PCA API's explained variance, substitute the full SVD $X=U\Sigma V^T$:

$$C=\frac{X^TX}{N-1} =V\frac{\Sigma^T\Sigma}{N-1}V^T.$$

Thus $v_i$ is a principal direction, with variance $\sigma_i^2/(N-1)$. For the first $k$ directions, scores are $Z=XV_k=U_k\Sigma_k$ and reconstruction is $ZV_k^T=X_k$. Add the mean back to recover original coordinates when the original data were not centered.

Our matrix $A$ has zero column sums, so it can serve directly as $X$. With $N=3$, the covariance is $C=\begin{bmatrix}4&2\\2&4\end{bmatrix}$ and component variances are 6 and 2. The first component's scores are $(\sqrt2,\sqrt2,-2\sqrt2)^T$. Multiplying by $v_1^T$ reconstructs $A_1$, with explained variance $6/(6+2)=75\%$.

Clustering and dimensionality reduction discusses scaling, fitting data, and downstream evaluation. The 75% measures variance retained in these coordinates; it does not establish that a classification task retains enough information.

Least squares: zero and small singular values

For $\min_x\|Ax-b\|_2^2$, MIT's pseudoinverse notes construct the pseudoinverse by reciprocating the nonzero singular values. The least-squares solution of minimum Euclidean norm is

$$x_*=A^+b=\sum_{i=1}^{r}\frac{u_i^Tb}{\sigma_i}v_i.$$

This projects $b$ into the column space and undoes the stretching direction by direction. Every other least-squares solution has the form $x_*+z$ with $Az=0$. These null-space directions are orthogonal to $x_*$, so $x_*$ has the smallest norm.

For our example with $b=(1,0,0)^T$, the coefficients along the two singular directions are $1/(6\sqrt2)$ and $1/(2\sqrt2)$. This yields

$$x_*=\begin{bmatrix}1/3\\-1/6\end{bmatrix},\quad Ax_*=\begin{bmatrix}2/3\\-1/3\\-1/3\end{bmatrix},\quad b-Ax_*=\frac13\begin{bmatrix}1\\1\\1\end{bmatrix}.$$

The residual lies in the left null space: $A^T(b-Ax_*)=0$. The full-column-rank $A$ gives unique coefficients. Replacing it with $A_1$ gives $A_1v_2=0$, so adding $tv_2$ to any solution leaves predictions unchanged. Identifiable coefficients require that the matrix erase no input direction completely.

A small nonzero singular value creates a different problem: division by it amplifies data perturbations. With $A$ fixed, $\delta b=\varepsilon u_i$ produces $\delta x=(\varepsilon/\sigma_i)v_i$. For example, with $D=\operatorname{diag}(1,0.001)$, a perturbation of $10^{-6}$ in the second output coordinate becomes a coefficient perturbation of $0.001$.

For full column rank, the spectral condition number is $\kappa_2(A)=\|A\|_2\|A^+\|_2=\sigma_1/\sigma_n$. MIT's numerical-methods notes, Lecture 5, relate these norms to the extreme eigenvalues of $A^TA$ and show that $\kappa_2(A^TA)=\kappa_2(A)^2$: the eigenvalues are squared singular values. Our example gives $\sqrt3$ and 3 respectively; $D$ has condition number 1000. Numerical analysis separates problem sensitivity from algorithmic stability; ordinary least squares and regularization develops solver choices and coefficient interpretation.

Truncating an SVD for dimensionality reduction controls matrix reconstruction error. Truncating a pseudoinverse for solving discards coefficient directions and changes the estimate. Choose a threshold using data error and the intended use; a small nonzero singular value is not a mathematical zero.

Checking the components with Python

This standard-library code checks the derived eigenvectors and singular components, then computes reconstruction, rank-one error, and least-squares coefficients. The helpers take lists and require matching vector lengths or matrix row lengths; incompatible dimensions raise ValueError.

from math import sqrt, isclose


def dot(x, y):
    if len(x) != len(y):
        raise ValueError("Vector lengths must match")
    return sum(a * b for a, b in zip(x, y))


def mv(matrix, vector):
    return [dot(row, vector) for row in matrix]


B = [[2, 1], [1, 2]]
vectors = [[1 / sqrt(2), 1 / sqrt(2)],
           [1 / sqrt(2), -1 / sqrt(2)]]
for eigenvalue, v in zip([3, 1], vectors):
    assert all(isclose(a, eigenvalue * b, abs_tol=1e-12)
               for a, b in zip(mv(B, v), v))

A = [[2, 0], [0, 2], [-2, -2]]
s = [sqrt(12), 2.0]
left = [[a / sigma for a in mv(A, v)]
        for sigma, v in zip(s, vectors)]
assert isclose(dot(left[0], left[1]), 0, abs_tol=1e-12)
assert all(isclose(dot(u, u), 1) for u in left)
components = [[[sigma * u[i] * v[j] for j in range(2)]
               for i in range(3)]
              for sigma, u, v in zip(s, left, vectors)]
reconstructed = [[sum(c[i][j] for c in components) for j in range(2)]
                 for i in range(3)]
assert all(isclose(reconstructed[i][j], A[i][j], abs_tol=1e-12)
           for i in range(3) for j in range(2))
A1 = components[0]
error = sqrt(sum((A[i][j] - A1[i][j]) ** 2
                 for i in range(3) for j in range(2)))
b = [1, 0, 0]
x = [sum(v[j] * dot(u, b) / sigma
         for sigma, u, v in zip(s, left, vectors)) for j in range(2)]
print("eigenvalues:", [3, 1])
print("singular values:", [round(t, 6) for t in s])
print("reconstructed:", [[round(t, 6) for t in row] for row in reconstructed])
print("rank-one:", [[round(t, 6) for t in row] for row in A1])
print(f"Frobenius error: {error:.6f}")
print(f"retained variance: {s[0] ** 2 / sum(t * t for t in s):.6f}")
print("least-squares x:", [round(t, 6) for t in x])

Output:

eigenvalues: [3, 1]
singular values: [3.464102, 2.0]
reconstructed: [[2.0, 0.0], [0.0, 2.0], [-2.0, -2.0]]
rank-one: [[1.0, 1.0], [1.0, 1.0], [-2.0, -2.0]]
Frobenius error: 2.000000
retained variance: 0.750000
least-squares x: [0.333333, -0.166667]

My reading · Opened notes are remembered in this browser.