Singular Value Decomposition

Why Singular Value Decomposition Matters

A matrix can contain many entries while most of its action comes from a few important directions. Singular value decomposition, or SVD, separates those directions from their strengths. Keeping the strongest terms gives a small approximation; examining weak terms reveals sensitivity; and using reciprocal strengths provides a principled least-squares solver.

Unlike a real eigenvalue decomposition, an SVD exists for every real rectangular matrix, including nonsquare, singular, and rank-deficient inputs. It writes the transformation as orthogonal input coordinates, nonnegative stretching, and orthogonal output coordinates. This makes it useful for data compression, principal component analysis, inverse problems, and understanding what a linear model can or cannot determine.

We will trace a one-sided Jacobi method, implement it in Python, C++, and Java, and check reconstruction and orthogonality independently. The programs support tall and wide matrices and make numerical truncation explicit. Their least-squares helper returns the minimum-norm solution for the retained numerical model, rather than attempting to divide by discarded singular values.

Core Concepts and Terminology

For $A\in\mathbb{R}^{m\times n}$, let $p=\min(m,n)$. A thin SVD is

$$ A=U\Sigma V^T, \qquad U\in\mathbb{R}^{m\times p}, \quad V\in\mathbb{R}^{n\times p}, \quad\Sigma=\operatorname{diag}(\sigma_1,\ldots,\sigma_p), $$

where $U^TU=V^TV=I_p$ and $\sigma_1\ge\cdots\ge\sigma_p\ge0$. Columns $u_i$ are left singular vectors, and columns $v_i$ are right singular vectors. Equivalently,

$$ A=\sum_{i=1}^{p}\sigma_i u_iv_i^T, \qquad Av_i=\sigma_i u_i, \qquad A^Tu_i=\sigma_i v_i. $$

These formulas explain the geometry: $V^T$ extracts input-direction coordinates, $\Sigma$ scales them, and $U$ assembles the output. Orthogonal factors may include reflections, so describing every factor as a rotation would be too restrictive.

A full SVD extends $U$ and $V$ to square orthogonal matrices and uses a rectangular diagonal $\Sigma$. A compact rank-$r$ SVD keeps only positive singular values and their vectors. The code returns thin factors with $p$ singular values, explicitly setting numerically discarded ones to zero and completing the corresponding vector columns orthogonally.

The exact rank is the number of positive singular values. Numerical rank depends on a chosen cutoff and floating-point accuracy. Our cutoff is $\tau_r\lVert A\rVert_F$, using a separate convergence tolerance $\tau$ for column orthogonality. This convention differs from the common cutoff $\tau_r\sigma_1$; neither is a universal noise model.

The LAPACK SVD overview documents the decomposition and singular-vector relations. Our real-valued lesson uses ordinary transposes; complex SVD uses conjugate transposes.

How Singular Value Decomposition Works, Step by Step

A one-sided Jacobi method mixes pairs of columns until they become mutually orthogonal. Their lengths then supply singular values. The corresponding normalized columns become left singular vectors, while the accumulated column transformations become right singular vectors.

For clarity, first assume $m\ge n$:

1. Validate and scale. Check nonempty rectangular shape, finite entries, finite tolerances in $(0,1)$, and a nonnegative rotation budget. Set $s=\max_{i,j}\lvert A_{ij}\rvert$, using $s=1$ for the zero matrix. Start $B=A/s$, $V=I_n$, and keep the fixed cutoff $\delta=\tau_r\lVert A/s\rVert_F$.
2. Choose a pivot pair. Recompute column lengths $\ell_j=\lVert b_j\rVert_2$. Ignore columns with $\ell_j\le\delta$. Among remaining pairs, find the largest absolute correlation $\lvert b_p^Tb_q\rvert/(\ell_p\ell_q)$, keeping the first on an exact tie. If it is at most $\tau$, the retained columns pass the orthogonality test.
3. Construct a stable plane rotation. Let $h=\max(\ell_p,\ell_q)$, $a=(\ell_p/h)^2$, $d=(\ell_q/h)^2$, and $g=(\ell_p/h)(\ell_q/h)\operatorname{corr}(b_p,b_q)$. These scaled quantities avoid explicitly squaring large column lengths. Set

$$ \zeta=\frac{d-a}{2g}, \qquad t=\frac{\operatorname{sgn}_+(\zeta)}{\lvert\zeta\rvert+\sqrt{1+\zeta^2}}, \qquad c=\frac1{\sqrt{1+t^2}},\quad s_g=ct, $$

where $\operatorname{sgn}_+(0)=1$. Norm evaluations in the code use hypot. A nonrepresentable rotation is reported as an error.

1. Rotate both matrices. Apply the same transformation to columns $p,q$ of $B$ and $V$:

$$ b_p' = cb_p-s_gb_q, \qquad b_q'=s_gb_p+cb_q. $$

Save both old entries before overwriting either. This is right multiplication by $J=\begin{pmatrix}c&s_g\\-s_g&c\end{pmatrix}$ embedded in the selected columns. It preserves $B=(A/s)V$.

1. Repeat and sort. Recheck correlations after every rotation. If the budget is exhausted before convergence, report failure. Sort final column lengths in descending order, moving the associated columns of both factors together.
2. Extract the thin factors. For retained length $\ell_j$, set $\sigma_j=s\ell_j$ and $u_j=b_j/\ell_j$. For a discarded column, report $\sigma_j=0$ and choose a unit vector orthogonal to previously formed left columns. The code selects a coordinate with the largest remaining projector-diagonal estimate and applies two Gram–Schmidt passes to complete this basis.

