Look-ahead Lanczos Quasi-Minimal Residual

Project.toml

@@ -3,6 +3,7 @@ uuid = "42fd0dbc-a981-5370-80f2-aaf504508153"
 version = "0.9.2"
 
 [deps]
+BlockDiagonals = "0a1fb500-61f7-11e9-3c65-f5ef3456f9f0"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
@@ -12,3 +13,4 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 [compat]
 RecipesBase = "0.6, 0.7, 0.8, 1.0"
 julia = "1.3"
+BlockDiagonals = "0.1"

docs/make.jl

@@ -22,9 +22,12 @@ makedocs(
 			"BiCGStab(l)" => "linear_systems/bicgstabl.md",
 			"IDR(s)" => "linear_systems/idrs.md",
 			"Restarted GMRES" => "linear_systems/gmres.md",
+			"QMR" => "linear_systems/qmr.md",
+			"LALQMR" => "linear_systems/lalqmr.md",
 			"LSMR" => "linear_systems/lsmr.md",
 			"LSQR" => "linear_systems/lsqr.md",
-			"Stationary methods" => "linear_systems/stationary.md"
+			"Stationary methods" => "linear_systems/stationary.md",
+			"LAL" => "linear_systems/lal.md"
 		],
 		"Eigenproblems" => [
 			"Power method" => "eigenproblems/power_method.md",

docs/src/linear_systems/lal.md

@@ -0,0 +1,17 @@
+# [Look-Ahead Lanczos Bi-Orthogonalization](@id LAL)
+
+The Look-ahead Lanczos Bi-orthogonalization process is a generalization of the Lanczos
+bi-orthoganilization (as used in, for instance, [QMR](@ref)) [^Freund1993], [^Freund1994]. The look-ahead process detects
+(near-)singularities during the construction of the Lanczos iterates, known as break-down,
+and skips over them. This is particularly advantageous for linear systems with large
+null-spaces or eigenvalues close to 0.
+
+We provide an iterator interface for the Lanczos decomposition (Freund, 1994),
+where the look-ahead Lanczos process is implemented as a two-term coupled recurrence.
+
+## Usage
+
+```@docs
+IterativeSolvers.LookAheadLanczosDecomp
+IterativeSolvers.LookAheadLanczosDecompLog
+```

docs/src/linear_systems/lalqmr.md

@@ -0,0 +1,30 @@
+# [LALQMR](@id LALQMR)
+
+Look-ahead Lanczos Quasi-minimal Residual (LALQMR) is the look-ahead variant of [QMR](@ref) for solving $Ax = b$ approximately for $x$ where $A$ is a linear operator and $b$ the right-hand side vector. $A$ may be non-symmetric [^Freund1990]. The Krylov subspace is generated via the [Look-ahead Lanczos process](@ref LAL)
+
+## Usage
+
+```@docs
+lalqmr
+lalqmr!
+```
+
+## Implementation details
+This implementation of LALQMR follows [^Freund1994] and [^Freund1993], where we generate the Krylov subspace via Lanczos bi-orthogonalization based on the matrix `A` and its transpose. In the regular Lanczos process, we may encounter singularities during the construction of the Krylov basis. The look-ahead process avoids this by building blocks instead, avoiding the singularities. Therefore, the LALQMR technique is well-suited towards problematic linear systems. Typically, most systems will not have blocks of size larger than 4 made, and most iterations are "regular", or blocks of size 1 [^Freund1994].
+
+For more detail on the implementation see the original paper [^Freund1994]. This implementation follows the paper, with the exception that if a block of maximum size is reached during the Lanczos process, the block is closed instead of restarted.
+
+For a solution vector of size `n`, the memory allocated will be `(12 + max_memory * 5) * n`.
+
+!!! tip
+    LALQMR can be used as an [iterator](@ref Iterators) via `lalqmr_iterable!`. This makes it possible to access the next, current, and previous Krylov basis vectors during the iteration.
+
+## References
+[^Freund1993]:
+Freund, R. W., Gutknecht, M. H., & Nachtigal, N. M. (1993). An Implementation of the Look-Ahead Lanczos Algorithm for Non-Hermitian Matrices. SIAM Journal on Scientific Computing, 14(1), 137–158. https://doi.org/10.1137/0914009
+
+[^Freund1994]:
+Freund, R. W., & Nachtigal, N. M. (1994). An Implementation of the QMR Method Based on Coupled Two-Term Recurrences. SIAM Journal on Scientific Computing, 15(2), 313–337. https://doi.org/10.1137/0915022
+[^Freund1990]:
+
+Freund, W. R., & Nachtigal, N. M. (1990). QMR : for a Quasi-Minimal Residual Linear Method Systems. (December).

docs/src/linear_systems/qmr.md

@@ -1,6 +1,6 @@
 # [QMR](@id QMR)
 
-QMR is a short-recurrence version of GMRES for solving $Ax = b$ approximately for $x$ where $A$ is a linear operator and $b$ the right-hand side vector. $A$ may be non-symmetric.
+QMR is a short-recurrence version of GMRES for solving $Ax = b$ approximately for $x$ where $A$ is a linear operator and $b$ the right-hand side vector. $A$ may be non-symmetric. This is a simplified version of [LALQMR](@ref)
 
 ## Usage
 
@@ -10,17 +10,9 @@ qmr!
 ```
 
 ## Implementation details
-QMR exploits the tridiagonal structure of the Hessenberg matrix. Although QMR is similar to GMRES, where instead of using the Arnoldi process, a pair of biorthogonal vector spaces $V$ and $W$ is constructed via the Lanczos process. It requires that the adjoint of $A$ `adjoint(A)` be available.
+QMR exploits the tridiagonal structure of the Hessenberg matrix. Although QMR is similar to GMRES, where instead of using the Arnoldi process, a pair of biorthogonal vector spaces $V$ and $W$ is constructed via the Lanczos process. It requires that the adjoint of $A$ `adjoint(A)` be available [^Saad2003], [^Freund1990].
 
 QMR enables the computation of $V$ and $W$ via a three-term recurrence. A three-term recurrence for the projection onto the solution vector can also be constructed from these values, using the portion of the last column of the Hessenberg matrix. Therefore we pre-allocate only eight vectors.
 
-For more detail on the implementation see the original paper [^Freund1990] or [^Saad2003].
-
 !!! tip
     QMR can be used as an [iterator](@ref Iterators) via `qmr_iterable!`. This makes it possible to access the next, current, and previous Krylov basis vectors during the iteration.
-
-[^Saad2003]:
-    Saad, Y. (2003). Interactive method for sparse linear system.
-[^Freund1990]:
-    Freund, W. R., & Nachtigal, N. M. (1990). QMR : for a Quasi-Minimal
-    Residual Linear Method Systems. (December).

src/IterativeSolvers.jl

@@ -13,6 +13,10 @@ include("history.jl")
 # Factorizations
 include("hessenberg.jl")
 
+# Krylov subspace methods
+include("limited_memory_matrices.jl")
+include("lal.jl")
+
 # Linear solvers
 include("stationary.jl")
 include("stationary_sparse.jl")
@@ -23,6 +27,7 @@ include("gmres.jl")
 include("chebyshev.jl")
 include("idrs.jl")
 include("qmr.jl")
+include("lalqmr.jl")
 
 # Eigensolvers
 include("simple.jl")