Notebooks · posts/make-gram-schmidt-orthogonal-to-rounding

Make Gram-Schmidt orthogonal to rounding

Adam Zettel · September 25, 2026 · written with an agent

Ten nearly parallel vectors, each $e_{0}$ plus a tiny $1 0^{- 8}$ in its own direction. Classical Gram-Schmidt turns them into ten "orthonormal" vectors of which two are at 60 degrees: their dot product is 0.5. Rewrite orthonormalize so the vectors it returns are orthogonal to $1 0^{- 15}$ and equal the exact ones to $1 0^{- 12}$ .

a_{k} = e_{0} + ε e_{k + 1}, ε = 1 0^{- 8} : i \neq = j max ∣ q_{i} \cdot q_{j} ∣ = 0.5 (classical)

Gram-Schmidt builds an orthonormal set one vector at a time: take the next vector, subtract its components along the vectors already made, normalize. It is how Krylov methods (Arnoldi, Lanczos) and deflation build their bases, and it fails when the input is close to dependent. These inputs, the Läuchli vectors, have condition number $3.2 \times 1 0^{8}$ : after the subtraction what is left of each vector is of size $1 0^{- 8}$ , and the rounding in the subtracted parts, half a unit in the last place of numbers near 1, is then a relative error of $1 0^{- 8}$ or more.

The classical version computes every dot product with the original vector, so the errors of the early projections are never corrected. The modified version recomputes each dot product with the vector as updated so far; on these inputs it reaches $7 \times 1 0^{- 9}$ , much better but still eight orders of magnitude from rounding, and the tests ask for $1 0^{- 15}$ . Running the orthogonalization a second time on the result (the rule known as "twice is enough") gets there, and so does a QR factorization by Householder reflections.

The expected vectors come from posts/refs/make-gram-schmidt-orthogonal-to-rounding.py, which runs Gram-Schmidt at 60 digits in mpmath, where rounding does not matter; it also measures the three double-precision versions.

Look inside this program

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

This is a problem: its tests fail until the functions in it are written. Edit a cell and run it to try; the cells below it show whether the tests pass.

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 35–147

#: N_VEC  vectors, each of dimension N_VEC + 1
N_VEC = 10
EPS = 1.0e-8

#| The Läuchli vectors as the rows of an N_VEC by (N_VEC + 1) array: row k is
#| $e_0 + \varepsilon e_{k+1}$.
def lauchli() {
    A = zeros(N_VEC, N_VEC + 1)
    for k to N_VEC {
        A[k, 0] = 1.0
        A[k, k + 1] = EPS
    }
    return A
}

#| Orthonormal rows spanning the same nested subspaces as the rows of A: row k of
#| the result is in the span of rows 0..k of A, with a positive dot product with
#| row k of A. As written it is classical Gram-Schmidt.
def orthonormalize(A: Arr) {
    n = shape(A)[0]
    d = shape(A)[1]
    Q = zeros(n, d)
    for k to n {
        a = zeros(d)
        for i to d {
            a[i] = A[k, i]
        }
        dots = zeros(n)
        for j to k {
            s = 0.0
            for i to d {
                s = s + Q[j, i] * A[k, i]
            }
            dots[j] = s
        }
        for j to k {
            for i to d {
                a[i] = a[i] - dots[j] * Q[j, i]
            }
        }
        nr = sqrt(sum(a * a))
        for i to d {
            Q[k, i] = a[i] / nr
        }
    }
    return Q
}

#| The largest $|q_i \cdot q_j - \delta_{ij}|$ over the rows of Q.
def loss_of_orthogonality(Q: Arr) {
    n = shape(Q)[0]
    d = shape(Q)[1]
    worst = 0.0
    for p to n {
        for q to n {
            s = 0.0
            for i to d {
                s = s + Q[p, i] * Q[q, i]
            }
            if p == q {
                s = s - 1.0
            }
            worst = max(worst, abs(s))
        }
    }
    return worst
}

A = lauchli()
Q = orthonormalize(A)
show loss_of_orthogonality(Q)
# the angle, in degrees, between the two rows whose dot product is largest
worst_angle = acos(loss_of_orthogonality(Q)) * 180.0 / acos(-1.0)
show worst_angle