For a wide matrix, apply the tall algorithm to $A^T$. If $A^T=U_t\Sigma V_t^T$, then $A=V_t\Sigma U_t^T$; exchange the left and right factors when returning.

The LAPACK DGESVJ reference describes a production Jacobi SVD with more extensive scaling, stopping, and sweep machinery. Our maximum-correlation search is chosen for a visible pivot rule, and has a higher search cost than efficient sweep implementations.

flowchart TD
    A["Validate shape, finite data, tolerances, and budget"] --> B{"Wide matrix?"}
    B -->|Yes| C["Transpose for tall computation"]
    B -->|No| D["Scale input; initialize B and right basis"]
    C --> D
    D --> E["Measure column lengths; choose largest active correlation"]
    E --> F{"Correlation within tolerance?"}
    F -->|No| G{"Rotation budget remains?"}
    G -->|No| H["Report nonconvergence"]
    G -->|Yes| I["Rotate pivot columns of B and right basis"]
    I --> E
    F -->|Yes| J["Sort strengths; normalize retained left columns"]
    J --> K["Set discarded strengths to zero; complete left basis"]
    K --> L["Swap returned factors if input was transposed"]
    L --> M["Check reconstruction and both orthogonality identities"]

Worked Example with a Full Trace

Take the rectangular matrix and measurement vector

$$ A=\begin{pmatrix}3&1\\1&3\\0&0\end{pmatrix}, \qquad b=\begin{pmatrix}4\\0\\3\end{pmatrix}. $$

We show intermediate columns in original units. The code scales by $s=3$, which leaves correlations and rotation angles unchanged. Its rank cutoff is $10^{-12}\sqrt{20}/3$ in normalized units, far below either final column length.

Initial state. Let $B_0=A$ and $V_0=I_2$. Both column lengths are $\sqrt{10}$; their dot product is $6$, so their correlation is $3/5$. There is only one pair, and $3/5\gt10^{-12}$, so it must be rotated.

Rotation parameters. Dividing the two lengths by their maximum gives $a=d=1$ and $g=3/5$. Therefore $\zeta=0$, $t=1$, and $c=s_g=1/\sqrt2$. The rotation is

$$ J=\frac1{\sqrt2}\begin{pmatrix}1&1\\-1&1\end{pmatrix}. $$

Apply the rotation once. The first new column is $(b_0-b_1)/\sqrt2=(\sqrt2,-\sqrt2,0)^T$. The second is $(b_0+b_1)/\sqrt2=(2\sqrt2,2\sqrt2,0)^T$. Consequently,

$$ B_1=AJ=\begin{pmatrix}\sqrt2&2\sqrt2\\-\sqrt2&2\sqrt2\\0&0\end{pmatrix}, \qquad V_1=J. $$

The new dot product is $4-4=0$. Their lengths are $2$ and $4$. The next correlation check passes, so the full rotation trace ends after exactly one rotation.

Sort and normalize. Move the length-$4$ column ahead of the length-$2$ column in both $B_1$ and $V_1$. Normalize those columns to obtain

$$ U=\frac1{\sqrt2}\begin{pmatrix}1&1\\1&-1\\0&0\end{pmatrix}, \qquad \Sigma=\begin{pmatrix}4&0\\0&2\end{pmatrix}, \qquad V=\frac1{\sqrt2}\begin{pmatrix}1&1\\1&-1\end{pmatrix}. $$

Both factors have orthonormal columns. Expanding the two rank-one terms verifies every reconstructed entry:

$$ 4u_1v_1^T=\begin{pmatrix}2&2\\2&2\\0&0\end{pmatrix}, \qquad 2u_2v_2^T=\begin{pmatrix}1&-1\\-1&1\\0&0\end{pmatrix}. $$

Their sum is $A$. The strongest term alone gives $A_1$, with $\lVert A-A_1\rVert_F=2$ and retained squared-Frobenius fraction $16/(16+4)=4/5$. The matrix condition number in the Euclidean norm is $4/2=2$.

Least squares. Project the measurements to get $U^Tb=(2\sqrt2,2\sqrt2)^T$. Divide by the positive singular values to obtain $(\sqrt2/2,\sqrt2)^T$, then multiply by $V$. The result is $x=(3/2,-1/2)^T$. Thus $Ax=(4,0,0)^T$ and the residual $Ax-b=(0,0,-3)^T$ has norm $3$. No coefficient choice can fit the third measurement because the third row of $A$ is zero.

For this full-column-rank example, the pseudoinverse is $A^+=\begin{pmatrix}3/8&-1/8&0\\-1/8&3/8&0\end{pmatrix}$. Multiplying it by $b$ independently confirms the solution.

Implementation in Python, C++, and Java

Each panel is a complete program using ordinary nested arrays, an explicit maximum-correlation pivot-pair search, and plane rotations. Run Python 3, use C++14 or newer, or save Java as Main.java. The public SVD routine supports any nonempty finite real rectangular matrix, including wide, rank-deficient, and zero inputs. It preserves the original matrix and returns thin $U$, a descending list of $p$ strengths, thin $V$ rather than $V^T$, and the number of rotations performed.

