LALQMR for complex numbers

Author: Pawel Latawiec

src/IterativeSolvers.jl

@@ -27,6 +27,7 @@ include("gmres.jl")
 include("chebyshev.jl")
 include("idrs.jl")
 include("qmr.jl")
+include("lalqmr.jl")
 
 # Eigensolvers
 include("simple.jl")

src/lalqmr.jl

@@ -0,0 +1,203 @@
+import Base: iterate
+
+export lalqmr, lalqmr!
+
+mutable struct LALQMRIterable{T, xT, rT, lanczosT}
+    x::xT
+
+    lanczos::lanczosT
+    resnorm::rT
+    tol::rT
+    maxiter::Int
+
+    t::Vector{T}
+    R::UpperTriangular{T, Matrix{T}}
+
+    G::Vector{Givens{T}}
+
+    d::Matrix{T}
+end
+
+function lalqmr_iterable!(
+    x, A, b;
+    abstol::Real = zero(real(eltype(b))),
+    reltol::Real = sqrt(eps(real(eltype(b)))),
+    maxiter::Int = size(A, 2),
+    initially_zero::Bool = false,
+    kwargs...
+)
+    T = eltype(x)
+    r = copy(b)
+    if !initially_zero
+        r -= A*x
+    end
+    resnorm = norm(r)
+    v = r/resnorm
+    w = copy(v)
+
+    lanczos = LookAheadLanczosDecomp(
+        A, v, w;
+        vw_normalized = true,
+        kwargs...
+    )
+
+    R = UpperTriangular(Matrix{T}(undef, 0, 0))
+    # Givens rotations
+    G = Vector{Givens{T}}()
+
+    # Projection operators
+    D = similar(x, size(x, 1), 0)
+
+    t = Vector{T}(undef, 1)
+    t[1] = resnorm
+
+    tolerance = max(reltol * resnorm, abstol)
+
+    return LALQMRIterable(
+        x,
+        lanczos, resnorm, tolerance,
+        maxiter,
+        t, R,
+        G, D
+    )
+end
+
+converged(q::LALQMRIterable) = q.resnorm ≤ q.tol
+start(::LALQMRIterable) = 1
+done(q::LALQMRIterable, iteration::Int) = iteration > q.maxiter || converged(q)
+
+function iterate(q::LALQMRIterable, n::Int=start(q))
+    # Check for termination first
+    if done(q, n)
+        return nothing
+    end
+
+    iterate(q.lanczos, n)
+
+    # Eq. 6.2, update QR factorization of L
+    Llastcol = q.lanczos.L[:, end]
+    for g in q.G
+        Llastcol = g*Llastcol
+    end
+    gend, r = givens(Llastcol, n, n+1)
+    push!(q.G, gend)
+    Llastcol[end-1] = r
+    Llastcol[end] = 0
+    q.R = UpperTriangular([[q.R; fill(0, 1, n-1)] Llastcol[1:end-1]])
+    @assert q.R[:, end] ≈ Llastcol[1:end-1]
+
+
+    # Eq. 6.2, update t
+    push!(q.t, 0)
+    q.t = gend * q.t
+  
+    # Eq. 6.3, calculate projection
+    # Dn = [Vn Un^-1 Rn^-1]
+    # => Dn Rn Un = Vn
+    RU = q.R*q.lanczos.U
+    d = q.lanczos.V[:, end-1]
+    for i in 1:size(RU, 1)-1
+        axpy!(-RU[i, end], q.d[:, i], d)
+    end
+    d = d / RU[end, end]
+    q.d = [q.d d]
+
+    # iterate x_n = x_n-1 + d_n τ_n
+    axpy!(q.t[end-1], d, q.x)
+
+    # Eq. 4.12, Freund 1990
+    q.resnorm = q.resnorm * abs(gend.s) * sqrt(n+1)/sqrt(n)
+
+    return q.resnorm, n + 1
+end
+
+"""
+    lalqmr(A, b; kwargs...) -> x, [history]
+
+Same as [`lalqmr!`](@ref), but allocates a solution vector `x` initialized with zeros.
+"""
+lalqmr(A, b; kwargs...) = lalqmr!(zerox(A, b), A, b; initially_zero = true, kwargs...)
+
+"""
+    lalqmr!(x, A, b; kwargs...) -> x, [history]
+
+Solves the problem ``Ax = b`` with the Quasi-Minimal Residual (QMR) method with Look-ahead. See [`LookAheadLanczosDecomp`](@ref).
+
+# Arguments
+- `x`: Initial guess, will be updated in-place;
+- `A`: linear operator;
+- `b`: right-hand side.
+
+## Keywords
+
+- `initally_zero::Bool`: If `true` assumes that `iszero(x)` so that one
+  matrix-vector product can be saved when computing the initial residual
+  vector;
+- `maxiter::Int = size(A, 2)`: maximum number of iterations;
+- `abstol::Real = zero(real(eltype(b)))`,
+  `reltol::Real = sqrt(eps(real(eltype(b))))`: absolute and relative
+  tolerance for the stopping condition
+  `|r_k| ≤ max(reltol * |r_0|, abstol)`, where `r_k = A * x_k - b`
+- `log::Bool`: keep track of the residual norm in each iteration;
+- `verbose::Bool`: print convergence information during the iteration.
+- `max_block_size`: maximum block size during look-ahead process
+- `max_memory`: maximum allowed memory for look-ahead process
+
+# Return values
+
+**if `log` is `false`**
+
+- `x`: approximate solution.
+
+**if `log` is `true`**
+
+- `x`: approximate solution;
+
+- `history`: convergence history.
+
+[^Freund1990]:
+    Freund, W. R., & Nachtigal, N. M. (1990). QMR : for a Quasi-Minimal
+    Residual Linear Method Systems. (December).
+"""
+function lalqmr!(
+    x, A, b;
+    abstol::Real = zero(real(eltype(b))),
+    reltol::Real = sqrt(eps(real(eltype(b)))),
+    maxiter::Int = size(A, 2),
+    log::Bool = false,
+    initially_zero::Bool = false,
+    verbose::Bool = false,
+    kwargs...
+)
+    history = ConvergenceHistory(partial = !log)
+    history[:abstol] = abstol
+    history[:reltol] = reltol
+    log && reserve!(history, :resnorm, maxiter)
+
+    iterable = lalqmr_iterable!(
+        x, A, b;
+        abstol=abstol,
+        reltol=reltol,
+        maxiter=maxiter,
+        initially_zero=initially_zero,
+        log=log,
+        verbose=verbose
+    )
+
+    verbose && @printf("=== qmr ===\n%4s\t%7s\n","iter","resnorm")
+
+    for (iteration, residual) = enumerate(iterable)
+        if log
+            nextiter!(history)
+            push!(history, :resnorm, residual)
+        end
+
+        verbose && @printf("%3d\t%1.2e\n", iteration, residual)
+    end
+
+    verbose && println()
+    log && setconv(history, converged(iterable))
+    log && shrink!(history)
+
+    log ? (x, history) : x
+end

