# Classical vs. modified Gram-Schmidt

It is often said that modified Gram-Schmidt is more robust with respect to rounding errors than classical Gram-Schmidt, but it is very hard to find a good explanation / example of why this is so. Can anyone provide such an explanation / example.

Consider what happens if you apply modified / classical Gram-Schmidt to the matrix $$\begin{pmatrix} 1 & 1 & 1 \\ 1 & 1 & 1 \\ 0 & \varepsilon & 0 \\ 0 & 0 & \varepsilon \\ \end{pmatrix}$$ where $$\varepsilon \ll 1$$.
The two algorithms agree on the first and second column and produce a $$Q$$ factor whose entries are given by $$\begin{pmatrix} 1 & \varepsilon_\mathrm{mach} \, \varepsilon^{-1} & * \\ 1 & \varepsilon_\mathrm{mach} \, \varepsilon^{-1} & * \\ 0 & 1 & * \\ 0 & 0 & * \\ \end{pmatrix}$$ if we assume finite machine arithmetic with precision $$\varepsilon_\mathrm{mach}$$ and restrict ourselves to order-of-magnitude estimates. The first two entries in the second column may be explained by noting that $$\tilde q_2 := a_2 - q_1^T a_2 \, q_1$$ would set the first two entries of $$\tilde q_2$$ to zero in exact arithmetic, but rounding errors will make them $$\mathcal{O}(\varepsilon_\mathrm{mach})$$ and normalising $$\tilde q_2$$ will lead to the values quoted above.
Classical Gram-Schmidt now computes the last column as follows: $$\begin{pmatrix} 1 & \varepsilon_\mathrm{mach} \, \varepsilon^{-1} & \varepsilon_\mathrm{mach} \, \varepsilon^{-1} \\ 1 & \varepsilon_\mathrm{mach} \, \varepsilon^{-1} & \varepsilon_\mathrm{mach} \, \varepsilon^{-1} \\ 0 & 1 & \varepsilon_\mathrm{mach} \, \varepsilon^{-2} \\ 0 & 0 & 1 \\ \end{pmatrix}$$ The first two entries in the last column can be derived as above, and the third entry follows from $$q_2^T \, a_3 = \mathcal{O}\bigl(\varepsilon_\mathrm{mach} \, \varepsilon^{-1}\bigr)$$ which yields $$(\tilde q_3)_3 = \mathcal{O}\bigl(\varepsilon_\mathrm{mach} \, \varepsilon^{-1}\bigr)$$ and hence $$(\tilde q_3)_3 = \mathcal{O}\bigl(\varepsilon_\mathrm{mach} \, \varepsilon^{-2}\bigr)$$ after normalisation. We observe that if $$\varepsilon = 10^{-8}$$ and $$\varepsilon_\mathrm{mach} = 10^{-16}$$, then the first column of $$Q$$ is orthogonal to the other two up to an error of about $$\varepsilon_\mathrm{mach} \, \varepsilon^{-1} = 10^{-8}$$, while the second and third column are orthogonal up to an error of $$\varepsilon_\mathrm{mach} \, \varepsilon^{-2} = 1$$, i.e. they are not orthogonal at all.
In contrast, the modified Gram-Schmidt algorithm leads to $$\begin{pmatrix} 1 & \varepsilon_\mathrm{mach} \, \varepsilon^{-1} & \varepsilon_\mathrm{mach} \, \varepsilon^{-1} \\ 1 & \varepsilon_\mathrm{mach} \, \varepsilon^{-1} & \varepsilon_\mathrm{mach} \, \varepsilon^{-1} \\ 0 & 1 & \varepsilon_\mathrm{mach}^2 \, \varepsilon^{-2} \\ 0 & 0 & 1 \\ \end{pmatrix}$$ since after projecting out the component in the first direction, the temporary vector $$\tilde q_3^{(1)}$$ is given by $$\tilde q_3^{(1)} = \begin{pmatrix} \varepsilon_\mathrm{mach} \, \varepsilon^{-1} \\ \varepsilon_\mathrm{mach} \, \varepsilon^{-1} \\ 0 \\ 1 \end{pmatrix}$$ such that $$q_2^T \tilde q_3^{(1)} = \mathcal{O}\bigl(\varepsilon_\mathrm{mach}^2 \, \varepsilon^{-1}\bigr)$$. With the values given above, we now have that the second and third column are orthogonal up to an error $$\varepsilon_\mathrm{mach}^2 \, \varepsilon^{-2} = 10^{-16}$$.