Notebooks · lib/svd

The singular value decomposition, by one-sided Jacobi

A = U S V^T, with U and V orthonormal and S diagonal and non-negative. This is the decomposition that answers what a matrix does to the unit sphere, what its rank is when the answer is not a whole number, and which directions of its input it barely notices.

The temptation is to get it from eig_sym(A^T A): the eigenvectors are V and the eigenvalues are the squares of the singular values. That works, and it is wrong for the small ones, because squaring the matrix squares its condition number and each singular value comes back with about half the digits it should have. The last test measures it on a matrix built to have singular values 1, 10⁻³, 10⁻⁶, 10⁻⁹, 1e-12. Relative error in each, one-sided Jacobi against the squared route:

1       1.1e-15     1.2e-15
1e-3    7.4e-15     5.1e-12
1e-6    8.5e-12     2.0e-6
1e-9    3.0e-8      0.34
1e-12   1.8e-5      1.0

The squared route has lost a third of the fourth value and the whole of the fifth, which it returns as zero. Some of what Jacobi loses further down that column is the test matrix rather than the method: writing Q1 diag(s) Q2^T into double precision already perturbs a singular value of 1e-12 relative to a largest of 1.

One-sided Jacobi does it differently. Take pairs of columns of A and rotate each pair until the two are orthogonal, which is a 2 by 2 eigenvalue problem with a closed form. Sweep over all pairs, repeat, and the columns converge to mutually orthogonal vectors. Their lengths are the singular values and their directions are U; the accumulated rotations are V. Nothing is ever squared except inside a single 2 by 2, so the small singular values keep their digits.

A must have at least as many rows as columns. For a wide matrix, decompose its transpose and swap U and V.

Used with import svd.

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

Shift+Enter in a cell runs it and moves to the next. Running a cell runs every cell above it too, as one program.

[1]lines 43–310

from schur import schur
from linalg import test_matrix

SVD_SWEEPS = 60
SVD_TOL = 1.0e-15

#| Decompose `A`, which is `m` by `n` with `m >= n`. `A` is overwritten with
#| `U S`, so its columns come back scaled by the singular values; `s` gets the
#| singular values and `V` the right vectors, in whatever order the sweeps left
#| them. Returns how many sweeps it took.
def svd_jacobi(A, s, V, m, n) {
    for i to n {
        for j to n {
            V[i, j] = 0.0
        }
        V[i, i] = 1.0
    }
    sweeps = 0
    for sweep to SVD_SWEEPS {
        sweeps = sweeps + 1
        off = 0.0
        for p to n - 1 {
            for qq to n - 1 - p {
                q = p + 1 + qq
                alpha = 0.0
                beta = 0.0
                gamma = 0.0
                for i to m {
                    alpha = alpha + A[i, p] * A[i, p]
                    beta = beta + A[i, q] * A[i, q]
                    gamma = gamma + A[i, p] * A[i, q]
                }
                scale = sqrt(alpha * beta)
                if scale > 0.0 {
                    off = max(off, abs(gamma) / scale)
                    if abs(gamma) > SVD_TOL * scale {
                        # the rotation that makes these two columns orthogonal
                        zeta = (beta - alpha) / (2.0 * gamma)
                        t = 1.0 / (abs(zeta) + sqrt(1.0 + zeta * zeta))
                        if zeta < 0.0 {
                            t = -t
                        }
                        c = 1.0 / sqrt(1.0 + t * t)
                        sn = c * t
                        for i to m {
                            ap = A[i, p]
                            aq = A[i, q]
                            A[i, p] = c * ap - sn * aq
                            A[i, q] = sn * ap + c * aq
                        }
                        for i to n {
                            vp = V[i, p]
                            vq = V[i, q]
                            V[i, p] = c * vp - sn * vq
                            V[i, q] = sn * vp + c * vq
                        }
                    }
                }
            }
        }
        if off < SVD_TOL {
            break
        }
    }
    for j to n {
        total = 0.0
        for i to m {
            total = total + A[i, j] * A[i, j]
        }
        s[j] = sqrt(total)
    }
    return 1.0 * sweeps
}