Two tolerances serve different purposes. tol, default $10^{-12}$, controls normalized correlation among retained columns. rank_tol, also default $10^{-12}$, discards column lengths at most rank_tol times the original Frobenius norm. Once those small strengths are set to zero, the factors describe a numerical approximation $A_{\tau_r}$ to $A$. A small but mathematically positive singular value can therefore be deliberately removed. These thresholds do not classify exact rank or infer a noise distribution.

The least-squares helper uses only retained positive strengths and returns $A_{\tau_r}^+b$: the minimum-norm least-squares solution for that retained model, up to numerical errors. If only exact zeros were removed, this is the minimum-norm solution for original $A$. Discarding a positive strength can change the original problem's minimizer. The LAPACK DGELSD reference describes SVD-based minimum-norm solves and its own effective-rank cutoff convention.

Internal helpers assume compatible factors from a successful SVD call. reconstruct(factors, keep) expects $0\le\text{keep}\le p$. Reported diagnostics measure the largest reconstruction-entry error relative to $\max\lvert A_{ij}\rvert$, the largest entry of either $U^TU-I$ or $V^TV-I$, the rank-one Frobenius error, and the original-system fit residual. Printed zeros are rounded. Scaling and hypot help with range, but ordinary sums and the simple pseudoinverse helper can still overflow on extreme inputs; detected nonrepresentable results are rejected.

import math


def norm(v):
    s = 0.0
    for x in v:
        s = math.hypot(s, x)
    return s


def tall_svd(a, tol, rank_tol, limit):
    m, n = len(a), len(a[0])
    scale = max(abs(x) for row in a for x in row) or 1.0
    b = [[x / scale for x in row] for row in a]
    v = [[float(i == j) for j in range(n)] for i in range(n)]
    reference = norm([x for row in b for x in row])
    cutoff = rank_tol * reference
    for rotations in range(limit + 1):
        lengths = [norm([b[i][j] for i in range(m)]) for j in range(n)]
        best, pivot = 0.0, None
        for p in range(n):
            for q in range(p + 1, n):
                if lengths[p] <= cutoff or lengths[q] <= cutoff:
                    continue
                corr = sum((b[i][p] / lengths[p]) * (b[i][q] / lengths[q])
                           for i in range(m))
                if abs(corr) > best:
                    best, pivot = abs(corr), (p, q, corr)
        if best <= tol:
            break
        if rotations == limit:
            raise RuntimeError("Jacobi SVD did not converge")
        p, q, corr = pivot
        large = max(lengths[p], lengths[q])
        x, y = lengths[p] / large, lengths[q] / large
        g = x * y * corr
        if g == 0.0:
            raise ValueError("rotation is not representable")
        zeta = (y * y - x * x) / (2.0 * g)
        if not math.isfinite(zeta):
            raise ValueError("rotation is not representable")
        t = (1.0 if zeta >= 0.0 else -1.0) / (abs(zeta) + math.hypot(1.0, zeta))
        c = 1.0 / math.hypot(1.0, t)
        s = c * t
        for matrix in (b, v):
            for row in matrix:
                old_p, old_q = row[p], row[q]
                row[p] = c * old_p - s * old_q
                row[q] = s * old_p + c * old_q
    order = sorted(range(n), key=lambda j: -lengths[j])
    u = [[0.0] * n for _ in range(m)]
    right = [[v[i][j] for j in order] for i in range(n)]
    singular = []
    for j, source in enumerate(order):
        length = lengths[source]
        if length > cutoff:
            value = length * scale
            if value == 0.0 or not math.isfinite(value):
                raise ValueError("singular value is not representable")
            singular.append(value)
            for i in range(m):
                u[i][j] = b[i][source] / length
        else:
            singular.append(0.0)
            # Complete a left orthonormal basis for the zero singular value.
            axis = max(range(m), key=lambda i: 1.0 - sum(u[i][k] ** 2 for k in range(j)))
            w = [float(i == axis) for i in range(m)]
            for repeat in range(2):
                for k in range(j):
                    d = sum(u[i][k] * w[i] for i in range(m))
                    for i in range(m):
                        w[i] -= d * u[i][k]
            length_w = norm(w)
            if length_w <= 1e-14:
                raise ValueError("orthogonal completion failed")
            for i in range(m):
                u[i][j] = w[i] / length_w
    return u, singular, right, rotations


def svd(a, tol=1e-12, rank_tol=1e-12, limit=10000):
    m = len(a)
    n = len(a[0]) if m else 0
    if n == 0 or any(len(row) != n for row in a):
        raise ValueError("A must be nonempty and rectangular")
    if any(not math.isfinite(x) for row in a for x in row):
        raise ValueError("A must be finite")
    if any(not math.isfinite(t) or not 0.0 < t < 1.0 for t in (tol, rank_tol)):
        raise ValueError("tolerances must be finite and between 0 and 1")
    if not isinstance(limit, int) or limit < 0:
        raise ValueError("limit must be a nonnegative integer")
    if m >= n:
        return tall_svd(a, tol, rank_tol, limit)
    transposed = [[a[i][j] for i in range(m)] for j in range(n)]
    ut, singular, vt, rotations = tall_svd(transposed, tol, rank_tol, limit)
    return vt, singular, ut, rotations


def least_squares(factors, b):
    u, singular, v, rotations = factors
    if len(b) != len(u) or any(not math.isfinite(x) for x in b):
        raise ValueError("b must have m finite entries")
    z = []
    for j, sigma in enumerate(singular):
        value = sum(u[i][j] * b[i] for i in range(len(u))) / sigma if sigma > 0.0 else 0.0
        if not math.isfinite(value):
            raise ValueError("nonfinite pseudoinverse coordinate")
        z.append(value)
    x = [sum(row[j] * z[j] for j in range(len(z))) for row in v]
    if any(not math.isfinite(value) for value in x):
        raise ValueError("nonfinite solution")
    return x


