Numerical Instability in QR decomposition

As we mentioned in class, the Gram-Schmidt algorithm used to compute the QR decomposition can become numerically unstable when the vectors are close to being linearly dependent. Here’s an explicit example of the failure of the GS algorithm (and even the modified GS algorithm mentioned in the notes) for the “Hilbert matrix” H:

We will see shortly that there’s a stable way to perform QR decomposition, but it requires a certain amount of insight in order to set up. This demonstration is just intended to convince you of the need for a better algorithm. (Again, this is adapted from a Wolfram Demonstration.)

cgs[a_]:=Module[{u=a},
For[k=2,k≤Length[a],k++,
u[[k]]=a[[k]]-Sum[Projection[a[[k]],u[[j]]], {j,1,k-1}]];
N@Transpose[Normalize/@u]
];

mgs[a_]:=Module[{aprime=a,q=a},
For[j=1,j≤Length[a],j++,
For[k=1,k<j,k++,
aprime[[j]]=aprime[[j]]-Projection[aprime[[j]],aprime[[k]]]
];
q[[j]]=Normalize[aprime[[j]]];
];
N@Transpose[q]
];

m[n_]:=Table[If[j==1,1,If[i==1,1/1000,If[i==j,1/1000,0]]], {i,1,n},{j,1,n}];

orthotest[matrix_]:=Max@(IdentityMatrix[Length[matrix]]-Transpose[matrix].matrix);

Manipulate[
matrix=Switch[testmatrix, 1, m[size], 2, HilbertMatrix[size]];
Pane[Text@Grid[{
{"","result", "measure of orthogonality"},
{"unstable GS", MatrixForm[Quiet@N@(Round[#, 1/1000]&/@(cgs[SetPrecision[matrix,precision]]))], orthotest[cgs[SetPrecision[matrix,precision]]]},
{"more stable GS", MatrixForm[Quiet@N@Round[#, 1/1000]&/@mgs[SetPrecision[matrix,precision]]], Quiet@N@orthotest[mgs[SetPrecision[matrix,precision]]]},
{"stable algorithm", MatrixForm[Quiet@N@Round[#, 1/1000]&/@Normalize/@Transpose[First[QRDecomposition[Transpose[matrix]]]]], Quiet@N@orthotest[Normalize/@First[QRDecomposition[matrix]]]}
},
Dividers->{{2->True, 3->True}, {2->Thick,3->True, 4->True}}], {500,425}, Alignment->Center],
{{precision, 3, "digits of precision"}, 1,15,1, Appearance->"Labeled"},
{{size,3,"size of matrix"}, 1,6,1, Appearance->"Labeled"},
{{testmatrix, 1, "test vectors"}, {1->"diagonal thousandths", 2->"Hilbert matrix"}},
SaveDefinitions->True,TrackedSymbols->Manipulate]

CAPTION

The Gram–Schmidt process is an algorithm for converting a set of linearly independent vectors into a set of orthonormal vectors with the same span. The classical Gram–Schmidt algorithm is numerically unstable, which means that when implemented on a computer, round-off errors can cause the output vectors to be significantly non-orthogonal. This instability can be improved with a small adjustment to the algorithm. This Demonstration tests the two algorithms on two families of linearly independent vectors. Selecting "diagonal thousandths" tests the algorithm on the vectors that form the columns of the matrix:

(({{1, 1, 1, …, 1}, {1/1000, 1/1000, 0, …, 0}, {1/1000, 0, 1/1000, …, 0}, {⋮, ⋮, ⋮,  , ⋮}, {1/1000, 0, 0, …, 1/1000}})).

Selecting "Hilbert matrix" tests the algorithm on the vectors that form the columns of the n×n Hilbert matrix:

(({{1, 1/2, 1/3, …, 1/n}, {1/2, 1/3, 1/4, …, 1/(n+1)}, {1/3, 1/4, 1/5, …, 1/(n+2)}, {⋮, ⋮, ⋮,  , ⋮}, {1/n, 1/(n+1), 1/(n+2), …, 1/(2n-1)}})).

Both of these families of matrices have small but nonzero determinants, an indicator of sorts that their columns are linearly independent, but barely so. Let Q be the matrix whose columns are the vectors produced by the algorithm. This Demonstration measures the orthogonality of the resulting vectors via the number (Subscript[‖I-(Q^T) Q‖, ∞]). The closer this number is to zero, the closer the columns of Q are to being orthogonal.

The Mathematica function QRDecomposition produces an orthogonal collection of vectors from a linearly independent one in a highly stable way. The result of this function is shown at the bottom of the pane.

THIS NOTEBOOK IS THE SOURCE CODE FROM

"Numerical Instability in the Gram-Schmidt Algorithm" from the Wolfram Demonstrations Project
 http://demonstrations.wolfram.com/NumericalInstabilityInTheGramSchmidtAlgorithm/

Contributed by: Chris Boucher

A full-function Wolfram Mathematica system (Version 6 or higher) is required to edit this notebook.
GET WOLFRAM MATHEMATICA »

© Wolfram Demonstrations Project & Contributors   |   Terms of Use   |   Make a new version of this Demonstration »