src/qmr.jl

@@ -266,7 +266,6 @@ function qmr!(x, A, b;
               abstol::Real = zero(real(eltype(b))),
               reltol::Real = sqrt(eps(real(eltype(b)))),
               maxiter::Int = size(A, 2),
-              lookahead::Bool = false,
               log::Bool = false,
               initially_zero::Bool = false,
               verbose::Bool = false)

test/lalqmr.jl

@@ -10,17 +10,19 @@ using SparseArrays
 #LALQMR
 @testset "LALQMR" begin
 
-rng = Random.MersenneTwister(123)
+Random.seed!(1234321)
 n = 10
 
 @testset "Matrix{$T}" for T in (Float32, Float64, ComplexF32, ComplexF64)
-    A = rand(rng, T, n, n)
-    b = rand(rng, T, n)
+    A = rand(T, n, n)
+    b = rand(T, n)
     F = lu(A)
     reltol = √eps(real(T))
 
+    # Test optimality condition: residual should be non-increasing
     x, history = lalqmr(A, b, log=true, maxiter=10, reltol=reltol);
     @test isa(history, ConvergenceHistory)
+    @test all(diff(history[:resnorm]) .<= 0.0)
 
     # Left exact preconditioner
     #x, history = lalqmr(A, b, Pl=F, maxiter=1, reltol=reltol, log=true)
@@ -34,13 +36,14 @@ n = 10
 end
 
 @testset "SparseMatrixCSC{$T, $Ti}" for T in (Float64, ComplexF64), Ti in (Int64, Int32)
-    A = sprand(rng, T, n, n, 0.5) + n * I
-    b = rand(rng, T, n)
+    A = sprand(T, n, n, 0.5) + I
+    b = rand(T, n)
     F = lu(A)
     reltol = √eps(real(T))
 
+    # Test optimality condition: residual should be non-increasing
     x, history = lalqmr(A, b, log = true, maxiter = 10);
-    @test norm(A * x - b) / norm(b) ≤ reltol
+    @test all(diff(history[:resnorm]) .<= 0.0)
 
     # Left exact preconditioner
     #x, history = lalqmr(A, b, Pl=F, maxiter=1, log=true)
@@ -53,29 +56,10 @@ end
     #@test norm(A * x - b) / norm(b) ≤ reltol
 end
 
-@testset "Block Creation {$T}" for T in (Float32, Float64, ComplexF32, ComplexF64)
-    # Guaranteed to create blocks during Lanczos process
-    # This satisfies the condition that in the V-W sequence, the first
-    # iterates are orthogonal: <Av - v<A, v>, Atv - v<At, v>> under transpose inner product
-    # These do _not_ work for regular qmr
-    dl = fill(one(T), n-1)
-    du = fill(one(T), n-1)
-    d = fill(one(T), n)
-    dl[1] = -1
-    A = Tridiagonal(dl, d, du)
-    b = fill(zero(T), n)
-    b[2] = 1.0
-
-    reltol = √eps(real(T))
-    
-    x, history = lalqmr(A, b, log = true)
-    @test norm(A * x - b) / norm(b) ≤ reltol
-end
-
 @testset "Linear operator defined as a function" begin
-    A = LinearMap(identity, identity, 100; ismutating=false)
+    A = LinearMap(cumsum!, 100; ismutating=true)
     b = rand(100)
-    reltol = 1e-6
+    reltol = 1e-5
 
     x = lalqmr(A, b; reltol=reltol, maxiter=2000)
     @test norm(A * x - b) / norm(b) ≤ reltol

test/runtests.jl

@@ -30,6 +30,7 @@ include("qmr.jl")
 
 # Look-Ahead Lanczos
 include("lal.jl")
+include("lalqmr.jl")
 
 #Chebyshev
 include("chebyshev.jl")