#| The full decomposition, sorted: `U` is `m` by `n`, `s` descends, and
#| `A = U diag(s) V^T`. `A` is left alone.
def svd(A, U, s, V, m, n) {
    work = zeros(m, n)
    for i to m {
        for j to n {
            work[i, j] = A[i, j]
        }
    }
    junk = svd_jacobi(work, s, V, m, n)
    # sort descending, carrying the columns of both factors along
    for a to n {
        best = a
        for b to n - a {
            j = a + b
            if s[j] > s[best] {
                best = j
            }
        }
        if best != a {
            tmp = s[a]
            s[a] = s[best]
            s[best] = tmp
            for i to m {
                tw = work[i, a]
                work[i, a] = work[i, best]
                work[i, best] = tw
            }
            for i to n {
                tv = V[i, a]
                V[i, a] = V[i, best]
                V[i, best] = tv
            }
        }
    }
    # each column over its singular value; a zero one gets a unit vector
    U[:, :] = work / s
    for j to n {
        if not (s[j] > 0.0) {
            U[:, j] = 0.0
            U[j, j] = 1.0
        }
    }
    return 0.0
}

#| The largest entry of `A - U diag(s) V^T`.
def svd_residual(A: Arr, U: Arr, s: Arr, V: Arr, m, n) {
    return max(abs((U[0:m, 0:n] * s[0:n]) @ transpose(V[0:n, 0:n]) - A[0:m, 0:n]))
}

#| How far the `n` columns of an `m` by `n` matrix are from orthonormal.
def svd_orthonormal(M: Arr, m, n) {
    C = M[0:m, 0:n]
    return max(abs(transpose(C) @ C - eye(n)))
}

#| A matrix with the singular values you ask for: `Q1 diag(want) Q2^T`, with the
#| two orthogonal factors taken from the Schur vectors of two unrelated
#| matrices, so nothing about the construction lines up with the axes.
def svd_with_values(want: Arr, n, seed) {
    T = zeros(n, n)
    Q1 = zeros(n, n)
    Q2 = zeros(n, n)
    junk = schur(test_matrix(n, seed), T, Q1)
    junk = schur(test_matrix(n, seed + 1000.0), T, Q2)
    return (Q1 * want) @ transpose(Q2)
}

test "a diagonal matrix has its own entries as singular values" {
    n = 4
    A = diag([3.0, -5.0, 0.5, 2.0])
    U = zeros(n, n)
    s = zeros(n)
    V = zeros(n, n)
    junk = svd(A, U, s, V, n, n)
    assert_near(s[0], 5.0, 1.0e-13)
    assert_near(s[1], 3.0, 1.0e-13)
    assert_near(s[2], 2.0, 1.0e-13)
    assert_near(s[3], 0.5, 1.0e-13)
    assert_lt(svd_residual(A, U, s, V, n, n), 1.0e-13)
}

test "U S V^T puts a square matrix back, with both factors orthonormal" {
    worst_res = 0.0
    worst_u = 0.0
    worst_v = 0.0
    ncases = 10
    for si to ncases {
        n = 2 + mod(si, 7.0)
        A = test_matrix(n, 17.0 + 31.0 * si)
        U = zeros(n, n)
        s = zeros(n)
        V = zeros(n, n)
        junk = svd(A, U, s, V, n, n)
        worst_res = max(worst_res, svd_residual(A, U, s, V, n, n))
        worst_u = max(worst_u, svd_orthonormal(U, n, n))
        worst_v = max(worst_v, svd_orthonormal(V, n, n))
        for j to n - 1 {
            assert_lt(-1.0e-15, s[j] - s[j + 1])
        }
    }
    assert_lt(worst_res, 1.0e-13)
    assert_lt(worst_u, 1.0e-13)
    assert_lt(worst_v, 1.0e-13)
}

