Gram–Schmidt Orthogonalisation

Choose how close together the columns are, then read the two grids: they give the angle between every pair of vectors the two versions hand back.

Worked example — Gram–Schmidt Orthogonalisation

Condition number of A	1.7321 × 10⁸
What is left of column 3	1.4142 × 10⁻⁸
What q₂ takes out of column 3	0
Largest dot product between two of them	0.5000
Angle between q₂ and q₃	60.00°
Loss of orthogonality, classical	0.7071
Loss of orthogonality, modified	1.1547 × 10⁻⁸
Loss of orthogonality, classical run twice	2.6158 × 10⁻²⁴

These are the values the calculator shows when the page opens. Change any input above to recompute them.

Loading interactive simulation...

Classical Gram–Schmidt hands back two vectors 60° apart 🖖

Press ε = 10⁻⁸, classical. In the classical grid two cells read 60.00°, the pair q₂ and q₃ in both corners, and every other angle in both grids reads 90.00°. The dot product behind them is 0.5000: a half, to within one unit in the last place of a double.

The matrix is three columns, each of them the vector (1, 0, 0, 0) with a single 10⁻⁸ dropped into a different slot. Peter Läuchli used it in 1961 and it has carried his name since. Columns that nearly repeat are the awkward case for any least-squares fit, and they turn up whenever two of your measurements are almost the same measurement.

Modified Gram–Schmidt is handed the same three columns and every angle it returns is 90.00°, for the same count of multiplications.

Move the separation back one notch, to 10⁻⁷, and that cell reads 89.24°. One click of the slider is the difference between being wrong by three quarters of a degree and being wrong by thirty.

The two versions read one number differently 🖖

Set the walkthrough to column 3 and press Classical, then Modified. The row What q₂ takes out of column 3 goes from 0 to 7.0711 × 10⁻⁹. At three columns that is the only coefficient the two algorithms compute differently. The diagonal entry below it, and every angle in the grid, follows from it.

Classical evaluates that coefficient as q₂ᵀa₃, against the column it was handed. Modified evaluates q₂ᵀv, against whatever is left of the column once q₁'s share has come off. On paper the two are the same number, because v differs from a₃ only by a multiple of q₁ and q₂ is perpendicular to q₁. In arithmetic q₂ is not quite perpendicular to q₁, and a₃ still has 7.0711 × 10⁻⁹ of itself pointing along q₂. Classical never looks there, so q₃ keeps it, and that is the 60°.

Drag the slider from 10⁻¹ to 10⁻⁸ and watch the two curves on the lower panel come apart. Fitted over 113 points across that range, classical's loss of orthogonality rises as κ^1.98 on these three columns and modified's as κ^0.98. The panel refits whichever column count you select, so those two exponents move as you widen the matrix. One reordered subtraction is worth a whole power of the condition number.

The basis is ruined and the factorisation is exact 🖖

The last row of the results table is ‖A − QR‖ / ‖A‖. At ε = 10⁻⁸ under Classical it reads 0. Multiply the ruined Q by the R that came out beside it and the original matrix comes back, entry for entry, with nothing left over.

Swept over every state the controls reach — four column counts, ten separations and both algorithms — that relative residual never rises above 5.5511 × 10⁻¹⁷. The gap between 1 and the next double is 2.2204 × 10⁻¹⁶, so it stays a quarter of one unit in the last place.

Gram–Schmidt did write A as a product of the vectors it returned and the coefficients it recorded. What it stopped doing was returning vectors at right angles to each other. Which of the two properties you need is a question about your own code: projecting a vector with Qᵀb leans on the one that failed, rebuilding A from Q and R leans on the one that held.

Problem solved in full

Where classical Gram–Schmidt loses its right angle 8 steps
The tool is showing three columns 10⁻⁸ apart and classical Gram–Schmidt's answer for them. Work the algorithm by hand, find the one rounding that costs it thirty degrees, and then say how far either version can be pushed before it stops being worth running.

Open this problem in the tool
1. Each column is the vector (1, 0, 0, 0) with one ε dropped into a different slot. AᵀA is the all-ones 3 × 3 matrix plus ε²I, whose eigenvalues are 3 + ε² once and ε² twice, so the singular values are √(3 + ε²) and ε and the condition number is exact rather than estimated. The top row of the table reads 1.7321 × 10⁸.
  
  $A = 1 ε 00 10 ε 0 100 ε, ε = 1 0^{- 8}, κ (A) = \frac{3 + ε ^{2}}{ε} = 1.7321 \times 1 0^{8}$
