A Deep Reading of SVD

The Summit

Talking about SVD in linear algebra can feel like jumping straight to the top of a mountain. Why begin at the highest point of the theory before we know the path that leads there?

Do not worry. We are not going to memorize the summit or disappear into a maze of proofs. We will stand there for a moment, then take it apart until we can see what the mountain is made of.

From Composition to Decomposition

Since the story is about decomposition, begin in the opposite direction: composition.

Put three matrices in front of us, $A$, $B$, and $C$. Their dimensions fit, so we can connect them:

That is composition in its simplest form: we know the pieces, we combine them, and a new object appears. If a vector $x$ enters the chain, the action begins on the right:

so

Take a small numerical example:

Together they become

Composition in Python

1import numpy as np2 3A = np.array([[1., 2.],4              [0., 1.]])5 6B = np.array([[2., 0.],7              [0., 0.5]])8 9C = np.array([[0., -1.],10              [1.,  0.]])11 12W = A @ B @ C13print(W)

This is the forward direction: known pieces produce a final matrix.

Here you can touch composition instead of merely reading its definition. Change $x$, press the matrices, and watch the vector move through $C$, then $B$, then $A$.

Now reverse the question.

Instead of giving you the pieces, I place one matrix $D$ in front of you. We can inspect its rows and columns, calculate its rank, and measure everything visible from the outside. Then I say: open it. Show me the structure hiding inside.

Without a method, we could guess until the stars burn out. That is precisely why linear algebra developed decompositions: we do not want arbitrary factors; we want pieces with structure and jobs we can name.

The pieces may be triangular, as in LU; orthonormal directions plus triangular coefficients, as in QR; or one triangular factor mirrored by its transpose, as in Cholesky. Sometimes we change the viewpoint itself, as in spectral decomposition and Schur. Sometimes we separate rotation from stretching, as in the polar decomposition. Then we reach SVD, which takes the stage from here.

The Decomposition Atlas

There are many names, but one instinct connects them: do not fight the matrix in the form in which it arrived. Find a representation in which each piece has a clear job.

In the interactive atlas below, choose LU, QR, Cholesky, Spectral, or SVD. Open the factors, then rebuild the original matrix. Do not memorize five formulas; watch the same move repeat.

A Small Decomposition Atlas in Python

1import numpy as np2 3A = np.array([[4., 2.],4              [2., 3.]])5 6# QR7Q, R = np.linalg.qr(A)8print(np.allclose(A, Q @ R))9 10# SVD11U, s, Vt = np.linalg.svd(A)12Sigma = np.diag(s)13print(np.allclose(A, U @ Sigma @ Vt))14 15# Spectral decomposition for symmetric A16lam, E = np.linalg.eigh(A)17print(np.allclose(A, E @ np.diag(lam) @ E.T))18 19# Cholesky for positive-definite A20L = np.linalg.cholesky(A)21print(np.allclose(A, L @ L.T))

With SciPy:

1import numpy as np2from scipy.linalg import lu, schur, polar3 4A = np.array([[4., 2.],5              [2., 3.]])6 7P, L, U_lu = lu(A)8T, Z = schur(A)9Qp, H = polar(A)10 11print(np.allclose(A, P @ L @ U_lu))12print(np.allclose(A, Z @ T @ Z.T))13print(np.allclose(A, Qp @ H))

From here, SVD no longer feels like a jump into a different chapter. It is simply the decomposition we are about to examine more closely than the others.

The Essence, Nothing More

For any real matrix

we can write

In the full picture, $U$ is an $m\times m$ orthogonal matrix, $V$ is an $n\times n$ orthogonal matrix, and $\Sigma$ is an $m\times n$ rectangular diagonal matrix.

On the diagonal of $\Sigma$ sit the singular values:

The columns of $V$ are the right singular vectors; the columns of $U$ are the left singular vectors. For complex matrices, $T$ is replaced by the conjugate transpose $^*$.

That is the formula. The meaning appears when each factor gets a job.

Compute and Check an SVD

1import numpy as np2 3A = np.array([[3., 1.],4              [1., 2.]])5 6U, s, Vt = np.linalg.svd(A)7Sigma = np.diag(s)8 9print("U =\n", U)10print("singular values =", s)11print("V^T =\n", Vt)12 13print("U^T U =\n", U.T @ U)14print("V^T V =\n", Vt @ Vt.T)15 16A_reconstructed = U @ Sigma @ Vt17print(A_reconstructed)18print("error =", np.linalg.norm(A - A_reconstructed))

For a rectangular matrix, use

1import numpy as np2 3A = np.array([[3., 1.],4              [1., 2.]])5U, s, Vt = np.linalg.svd(A, full_matrices=False)

for the compact/economy form when that is all you need.

The Story Begins on the Right

Feed in a vector $x$:

