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

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.

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

Mathematician and
Software Engineer

Researcher, Engineer, Tinkerer, Scholar, and Philosopher.

Related