3.8.1 The singular value decomposition

Turn, stretch, turn

Lesson 3.6.3 found that a symmetric matrix stretches the plane along two perpendicular axes: $A = Q\Lambda Q^\top$ . Most matrices are not symmetric, and many have no real eigenvectors at all. Yet every matrix, of any shape, still stretches along perpendicular axes.

The difference is that the axes it stretches from are not the axes it stretches onto. Take

A = \begin{pmatrix} 2 & 2 \\ -1 & 2 \end{pmatrix}

Its eigenvalues are $2 \pm 1.4142i$ , not real (lesson 3.6.1), so no direction is kept. But it sends the unit vector $\mathbf{v}_1 = \frac{1}{\sqrt 5}(1, 2)$ to $\frac{1}{\sqrt 5}(6, 3)$ , which is 3 times the unit vector $\mathbf{u}_1 = \frac{1}{\sqrt 5}(2, 1)$ . It sends the perpendicular $\mathbf{v}_2 = \frac{1}{\sqrt 5}(-2, 1)$ to $\frac{1}{\sqrt 5}(-2, 4)$ , twice $\mathbf{u}_2 = \frac{1}{\sqrt 5}(-1, 2)$ .

So two perpendicular directions go to two perpendicular directions, stretched by 3 and by 2. Collect the $\mathbf{v}$ 's as the columns of $V$ and the $\mathbf{u}$ 's as the columns of $U$ , and put 3 and 2 on the diagonal of $S$ :

\begin{gathered} U = \tfrac{1}{\sqrt 5}\begin{pmatrix} 2 & -1 \\ 1 & 2 \end{pmatrix} \\[6pt] S = \begin{pmatrix} 3 & 0 \\ 0 & 2 \end{pmatrix} \\[6pt] V = \tfrac{1}{\sqrt 5}\begin{pmatrix} 1 & -2 \\ 2 & 1 \end{pmatrix} \end{gathered}

Then $A = USV^\top$ . Read it right to left, as in lesson 3.2.2: $V^\top$ turns the plane, $S$ stretches it along the axes, and $U$ turns it again.

Rule: The singular value decomposition

\begin{gathered} A = U S V^\top \\[4pt] A\mathbf{v}_i = \sigma_i\,\mathbf{u}_i \end{gathered}

Every $m \times n$ matrix has one (stated here without proof). The columns $\mathbf{v}_i$ of $V$ and $\mathbf{u}_i$ of $U$ are orthonormal: unit length and at right angles. $S$ holds the singular values $\sigma_1 \ge \sigma_2 \ge \cdots \ge 0$ on its diagonal and zeros elsewhere.

The $\mathbf{v}_i$ are the right singular vectors, directions in the input. The $\mathbf{u}_i$ are the left singular vectors, directions in the output. A square matrix with orthonormal columns is orthogonal (lessons 3.5.2 and 3.6.3): $Q^\top Q = I$ , so $Q^{-1} = Q^\top$ , and it keeps every length and angle.

An orthogonal $2 \times 2$ matrix is a rotation, or a rotation and a reflection. Here $V$ is the rotation through 63.4° and $U$ the rotation through 26.6° (lesson 3.2.2), so $V^\top$ turns back through 63.4°. A singular value has the units of the entries of $A$ : it is a stretch factor.

The picture: a circle becomes an ellipse

The unit circle holds every unit vector $\mathbf{x}$ . Follow the three steps of $A = USV^\top$ . $V^\top$ turns the circle, which leaves it a circle, with $\mathbf{v}_1$ and $\mathbf{v}_2$ now on the axes.

$S$ stretches it into an ellipse with semi-axes 3 and 2 along the axes. $U$ turns that ellipse so that its semi-axes point along $\mathbf{u}_1$ and $\mathbf{u}_2$ .

$A = USV^\top$ in three steps. The unit circle (dashed in the later panels) ends as an ellipse whose semi-axes, $3\mathbf{u}_1$ and $2\mathbf{u}_2$ , are the images of $\mathbf{v}_1$ and $\mathbf{v}_2$ .

So the singular values are the longest and shortest stretches. A table of unit vectors at different angles shows it:

angle of $\mathbf{x}$	0°	30°	63.4°	90°	120°	153.4°
length of $A\mathbf{x}$	2.2361	2.7353	3	2.8284	2.349	2

