# QR Decompositions with Reorthogonalisation

## 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. If this is the case we simply repeat the orthogonalisation (the result can’t get worse as we are already more orthogonal than before and the improved orthogonality makes the whole problem well-behaved). The only downside is a higher computational cost. As mentioned in [1], the idea goes back to Heinz Rutishauser.

We need a criterium which tells us when a reorthogonalisation becomes necessary. Here we use the criteria of Rutishauser as described in [1]. Cancellation has occurred if the following condition holds.

$$ \lVert A_{*,k} - \sum_{i=1}^{k-1} R_{i,k}Q_{*,i}\rVert \ll \lVert A_{*,k} \rVert $$

As a rule of thumb one may say that we have lost at least one decimal digit if the lefthand side is smaller than the righthand side by more than a factor 10. Adapting our previously developed code is now straightforward.

```
using LinearAlgebra
function qr_reorth_gram_schmidt(matrix)
# perform a QR decomposition using classical Gram Schmidt with reorthogonalisation
# We use the rule of thumb suggested by Rutishauser to decide when to reorthogonalise
# Note that this version does not work for rank deficient setups!
num_vectors = size(matrix)[2]
orth_matrix = copy(matrix)
r_matrix = zeros(num_vectors, num_vectors)
for vec_idx = 1:num_vectors
norm_orth_column = 0
perform_reorthogonalisation = true
while (perform_reorthogonalisation)
norm_current_column = norm(orth_matrix[:, vec_idx])
for span_base_idx = 1:(vec_idx-1) # orthogonalisation
projection_length = dot(orth_matrix[:, span_base_idx], orth_matrix[:, vec_idx])
if (norm_orth_column == 0)
# do not update R when performing reorthogonalisation
# norm_orth_column is exactly 0 on the first pass thgough
# this ensures that R gets updated properly
r_matrix[span_base_idx, vec_idx] = projection_length
end
orth_matrix[:, vec_idx] -= projection_length*orth_matrix[:, span_base_idx]
end
norm_orth_column = norm(orth_matrix[:, vec_idx])
# check whether reorthogonalisation is necessary.
perform_reorthogonalisation = (norm_orth_column < norm_current_column/10)
end
# normalise vector
orth_matrix[:, vec_idx] = orth_matrix[:, vec_idx]/norm(orth_matrix[:, vec_idx])
r_matrix[vec_idx, vec_idx] = dot(orth_matrix[:, vec_idx], matrix[:, vec_idx])
end
return (orth_matrix, r_matrix)
end
```

We use the same benchmark as for QR decompositions, that is, different size of the Hilbert matrix, and we compare on the one side the orthogonality of the $Q$ matrix as well as how good the product $QR$ approximates the input matrix.

```
error_qr_cgs_1 = []
error_qr_mgs_1 = []
error_qr_cgs_2 = []
error_qr_mgs_2 = []
error_qr_reorth_1 = []
error_qr_reorth_2 = []
for n = 1:10
H = [1.0/(i+j-1) for i=1:2^n, j=1:2^n] + 0.00001 * I # Hilbert Matrix
eye = Matrix{Float64}(I, 2^n, 2^n) # Identity Matrix
(Q_cgs, R_cgs) = qr_classical_gram_schmidt(H)
(Q_mgs, R_mgs) = qr_modified_gram_schmidt(H)
(Q_cgs_reorth, R_cgs_reorth) = qr_reorth_gram_schmidt(H)
append!(error_qr_cgs_1, norm(eye - transpose(Q_cgs)*Q_cgs, Inf))
append!(error_qr_mgs_1, norm(eye - transpose(Q_mgs)*Q_mgs, Inf))
append!(error_qr_reorth_1, norm(eye - transpose(Q_cgs_reorth)*Q_cgs_reorth, Inf))
append!(error_qr_cgs_2, norm(H - Q_cgs*R_cgs, Inf))
append!(error_qr_mgs_2, norm(H - Q_mgs*R_mgs, Inf))
append!(error_qr_reorth_2, norm(H - Q_cgs_reorth*R_cgs_reorth, Inf))
end
```

For the orthogonality we get the plot below.

For the accuracy of the reconstruction we get the following plot.

In both cases we have achieved an accuracy which remains stable slightly above $M_\varepsilon$.

## References

[1] W. Gander, Algorithms for the QR-Decomposition, Eidgenössische Technische Hochschule, Research Report No. 80-02, 1980. Retyped 2003

## Notes

All the code was evaluated with Julia version 1.5.1 (2020-08-25) using the official Julia Docker image.

This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License