Notebook

In [1]:

using LinearAlgebra

Classical vs. Modified Gram–Schmidt¶

This is an example adapted from Trefethen and Bau, Numerical Linear Algebra, lecture 9. We construct a square matrix $A$ with exponentially graded singular values between $2^{-1}$ and $2^{-80}$, perform QR factorization in various ways, and check the accuracy and orthogonality of the results.

In [2]:

U = qr(randn(80,80)).Q # U = random 80×80 orthogonal matrix
V = qr(randn(80,80)).Q # V = random 80×80 orthogonal matrix
Σ = Diagonal(2.0 .^ (-1:-1:-80)) # Σ = diagonal matrix with entries 2^-1, 2^-2, ..., 2^-80
A = U * Σ * V'

Out[2]:

80×80 Array{Float64,2}:
  0.00231288   -0.00188713    0.00115339   …  -0.000777549  -0.00196706 
  0.00447612    0.0119956     0.00404125       0.00753016    0.00083697 
  0.000730346   0.00228595    0.00128298       0.00071391    0.000510262
 -0.00165406    0.00106365    0.000736025      4.42455e-5    0.00118117 
  0.003155     -0.00739048    0.000132767     -0.0014877    -0.00334916 
  0.000346554  -0.00636389   -0.00172838   …  -0.00366515   -0.000830444
 -0.00235921   -0.00687301   -0.000762344     -0.00675667    0.000155346
  0.000472659   0.00683781    0.00318707       0.00217628    0.00242981 
 -0.00434494    0.00826401    0.00196047       0.00511342    0.00391427 
 -0.00131739    0.0088279     0.0003912        0.0055408     0.00253646 
 -0.000332517  -0.00588532   -0.000676342  …  -0.00485406    0.000115519
 -0.000461368  -0.000437043  -0.00198776       0.000728049   0.00106017 
  0.00210663    0.00154606    9.68214e-5       0.00219168   -0.000360623
  ⋮                                        ⋱                            
 -0.000297086  -0.000820249  -0.000581252     -0.00193116    0.000442198
 -0.00257615    0.0033895    -0.000220947      0.000960119   0.00199674 
  0.00276878    0.00638691    0.00230822   …  -0.00269343    0.00112297 
 -0.00217747    0.00454191    0.000260071      0.00577961    0.00127744 
 -0.00365117   -0.00727022   -0.00288899      -0.00559812    2.08879e-5 
 -0.00160264    0.00883925    0.000150251      0.00276509    0.00288517 
 -0.000872615   0.00200195    0.000797861     -0.00144595    0.00204208 
 -0.00141865   -0.00393479   -0.00054564   …  -0.00428038   -0.00132189 
 -0.000173487   0.00315454   -0.00112636       0.00162983    0.000267716
  0.00138451    0.00430689    0.000437716      0.000179362   0.000557466
 -0.000243029  -0.0106028    -0.00260152      -0.00219126   -0.0019942  
 -0.00246574   -0.00295037   -0.000836444     -0.00417088    0.000334032

Now, let's implement classical and modified Gram–Schmidt algorithms to return the (reduced) QR factorizations of a matrix A:

In [3]:

# Classical Gram–Schmidt (Trefethen algorithm 7.1), implemented in the simplest way
# (We could make it faster by unrolling loops to avoid temporaries arrays etc.)
function clgs(A)
    m,n = size(A)
    T = float(eltype(A))
    Q = similar(A, T)
    R = zeros(T,n,n)
    @views for j = 1:n
        aⱼ = A[:,j]
        vⱼ = copy(aⱼ) # use copy so that modifying vⱼ doesn't change aⱼ
        for i = 1:j-1
            qᵢ = Q[:,i]
            R[i,j] = qᵢ'aⱼ
            vⱼ .-= R[i,j] .* qᵢ
        end
        R[j,j] = norm(vⱼ)
        Q[:,j] .= vⱼ ./ R[j,j]
    end
    return Q, R
end

# Modified Gram–Schmidt (Trefethen algorithm 8.1)
function mgs(A)
    m,n = size(A)
    T = float(eltype(A))
    Q = similar(A, T)
    R = zeros(T,n,n)
    @views for j = 1:n
        aⱼ = A[:,j]
        vⱼ = copy(aⱼ)
        for i = 1:j-1
            qᵢ = Q[:,i]
            R[i,j] = qᵢ'vⱼ # ⟵ NOTICE: mgs has vⱼ, clgs has aⱼ
            vⱼ .-= R[i,j] .* qᵢ
        end
        R[j,j] = norm(vⱼ)
        Q[:,j] .= vⱼ ./ R[j,j]
    end
    return Q, R