test "a tall matrix works too, and its singular values are the square roots of eig(A^T A)" {
    m = 9
    n = 4
    A = zeros(m, n)
    B = test_matrix(m, 71.0)
    for i to m {
        for j to n {
            A[i, j] = B[i, j]
        }
    }
    U = zeros(m, n)
    s = zeros(n)
    V = zeros(n, n)
    junk = svd(A, U, s, V, m, n)
    assert_lt(svd_residual(A, U, s, V, m, n), 1.0e-13)
    assert_lt(svd_orthonormal(U, m, n), 1.0e-13)
    G = transpose(A) @ A
    vals = zeros(n)
    vecs = zeros(n, n)
    eig_sym(G, vals, vecs)
    # eig_sym returns them ascending, svd descending
    for j to n {
        assert_near(s[j], sqrt(vals[n - 1 - j]), 1.0e-11)
    }
}

test "the small singular values survive, where going through A^T A loses them" {
    n = 5
    want = zeros(n)
    want[0] = 1.0
    want[1] = 1.0e-3
    want[2] = 1.0e-6
    want[3] = 1.0e-9
    want[4] = 1.0e-12
    A = svd_with_values(want, n, 23.0)
    U = zeros(n, n)
    s = zeros(n)
    V = zeros(n, n)
    junk = svd(A, U, s, V, n, n)
    # the squared route: eigenvalues of A^T A, whose condition number is squared
    G = transpose(A) @ A
    vals = zeros(n)
    vecs = zeros(n, n)
    eig_sym(G, vals, vecs)
    jac = zeros(n)
    sq = zeros(n)
    for j to n {
        jac[j] = abs(s[j] - want[j]) / want[j]
        sq[j] = abs(sqrt(max(vals[n - 1 - j], 0.0)) - want[j]) / want[j]
    }
    assert_lt(jac[0], 1.0e-14)
    assert_lt(jac[1], 1.0e-13)
    assert_lt(jac[2], 1.0e-10)
    assert_lt(jac[3], 1.0e-6)
    assert_lt(jac[4], 1.0e-3)
    # the squared route has thrown away the fourth and the fifth entirely
    assert_lt(0.1, sq[3])
    assert_lt(0.5, sq[4])
    # and is worse at every level below the largest
    for j to n - 1 {
        assert_lt(jac[j + 1], sq[j + 1])
    }
}

test "a rank-deficient matrix reports the rank it has" {
    n = 6
    want = zeros(n)
    want[0] = 4.0
    want[1] = 2.0
    want[2] = 0.5
    A = svd_with_values(want, n, 137.0)
    U = zeros(n, n)
    s = zeros(n)
    V = zeros(n, n)
    junk = svd(A, U, s, V, n, n)
    assert_near(s[0], 4.0, 1.0e-12)
    assert_near(s[1], 2.0, 1.0e-12)
    assert_near(s[2], 0.5, 1.0e-12)
    rank = 3
    for j to n - rank {
        assert_lt(s[rank + j], 1.0e-14)
    }
    assert_lt(svd_residual(A, U, s, V, n, n), 1.0e-13)
}

# SPDX-License-Identifier: MIT OR Apache-2.0
# SPDX-FileCopyrightText: 2026 The Viv authors

from schur import schur
from linalg import test_matrix

SVD_SWEEPS = 60
SVD_TOL = 1.0e-15