def reconstruct(factors, keep):
    u, singular, v, rotations = factors
    return [[sum(u[i][k] * singular[k] * v[j][k] for k in range(keep))
             for j in range(len(v))] for i in range(len(u))]


def orthogonality(q):
    n = len(q[0])
    return max(abs(sum(row[i] * row[j] for row in q) - float(i == j))
               for i in range(n) for j in range(n))


def vector_text(x):
    return "[" + ", ".join(f"{(0.0 if v == 0.0 else v):.6f}" for v in x) + "]"


def main():
    a = [[3.0, 1.0], [1.0, 3.0], [0.0, 0.0]]
    b = [4.0, 0.0, 3.0]
    f = svd(a)
    u, singular, v, rotations = f
    full, low = reconstruct(f, len(singular)), reconstruct(f, 1)
    scale = max(abs(x) for row in a for x in row)
    error = max(abs(a[i][j] - full[i][j]) for i in range(3) for j in range(2)) / scale
    low_error = norm([a[i][j] - low[i][j] for i in range(3) for j in range(2)])
    x = least_squares(f, b)
    residual = norm([sum(a[i][j] * x[j] for j in range(2)) - b[i] for i in range(3)])
    orth = max(orthogonality(u), orthogonality(v))
    print("singular values = " + vector_text(singular))
    print(f"rotations = {rotations}")
    print("least-squares x = " + vector_text(x))
    print(f"fit residual norm = {residual:.6f}")
    print(f"rank-one error = {low_error:.6f}")
    print(f"relative reconstruction error = {error:.6f}")
    print(f"orthogonality error = {orth:.6f}")
    good = max(error, orth, abs(low_error - 2.0), abs(residual - 3.0)) < 1e-10
    print("check = " + ("PASS" if good else "FAIL"))


if __name__ == "__main__":
    main()

#include <algorithm>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <numeric>
#include <stdexcept>
#include <vector>

using Vec = std::vector<double>;
using Matrix = std::vector<Vec>;
struct SVD { Matrix u; Vec singular; Matrix v; int rotations; };

