problemAINumerics
Make Gram-Schmidt orthogonal to rounding
Ten nearly parallel vectors, each plus a tiny 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 and equal the exact ones to .
Open in the workbenchDownload .vivOutputs below are from the site build’s run; Run all runs it again in your browser.
Look inside this program
- Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.
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 : after the subtraction what is left of each vector is of size , and the rounding in the subtracted parts, half a unit in the last place of numbers near 1, is then a relative error of 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 , much better but still eight orders of magnitude from rounding, and the tests ask for . 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.
Your task
Write these 3 functions:
lauchli()- 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}.
orthonormalize(A: Arr)- 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.
loss_of_orthogonality(Q: Arr)- The largest |q_i \cdot q_j - \delta_{ij}| over the rows of Q.
Make these tests pass 4 of 4 fail
- the rows are orthonormal to 1e-15 (classical Gram-Schmidt gives 0.5, modified 7.1e-9) assert_lt: 0.5 >= 1e-15 (fails)
- the second vector is mpmath's (7.07e-9, -1/sqrt 2, 1/sqrt 2, 0, ...) assert_near: |0 - 7.07107e-09| = 7.07107e-09 > 1e-12 (fails)
- the last vector is mpmath's (1.05e-9, -0.10541 nine times, 0.94868) assert_near: |0 - 1.05409e-09| = 1.05409e-09 > 1e-12 (and 10 more) (fails)
- each input row is rebuilt from the rows before it: a_k = sum over j <= k of (q_j . a_k) q_j assert_near: |-5e-09 - 0| = 5e-09 > 1e-15 (and 17 more) (fails)
Printed