๐Ÿ” Full SVD: Pseudoinverse, Truncation and Rank Control

By Better Wise

QR and Cholesky both assume you can solve the system you were handed. The SVD tells you whether that assumption is worth anything, and what to do when it is not. If $A = QR$ fails because $A$ is rank-deficient or so nearly so that Cholesky's pivots underflow, the SVD is where you go to find out which directions actually carry information.

After this lesson you can:

  • Write the full and thin SVD and say what each factor contributes.
  • Build the Moore-Penrose pseudoinverse and state exactly which system it solves.
  • Apply truncated SVD to a rank-deficient least-squares problem by hand.
  • Read a singular value spectrum and defend a truncation rank choice.

Recall from the previous lesson: QR gives you an orthogonal basis for the column space via Householder reflections, and Cholesky gives $A = LL^\top$ for SPD $A$. Both need full column rank. Neither tells you how close to rank-deficient you are. SVD does both jobs in one factorization.

The full SVD and what each factor does

For $A \in \mathbb{R}^{m \times n}$ with $m \ge n$:

$$ A = U \Sigma V^\top $$

$U \in \mathbb{R}^{m \times m}$ is orthogonal, $V \in \mathbb{R}^{n \times n}$ is orthogonal, and $\Sigma \in \mathbb{R}^{m \times n}$ is diagonal with $\sigma_1 \ge \sigma_2 \ge \cdots \ge \sigma_n \ge 0$.

Three readings of the same identity, all useful:

  • Geometry. $V^\top$ rotates the input into a coordinate frame where $A$ acts as pure axis scaling. $\Sigma$ scales. $U$ rotates into the output frame.
  • Sum of rank-1 terms. $A = \sum_{i=1}^{n} \sigma_i u_i v_i^\top$. Each term is an outer product; each $\sigma_i$ says how much that term contributes to $A$.
  • Column and row spaces. The columns $u_1 \dots u_r$ span the column space, the $v_1 \dots v_n$ give the input directions, and $r$ is the number of nonzero singular values = $\operatorname{rank}(A)$.

The thin (economy) SVD keeps only the first $n$ columns of $U$ and the top $n \times n$ block of $\Sigma$. Memory $O(mn)$ instead of $O(m^2)$, same information for solving. Use thin in code, think full when reasoning.

Note: Singular values are not eigenvalues. For SPD matrices the two coincide. In general $\sigma_i(A) = \sqrt{\lambda_i(A^\top A)}$, which squares the condition number. That squaring is why normal equations are dangerous and SVD is not.

The pseudoinverse: the solution that always exists

If $A$ is square and invertible, $A^{-1}$ exists and solves everything. If it is not, the equation $Ax = b$ may have no solution (overdetermined, inconsistent) or infinitely many (underdetermined or rank-deficient). The Moore-Penrose pseudoinverse handles all cases uniformly:

$$ A^{+} = V \Sigma^{+} U^\top, \qquad \Sigma^{+}_{ii} = \begin{cases} 1/\sigma_i & \sigma_i > 0 \ 0 & \sigma_i = 0 \end{cases} $$

So $A^{+}$ inverts every nonzero singular value and sends exact zeros to zero. Because $\sigma_i > 0$ is a test on a floating-point number, in practice you compare against a tolerance $\tau$, not against zero (more below).

What $x = A^{+}b$ actually is. Not "the solution", but the minimum-norm least-squares solution: it minimizes $|Ax - b|_2$ and, among all minimizers, minimizes $|x|_2$. That single sentence covers four classical cases: unique solution, inconsistent system (returns the projection of $b$ onto the column space), underdetermined system (returns the smallest-norm answer), and rank-deficient system.

Worked example 1: a rank-1 system, by hand

Let

$$ A = \begin{pmatrix} 1 & 1 \ 1 & 1 \end{pmatrix}, \qquad b = \begin{pmatrix} 2 \ 1 \end{pmatrix} $$

The system is inconsistent and $\operatorname{rank}(A) = 1$. Compute $A^\top A = \begin{pmatrix} 2 & 2 \ 2 & 2 \end{pmatrix}$, eigenvalues $4$ and $0$, so $\sigma_1 = 2$, $\sigma_2 = 0$. The unit eigenvectors give $v_1 = \tfrac{1}{\sqrt{2}}(1,1)^\top$, $v_2 = \tfrac{1}{\sqrt{2}}(1,-1)^\top$. And $u_1 = A v_1 / \sigma_1 = \tfrac{1}{\sqrt{2}}(1,1)^\top$.

Now the pseudoinverse, treating $\sigma_2 = 0$ correctly:

$$ A^{+} = \frac{1}{2} u_1 v_1^\top = \frac{1}{2} \cdot \frac{1}{2}\begin{pmatrix} 1 \ 1 \end{pmatrix}\begin{pmatrix} 1 & 1 \end{pmatrix} = \begin{pmatrix} 0.25 & 0.25 \ 0.25 & 0.25 \end{pmatrix} $$

$$ x = A^{+} b = \begin{pmatrix} 0.75 \ 0.75 \end{pmatrix} $$

Check: $Ax = (1.5, 1.5)^\top$, residual $(0.5, -0.5)^\top$, which is orthogonal to the column space direction $(1,1)^\top$. Good, and $|x|$ is as small as possible among least-squares solutions. Note what a naive solver would have done: forming $A^\top A$ and inverting gives a singular matrix, and adding $10^{-12} I$ "for safety" gives $x \approx (6.25, 6.25)$, nonsense of magnitude far larger than the data. The zero singular value told us to throw that direction away, not to nudge it.

Worked example 2: truncation on a nearly-rank-1 system

Now the realistic case. Take the same picture but with a tiny perturbation:

$$ A = \begin{pmatrix} 1 & 1 \ 1 & 1 + \epsilon \end{pmatrix}, \qquad \epsilon = 10^{-8} $$

Now $\sigma_2 \approx \epsilon / \sqrt{2} \approx 7 \times 10^{-9}$: nonzero, so exact pseudoinverse inverts it and produces a coefficient of order $10^8$. Concretely, in float64 the computed $\sigma_2$ carries relative error around $10^{-16}/\epsilon = 10^{-8}$, so the recovered $1/\sigma_2$ term has meaning only in its first eight or nine digits and, more importantly, an irrelevant direction has been amplified to dominate the answer. In float32 ($\varepsilon \approx 1.2 \times 10^{-7}$) the computed $\sigma_2$ may be entirely noise.

The truncated pseudoinverse replaces the tail with zeros:

$$ A_k^{+} = \sum_{i=1}^{k} \frac{1}{\sigma_i} v_i u_i^\top, \qquad k = #{i : \sigma_i > \tau} $$

Choosing $k=1$ here returns the minimum-norm least-squares solution of the nearest rank-1 problem, which is exactly what a stable answer for this data looks like. This is the central trade: you are choosing to solve a slightly different, well-posed problem rather than the literal one.

Reading a spectrum to choose the rank

Plot $\sigma_i$ on a log scale against index $i$. The shape of that curve is your diagnostic, and there are three characteristic patterns:

flowchart TD
    S["๐Ÿ“‰ Log singular value spectrum"] --> A["Gap: flat then cliff<br/>ฯƒ drops by orders of magnitude"]
    S --> B["Smooth decay: ฯƒ_i ~ i^-ฮฑ<br/>power-law, no gap"]
    S --> C["Plateau at noise floor<br/>ฯƒ_i โ‰ˆ ฮต ยท ฯƒ_max"]
    A --> A1["Truncate at the cliff<br/>numerical rank is unambiguous"]
    B --> B1["No natural rank<br/>choose k by validation or budget"]
    C --> C1["Everything below plateau<br/>is discretization noise"]
  • Gap. If $\sigma_k / \sigma_{k+1}$ is large (say $10^3$ or more), truncate at $k$. This is the easy case: the data was generated by a low-rank process plus small noise.
  • Smooth decay. Power-law spectra (image operators, kernel matrices, discretized inverse problems) have no gap. There is no correct $k$; you pick one by held-out error, by a variance-explained target, or by compute budget, and you report that you did so.
  • Noise floor. Any $\sigma_i$ below roughly $\varepsilon \cdot \sigma_{\max} \cdot \sqrt{\max(m,n)}$ is not a property of your matrix, it is a property of your float format. Truncating there is not regularization, it is bug avoidance.

Practical tolerance. A defensible default is

$$ \tau = \sigma_{\max} \cdot \max(m, n) \cdot \varepsilon $$

with $\varepsilon$ the machine epsilon of your working precision. This is the standard rank test used by LAPACK-style matrix_rank. It is conservative, so treat it as the smallest truncation you should consider, not a recommendation.

Rules of thumb for choosing $k$:

  1. Report the condition number of the truncated operator, $\sigma_1/\sigma_k$, not of $A$. Truncation is a conditioning fix, and you should state what condition number you bought.
  2. Add singular values in decreasing order and track the reconstruction error $|A - A_k|F^2 = \sum{i>k}\sigma_i^2$. Stop when the error is inside your noise budget. This is the Eckart-Young guarantee: truncated SVD is the optimal rank-$k$ approximation in both Frobenius and spectral norm, so no other rank-$k$ matrix does better.
  3. In a solve, check the residual you actually got against the residual the truncation predicted, $|b - Ax| \approx |(I - U_kU_k^\top)b|$. If the real residual is much larger, your model $A$ is wrong, not your solver.
  4. Do not tune $k$ on the same data you evaluate on. A too-large $k$ fits noise and looks great in-sample.