The longest image, 3, comes from $\mathbf{v}_1$ at 63.4°; the shortest, 2, from $\mathbf{v}_2$ at 153.4°, a right angle further on. Unit vectors between them land between 2 and 3. Turning $\mathbf{x}$ round to $-\mathbf{x}$ gives the same length, so half the circle tells the whole story.

The ellipse also shows the area rule of lesson 3.4.1. The circle has area $\pi$ and the ellipse $\pi \times 3 \times 2$ , so $A$ multiplies areas by $3 \times 2 = 6$ . That is $\det A = 2 \cdot 2 - 2 \cdot (-1) = 6$ .

\lVert A\mathbf{x} \rVert

import numpy as np A = np.array([[2, 2], [-1, 2]]) v1_angle = np.degrees(np.arctan2(2, 1)) for deg in [0, 30, v1_angle, 90, 120, v1_angle + 90]: t = np.radians(deg) x = np.array([np.cos(t), np.sin(t)]) print(round(deg, 1), np.round(A @ x, 4), round(float(np.linalg.norm(A @ x)), 4)) t = np.linspace(0, np.pi, 1801) lengths = np.linalg.norm(A @ np.vstack([np.cos(t), np.sin(t)]), axis=0) print(round(lengths.max(), 4), round(np.degrees(t[lengths.argmax()]), 1)) print(round(lengths.min(), 4), round(np.degrees(t[lengths.argmin()]), 1))

⌘/Ctrl + Enter runs the cell

U

import numpy as np A = np.array([[2, 2], [-1, 2]]) U, s, Vt = np.linalg.svd(A) flip = np.where(Vt[:, -1] < 0, -1, 1) U, Vt = U * flip, Vt * flip[:, None] print(np.round(s, 4)) print(np.round(U * np.sqrt(5), 4)) print(np.round(Vt * np.sqrt(5), 4)) print(np.allclose(U @ np.diag(s) @ Vt, A)) print(np.round(A @ Vt[0], 4), np.round(s[0] * U[:, 0], 4))

⌘/Ctrl + Enter runs the cell

Finding them: the eigenvalues of AᵀA

The squared length of an image is $\lVert A\mathbf{x} \rVert^2 = (A\mathbf{x})^\top(A\mathbf{x}) = \mathbf{x}^\top(A^\top A)\mathbf{x}$ , by the transpose rule of lesson 3.2.3. So the stretches of $A$ are governed by $A^\top A$ , which is symmetric. By the spectral theorem (lesson 3.6.3) it has real eigenvalues and perpendicular eigenvectors.

Those eigenvectors are the $\mathbf{v}_i$ , and the eigenvalues are the $\sigma_i^2$ . In symbols, $A^\top A = VS^\top U^\top U S V^\top = V(S^\top S)V^\top$ , since $U^\top U = I$ .

Rule: Singular values from AᵀA

\begin{gathered} A^\top A = V\,(S^\top S)\,V^\top \\[4pt] \sigma_i = \sqrt{\lambda_i(A^\top A)} \\[4pt] \mathbf{u}_i = \frac{A\mathbf{v}_i}{\sigma_i} \end{gathered}

$S^\top S$ is the diagonal matrix of the $\sigma_i^2$ , and $\lambda_i(A^\top A)$ is the $i$ th eigenvalue of $A^\top A$ , largest first. Its eigenvectors are the $\mathbf{v}_i$ ; each $\mathbf{u}_i$ with $\sigma_i > 0$ follows from $A\mathbf{v}_i$ . In the same way $AA^\top = U(SS^\top)U^\top$ , so the $\mathbf{u}_i$ are its eigenvectors.

1
$A^\top A = \begin{pmatrix} 5 & 2 \\ 2 & 8 \end{pmatrix}$
2
$\lambda^2 - 13\lambda + 36 = (\lambda - 9)(\lambda - 4)$
3
$\sigma_1 = 3, \quad \sigma_2 = 2$
4
$\begin{gathered} \mathbf{v}_1 = \tfrac{1}{\sqrt 5}(1, 2) \\[4pt] \mathbf{v}_2 = \tfrac{1}{\sqrt 5}(-2, 1) \end{gathered}$
5
$\begin{gathered} \mathbf{u}_1 = \tfrac{1}{3\sqrt 5}(6, 3) = \tfrac{1}{\sqrt 5}(2, 1) \\[4pt] \mathbf{u}_2 = \tfrac{1}{2\sqrt 5}(-2, 4) = \tfrac{1}{\sqrt 5}(-1, 2) \end{gathered}$

