1
$\begingroup$

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.

$\endgroup$
4
$\begingroup$

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}$.

$\endgroup$
  • $\begingroup$ You can select your answer $\endgroup$ – nicoguaro Aug 25 at 17:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.