Which singular values are trustworthy is a claim about $\sigma_i$ of the form $\sigma_i$: _____ in float arithmetic. And note the relationship you should never forget: {{c2::the pseudoinverse inverts nonzero singular values and sends zero ones to zero::A^{+} = V \Sigma^{+} U^T}. The truncated pseudoinverse with rank $k$ is written {{c3::A_k^+ = sum_{i=1}^{k} (1/sigma_i) v_i u_i^T}}.

Practitioner notes on SVD-based solving

  • Never form $A^\top A$ yourself. Its condition number is $\kappa(A)^2$, and you lose half your mantissa digits before the solver starts. Call svd or lstsq on $A$ directly.
  • torch.linalg.svd returns $V^\top$ (or Vh), not $V$. The single most common bug in pseudoinverse code is transposing the wrong factor. Check shapes and check $U\Sigma V^\top \approx A$ before trusting anything downstream.
  • Sort order. Library SVDs return singular values in descending order, but if you build $\Sigma^+$ yourself, a mismatched order silently pairs $1/\sigma_1$ with the smallest direction. Verify with a tiny matrix where you know the answer.
  • Differentiability. Gradients of the SVD blow up as $\sigma_i \to \sigma_{i+1}$; truncation also removes those terms from the backward pass, which is usually a feature (it stabilizes training) and occasionally a silent zero-gradient bug.
  • Rank control is a modeling decision. Truncation changes which problem you solve. Say in your report that you solved the rank-$k$ problem, and why $k$ was chosen.

Wrap-up

The SVD gives you the whole picture of a linear map in one factorization: which directions exist, how much each is amplified, and which are numerically fictional. The pseudoinverse is the universal solver, truncation is the knob that keeps it finite, and the singular value spectrum is the evidence for where to set that knob. Optimality of the truncation comes free from Eckart-Young; defensibility of $k$ does not.

Explicit truncation is a hard cutoff, and hard cutoffs have a soft cousin: adding a penalty that shrinks all singular directions smoothly, with the same limiting behavior. The next lesson makes that precise with conditioning numbers and Tikhonov regularization, and shows when the smooth version is the better tool.

Front: What is the rank of a matrix in terms of its SVD? Back: The number of strictly positive singular values, sigma_i > 0.

Front: What does the Moore-Penrose pseudoinverse return for an inconsistent overdetermined system? Back: The minimum-norm least-squares solution, x = A^+ b.

Front: How is Sigma^+ defined? Back: Transpose of Sigma with each nonzero sigma_i replaced by 1/sigma_i and zeros kept as zeros.

Front: Why is forming A^T A before solving a bad idea? Back: It squares the condition number, so you lose roughly half your significant digits.

Front: What is the Eckart-Young theorem? Back: The truncated SVD A_k is the optimal rank-k approximation in both Frobenius and spectral norm.

Front: How do you compute the Frobenius reconstruction error of a rank-k truncation? Back: The sum of squares of the discarded singular values, sum over i>k of sigma_i^2.

Front: What does a singular value plateau near machine epsilon mean? Back: Those directions are float noise, not structure; truncating there is bug avoidance, not regularization.

Front: What does a smooth power-law spectrum tell you about the right truncation rank? Back: There is no sharp correct rank; choose k by held-out error or a variance budget and report it.

Front: When solving with rank-k truncation, which condition number should you report? Back: The truncated one, sigma_1 / sigma_k, not sigma_1 / sigma_n of the original matrix.

Front: What tolerance is a defensible default for testing sigma_i > 0? Back: tau = sigma_max * max(m, n) * machine epsilon, the standard matrix_rank test.

Front: What does a least-squares residual look like geometrically? Back: It is orthogonal to the column space, i.e. to the span of the leading left singular vectors.

Front: Why does truncation stabilize an ill-posed solve? Back: It stops tiny singular values from contributing amplification factors of order 1/sigma_i.

The truncated pseudoinverse keeps only the first _____ singular triplets, where k is the count of singular values above _____. This is equivalent to solving the _____ problem instead of the literal one. For an underdetermined consistent system, x = A^+ b is the solution of ____. A gap in the spectrum, with sigma_k / sigma{k+1} around _____ or larger, makes the truncation rank unambiguous.