A tall example

The method works for any shape. Take the $3 \times 2$ matrix

B = \begin{pmatrix} 2 & 0 \\ 1 & 1 \\ 0 & 2 \end{pmatrix}

Then $B^\top B = \begin{pmatrix} 5 & 1 \\ 1 & 5 \end{pmatrix}$ , with eigenvalues 6 and 4. So $\sigma_1 = \sqrt 6 \approx 2.4495$ and $\sigma_2 = 2$ , with $\mathbf{v}_1 = \frac{1}{\sqrt 2}(1, 1)$ and $\mathbf{v}_2 = \frac{1}{\sqrt 2}(-1, 1)$ .

The left singular vectors live in $\mathbb{R}^3$ . $B\mathbf{v}_1 = \frac{1}{\sqrt 2}(2, 2, 2)$ , so $\mathbf{u}_1 = \frac{1}{\sqrt 3}(1, 1, 1)$ ; $B\mathbf{v}_2 = \frac{1}{\sqrt 2}(-2, 0, 2)$ , so $\mathbf{u}_2 = \frac{1}{\sqrt 2}(-1, 0, 1)$ . $B$ maps the unit circle onto an ellipse with semi-axes $\sqrt 6$ and 2, lying flat in a plane of $\mathbb{R}^3$ .

$BB^\top$ is $3 \times 3$ , with eigenvalues 6, 4 and 0. The extra 0 belongs to $\mathbf{u}_3 = \frac{1}{\sqrt 6}(1, -2, 1)$ , perpendicular to that plane: a direction $B$ never reaches.

B

import numpy as np A = np.array([[2, 2], [-1, 2]]) lam, V = np.linalg.eigh(A.T @ A) print(A.T @ A, np.round(lam, 4)) print(np.round(np.sqrt(lam[::-1]), 4), np.round(np.linalg.svd(A, compute_uv=False), 4)) v1 = V[:, 1] * np.sign(V[-1, 1]) print(np.round(v1 * np.sqrt(5), 4), np.round(A @ v1 / 3 * np.sqrt(5), 4)) B = np.array([[2, 0], [1, 1], [0, 2]]) print(B.T @ B, np.round(np.linalg.eigvalsh(B.T @ B), 4)) print(np.round(np.linalg.eigvalsh(B @ B.T), 4) + 0.0) print(np.round(np.linalg.svd(B, compute_uv=False), 4))

⌘/Ctrl + Enter runs the cell

This is the method on paper, but a computer avoids forming $A^\top A$ . Squaring turns a singular value of $10^{-8}$ , beside one of 1, into an eigenvalue of $10^{-16}$ : the size of rounding errors, so its digits are lost. np.linalg.svd works on $A$ itself, as lesson 3.5.3's lstsq does; the last thing lstsq returns is the list of singular values of $A$ .

Rank, determinant and nearly flat ellipses

Lesson 3.2.3's singular matrix

Z = \begin{pmatrix} 1 & 2 \\ 0.5 & 1 \end{pmatrix}

has singular values 2.5 and 0. The ellipse has collapsed to a segment from $-2.5\mathbf{u}_1$ to $2.5\mathbf{u}_1$ , where $\mathbf{u}_1 = \frac{1}{\sqrt 5}(2, 1)$ .

The segment lies on the line $x_2 = x_1/2$ that lesson 3.2.3 found. The direction $\mathbf{v}_2 = \frac{1}{\sqrt 5}(-2, 1)$ is sent to $\mathbf{0}$ .

Each zero singular value is a direction squashed flat. So the rank, the number of pivots in lesson 3.3.2, is also the number of singular values that are not 0.

Rule: Rank and determinant

\begin{gathered} \operatorname{rank} A = \text{number of } \sigma_i > 0 \\[4pt] \lvert\det A\rvert = \sigma_1\sigma_2\cdots\sigma_n \end{gathered}