double norm(const Vec& v) {
    double s = 0.0; for (double x : v) s = std::hypot(s, x); return s;
}
double column_norm(const Matrix& b, std::size_t j) {
    double s = 0.0; for (const auto& row : b) s = std::hypot(s, row[j]); return s;
}
SVD tall_svd(const Matrix& a, double tol, double rank_tol, int limit) {
    const std::size_t m = a.size(), n = a[0].size();
    double scale = 0.0;
    for (const auto& row : a) for (double x : row) scale = std::max(scale, std::abs(x));
    if (scale == 0.0) scale = 1.0;
    Matrix b = a, v(n, Vec(n, 0.0));
    double reference = 0.0;
    for (auto& row : b) for (auto& x : row) { x /= scale; reference = std::hypot(reference, x); }
    for (std::size_t i = 0; i < n; ++i) v[i][i] = 1.0;
    const double cutoff = rank_tol * reference;
    Vec lengths(n);
    int rotations = 0;
    for (;;) {
        for (std::size_t j = 0; j < n; ++j) lengths[j] = column_norm(b, j);
        double best = 0.0, chosen_corr = 0.0;
        std::size_t p = 0, q = 0;
        for (std::size_t i = 0; i < n; ++i) for (std::size_t j = i + 1; j < n; ++j) {
            if (lengths[i] <= cutoff || lengths[j] <= cutoff) continue;
            double corr = 0.0;
            for (std::size_t k = 0; k < m; ++k) corr += (b[k][i] / lengths[i]) * (b[k][j] / lengths[j]);
            if (std::abs(corr) > best) { best = std::abs(corr); p = i; q = j; chosen_corr = corr; }
        }
        if (best <= tol) break;
        if (rotations == limit) throw std::runtime_error("Jacobi SVD did not converge");
        const double large = std::max(lengths[p], lengths[q]);
        const double x = lengths[p] / large, y = lengths[q] / large, g = x * y * chosen_corr;
        if (g == 0.0) throw std::runtime_error("rotation is not representable");
        const double zeta = (y * y - x * x) / (2.0 * g);
        if (!std::isfinite(zeta)) throw std::runtime_error("rotation is not representable");
        const double t = (zeta >= 0.0 ? 1.0 : -1.0) / (std::abs(zeta) + std::hypot(1.0, zeta));
        const double c = 1.0 / std::hypot(1.0, t), s = c * t;
        for (Matrix* matrix : {&b, &v}) for (auto& row : *matrix) {
            const double old_p = row[p], old_q = row[q];
            row[p] = c * old_p - s * old_q; row[q] = s * old_p + c * old_q;
        }
        ++rotations;
    }
    std::vector<std::size_t> order(n); std::iota(order.begin(), order.end(), 0);
    std::stable_sort(order.begin(), order.end(), [&](std::size_t i, std::size_t j) { return lengths[i] > lengths[j]; });
    Matrix u(m, Vec(n, 0.0)), right(n, Vec(n)); Vec singular(n);
    for (std::size_t i = 0; i < n; ++i) for (std::size_t j = 0; j < n; ++j) right[i][j] = v[i][order[j]];
    for (std::size_t j = 0; j < n; ++j) {
        const std::size_t source = order[j]; const double length = lengths[source];
        if (length > cutoff) {
            singular[j] = length * scale;
            if (singular[j] == 0.0 || !std::isfinite(singular[j])) throw std::runtime_error("singular value is not representable");
            for (std::size_t i = 0; i < m; ++i) u[i][j] = b[i][source] / length;
        } else {
            std::size_t axis = 0; double best = -1.0;
            for (std::size_t i = 0; i < m; ++i) {
                double score = 1.0; for (std::size_t k = 0; k < j; ++k) score -= u[i][k] * u[i][k];
                if (score > best) { best = score; axis = i; }
            }
            Vec w(m, 0.0); w[axis] = 1.0;
            for (int repeat = 0; repeat < 2; ++repeat) for (std::size_t k = 0; k < j; ++k) {
                double d = 0.0; for (std::size_t i = 0; i < m; ++i) d += u[i][k] * w[i];
                for (std::size_t i = 0; i < m; ++i) w[i] -= d * u[i][k];
            }
            const double length_w = norm(w);
            if (length_w <= 1e-14) throw std::runtime_error("orthogonal completion failed");
            for (std::size_t i = 0; i < m; ++i) u[i][j] = w[i] / length_w;
        }
    }
    return {u, singular, right, rotations};
}
SVD svd(const Matrix& a, double tol = 1e-12, double rank_tol = 1e-12, int limit = 10000) {
    const std::size_t m = a.size(), n = m ? a[0].size() : 0;
    if (n == 0) throw std::invalid_argument("A must be nonempty");
    for (const auto& row : a) {
        if (row.size() != n) throw std::invalid_argument("A must be rectangular");
        for (double x : row) if (!std::isfinite(x)) throw std::invalid_argument("A must be finite");
    }
    if (!std::isfinite(tol) || !std::isfinite(rank_tol) || !(tol > 0 && tol < 1) || !(rank_tol > 0 && rank_tol < 1) || limit < 0)
        throw std::invalid_argument("invalid tolerances or rotation limit");
    if (m >= n) return tall_svd(a, tol, rank_tol, limit);
    Matrix at(n, Vec(m));
    for (std::size_t i = 0; i < m; ++i) for (std::size_t j = 0; j < n; ++j) at[j][i] = a[i][j];
    auto result = tall_svd(at, tol, rank_tol, limit);
    std::swap(result.u, result.v); return result;
}
Vec least_squares(const SVD& f, const Vec& b) {
    if (b.size() != f.u.size()) throw std::invalid_argument("b must have m entries");
    for (double x : b) if (!std::isfinite(x)) throw std::invalid_argument("b must be finite");
    Vec z(f.singular.size(), 0.0);
    for (std::size_t j = 0; j < z.size(); ++j) if (f.singular[j] > 0.0) {
        for (std::size_t i = 0; i < b.size(); ++i) z[j] += f.u[i][j] * b[i];
        z[j] /= f.singular[j];
        if (!std::isfinite(z[j])) throw std::runtime_error("nonfinite pseudoinverse coordinate");
    }
    Vec x(f.v.size(), 0.0);
    for (std::size_t i = 0; i < x.size(); ++i) {
        for (std::size_t j = 0; j < z.size(); ++j) x[i] += f.v[i][j] * z[j];
        if (!std::isfinite(x[i])) throw std::runtime_error("nonfinite solution");
    }
    return x;
}
Matrix reconstruct(const SVD& f, std::size_t keep) {
    Matrix a(f.u.size(), Vec(f.v.size(), 0.0));
    for (std::size_t i = 0; i < a.size(); ++i) for (std::size_t j = 0; j < a[0].size(); ++j)
        for (std::size_t k = 0; k < keep; ++k) a[i][j] += f.u[i][k] * f.singular[k] * f.v[j][k];
    return a;
}
double orthogonality(const Matrix& q) {
    double error = 0.0;
    for (std::size_t i = 0; i < q[0].size(); ++i) for (std::size_t j = 0; j < q[0].size(); ++j) {
        double s = 0.0; for (const auto& row : q) s += row[i] * row[j];
        error = std::max(error, std::abs(s - (i == j ? 1.0 : 0.0)));
    }
    return error;
}
void print_vector(const Vec& x) {
    std::cout << '['; for (std::size_t i = 0; i < x.size(); ++i) {
        if (i) std::cout << ", ";
        std::cout << (x[i] == 0.0 ? 0.0 : x[i]);
    } std::cout << ']';
}
int main() {
    const Matrix a = { {3.0, 1.0}, {1.0, 3.0}, {0.0, 0.0} }; const Vec b = {4.0, 0.0, 3.0};
    const auto f = svd(a); const auto full = reconstruct(f, f.singular.size()), low = reconstruct(f, 1);
    double error = 0.0; Vec differences;
    for (std::size_t i = 0; i < 3; ++i) for (std::size_t j = 0; j < 2; ++j) {
        error = std::max(error, std::abs(a[i][j] - full[i][j]) / 3.0); differences.push_back(a[i][j] - low[i][j]);
    }
    const double low_error = norm(differences); const auto x = least_squares(f, b); Vec residual;
    for (std::size_t i = 0; i < 3; ++i) residual.push_back(a[i][0] * x[0] + a[i][1] * x[1] - b[i]);
    const double re = norm(residual), orth = std::max(orthogonality(f.u), orthogonality(f.v));
    std::cout << std::fixed << std::setprecision(6) << "singular values = "; print_vector(f.singular);
    std::cout << "\nrotations = " << f.rotations << "\nleast-squares x = "; print_vector(x);
    std::cout << "\nfit residual norm = " << re << "\nrank-one error = " << low_error
              << "\nrelative reconstruction error = " << error << "\northogonality error = " << orth
              << "\ncheck = " << (std::max({error, orth, std::abs(low_error - 2), std::abs(re - 3)}) < 1e-10 ? "PASS" : "FAIL") << '\n';
}