#| Decompose `A`, which is `m` by `n` with `m >= n`. `A` is overwritten with
#| `U S`, so its columns come back scaled by the singular values; `s` gets the
#| singular values and `V` the right vectors, in whatever order the sweeps left
#| them. Returns how many sweeps it took.
def svd_jacobi(A, s, V, m, n) {
    for i to n {
        for j to n {
            V[i, j] = 0.0
        }
        V[i, i] = 1.0
    }
    sweeps = 0
    for sweep to SVD_SWEEPS {
        sweeps = sweeps + 1
        off = 0.0
        for p to n - 1 {
            for qq to n - 1 - p {
                q = p + 1 + qq
                alpha = 0.0
                beta = 0.0
                gamma = 0.0
                for i to m {
                    alpha = alpha + A[i, p] * A[i, p]
                    beta = beta + A[i, q] * A[i, q]
                    gamma = gamma + A[i, p] * A[i, q]
                }
                scale = sqrt(alpha * beta)
                if scale > 0.0 {
                    off = max(off, abs(gamma) / scale)
                    if abs(gamma) > SVD_TOL * scale {
                        # the rotation that makes these two columns orthogonal
                        zeta = (beta - alpha) / (2.0 * gamma)
                        t = 1.0 / (abs(zeta) + sqrt(1.0 + zeta * zeta))
                        if zeta < 0.0 {
                            t = -t
                        }
                        c = 1.0 / sqrt(1.0 + t * t)
                        sn = c * t
                        for i to m {
                            ap = A[i, p]
                            aq = A[i, q]
                            A[i, p] = c * ap - sn * aq
                            A[i, q] = sn * ap + c * aq
                        }
                        for i to n {
                            vp = V[i, p]
                            vq = V[i, q]
                            V[i, p] = c * vp - sn * vq
                            V[i, q] = sn * vp + c * vq
                        }
                    }
                }
            }
        }
        if off < SVD_TOL {
            break
        }
    }
    for j to n {
        total = 0.0
        for i to m {
            total = total + A[i, j] * A[i, j]
        }
        s[j] = sqrt(total)
    }
    return 1.0 * sweeps
}

#| The full decomposition, sorted: `U` is `m` by `n`, `s` descends, and
#| `A = U diag(s) V^T`. `A` is left alone.
def svd(A, U, s, V, m, n) {
    work = zeros(m, n)
    for i to m {
        for j to n {
            work[i, j] = A[i, j]
        }
    }
    junk = svd_jacobi(work, s, V, m, n)
    # sort descending, carrying the columns of both factors along
    for a to n {
        best = a
        for b to n - a {
            j = a + b
            if s[j] > s[best] {
                best = j
            }
        }
        if best != a {
            tmp = s[a]
            s[a] = s[best]
            s[best] = tmp
            for i to m {
                tw = work[i, a]
                work[i, a] = work[i, best]
                work[i, best] = tw
            }
            for i to n {
                tv = V[i, a]
                V[i, a] = V[i, best]
                V[i, best] = tv
            }
        }
    }
    # each column over its singular value; a zero one gets a unit vector
    U[:, :] = work / s
    for j to n {
        if not (s[j] > 0.0) {
            U[:, j] = 0.0
            U[j, j] = 1.0
        }
    }
    return 0.0
}

#| The largest entry of `A - U diag(s) V^T`.
def svd_residual(A: Arr, U: Arr, s: Arr, V: Arr, m, n) {
    return max(abs((U[0:m, 0:n] * s[0:n]) @ transpose(V[0:n, 0:n]) - A[0:m, 0:n]))
}

#| How far the `n` columns of an `m` by `n` matrix are from orthonormal.
def svd_orthonormal(M: Arr, m, n) {
    C = M[0:m, 0:n]
    return max(abs(transpose(C) @ C - eye(n)))
}

#| A matrix with the singular values you ask for: `Q1 diag(want) Q2^T`, with the
#| two orthogonal factors taken from the Schur vectors of two unrelated
#| matrices, so nothing about the construction lines up with the axes.
def svd_with_values(want: Arr, n, seed) {
    T = zeros(n, n)
    Q1 = zeros(n, n)
    Q2 = zeros(n, n)
    junk = schur(test_matrix(n, seed), T, Q1)
    junk = schur(test_matrix(n, seed + 1000.0), T, Q2)
    return (Q1 * want) @ transpose(Q2)
}

test "a diagonal matrix has its own entries as singular values" {
    n = 4
    A = diag([3.0, -5.0, 0.5, 2.0])
    U = zeros(n, n)
    s = zeros(n)
    V = zeros(n, n)
    junk = svd(A, U, s, V, n, n)
    assert_near(s[0], 5.0, 1.0e-13)
    assert_near(s[1], 3.0, 1.0e-13)
    assert_near(s[2], 2.0, 1.0e-13)
    assert_near(s[3], 0.5, 1.0e-13)
    assert_lt(svd_residual(A, U, s, V, n, n), 1.0e-13)
}