The second line is for a square $n \times n$ matrix: the ellipse's area (or the box's volume) over the circle's, as in lesson 3.4.1. One zero singular value makes the determinant 0.

On a computer, $Z$ 's second singular value prints as a tiny number such as $10^{-17}$ , not 0. np.linalg.matrix_rank counts the singular values above a tolerance: a few times $10^{-16}$ , times the largest one. That is how it allowed for rounding in lesson 3.3.2, where $\begin{pmatrix} 0.1 & 0.7 \\ 0.3 & 2.1 \end{pmatrix}$ came out with rank 1.

A singular value that is small but not 0 means a nearly flat ellipse. Lesson 3.3.2's $\begin{pmatrix} 1 & 1 \\ 1 & 1.001 \end{pmatrix}$ has singular values 2.0005 and 0.0005. Their ratio $\sigma_1/\sigma_n$ is the condition number that np.linalg.cond printed there: 4002.

Singular values against eigenvalues

For most matrices the two lists differ. The table compares them for $A$ , two symmetric matrices, and maps from lessons 3.2.2 and 3.2.3:

matrix	singular values	eigenvalues	$\lvert\det\rvert$
$A$	3, 2	none real	6
$\begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix}$	3, 1	3, 1	3
$\begin{pmatrix} 1 & 2 \\ 2 & 1 \end{pmatrix}$	3, 1	3, $-1$	3
quarter turn	1, 1	none real	1
shear $\begin{pmatrix} 1 & 1 \\ 0 & 1 \end{pmatrix}$	1.618, 0.618	1, 1	1
$Z$	2.5, 0	2, 0	0

Symmetric with eigenvalues at least 0: the SVD is the spectral decomposition, $U = V = Q$ and $S = \Lambda$ .
Symmetric with a negative eigenvalue: singular values are never negative. They are the sizes of the eigenvalues, and the minus sign moves into $U$ .
A rotation keeps every length, so the circle stays a circle: both singular values are 1.
The shear keeps only one direction, with eigenvalue 1 twice (lesson 3.6.1), yet it stretches by up to 1.618.

2 \times 2

import numpy as np def eig2(M): t = float(M[0, 0] + M[1, 1]) d = float(M[0, 0] * M[1, 1] - M[0, 1] * M[1, 0]) disc = t * t / 4 - d if disc < 0: return "none real" return round(t / 2 + disc**0.5, 4), round(t / 2 - disc**0.5, 4) examples = { "A": [[2, 2], [-1, 2]], "symmetric": [[2, 1], [1, 2]], "symmetric, negative eigenvalue": [[1, 2], [2, 1]], "quarter turn": [[0, -1], [1, 0]], "shear": [[1, 1], [0, 1]], "Z": [[1, 2], [0.5, 1]], } for name, M in examples.items(): M = np.array(M) s = np.linalg.svd(M, compute_uv=False) print(name, np.round(s, 4) + 0.0, eig2(M), round(abs(np.linalg.det(M)), 4), round(float(np.prod(s)), 4)) Z = np.array([[1, 2], [0.5, 1]]) s = np.linalg.svd(Z, compute_uv=False) tol = s.max() * max(Z.shape) * np.finfo(float).eps print(round(s[0], 4), s[1] < 1e-12, int(np.sum(s > tol)), np.linalg.matrix_rank(Z)) C = np.array([[1, 1], [1, 1.001]]) s = np.linalg.svd(C, compute_uv=False) print(round(float(s[0]), 4), round(float(s[1]), 4), round(s[0] / s[1]), round(np.linalg.cond(C)))

⌘/Ctrl + Enter runs the cell

A sum of rank-one pieces

Multiply out $USV^\top$ column by row, and the matrix splits into one piece per singular value. Each piece $\mathbf{u}_i\mathbf{v}_i^\top$ is a column times a row, a matrix (lesson 3.2.3). Every column of it is a multiple of $\mathbf{u}_i$ , so its rank is 1.

Rule: Rank-one pieces

\begin{gathered} A = \sigma_1\mathbf{u}_1\mathbf{v}_1^\top + \sigma_2\mathbf{u}_2\mathbf{v}_2^\top \\ + \cdots + \sigma_r\mathbf{u}_r\mathbf{v}_r^\top \end{gathered}

