Computing matrix-function-vector products using Krylov

Krylov Funmv.md

include("Phi.jl")
using Phi: phim_alt, phimc

Arnoldi iteration

The Arnoldi iteratrion algorithm computes a basis for the Krylov subspace

$$ K_m(A,b) = \mathrm{span}{b,Ab,A^2b,\cdots,A^{m-1}b} $$

For some linear operator $A$ and vector $b$. Typically $m$ is much smaller compared to $n$, the dimension of $A$.

The output is the $n \times (m+1)$ (extended) unitary basis vectors $V$ and the $(m+1) \times m$ (extended) upper Heisenberg matrix $H$, which are related by the recurrence formula

$$v_1=b,\quad Av_j = \sum_{i=1}^{j+1}h_{ij}v_i\quad(j = 1,2,\ldots,m)$$.

TODO

1. Add "happy-breakdown" test.

2. More efficient in-place operations.

3. Lanczos for symmetric/Hermitian case.

function arnoldi(A, b::AbstractVector{T}, m::Int) where {T <: Number}
    @assert size(A,1) == size(A,2) == length(b)
    b = b / norm(b)
    
    # Allocate output
    H = zeros(T, m+1, m)
    V = zeros(T, size(A,1), m+1)
    V[:,1] = b
    
    for j = 1:m
        w = A * V[:,j]
        for i = 1:j
            H[i,j] = V[:,i]' * w
            w -= H[i,j] * V[:,i]
        end
        H[j+1,j] = norm(w)
        V[:,j+1] = w / H[j+1,j]
    end
    return V, H
end;

n = 20; m = 5
A = randn(n,n)
b = randn(n)
V,H = arnoldi(A, b, m);

# Check orthogonality
V' * V

6×6 Array{Float64,2}:
  1.0           7.22784e-18  -2.63102e-17  …  -5.63211e-17   8.58313e-17
  7.22784e-18   1.0           2.64428e-19      8.18148e-17  -6.3933e-17 
 -2.63102e-17   2.64428e-19   1.0              5.3139e-17   -7.53746e-17
 -1.94614e-17  -1.83383e-17  -2.7792e-17       1.22171e-17  -2.41003e-17
 -5.63211e-17   8.18148e-17   5.3139e-17       1.0          -3.66553e-17
  8.58313e-17  -6.3933e-17   -7.53746e-17  …  -3.66553e-17   1.0        

# Check recurrence relation
A*V[:,1:m] ≈ V * H

true

Compute matrix-function-vector products using Krylov

The approximation formula is

$$ f(A)b \approx |b| Vf(H)e_1 $$.

where $f$ is an analytic scalar function and $V$ and $H$ are computed using Arnoldi (the non-extended version). Note we still need to evaluate $f$ on the matrix $H$, so it will be up to the user to provide such an implementation (e.g. expm for matrix exponentiation).

function funmv(f, A, b, m)
    V, H = arnoldi(A, b, m)
    return norm(b) * V[:,1:m] * (f(H[1:m,:]))[:,1]
end;

n = 100
A = 1e-3 * randn(n,n)
b = randn(n)
w = expm(A) * b
w_approx = funmv(expm, A, b, m)
err = norm((w - w_approx) ./ w, Inf)

3.0295386161892574e-9

Special handling of phimv using the Sidje formula

"""
    phimv(A,b,k,m) -> [phi_0(A)*b phi_1(A)*b ... phi_k(A)*b]
"""
function phimv(A, b, k, m)
    V, H = arnoldi(A, b, m)
    e1 = zeros(m); e1[1] = 1.0
    F = phimc(H[1:m, :], e1, k) # [phi_0(H)*e1 phi_1(H)*e1 ... phi_k(H)*e1]
    return norm(b) * V[:, 1:m] * F
end;

K = 4
P = phim_alt(A, K)
W = zeros(n, K+1)
for i = 1:K+1
    W[:,i] = P[i] * b
end
W_approx = phimv(A, b, K, m)
err = norm((W - W_approx) ./ W, Inf)

3.05389736883273e-9