test "U S V^T puts a square matrix back, with both factors orthonormal" {
    worst_res = 0.0
    worst_u = 0.0
    worst_v = 0.0
    ncases = 10
    for si to ncases {
        n = 2 + mod(si, 7.0)
        A = test_matrix(n, 17.0 + 31.0 * si)
        U = zeros(n, n)
        s = zeros(n)
        V = zeros(n, n)
        junk = svd(A, U, s, V, n, n)
        worst_res = max(worst_res, svd_residual(A, U, s, V, n, n))
        worst_u = max(worst_u, svd_orthonormal(U, n, n))
        worst_v = max(worst_v, svd_orthonormal(V, n, n))
        for j to n - 1 {
            assert_lt(-1.0e-15, s[j] - s[j + 1])
        }
    }
    assert_lt(worst_res, 1.0e-13)
    assert_lt(worst_u, 1.0e-13)
    assert_lt(worst_v, 1.0e-13)
}

test "a tall matrix works too, and its singular values are the square roots of eig(A^T A)" {
    m = 9
    n = 4
    A = zeros(m, n)
    B = test_matrix(m, 71.0)
    for i to m {
        for j to n {
            A[i, j] = B[i, j]
        }
    }
    U = zeros(m, n)
    s = zeros(n)
    V = zeros(n, n)
    junk = svd(A, U, s, V, m, n)
    assert_lt(svd_residual(A, U, s, V, m, n), 1.0e-13)
    assert_lt(svd_orthonormal(U, m, n), 1.0e-13)
    G = transpose(A) @ A
    vals = zeros(n)
    vecs = zeros(n, n)
    eig_sym(G, vals, vecs)
    # eig_sym returns them ascending, svd descending
    for j to n {
        assert_near(s[j], sqrt(vals[n - 1 - j]), 1.0e-11)
    }
}

test "the small singular values survive, where going through A^T A loses them" {
    n = 5
    want = zeros(n)
    want[0] = 1.0
    want[1] = 1.0e-3
    want[2] = 1.0e-6
    want[3] = 1.0e-9
    want[4] = 1.0e-12
    A = svd_with_values(want, n, 23.0)
    U = zeros(n, n)
    s = zeros(n)
    V = zeros(n, n)
    junk = svd(A, U, s, V, n, n)
    # the squared route: eigenvalues of A^T A, whose condition number is squared
    G = transpose(A) @ A
    vals = zeros(n)
    vecs = zeros(n, n)
    eig_sym(G, vals, vecs)
    jac = zeros(n)
    sq = zeros(n)
    for j to n {
        jac[j] = abs(s[j] - want[j]) / want[j]
        sq[j] = abs(sqrt(max(vals[n - 1 - j], 0.0)) - want[j]) / want[j]
    }
    assert_lt(jac[0], 1.0e-14)
    assert_lt(jac[1], 1.0e-13)
    assert_lt(jac[2], 1.0e-10)
    assert_lt(jac[3], 1.0e-6)
    assert_lt(jac[4], 1.0e-3)
    # the squared route has thrown away the fourth and the fifth entirely
    assert_lt(0.1, sq[3])
    assert_lt(0.5, sq[4])
    # and is worse at every level below the largest
    for j to n - 1 {
        assert_lt(jac[j + 1], sq[j + 1])
    }
}

test "a rank-deficient matrix reports the rank it has" {
    n = 6
    want = zeros(n)
    want[0] = 4.0
    want[1] = 2.0
    want[2] = 0.5
    A = svd_with_values(want, n, 137.0)
    U = zeros(n, n)
    s = zeros(n)
    V = zeros(n, n)
    junk = svd(A, U, s, V, n, n)
    assert_near(s[0], 4.0, 1.0e-12)
    assert_near(s[1], 2.0, 1.0e-12)
    assert_near(s[2], 0.5, 1.0e-12)
    rank = 3
    for j to n - rank {
        assert_lt(s[rank + j], 1.0e-14)
    }
    assert_lt(svd_residual(A, U, s, V, n, n), 1.0e-13)
}

# SPDX-License-Identifier: MIT OR Apache-2.0
# SPDX-FileCopyrightText: 2026 The Viv authors