Orthogonalisation

We discuss the necessary changes to our QR decomposition algorithms to handle matrices which do not have full rank.

We developed a QR decomposition algorithm, based on the orthogonalisation process of Gram-Schmidt in a series of posts here, here, here, and here. Let’s have a look how good this algorithm performs against built-in implementations from julia and other programming languages. Test criteria For a given matrix $A$ we get a decomposition $A=QR$ where $Q$ is orthogonal and where $R$ is upper triangular. If $A$ does not have full column rank, then some columns in our matrix $Q$ will be 0 (at least in our implementation).

Problem Formulation We already discussed QR decompositions and showed that using the modified formulation of Gram-Schmidt significantly improves the accuracy of the results. However, there is still an error of about $10^3 M_\varepsilon$ (where $M_\varepsilon$ is the machine epsilon) when using the modified Gram Schmidt as base algorithm for the orthogonalisation. These errors are due to cancellations due to the limited precision in our floating point representation. Reorthogonalisations We consider the following idea: each time we compute the next orthogonal vector, we check whether cancellations occurred and whether our result might be inaccurate.

We consider the necessary changes to the Gram-Schmidt orthogonalisation to obtain a QR Decomposition

We compare the accuracy of the classical Gram-Schmidt algorithm to the modified Gram-Schmidt algorithm.

Floating point computations on computers may behave differently than one might expect. Every software developer should be aware of these since computed results may be off by orders of magnitude in the worst case. My first encounter with this subject was when I implemented the algorithms for the inpainting mask optimization. Recently I investigated the accuracy of QR decompositions that use Gram-Schmidt for the orthogonalisation process. The findings have been posted on this site in a series of posts in 2020.