2. ‖a₁‖² is 1 + ε², and ε² is 10⁻¹⁶. A double holds 1, and then the next number it can hold is 2.2204 × 10⁻¹⁶ further along, so anything below half that gap, 1.1102 × 10⁻¹⁶, disappears when it is added to 1. Ten to the minus sixteen is below it. The square root is 1, dividing by it changes nothing, and q₁ comes out as the column itself with the ε² struck from the record. Every other rounding in this problem is worth one unit in the last place. This one is worth thirty degrees.
  
  $∥ a_{1} ∥^{2} = 1 + ε^{2} = 1 + 1 0^{- 16}, 1 0^{- 16} < 1.1102 \times 1 0^{- 16} ⟹ fl (1 + ε^{2}) = 1, q_{1} = a_{1}$
3. Column 2 next, and both versions agree here because there is only one vector behind it. The projection coefficient is q₁ᵀa₂ = 1, subtracting it leaves (0, −ε, ε, 0), and that has length ε√2 = 1.4142 × 10⁻⁸. Column 2 arrived with length 1 as far as a double can tell, and all but 1.4142 × 10⁻⁸ of it was already covered by q₁.
  
  $r_{12} = q_{1}^{T} a_{2} = 1, v_{2} = a_{2} - q_{1} = (0, - ε, ε, 0)^{T}, ∥ v_{2} ∥ = ε 2 = 1.4142 \times 1 0^{- 8}$
4. Column 3, classically. The coefficient against q₁ is 1 again. The coefficient against q₂ is q₂ᵀa₃, and q₂ lives entirely in slots 2 and 3 while a₃ keeps its ε in slot 4, so that dot product is 0 — which is what the What q₂ takes out of column 3 row says. Subtracting q₁ alone leaves (0, −ε, 0, ε), of length ε√2 again, so What is left of column 3 also reads 1.4142 × 10⁻⁸.
  
  $r_{13} = q_{1}^{T} a_{3} = 1, r_{23} = q_{2}^{T} a_{3} = 0, v_{3} = a_{3} - q_{1} = (0, - ε, 0, ε)^{T}, q_{3} = \frac{1}{2} (0, - 1, 0, 1)^{T}$
5. Both q₂ and q₃ were built by taking q₁ away from a column that shares its leading 1, so both carry −1/√2 in slot 2 and nothing else in common. Their dot product is a half. The tool prints 0.5000 and 60.00°, and the value it stores is 0.4999999999999999, one unit in the last place away from the half the algebra gives, because ±1/√2 is not a binary fraction. That is the size of error this arithmetic is capable of. The thirty degrees are not it. The classical loss row reads 0.7071, which is √(2 × 0.5²): the same half counted once in each corner of QᵀQ.
  
  $q_{2}^{T} q_{3} = \frac{1}{2} [(- 1) (- 1) + (1) (0)] = \frac{1}{2}, arccos \frac{1}{2} = 6 0^{\circ}, 2 \cdot 0 . 5^{2} = 0.7071$
6. Modified, on the same column. After q₁'s share comes off, the remainder is (0, −ε, 0, ε), and this version measures against that instead of against a₃. Dotting q₂ into it gives ε/√2 = 7.0711 × 10⁻⁹, not zero. Take that away and what is left is (0, −ε/2, −ε/2, ε), of length ε√(3/2) = 1.2247 × 10⁻⁸, and q₂ᵀq₃ comes out at 1.1102 × 10⁻¹⁶ instead of a half. The modified loss row reads 1.1547 × 10⁻⁸, and all that is left in it is what q₁ could not take away: the subtraction leaves −ε² lying along q₁, and dividing by a remainder only ε√2 long turns that into ε/√2 again.
  
  $r_{23} = q_{2}^{T} v_{3} = \frac{ε}{2} = 7.0711 \times 1 0^{- 9}, v_{3} \mapsto (0, - \frac{ε}{2}, - \frac{ε}{2}, ε)^{T}, ∥ v_{3} ∥ = ε \frac{3}{2} = 1.2247 \times 1 0^{- 8}, q_{2}^{T} q_{3} = 0$