end

Out[3]:

mgs (generic function with 1 method)

Next, use these on our matrix $A$ from above:

In [4]:

QC, RC = clgs(A);
QM, RM = mgs(A);

Finally, we plot the diagonal elements of $R$, because $r_{jj} = \Vert v_j \Vert$, the norm of the projection of $a_j$ on to the space orthogonal to the previous columns.

In [5]:

using PyPlot

n = size(A,2)
semilogy(diag(RC), "bo")
semilogy(diag(RM), "rx")
semilogy(2.0 .^ -(1:n), "k-")
semilogy(fill(sqrt(eps()), n), "b--")
semilogy(fill(eps(), n), "r--")
legend(["classical","modified",L"2^{-j}", L"\sqrt{\epsilon_\mathrm{mach}}", L"\epsilon_\mathrm{mach}"], loc="lower left")
ylabel(L"r_{jj}")
xlabel(L"j")

Out[5]:

PyObject Text(0.5, 28.0, '$j$')

Note that the errors in classical Gram–Schmidt are at the level of $\sqrt{\epsilon_\mathrm{mach}}$, whereas they are at the level of $\epsilon_\mathrm{mach}$ for modified Gram–Schmidt. Clearly, the errors are accumulating too fast for stability in the former algorithm.

Arbitrary-precision Gram–Schmidt¶

It turns out that the "errors" in the modifed Gram–Schmidt $r_{jj}$ are mainly due to errors in constructing $A$ in the first place (which make A less singular than we thought: its actual singular values don't go down to $2^{-80}$), not errors during the Gram–Schmidt process. We can see this by comparing to a computation in arbitrary precision arithmetic via big in Julia. (The default is about 80 decimal places, which is more than sufficient for us here.) In arbitrary precision, even classical Gram–Schmidt is fine.

In [6]:

QCbig, RCbig = clgs(big.(A))

semilogy(Float64.(diag(RCbig)), "o", mfc="none") # convert back to Float64 for plotting, since PyPlot doesn't understand BigFloat
semilogy(diag(RM), "rx")
legend(["classical GS in arbitrary precision","modified GS"], loc="upper right")
ylabel(L"r_{jj}")
xlabel(L"j")

Out[6]:

PyObject Text(0.5, 28.0, '$j$')

Orthogonality of $Q$¶

We can also look at how close $Q$ is to being orthogonal by computing $\Vert Q^*Q - I\Vert_1$:

In [7]:

opnorm(QC'*QC - I, 1), opnorm(QM'*QM - I, 1), opnorm(QCbig'*QCbig - I, 1)

Out[7]:

(52.22409280066077, 4.203714197941792, 3.363088553594440047513746051409831884824413129653818425769182742811618036868859e-37)

Unfortunately, both classical and modified Gram–Schmidt can yield very non-orthogonal $Q$ matrices! In general, they are both unstable as a way to compute $Q$, but the modified Gram–Schmidt algorithm is still useful for least-squares algorithms because $R$ is computed in a backwards-stable way, and the algorithm can be tweaked to solve the least-squares problem $Ax \approx b$ as explained in Trefethen chapter 19. (The key is to not compute $Q^* b$ via $Q$, but rather by adding $b$ as an additional column of $A$...)

Householder QR¶

The qr function in Julia (and Matlab and SciPy), in contrast, defaults to using a Householder QR algorithm. This gives an orthogonal $Q$ matrix and an accurate $R$:

In [8]:

QRH = qr(A)
QH, RH = QRH.Q, QRH.R
semilogy(abs.(diag(RH)), "s", mfc="none") # need abs because householder has arbitrary phase factor in R
semilogy(diag(RM), "rx")
legend(["Householder QR","modified GS"], loc="upper right")
ylabel(L"r_{jj}")
xlabel(L"j")

println("L1 NORM of Householder Q'Q - I: ", opnorm(QH'*QH - I, 1))

L1 NORM of Householder Q'Q - I: 7.997590220285011e-15

On top of that, Householder is faster than Gram–Schmidt, so what's not to love? But modified Gram–Schmidt is still useful to know since it is the basis of the Arnoldi and GMRES iterative algorithms for large (usually sparse) systems.