0

Can anyone explain the stability of Gram-Schmidt and how to use it?

For example, how can I use it to find the $QR$ factorization of the coefficient matrix $A$ and solve this system:

(I know how to use Gram-Schmidt, just not stable Gram-Schmidt)

$ A=\left( \begin{array}{ccc} 2 & 1 & -1 \\ 1 & 0 & 2 \\ 2 & -1 & 3 \end{array} \right)$

$ \left( \begin{array}{ccc} 2 & 1 & -1 \\ 1 & 0 & 2 \\ 2 & -1 & 3 \end{array} \right)$ $ \left( \begin{array}{ccc} x \\ y \\ z \end{array} \right)$ $ =\left( \begin{array}{ccc} 2 \\ -1 \\ 0 \end{array} \right)$

hawk2015
  • 127

1 Answers1

2

"Classical" Gram-Schmidt, in which you subtract off the projections of the $(k+1)$th vector onto the first $k$ vectors, is quite unstable, especially in high dimensions, because you essentially ensure that your new vector is orthogonal to the input vector in question but fail to ensure that the vectors you get at the end of the process are orthogonal to each other. Combine that with the fact that you can wind up subtracting nearly equal numbers and you get a bad situation.

"Modified" Gram-Schmidt, in which you subtract off the projection onto the first vector, then subtract off the projection of that result onto the second vector, etc., is stable, albeit slightly less so than the Householder transformation approach.

The two are exactly the same if your matrix only has two columns, so the smallest interesting example is $3 \times 3$. An example:

myeps=1e-12;
x1=rand(3,1);
x2=rand(3,1);
x3=x1+x2+myeps*rand(3,1);
y1=x1;
q1=y1/norm(y1);
y2=x2-(q1'*x2)*q1;
q2=y2/norm(y2);
y3=x3-(q1'*x3)*q1;
y3=y3-(q2'*y3)*q2;
q3=y3/norm(y3);

For classical Gram-Schmidt, the second to last line would be

y3=y3-(q2'*x3)*q2;

instead, but that would be the only difference (in the 3x3 setting). The stability issue is that you do arithmetic with y3, introducing some rounding errors in those digits, and then you multiply them by roughly $10^{12}$ when you normalize to make q3. This actually happens either way, but it is worse with classical GS.

Ian
  • 104,572