7. The same argument at every column count the tool offers. Classical builds every later q by taking q₁ away from a column that shares its leading 1, so every pair among them comes out at 60°: (n − 1)(n − 2)/2 pairs, and a loss of orthogonality equal to the square root of half that count. With two columns there are no such pairs, and what is left is the ε-sized tilt both versions share.
  
  $n 2345 κ (A) 1.4142 \times 1 0^{8} 1.7321 \times 1 0^{8} 2.0000 \times 1 0^{8} 2.2361 \times 1 0^{8} # {6 0^{\circ}} 0136 ∥ Q_{C}^{T} Q_{C} - I ∥_{F} 1.0000 \times 1 0^{- 8} 0.7071 1.2247 1.7321 ∥ Q_{M}^{T} Q_{M} - I ∥_{F} 1.0000 \times 1 0^{- 8} 1.1547 \times 1 0^{- 8} 1.2247 \times 1 0^{- 8} 1.2649 \times 1 0^{- 8}$
8. Classical's fitted line reaches a loss of 1 at κ = 2.35 × 10⁸, and the curve itself has stopped rising well before then: past 10⁻⁸ every pair after the first sits at exactly 60°, and the loss can go no higher than the 0.7071 step 7 derived. Modified gets a measured bound instead of an extrapolated one, because its own curve turns over past 10⁻⁸ and a fitted line carried out that far would be refuted by the panel underneath it: across 724 points of the same curve, 181 separations at each column count, its loss never exceeded 1.25 machine epsilons times κ, and that bound reaches 1 only at κ = 3.60 × 10¹⁵. One over the machine epsilon is 4.50 × 10¹⁵.
  
  $lo g_{10} κ = \frac{16 . 5755}{1 . 9801} 1.25 \cdot (2.2204 \times 1 0^{- 16}) κ (2.2204 \times 1 0^{- 16})^{- 1} = = = 8.3710 ⟹ κ = 2.35 \times 1 0^{8} 1 ⟹ κ = 3.60 \times 1 0^{15} 4.50 \times 1 0^{15}$
Answer

The thirty degrees come out of one line. Classical reads q₂ᵀa₃ = 0 where modified reads q₂ᵀv₃ = 7.0711 × 10⁻⁹, and on paper those are the same number.

The separation of 10⁻⁸ is the first notch at which ε² vanishes from 1 + ε², since √(half the machine epsilon) is 1.0537 × 10⁻⁸. From there down the cell reads 60.00° however much further you drag; one notch back, at 10⁻⁷, it reads 89.24° and classical has already lost 1.88 × 10⁻² of orthogonality. The collapse starts well before the slider reaches the notch that makes it clean.

Åke Björck published the analysis of the modified version in 1967, and the reordering has been the default ever since. The two run the same number of multiplications, and at this separation the loss of orthogonality between them differs by a factor of 6.1 × 10⁷. Modified keeps its right angles until about the point where a double stops being able to tell the columns apart at all; classical gives them up seven orders of magnitude earlier.

$q_{2}^{T} a_{3} = 0 \neq = q_{2}^{T} v_{3} = \frac{ε}{2} = 7.0711 \times 1 0^{- 9}$

References (5)

The matrix the tool is built on, and the name it carries: P. Läuchli, "Jordan-Elimination und Ausgleichung nach kleinsten Quadraten", Numerische Mathematik 3(1), 226–240, 1961. doi:10.1007/BF01386022
The analysis of the modified variant that made it the default: Åke Björck, "Solving linear least squares problems by Gram-Schmidt orthogonalization", BIT 7(1), 1–21, 1967. doi:10.1007/BF01934122
The corresponding rounding-error analysis of the classical variant: Luc Giraud, Julien Langou, Miroslav Rozložník and Jasper van den Eshof, "Rounding error analysis of the classical Gram-Schmidt orthogonalization process", Numerische Mathematik 101(1), 87–100, 2005. doi:10.1007/s00211-005-0615-4
Gram's paper of 1883, where the process is set out for least squares: J. P. Gram, "Ueber die Entwickelung reeller Functionen in Reihen mittelst der Methode der kleinsten Quadrate", Journal für die reine und angewandte Mathematik 1883(94), 41–73. doi:10.1515/crll.1883.94.41
And Schmidt's of 1907, where the form used today appears: Erhard Schmidt, "Zur Theorie der linearen und nichtlinearen Integralgleichungen", Mathematische Annalen 63, 433–476, 1907. doi:10.1007/BF01449770

By the 42Math team · Reviewed against published sources. Corrections welcome.