Editing Gram–Schmidt process (section)

== Numerical stability ==
When this process is implemented on a computer, the vectors <math>\mathbf{u}_k</math> are often not quite orthogonal, due to [[round-off error|rounding errors]]. For the Gram–Schmidt process as described above (sometimes referred to as "classical Gram–Schmidt") this loss of orthogonality is particularly bad; therefore, it is said that the (classical) Gram–Schmidt process is [[numerical stability|numerically unstable]].

The Gram–Schmidt process can be stabilized by a small modification; this version is sometimes referred to as '''modified Gram-Schmidt''' or MGS. This approach gives the same result as the original formula in exact arithmetic and introduces smaller errors in finite-precision arithmetic.

Instead of computing the vector {{math|'''u'''<sub>''k''</sub>}} as
<math display="block"> \mathbf{u}_k = \mathbf{v}_k - \operatorname{proj}_{\mathbf{u}_1} (\mathbf{v}_k) - \operatorname{proj}_{\mathbf{u}_2} (\mathbf{v}_k) - \cdots - \operatorname{proj}_{\mathbf{u}_{k-1}} (\mathbf{v}_k), </math>
it is computed as
<math display="block"> \begin{align}
\mathbf{u}_k^{(1)} &= \mathbf{v}_k - \operatorname{proj}_{\mathbf{u}_1} (\mathbf{v}_k), \\
\mathbf{u}_k^{(2)} &= \mathbf{u}_k^{(1)} - \operatorname{proj}_{\mathbf{u}_2} \left(\mathbf{u}_k^{(1)}\right), \\
& \;\; \vdots \\
\mathbf{u}_k^{(k-2)} &= \mathbf{u}_k^{(k-3)} - \operatorname{proj}_{\mathbf{u}_{k-2}} \left(\mathbf{u}_k^{(k-3)}\right), \\
\mathbf{u}_k^{(k-1)} &= \mathbf{u}_k^{(k-2)} - \operatorname{proj}_{\mathbf{u}_{k-1}} \left(\mathbf{u}_k^{(k-2)}\right), \\
\mathbf{e}_k &=  \frac{\mathbf{u}_k^{(k-1)}}{\left\|\mathbf{u}_k^{(k-1)}\right\|}
\end{align} </math>

This method is used in the previous animation, when the intermediate <math>\mathbf{v}'_3</math> vector is used when orthogonalizing the blue vector <math>\mathbf{v}_3</math>.

Here is another description of the modified algorithm. Given the vectors <math>\mathbf{v}_1, \mathbf{v}_2, \dots, \mathbf{v}_n</math>, in our first step we produce vectors <math>\mathbf{v}_1, \mathbf{v}_2^{(1)}, \dots, \mathbf{v}_n^{(1)}</math>by removing components along the direction of <math>\mathbf{v}_1</math>. In formulas, <math>\mathbf{v}_k^{(1)} := \mathbf{v}_k - \frac{\langle \mathbf{v}_k, \mathbf{v}_1 \rangle}{\langle \mathbf{v}_1, \mathbf{v}_1 \rangle} \mathbf{v}_1</math>. After this step we already have two of our desired orthogonal vectors <math>\mathbf{u}_1, \dots, \mathbf{u}_n</math>, namely <math>\mathbf{u}_1 = \mathbf{v}_1, \mathbf{u}_2 = \mathbf{v}_2^{(1)}</math>, but we also made <math>\mathbf{v}_3^{(1)}, \dots, \mathbf{v}_n^{(1)}</math> already orthogonal to <math>\mathbf{u}_1</math>. Next, we orthogonalize those remaining vectors against <math>\mathbf{u}_2 = \mathbf{v}_2^{(1)}</math>. This means we compute <math>\mathbf{v}_3^{(2)}, \mathbf{v}_4^{(2)}, \dots, \mathbf{v}_n^{(2)}</math> by subtraction <math>\mathbf{v}_k^{(2)} := \mathbf{v}_k^{(1)} - \frac{\langle \mathbf{v}_k^{(1)}, \mathbf{u}_2 \rangle}{\langle \mathbf{u}_2, \mathbf{u}_2 \rangle} \mathbf{u}_2</math>. Now we have stored the vectors <math>\mathbf{v}_1, \mathbf{v}_2^{(1)}, \mathbf{v}_3^{(2)}, \mathbf{v}_4^{(2)}, \dots, \mathbf{v}_n^{(2)}</math> where the first three vectors are already <math>\mathbf{u}_1, \mathbf{u}_2, \mathbf{u}_3</math> and the remaining vectors are already orthogonal to <math>\mathbf{u}_1, \mathbf{u}_2</math>. As should be clear now, the next step orthogonalizes <math>\mathbf{v}_4^{(2)}, \dots, \mathbf{v}_n^{(2)}</math> against <math>\mathbf{u}_3 = \mathbf{v}_3^{(2)}</math>. Proceeding in this manner we find the full set of orthogonal vectors <math>\mathbf{u}_1, \dots, \mathbf{u}_n</math>. If orthonormal vectors are desired, then we normalize as we go, so that the denominators in the subtraction formulas turn into ones.