There are $r$ pieces, one for each nonzero singular value, where $r$ is the rank. The pieces come in order of size, largest first. Lesson 3.6.3 wrote a symmetric matrix the same way, as $\sum \lambda_i\mathbf{q}_i\mathbf{q}_i^\top$ .

For $B$ the pieces are whole numbers:

\begin{gathered} \sqrt 6\,\mathbf{u}_1\mathbf{v}_1^\top = \begin{pmatrix} 1 & 1 \\ 1 & 1 \\ 1 & 1 \end{pmatrix} \\[6pt] 2\,\mathbf{u}_2\mathbf{v}_2^\top = \begin{pmatrix} 1 & -1 \\ 0 & 0 \\ -1 & 1 \end{pmatrix} \end{gathered}

Their sum is $B$ . The first piece, $\sqrt 6 \cdot \frac{1}{\sqrt 3} \cdot \frac{1}{\sqrt 2} = 1$ times a table of ones, is the part of $B$ that points along $(1, 1, 1)$ . For $A$ the pieces are $\begin{pmatrix} 1.2 & 2.4 \\ 0.6 & 1.2 \end{pmatrix}$ and $\begin{pmatrix} 0.8 & -0.4 \\ -1.6 & 0.8 \end{pmatrix}$ .

Lesson 3.8.2 keeps only the first few pieces, to make the closest copy of a matrix with a lower rank. Case study 3.8.4 does it with a table of asset returns.

\mathbf{u}\mathbf{v}^\top

import numpy as np def pieces(M): U, s, Vt = np.linalg.svd(M, full_matrices=False) flip = np.where(Vt[:, -1] < 0, -1, 1) U, Vt = U * flip, Vt * flip[:, None] return [s[i] * np.outer(U[:, i], Vt[i]) for i in range(len(s))] B = np.array([[2, 0], [1, 1], [0, 2]]) for p in pieces(B): print(np.round(p, 4) + 0.0) print(np.round(sum(pieces(B)), 4) + 0.0) print([int(np.linalg.matrix_rank(p)) for p in pieces(B)]) A = np.array([[2, 2], [-1, 2]]) for p in pieces(A): print(np.round(p, 4) + 0.0)

⌘/Ctrl + Enter runs the cell

The SVD in NumPy

U, s, Vt = np.linalg.svd(A) returns three arrays:

s is a 1-D array of the singular values, largest first, not the matrix $S$ . For an $m \times n$ matrix it holds $\min(m, n)$ of them.
Vt is $V^\top$ , so its rows are the $\mathbf{v}_i$ .
U holds the $\mathbf{u}_i$ as its columns.

By default $U$ is $m \times m$ and $V^\top$ is $n \times n$ . For the tall $B$ , $U$ is $3 \times 3$ : it includes $\mathbf{u}_3$ , which no singular value uses. With full_matrices=False, $U$ keeps only the first $\min(m, n)$ columns, here $3 \times 2$ , and U @ np.diag(s) @ Vt rebuilds $B$ .

Rule: Shapes from np.linalg.svd

\begin{gathered} A:\ m \times n, \quad k = \min(m, n) \\[4pt] U:\ m \times k, \quad s:\ k, \quad V^\top:\ k \times n \end{gathered}

These are the shapes with full_matrices=False, the form to use for rebuilding $A$ or keeping pieces. A wide matrix works the same way: $B^\top$ has the same singular values as $B$ , with $U$ and $V$ swapped.

Signs. Flip both $\mathbf{u}_i$ and $\mathbf{v}_i$ and the piece $\sigma_i\mathbf{u}_i\mathbf{v}_i^\top$ is unchanged, so NumPy may return either sign, and two computers can differ. Fix a rule before printing or comparing; this lesson makes the last entry of each $\mathbf{v}_i$ positive.

B

import numpy as np B = np.array([[2, 0], [1, 1], [0, 2]]) U, s, Vt = np.linalg.svd(B) print(U.shape, s.shape, Vt.shape) Ur, s, Vt = np.linalg.svd(B, full_matrices=False) print(Ur.shape, np.allclose(Ur @ np.diag(s) @ Vt, B), np.allclose((Ur * s) @ Vt, B)) print(np.allclose(U.T @ U, np.eye(3)), np.allclose(Vt @ Vt.T, np.eye(2))) u3 = U[:, 2] * np.sign(U[0, 2]) print(np.round(u3 * np.sqrt(6), 4), np.round(B.T @ u3, 4) + 0.0) print(np.round(np.linalg.svd(B.T, compute_uv=False), 4))