import java.util.Arrays;
import java.util.Locale;

public class Main {
    static class SVD {
        final double[][] u, v; final double[] singular; final int rotations;
        SVD(double[][] u,double[] singular,double[][] v,int rotations) {
            this.u=u;this.singular=singular;this.v=v;this.rotations=rotations;
        }
    }
    static double norm(double[] v) {
        double s=0.0;for(double x:v)s=Math.hypot(s,x);return s;
    }
    static double columnNorm(double[][] b,int j) {
        double s=0.0;for(double[]row:b)s=Math.hypot(s,row[j]);return s;
    }
    static SVD tallSVD(double[][] a,double tol,double rankTol,int limit) {
        int m=a.length,n=a[0].length;double scale=0.0;
        for(double[]row:a)for(double x:row)scale=Math.max(scale,Math.abs(x));if(scale==0.0)scale=1.0;
        double[][] b=new double[m][n],v=new double[n][n];double reference=0.0;
        for(int i=0;i<m;++i)for(int j=0;j<n;++j){b[i][j]=a[i][j]/scale;reference=Math.hypot(reference,b[i][j]);}
        for(int i=0;i<n;++i)v[i][i]=1.0;
        double cutoff=rankTol*reference;double[] lengths=new double[n];int rotations=0;
        for(;;) {
            for(int j=0;j<n;++j)lengths[j]=columnNorm(b,j);
            double best=0.0,chosenCorr=0.0;int p=0,q=0;
            for(int i=0;i<n;++i)for(int j=i+1;j<n;++j) {
                if(lengths[i]<=cutoff||lengths[j]<=cutoff)continue;
                double corr=0.0;for(int k=0;k<m;++k)corr+=(b[k][i]/lengths[i])*(b[k][j]/lengths[j]);
                if(Math.abs(corr)>best){best=Math.abs(corr);p=i;q=j;chosenCorr=corr;}
            }
            if(best<=tol)break;
            if(rotations==limit)throw new IllegalStateException("Jacobi SVD did not converge");
            double large=Math.max(lengths[p],lengths[q]),x=lengths[p]/large,y=lengths[q]/large,g=x*y*chosenCorr;
            if(g==0.0)throw new IllegalArgumentException("rotation is not representable");
            double zeta=(y*y-x*x)/(2.0*g);
            if(!Double.isFinite(zeta))throw new IllegalArgumentException("rotation is not representable");
            double t=(zeta>=0.0?1.0:-1.0)/(Math.abs(zeta)+Math.hypot(1.0,zeta));
            double c=1.0/Math.hypot(1.0,t),s=c*t;
            for(double[][]matrix:new double[][][]{b,v})for(double[]row:matrix) {
                double oldP=row[p],oldQ=row[q];row[p]=c*oldP-s*oldQ;row[q]=s*oldP+c*oldQ;
            }
            ++rotations;
        }
        Integer[] order=new Integer[n];for(int j=0;j<n;++j)order[j]=j;
        Arrays.sort(order,(i,j)->Double.compare(lengths[j],lengths[i]));
        double[][] u=new double[m][n],right=new double[n][n];double[] singular=new double[n];
        for(int i=0;i<n;++i)for(int j=0;j<n;++j)right[i][j]=v[i][order[j]];
        for(int j=0;j<n;++j) {
            int source=order[j];double length=lengths[source];
            if(length>cutoff) {
                singular[j]=length*scale;
                if(singular[j]==0.0||!Double.isFinite(singular[j]))throw new IllegalArgumentException("singular value is not representable");
                for(int i=0;i<m;++i)u[i][j]=b[i][source]/length;
            } else {
                int axis=0;double best=-1.0;
                for(int i=0;i<m;++i) {
                    double score=1.0;for(int k=0;k<j;++k)score-=u[i][k]*u[i][k];
                    if(score>best){best=score;axis=i;}
                }
                double[]w=new double[m];w[axis]=1.0;
                for(int repeat=0;repeat<2;++repeat)for(int k=0;k<j;++k) {
                    double d=0.0;for(int i=0;i<m;++i)d+=u[i][k]*w[i];
                    for(int i=0;i<m;++i)w[i]-=d*u[i][k];
                }
                double lengthW=norm(w);
                if(lengthW<=1e-14)throw new IllegalArgumentException("orthogonal completion failed");
                for(int i=0;i<m;++i)u[i][j]=w[i]/lengthW;
            }
        }
        return new SVD(u,singular,right,rotations);
    }
    static SVD svd(double[][] a,double tol,double rankTol,int limit) {
        int m=a.length,n=m>0&&a[0]!=null?a[0].length:0;
        if(n==0)throw new IllegalArgumentException("A must be nonempty");
        for(double[]row:a) {
            if(row==null||row.length!=n)throw new IllegalArgumentException("A must be rectangular");
            for(double x:row)if(!Double.isFinite(x))throw new IllegalArgumentException("A must be finite");
        }
        if(!Double.isFinite(tol)||!Double.isFinite(rankTol)||!(tol>0&&tol<1)||!(rankTol>0&&rankTol<1)||limit<0)
            throw new IllegalArgumentException("invalid tolerances or rotation limit");
        if(m>=n)return tallSVD(a,tol,rankTol,limit);
        double[][]at=new double[n][m];for(int i=0;i<m;++i)for(int j=0;j<n;++j)at[j][i]=a[i][j];
        SVD result=tallSVD(at,tol,rankTol,limit);return new SVD(result.v,result.singular,result.u,result.rotations);
    }
    static double[] leastSquares(SVD f,double[] b) {
        if(b.length!=f.u.length)throw new IllegalArgumentException("b must have m entries");
        for(double x:b)if(!Double.isFinite(x))throw new IllegalArgumentException("b must be finite");
        double[]z=new double[f.singular.length];
        for(int j=0;j<z.length;++j)if(f.singular[j]>0.0) {
            for(int i=0;i<b.length;++i)z[j]+=f.u[i][j]*b[i];z[j]/=f.singular[j];
            if(!Double.isFinite(z[j]))throw new IllegalArgumentException("nonfinite pseudoinverse coordinate");
        }
        double[]x=new double[f.v.length];
        for(int i=0;i<x.length;++i) {
            for(int j=0;j<z.length;++j)x[i]+=f.v[i][j]*z[j];
            if(!Double.isFinite(x[i]))throw new IllegalArgumentException("nonfinite solution");
        }
        return x;
    }
    static double[][] reconstruct(SVD f,int keep) {
        double[][]a=new double[f.u.length][f.v.length];
        for(int i=0;i<a.length;++i)for(int j=0;j<a[0].length;++j)for(int k=0;k<keep;++k)
            a[i][j]+=f.u[i][k]*f.singular[k]*f.v[j][k];return a;
    }
    static double orthogonality(double[][]q) {
        double error=0.0;
        for(int i=0;i<q[0].length;++i)for(int j=0;j<q[0].length;++j) {
            double s=0.0;for(double[]row:q)s+=row[i]*row[j];error=Math.max(error,Math.abs(s-(i==j?1.0:0.0)));
        }
        return error;
    }
    static String vectorText(double[]x) {
        StringBuilder s=new StringBuilder("[");for(int i=0;i<x.length;++i) {
            if(i>0)s.append(", ");s.append(String.format(Locale.ROOT,"%.6f",x[i]==0.0?0.0:x[i]));
        }return s.append(']').toString();
    }
    public static void main(String[]args) {
        double[][]a={ {3.0,1.0},{1.0,3.0},{0.0,0.0} };double[]b={4.0,0.0,3.0};
        SVD f=svd(a,1e-12,1e-12,10000);double[][]full=reconstruct(f,f.singular.length),low=reconstruct(f,1);
        double error=0.0;double[]differences=new double[6];int pos=0;
        for(int i=0;i<3;++i)for(int j=0;j<2;++j) {
            error=Math.max(error,Math.abs(a[i][j]-full[i][j])/3.0);differences[pos++]=a[i][j]-low[i][j];
        }
        double lowError=norm(differences);double[]x=leastSquares(f,b),residual=new double[3];
        for(int i=0;i<3;++i)residual[i]=a[i][0]*x[0]+a[i][1]*x[1]-b[i];
        double re=norm(residual),orth=Math.max(orthogonality(f.u),orthogonality(f.v));
        System.out.println("singular values = "+vectorText(f.singular));
        System.out.println("rotations = "+f.rotations);System.out.println("least-squares x = "+vectorText(x));
        System.out.printf(Locale.ROOT,"fit residual norm = %.6f%nrank-one error = %.6f%nrelative reconstruction error = %.6f%northogonality error = %.6f%n",re,lowError,error,orth);
        boolean good=Math.max(Math.max(error,orth),Math.max(Math.abs(lowError-2),Math.abs(re-3)))<1e-10;
        System.out.println("check = "+(good?"PASS":"FAIL"));
    }
}