Begin with $V^T$. It asks how much of $x$ lies along each special direction $v_i$. If

then

So $V^T$ rewrites the input in coordinates designed specifically for this matrix.

Then $\Sigma$ acts. Nothing complicated is mixed together anymore; each coordinate is multiplied by one number:

This is the quiet center of SVD: in the right coordinates, a complicated transformation becomes independent stretching or shrinking along separate axes.

Finally, $U$ places those scaled components into their output directions.

Remember the story rather than the letters:

And one equation carries the whole idea:

Send in the special direction $v_i$, and $A$ sends it out along $u_i$, scaled by $\sigma_i$.

The Secret Behind the Beauty

Eigenvalue language is naturally square:

requires $Av$ and $v$ to live in the same vector space.

If $A:\mathbb R^n\rightarrow\mathbb R^m$ with $m\neq n$, that is not true in general.

SVD avoids the problem by using two different orthonormal coordinate systems:

Here $V$ belongs to the input space $\mathbb R^n$, while $U$ belongs to the output space $\mathbb R^m$.

The right singular vector $v_i$ lives where the input lives.

The left singular vector $u_i$ lives where the output lives.

They are connected by the equation

This equation is the heart of SVD.

Read it literally:

Feed the special input direction $v_i$ into $A$. The matrix does not throw it into some arbitrary mess. It sends it exactly into the special output direction $u_i$, scaled by $\sigma_i$.

There is a companion equation

for real matrices, or

for complex matrices.

When the Circle Changes Shape

Take the unit sphere in the input space.

The first orthogonal factor $V^T$ cannot deform it; rotations and reflections preserve a sphere.

Then $\Sigma$ stretches the sphere by different amounts along orthogonal axes. The sphere becomes an ellipsoid.

Finally $U$ rotates or reflects that ellipsoid into its final orientation in the output space.

The lengths of the principal semiaxes of the ellipsoid are the singular values.

That is why a large $\sigma_i$ means that the corresponding direction has a strong effect, while a tiny $\sigma_i$ means that direction barely survives the transformation.

The companion code can still generate static checkpoints of these stages, but the interactive version is more faithful to the idea: move through the stages and watch the same object change rather than comparing four disconnected pictures.

Let the Geometry Move

1import numpy as np2import matplotlib.pyplot as plt3 4A = np.array([[3., 1.],5              [1., 2.]])6 7U, s, Vt = np.linalg.svd(A)8Sigma = np.diag(s)9 10theta = np.linspace(0, 2*np.pi, 500)11circle = np.vstack([np.cos(theta), np.sin(theta)])12 13after_vt = Vt @ circle14after_sigma = Sigma @ after_vt15after_u = U @ after_sigma16 17plt.figure(figsize=(6, 6))18plt.plot(after_u[0], after_u[1])19plt.axhline(0)20plt.axvline(0)21plt.gca().set_aspect("equal", adjustable="box")22plt.title("Final image of the unit circle under A")23plt.show()

For teaching, it is even better to plot circle, after_vt, after_sigma, and after_u in separate figures so that the transformation can be watched one stage at a time.

Where Do These Pieces Come From?

The formula is beautiful, but we should not accept mysterious matrices that appear from smoke.

The first doorway is

For any real $A$, the matrix $A^TA$ is square and symmetric:

It is also positive semidefinite because for every vector $x$,

Therefore its eigenvalues are real and nonnegative, and it has an orthonormal eigenbasis.

Suppose

Then

So

Since $U^TU=I$,

But $\Sigma^T\Sigma$ is diagonal, with entries

Therefore:

and

So

Similarly,

so the left singular vectors are eigenvectors of $AA^T$.

This is the computational bridge that makes SVD feel less magical.

First Clue: $A^TA$

For a real matrix, an educational construction is easier to remember as a journey than as a recipe list. Begin with $A^TA$. Its orthonormal eigenvectors give the directions $v_i$, and its nonnegative eigenvalues $\lambda_i$ reveal the scales through

Whenever $\sigma_i\neq0$, send $v_i$ through $A$ and normalize the result:

So the three parts are not pulled from a hat: $A^TA$ chooses the important input directions, the eigenvalues determine how strongly they are stretched, and $A$ itself shows us where those directions land. Then $Av_i=\sigma_i u_i$.

Why are the $u_i$ unit vectors?

Why are different $u_i$ perpendicular?

For $i\neq j$,

because $v_i\perp v_j$.

So the nonzero singular directions already give orthonormal columns in both spaces.

This route is extremely useful for understanding and for hand calculations.

For the full existence proof, we now use the Hermitian block construction from the study notes. It looks like a trick at first; in a moment, it will feel inevitable.

Build SVD from $A^TA$

This is useful for understanding the theory. For numerical production work, prefer np.linalg.svd.