test "the rows are orthonormal to 1e-15 (classical Gram-Schmidt gives 0.5, modified 7.1e-9)" {
    assert_lt(loss_of_orthogonality(orthonormalize(lauchli())), 1.0e-15)
}
test "the second vector is mpmath's (7.07e-9, -1/sqrt 2, 1/sqrt 2, 0, ...)" {
    Q = orthonormalize(lauchli())
    assert_near(Q[1, 0], 7.0710678118654749e-9, 1.0e-12)
    assert_near(Q[1, 1], -0.70710678118654747, 1.0e-12)
    assert_near(Q[1, 2], 0.70710678118654754, 1.0e-12)
    assert_near(Q[1, 3], 0.0, 1.0e-12)
}
test "the last vector is mpmath's (1.05e-9, -0.10541 nine times, 0.94868)" {
    Q = orthonormalize(lauchli())
    assert_near(Q[N_VEC - 1, 0], 1.0540925533894598e-9, 1.0e-12)
    for i from 1 to N_VEC {
        assert_near(Q[N_VEC - 1, i], -0.10540925533894598, 1.0e-12)
    }
    assert_near(Q[N_VEC - 1, N_VEC], 0.9486832980505138, 1.0e-12)
}
test "each input row is rebuilt from the rows before it: a_k = sum over j <= k of (q_j . a_k) q_j" {
    A = lauchli()
    Q = orthonormalize(A)
    for k to N_VEC {
        for i to N_VEC + 1 {
            r = A[k, i]
            for j from 0 to k + 1 {
                c = 0.0
                for m to N_VEC + 1 {
                    c = c + Q[j, m] * A[k, m]
                }
                r = r - c * Q[j, i]
            }
            assert_near(r, 0.0, 1.0e-15)
        }
    }
}

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

#: N_VEC  vectors, each of dimension N_VEC + 1
N_VEC = 10
EPS = 1.0e-8

#| The Läuchli vectors as the rows of an N_VEC by (N_VEC + 1) array: row k is
#| $e_0 + \varepsilon e_{k+1}$.
def lauchli() {
    A = zeros(N_VEC, N_VEC + 1)
    for k to N_VEC {
        A[k, 0] = 1.0
        A[k, k + 1] = EPS
    }
    return A
}

#| Orthonormal rows spanning the same nested subspaces as the rows of A: row k of
#| the result is in the span of rows 0..k of A, with a positive dot product with
#| row k of A. As written it is classical Gram-Schmidt.
def orthonormalize(A: Arr) {
    n = shape(A)[0]
    d = shape(A)[1]
    Q = zeros(n, d)
    for k to n {
        a = zeros(d)
        for i to d {
            a[i] = A[k, i]
        }
        dots = zeros(n)
        for j to k {
            s = 0.0
            for i to d {
                s = s + Q[j, i] * A[k, i]
            }
            dots[j] = s
        }
        for j to k {
            for i to d {
                a[i] = a[i] - dots[j] * Q[j, i]
            }
        }
        nr = sqrt(sum(a * a))
        for i to d {
            Q[k, i] = a[i] / nr
        }
    }
    return Q
}

#| The largest $|q_i \cdot q_j - \delta_{ij}|$ over the rows of Q.
def loss_of_orthogonality(Q: Arr) {
    n = shape(Q)[0]
    d = shape(Q)[1]
    worst = 0.0
    for p to n {
        for q to n {
            s = 0.0
            for i to d {
                s = s + Q[p, i] * Q[q, i]
            }
            if p == q {
                s = s - 1.0
            }
            worst = max(worst, abs(s))
        }
    }
    return worst
}

A = lauchli()
Q = orthonormalize(A)
show loss_of_orthogonality(Q)
# the angle, in degrees, between the two rows whose dot product is largest
worst_angle = acos(loss_of_orthogonality(Q)) * 180.0 / acos(-1.0)
show worst_angle

test "the rows are orthonormal to 1e-15 (classical Gram-Schmidt gives 0.5, modified 7.1e-9)" {
    assert_lt(loss_of_orthogonality(orthonormalize(lauchli())), 1.0e-15)
}
test "the second vector is mpmath's (7.07e-9, -1/sqrt 2, 1/sqrt 2, 0, ...)" {
    Q = orthonormalize(lauchli())
    assert_near(Q[1, 0], 7.0710678118654749e-9, 1.0e-12)
    assert_near(Q[1, 1], -0.70710678118654747, 1.0e-12)
    assert_near(Q[1, 2], 0.70710678118654754, 1.0e-12)
    assert_near(Q[1, 3], 0.0, 1.0e-12)
}
test "the last vector is mpmath's (1.05e-9, -0.10541 nine times, 0.94868)" {
    Q = orthonormalize(lauchli())
    assert_near(Q[N_VEC - 1, 0], 1.0540925533894598e-9, 1.0e-12)
    for i from 1 to N_VEC {
        assert_near(Q[N_VEC - 1, i], -0.10540925533894598, 1.0e-12)
    }
    assert_near(Q[N_VEC - 1, N_VEC], 0.9486832980505138, 1.0e-12)
}
test "each input row is rebuilt from the rows before it: a_k = sum over j <= k of (q_j . a_k) q_j" {
    A = lauchli()
    Q = orthonormalize(A)
    for k to N_VEC {
        for i to N_VEC + 1 {
            r = A[k, i]
            for j from 0 to k + 1 {
                c = 0.0
                for m to N_VEC + 1 {
                    c = c + Q[j, m] * A[k, m]
                }
                r = r - c * Q[j, i]
            }
            assert_near(r, 0.0, 1.0e-15)
        }
    }
}

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