All three programs print the same eight lines:

• singular values = [4.000000, 2.000000]
• rotations = 1
• least-squares x = [1.500000, -0.500000]
• fit residual norm = 3.000000
• rank-one error = 2.000000
• relative reconstruction error = 0.000000
• orthogonality error = 0.000000
• check = PASS

Complexity, Correctness, and Invariants

Let $M=\max(m,n)$ and $p=\min(m,n)$. Each global pivot search checks up to $p(p-1)/2$ correlations of length $M$, costing $O(Mp^2)$. A rotation itself updates two data columns and two basis columns in $O(M+p)$ time. After $K$ rotations, total work is $O((K+1)Mp^2)$, including correlation checks and $O(Mp^2)$ basis-completion work. This search-based teaching implementation is slower than efficient cyclic or blocked Jacobi methods.

There is no fixed successful rotation count for every input and tolerance. An already orthogonal input can finish without rotations, while a zero budget rejects a nonorthogonal active pair. With a finite rotation budget, the routine returns a converged retained factorization or reports failure. Typical cost depends on the matrix, requested accuracy, and pivot sequence; a single universal average-case iteration count would be misleading.

Thin factors and working arrays use $O(Mp+p^2)=O((m+n)p)$ storage. A pseudoinverse solve with stored factors costs $O((m+n)p)$ time and $O(n+p)$ working/output space. Reconstructing all entries from $k$ terms takes $O(mnk)$ time; storing just those truncated factors takes $O((m+n)k)$ space. Full square orthogonal factors would require extra storage and are not built here.