1import numpy as np2 3def svd_from_ata(A, tol=1e-12):4    A = np.asarray(A, dtype=float)5 6    # 1. Symmetric positive-semidefinite Gram matrix7    G = A.T @ A8 9    # 2. Its orthonormal eigenvectors become right singular vectors10    eigenvalues, V = np.linalg.eigh(G)11 12    # 3. Sort from largest to smallest13    order = np.argsort(eigenvalues)[::-1]14    eigenvalues = np.maximum(eigenvalues[order], 0.0)15    V = V[:, order]16 17    # 4. Singular values are square roots18    s = np.sqrt(eigenvalues)19 20    # 5. Keep the nonzero part for the condensed SVD21    keep = s > tol22    s = s[keep]23    V = V[:, keep]24 25    # 6. Av_i = sigma_i u_i26    U = A @ V / s27 28    return U, s, V.T

The Existence Proof

What Are We Actually Trying to Prove?

Let

have rank

We want to prove the existence of a condensed SVD

where

have orthonormal columns,

and

For real matrices, replace $^*$ with $^T$.

The proof has one main trick:

The Tool Behind the Door

A matrix $W$ is Hermitian if

The spectral theorem says that a Hermitian matrix has an orthonormal basis of eigenvectors and real eigenvalues. Equivalently,

where $Z$ is unitary and $\Lambda$ is real diagonal.

This theorem is the engine of the proof.

But $A$ may be rectangular, so $A$ itself may not even possess an ordinary eigendecomposition.

We need to build a square Hermitian matrix from it.

The Trick That Opens the Door

Define

Check the Block-Matrix Proof Numerically

1import numpy as np2 3A = np.array([[3., 1.],4              [0., 2.],5              [2., 2.]])6 7U, s, Vt = np.linalg.svd(A, full_matrices=False)8 9m, n = A.shape10W = np.block([11    [np.zeros((m, m)), A],12    [A.T, np.zeros((n, n))]13])14 15# Hermitian/symmetric, so eigvalsh is appropriate16w_eigenvalues = np.linalg.eigvalsh(W)17 18expected = np.sort(19    np.concatenate([20        -s,21        np.zeros(m + n - 2*len(s)),22        s23    ])24)25 26print(w_eigenvalues)27print(expected)28print(np.allclose(w_eigenvalues, expected))

This numerical experiment displays exactly what the proof predicts:

The upper-left zero block is $m\times m$.

The lower-right zero block is $n\times n$.

Therefore

so $W$ is square.

Now take the conjugate transpose:

Thus

So $W$ is Hermitian, and the spectral theorem applies.

That single construction is the door through which the SVD enters.

Split the Hidden Vector

Take an eigenvector $z$ of $W$ with eigenvalue $\sigma$:

Because $z\in\mathbb C^{m+n}$, write it in two blocks:

with

Substitute:

Block multiplication gives

Therefore

and

Stop here for a moment.

This is already the SVD relationship.

Compare

with

So the upper block $x$ behaves like a left singular vector, the lower block $y$ behaves like a right singular vector, and the eigenvalue $\sigma$ behaves like a singular value.

The singular structure of $A$ was hiding inside the eigenstructure of $W$.

Why the Values Come in Pairs

Suppose

has eigenvalue $+\sigma$. Then

Now consider

Applying $W$,

Hence

The nonzero eigenvalues of $W$ therefore occur in opposite pairs.

If $r=\operatorname{rank}(A)$, the nonzero part is organized as

The remaining eigenvalues are zero.

Give the Pieces Unit Length

For a positive eigenvalue $\sigma$, use the pair

Because $W$ is Hermitian and $+\sigma\neq-\sigma$ when $\sigma>0$, the eigenvectors belonging to those distinct eigenvalues are orthogonal:

Expanding,

So

Choose the block vector so that

Then

Together with equality of the two norms,

So each candidate singular vector has unit length.

Collect the Pieces

For the positive values

obtain vectors

and

such that

Define

and

The associated normalized eigenvectors of $W$ can be arranged as

and their eigenvalues as

The zero-eigenvalue eigenvectors do not contribute to $W=Z\Lambda Z^*$, because their entries in $\Lambda$ are zero.

Therefore the nonzero part alone gives

Put the Matrix Back Together

Substitute the block matrices:

The first two factors give

Multiplying the final block matrix gives

But by definition

Corresponding blocks must be equal. Hence

We have recovered the condensed SVD formula.

Why the Directions Stay Clean

Unit length has already been established.

We still need different columns to be perpendicular.

For two indices $i\neq j$, choose the positive-eigenvalue block vectors orthogonally:

Thus

Also compare the positive vector for index $i$ with the negative partner of index $j$:

which gives

Add the equations:

so

Subtract them:

so

Therefore

and

All requirements of the condensed SVD are satisfied.

Thus every matrix possesses a singular value decomposition.

