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.

Error evaluation

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

Error evaluation

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

Update (2021/06/05): I was recently made aware of the publication [2]. Its author discusses several key aspects of QR algorithms such as the column pivoting and the need for reorthogonalisation. The paper contains an interesting anaylsis on the number of reorthogonalisations that might be required. It states that two reorthogonalisations might not always be enough but that more than 4 reorthogonalisations are highly unlikely. In addition, there’s a good list of references that show that these ideas were studied intensively in the 1970s.

References

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

[2] A. Dax, A modified Gram–Schmidt algorithm with iterative orthogonalization and column pivoting, Linear Algebra and its Applications (310), 25–42, 2000

Notes

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

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

Changelog

  • 2021/06/05: Added a reference.
Mathematician and
Software Engineer

Researcher, Engineer, Tinkerer, Scholar, Philosopher, and Hyrox afficionado