The main Jacobi invariant is $B=(A/s)V$ with $V^TV=I$. Initially it holds because $V=I$. If both matrices are right-multiplied by the same orthogonal plane rotation $J$, then $BJ=(A/s)(VJ)$, so the relation and the Frobenius norm are preserved.

For a chosen pair, its two-column Gram block, after a harmless common scaling, is $\begin{pmatrix}a&g\\g&d\end{pmatrix}$. The formula for $t$ makes the off-diagonal entry of $J^T\begin{pmatrix}a&g\\g&d\end{pmatrix}J$ zero in exact arithmetic. The code computes the necessary scalar quantities directly from the two columns rather than storing the whole Gram matrix.

If all nonzero columns are orthogonal, their normalized versions form $U$ and their lengths form $\Sigma/s$. Thus $B=U(\Sigma/s)$ and $A=U\Sigma V^T$. Sorting matched columns preserves that product. Exact zero columns can use any left orthonormal completion because their strengths are zero. In floating-point arithmetic, retained columns are only approximately orthogonal; numerical truncation and rounding additionally perturb reconstruction.

For least squares, orthogonal coordinates separate the model into independent scalar problems. Positive-strength directions use $z_i=(u_i^Tb)/\sigma_i$, while zero-strength right directions do not affect the fit. Setting the latter coordinates to zero minimizes the solution norm. This explains the pseudoinverse and why it avoids arbitrary division by zero.

For exact SVD terms sorted by strength, the Eckart–Young–Mirsky theorem gives the best rank-at-most-$k$ approximation:

$$ A_k=\sum_{i=1}^{k}\sigma_i u_iv_i^T, \qquad \lVert A-A_k\rVert_2=\sigma_{k+1}, \qquad \lVert A-A_k\rVert_F=\sqrt{\sum_{i=k+1}^{p}\sigma_i^2} \quad(k\lt p). $$

The Frobenius error follows because the rank-one terms are orthonormal under the matrix Frobenius inner product; optimality is the additional theorem. Numerical factors achieve these statements approximately, and a prior rank cutoff changes the retained matrix being approximated.

Variants, Pitfalls, and Optimizations

Do not turn a theoretical identity into an unstable shortcut. In exact arithmetic, $A^TA$ has eigenvalues $\sigma_i^2$. But for $A=\begin{pmatrix}1&1\\\varepsilon&0\\0&\varepsilon\end{pmatrix}$ with $\varepsilon=10^{-10}$, the exact Gram matrix is $\begin{pmatrix}1+\varepsilon^2&1\\1&1+\varepsilon^2\end{pmatrix}$. Ordinary double precision rounds $1+10^{-20}$ to $1$, losing the weak singular value $10^{-10}$. A direct column rotation retains the small entries instead of adding their squares to a unit diagonal.

Sorting and signs are coupled choices. Sorting singular values alone breaks the decomposition. Move the matching left and right columns too. Flipping both $u_i$ and $v_i$ preserves their rank-one term; flipping only one does not. Within a repeated-singular-value subspace, different orthonormal bases can all be correct. Compare reconstruction, orthogonality, and subspaces instead of demanding identical signed vectors.

Zero strength and small strength need different interpretations. Exact zero directions carry no information, while small positive directions may be meaningful or noise-sensitive. The code's cutoff creates a truncated model and is intentionally visible. Inverting a tiny retained strength amplifies measurement error; discarding it regularizes the inverse problem by changing the model. A rank-deficient matrix's overall invertible-system condition number is infinite; a ratio restricted to retained strengths measures only that retained subspace.

Thin factors do not always list the entire nullspace. The $p$ thin columns suffice for reconstruction. If $n\gt m$, they do not include all $n-m$ additional right-nullspace directions of a full SVD. Minimum-norm solving still works: omitted right directions receive zero coefficient. Similarly, tall thin $U$ need not list every left-nullspace direction.

Choose the algorithm for the needed output. Bidiagonal reduction, divide-and-conquer, and optimized Jacobi solvers are established dense approaches. A QR reduction can help for very tall matrices. Lanczos bidiagonalization and randomized methods target a few singular triplets for large sparse or implicit matrices. Computing only singular values or only a truncated basis can save substantial work and storage.

Interpret applications carefully. For PCA, center the data first when covariance variation is the goal; singular values of centered data relate to variances after the chosen sample normalization. Low-rank compression can suppress useful detail as well as noise. Retained Frobenius energy is a useful numeric criterion, not a guarantee of perceptual or predictive quality. Reuse factors for multiple right-hand sides while the matrix is fixed.

Practice Problems and Real-World Use

Start with a known orthonormal pair of factor matrices and a chosen diagonal spectrum, then construct $A=U\Sigma V^T$. This gives an independent singular-value oracle for tall, wide, rank-deficient, repeated-strength, and zero cases. Check returned shapes, nonnegative descending strengths, orthogonality, reconstruction, and the effects of the declared rank cutoff.

For minimum-norm solving, test both the fit and the absence of unnecessary nullspace components. A wide or rank-deficient system may have infinitely many equally good fits; agreement with the data alone does not prove minimum norm. Compare low-rank errors with the discarded singular-value tail, and repeat well-scaled cases after uniform scaling.