⌘/Ctrl + Enter runs the cell

Common mistakes

Vt is $V^\top$ . Its rows are the $\mathbf{v}_i$ . Taking its columns gives the wrong vectors.
Flip signs in pairs. Changing the sign of $\mathbf{u}_i$ alone changes the matrix; flip $\mathbf{u}_i$ and $\mathbf{v}_i$ together, and fix a sign rule before comparing results.
Singular values are not eigenvalues. They agree only for symmetric matrices with no negative eigenvalues. The shear has eigenvalues 1 and 1 but singular values 1.618 and 0.618.
Singular values are never negative, and come largest first. A minus sign belongs in $U$ or $V$ .
The eigenvalues of $A^\top A$ are $\sigma_i^2$ , not $\sigma_i$ . Take square roots.
A tiny singular value prints as $10^{-17}$ , not 0. Count the rank with a tolerance, as matrix_rank does.
This $S$ is not a covariance matrix. Books call it $\Sigma$ ; here $\Sigma$ stays the covariance matrix, and $\sigma_i$ are not volatilities.

The singular value decomposition

0 of 12 correct

1.
$A^\top A$ has eigenvalues 16 and 9. What is $\sigma_1$ , the largest singular value of $A$ ?
2.
A $2 \times 2$ matrix has singular values 5 and 2. What is the shortest length of $A\mathbf{x}$ for a unit vector $\mathbf{x}$ ?
3.
The same matrix maps the unit circle to an ellipse. What is the ellipse's area? (type pi for π)
4.
The same matrix: what is $\lvert\det A\rvert$ ?
5.
Find $\sigma_1$ for the diagonal matrix $\begin{pmatrix} 3 & 0 \\ 0 & -5 \end{pmatrix}$
6.
A symmetric matrix has eigenvalues 3 and $-2$ . What is its smaller singular value?
7.
Find the only singular value of the $1 \times 2$ matrix $\begin{pmatrix} 3 & 4 \end{pmatrix}$
8.
A $4 \times 4$ matrix has singular values 6, 2, 0.5 and 0. What is its rank?
9.
A matrix has singular values 10 and 0.01. What is its condition number $\sigma_1 / \sigma_2$ ?
10.
$A = 3\mathbf{u}_1\mathbf{v}_1^\top + \mathbf{u}_2\mathbf{v}_2^\top$ with $\mathbf{u}_1 = (1, 0)$ , $\mathbf{v}_1 = (0.6, 0.8)$ , $\mathbf{u}_2 = (0, 1)$ , $\mathbf{v}_2 = (-0.8, 0.6)$ . Find $a_{11}$
11.
$A$ is $6 \times 4$ . With full_matrices=False, how many columns does $U$ have?
12.
$A$ is $6 \times 4$ . With the default full_matrices=True, how many columns does $U$ have?

How to type answers

Powers with ^: x^7, 2^(1/3), a^-2.
Multiply with * or by writing side by side: 3x^2, 2(x+1).
ln(x) is the natural log, log(x, 2) is log base 2, e is Euler's number, sqrt(x) is a square root.
sin(x), cos(x) and tan(x) take radians; write pi for π and sin(x)^2 for sin²x.
Your first answer to each question is saved: it counts towards your score and brings the question back for review.
Exact answers like ln(20)/ln(3) are best. Decimals need about 4 significant figures.

m \times n

The singular value decomposition

Turn, stretch, turn

Rule: The singular value decomposition

The picture: a circle becomes an ellipse

Finding them: the eigenvalues of AᵀA

Rule: Singular values from AᵀA

A tall example

Rank, determinant and nearly flat ellipses

Rule: Rank and determinant

Singular values against eigenvalues

A sum of rank-one pieces

Rule: Rank-one pieces

The SVD in NumPy

Rule: Shapes from np.linalg.svd

Common mistakes

The singular value decomposition

Summary card

Free further reading