The Proof in One Breath

If the details begin to blur, remember the argument. Start with $A$ and place it inside the Hermitian block matrix

The spectral theorem gives $W$ an orthonormal eigenbasis. Split each eigenvector into two parts, $x_i$ and $y_i$; the eigenvalue equations then give

Collect the $x_i$ and $y_i$ as columns of $X$ and $Y$, and place the singular values in $\Sigma_r$. Together, these pieces recover

The proof does not invent singular vectors from nothing. It finds them inside the eigenvectors of a carefully constructed Hermitian matrix.

From Product to Layers

Now comes a second interpretation of SVD, and for data science this one is just as important as the geometric picture.

Write

Because $\Sigma$ is diagonal,

Each outer product $u_i v_i^T$ has rank 1.

Why?

Every column of $u_iv_i^T$ is a scalar multiple of $u_i$. So its column space contains only one independent direction.

Therefore SVD says:

A matrix is an ordered stack of rank-1 layers.

The number $\sigma_i$ tells us how strong layer $i$ is.

Large singular value: strong layer.

Small singular value: weak layer.

Zero singular value: no layer at all.

Hence

The companion visualization script extracts the first rank-1 layers separately:

Looking at a Matrix

A very simple matrix heatmap is often enough to make structure visible:

1import numpy as np2import matplotlib.pyplot as plt3 4A = np.array([5    [5., 4., 3., 2., 1.],6    [4., 3.2, 2.4, 1.6, 0.8],7    [3., 2.4, 1.8, 1.2, 0.6],8    [2., 1.6, 1.2, 0.8, 0.4],9    [1., 0.8, 0.6, 0.4, 0.2],10])11 12plt.figure(figsize=(7, 5))13im = plt.imshow(A, aspect="auto")14plt.colorbar(im)15plt.title("Matrix values")16plt.xlabel("column")17plt.ylabel("row")18plt.show()

Use the same visualization as a sequence: begin with the original matrix $A$, isolate one rank-1 layer, rebuild a truncated approximation $A_k$, and finally inspect the residual $A-A_k$. That makes matrix decomposition feel like an object being opened rather than four unrelated pictures.

Pull Out the Rank-1 Layers

1import numpy as np2 3A = np.array([4    [5., 4., 3., 2., 1.],5    [4., 3.2, 2.4, 1.6, 0.8],6    [3., 2.4, 1.8, 1.2, 0.6],7    [2., 1.6, 1.2, 0.8, 0.4],8    [1., 0.8, 0.6, 0.4, 0.2],9])10U, s, Vt = np.linalg.svd(A, full_matrices=False)11layers = [s[i] * np.outer(U[:, i], Vt[i, :]) for i in range(len(s))]12 13# Exact reconstruction14A_again = sum(layers)15print(np.allclose(A, A_again))

Each layer_i has rank at most 1.

To visualize a layer:

1import matplotlib.pyplot as plt2import numpy as np3 4A = np.array([5    [5., 4., 3., 2., 1.],6    [4., 3.2, 2.4, 1.6, 0.8],7    [3., 2.4, 1.8, 1.2, 0.6],8    [2., 1.6, 1.2, 0.8, 0.4],9    [1., 0.8, 0.6, 0.4, 0.2],10])11U, s, Vt = np.linalg.svd(A, full_matrices=False)12first_layer = s[0] * np.outer(U[:, 0], Vt[0, :])13 14plt.figure()15plt.imshow(first_layer, aspect="auto")16plt.colorbar()17plt.title("First rank-1 SVD layer")18plt.show()

How Many Pieces Do We Need?

This connects SVD with the deeper meaning of matrix rank.

If $\operatorname{rank}(A)=r$, then $A$ can be written as a sum of $r$ rank-1 matrices:

SVD provides one such representation:

Why can we not do it with fewer than $r$ rank-1 pieces?

Because rank is subadditive:

If $A$ were a sum of only $r-1$ rank-1 matrices, then $\operatorname{rank}(A)\le r-1$, contradicting $\operatorname{rank}(A)=r$.

Therefore

This is one of the cleanest bridges from abstract rank to something you can almost touch.

Rank Factorization Meets SVD

A rank-$r$ matrix also admits a rank factorization

with

The columns of $Y$ provide basic directions, while rows of $Z^T$ provide coefficients for mixing them.

The inner dimension cannot be smaller than $r$. If

with only $r'<r$ inner columns, then every column of $A$ would lie in the span of the $r'$ columns of $Y'$, implying

which is impossible.

SVD is a special, highly structured rank factorization. In condensed form,

One can group the first two factors:

so

What SVD adds is orthonormality and a natural ordering by importance.

When We Agree to Lose a Little

Exact decomposition is not always the goal.

Suppose

with

If the later singular values are small, their layers contribute little.

So keep only the first $k$ layers: