Add InverseMap for lazy inverse of a linear map.

Author: fredrikekre

Project.toml

@@ -1,6 +1,6 @@
 name = "LinearMaps"
 uuid = "7a12625a-238d-50fd-b39a-03d52299707e"
-version = "3.7.0"
+version = "3.8.0"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

docs/src/history.md

@@ -1,5 +1,14 @@
 # Version history
 
+## What's new in v3.8
+
+* A new map called [`InverseMap`](@ref) is introduced. Letting an `InverseMap` act on a
+  vector is equivalent to solving the linear system, i.e. `InverseMap(A) * b` is the same as
+  `A \ b`. The default solver is `ldiv!`, but can be specified with the `solver` keyword
+  argument to the constructor (see the docstring for details). Note that `A` must be
+  compatible with the solver: `A` can, for example, be a factorization, or another
+  `LinearMap` in combination with an iterative solver.
+
 ## What's new in v3.7
 
 * `mul!(M::AbstractMatrix, A::LinearMap, s::Number, a, b)` methods are provided, mimicking

docs/src/index.md

@@ -111,6 +111,51 @@ KrylovKit.eigsolve(Δ, size(Δ, 1), 3, :LR)
 
 leveraging the fact that objects of type `L <: LinearMap` are callable.
 
+
+The `InverseMap` type can be used to lazily represent the inverse of an operator.
+When this map acts on a vector the linear system is solved. This can be used to solve a
+system of the form ``Sx = (C A^{-1} B) x = b`` without explicitly computing ``A^{-1}``
+(see for example solving a linear system using the
+[Schur complement](https://en.wikipedia.org/wiki/Schur_complement#Application_to_solving_linear_equations)).
+
+```julia
+using LinearMaps, IterativeSolvers
+
+A = [2.0 1.5 0.0
+     1.5 3.0 0.0
+     0.0 0.0 4.0]
+B = [2.0 0.0
+     0.0 1.0
+     0.0 0.0]
+C = B'
+b = [2.0, 3.0]
+
+# Use IterativeSolvers.cg! to solve the system with 0 as the initial guess
+linsolve = (x, A, b) -> IterativeSolvers.cg!(fill!(x, 0), A, b)
+
+# Construct the linear map S
+S = C * InverseMap(A; solver=linsolve) * B
+
+# Solve the system
+IterativeSolvers.cg(S, b)
+```
+
+In every CG iteration `S` acts on the vector `b`. Since `S` is a composed linear map,
+`S * b` is roughly equivalent to
+
+```julia
+# Apply first linear map B to b
+tmp1 = B * b
+# Apply second linear map: solve linear system with vector tmp1 as RHS
+tmp2 = A \ tmp1
+# Apply third linear map C to tmp2
+result = C * tmp2
+```
+
+i.e. inside the CG solver for solving `Sx = b` we use CG to solve another inner linear
+system.
+
+
 ## Philosophy
 
 Several iterative linear algebra methods such as linear solvers or eigensolvers

docs/src/types.md

@@ -96,6 +96,15 @@ LinearMaps.FillMap
 
 Type for representing linear maps that are embedded in larger zero maps.
 
+
+### `InverseMap`
+
+Type for lazy inverse of another linear map.
+
+```@docs
+LinearMaps.InverseMap
+```
+
 ## Methods
 
 ### Multiplication methods

src/LinearMaps.jl

@@ -3,6 +3,7 @@ module LinearMaps
 export LinearMap
 export ⊗, kronsum, ⊕
 export FillMap
+export InverseMap
 
 using LinearAlgebra
 import LinearAlgebra: mul!
@@ -344,6 +345,7 @@ include("embeddedmap.jl") # embedded linear maps
 include("conversion.jl") # conversion of linear maps to matrices
 include("show.jl") # show methods for LinearMap objects
 include("getindex.jl") # getindex functionality
+include("inversemap.jl")
 
 """
     LinearMap(A::LinearMap; kwargs...)::WrappedMap

src/inversemap.jl

@@ -0,0 +1,54 @@
+struct InverseMap{T, F, S} <: LinearMap{T}
+    A::F
+    ldiv!::S
+end
+
+"""
+    InverseMap(A; solver = ldiv!)
+
+Lazy inverse of `A` such that `InverseMap(A) * x` is the same as `A \\ x`. Letting an
+`InverseMap` act on a vector thus requires solving a linear system.
+
+A solver function can be passed with the `solver` keyword argument. The solver should
+be of the form `f(y, A, x)` where `A` is the wrapped map, `x` the right hand side, and
+`y` a preallocated output vector in which the result should be stored. The default solver
+is `LinearAlgebra.ldiv!`.
+
+Note that `A` must be compatible with the solver function. `A` can, for example`, be a
+factorization of a matrix, or another `LinearMap` (in combination with an iterative solver
+such as conjugate gradient).
+
+# Examples
+```julia
+julia> using LinearMaps, LinearAlgebra
+
+julia> A = rand(2, 2); b = rand(2);
+
+julia> InverseMap(lu(A)) * b
+2-element Vector{Float64}:
+  1.0531895201271027
+ -0.4718540250893251
+
+julia> A \\ b
+2-element Vector{Float64}:
+  1.0531895201271027
+ -0.4718540250893251
+```
+"""
+function InverseMap(A::F; solver::S=LinearAlgebra.ldiv!) where {F, S}
+    T = eltype(A)
+    InverseMap{T,F,S}(A, solver)
+end
+
+Base.size(imap::InverseMap) = size(imap.A)
+Base.transpose(imap::InverseMap) = InverseMap(transpose(imap.A); solver=imap.ldiv!)
+Base.adjoint(imap::InverseMap) = InverseMap(adjoint(imap.A); solver=imap.ldiv!)
+
+LinearAlgebra.issymmetric(imap::InverseMap) = issymmetric(imap.A)
+LinearAlgebra.ishermitian(imap::InverseMap) = ishermitian(imap.A)
+LinearAlgebra.isposdef(imap::InverseMap) = isposdef(imap.A)
+
+function LinearMaps._unsafe_mul!(y::AbstractVecOrMat, imap::InverseMap, x::AbstractVecOrMat)
+    imap.ldiv!(y, imap.A, x)
+    return y
+end

test/Project.toml

@@ -2,6 +2,7 @@
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
+IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Quaternions = "94ee1d12-ae83-5a48-8b1c-48b8ff168ae0"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

test/inversemap.jl

@@ -0,0 +1,102 @@
+using Test, LinearMaps, LinearAlgebra, SparseArrays, IterativeSolvers
+
+@testset "inversemap" begin
+    # First argument to cg! doubles as the initial guess so we make sure it is 0 instead
+    # of potential garbage since we don't control the allocation of the output vector.
+    cgz! = (x, A, b) -> IterativeSolvers.cg!(fill!(x), A, b)
+    # Dense test data
+    A = rand(10, 10) + 5I; A = A'A
+    B = rand(10, 10)
+    b = rand(10)
+    # LU
+    @test A \ b ≈ InverseMap(lu(A)) * b
+    # Cholesky
+    @test A \ b ≈ InverseMap(cholesky(A)) * b
+    # Specify solver
+    @test A \ b ≈ InverseMap(A; solver=cgz!) * b atol=1e-8
+    # Composition with other maps
+    @test A \ B * b ≈ InverseMap(lu(A)) * B * b
+    # Composition: make sure solvers called with vector B * b and not matrix B
+    my_ldiv! = (y, A, x) -> begin
+        @test x isa AbstractVector
+        @test x ≈ B * b
+        return ldiv!(y, lu(A), x)
+    end
+    @test A \ B * b ≈ InverseMap(A; solver=my_ldiv!) * B * b
+
+    # Sparse test data
+    A = sprand(10, 10, 0.2) + 5I; A = A'A
+    B = sprand(10, 10, 0.5)
+    @test A \ b ≈ InverseMap(lu(A)) * b
+    # Cholesky (CHOLMOD doesn't support inplace ldiv!)
+    my_ldiv! = (y, A, x) -> copy!(y, A \ x)
+    @test A \ b ≈ InverseMap(cholesky(A); solver=my_ldiv!) * b
+    # Specify solver
+    @test A \ b ≈ InverseMap(A; solver=cgz!) * b atol=1e-8
+    # Composition with other maps
+    @test A \ (B * b) ≈ InverseMap(lu(A)) * B * b
+    # Composition: make sure solver is called with vector B * b and not matrix B
+    my_ldiv! = (y, A, x) -> begin
+        @test x isa AbstractVector
+        @test x ≈ B * b
+        return ldiv!(y, lu(A), x)
+    end
+    @test A \ (B * b) ≈ InverseMap(A; solver=my_ldiv!) * B * b
+
+    # Interface testing: note that not all combinations of factorization and
+    # is(symmetric|hermitian|posdef) and transpose/adjoint are supported, so some tests are
+    # run with just wrapping the matrix and specifying a solver that factorizes the matrix.
+    ## Real symmetric (and Hermitian)
+    A = Float64[3 2; 2 5]; x = rand(2)
+    iA = InverseMap(A; solver=(y, A, x)->ldiv!(y, cholesky(A), x))
+    @test ishermitian(A) == ishermitian(iA) == true
+    @test issymmetric(A) == issymmetric(iA) == true
+    @test transpose(A) \ x ≈ transpose(iA) * x
+    iA = InverseMap(cholesky(A))
+    @test isposdef(A) == isposdef(iA) == true
+    @test adjoint(A) \ x ≈ adjoint(iA) * x
+    ## Real non-symmetric
+    A = Float64[3 2; -2 5]; x = rand(2)
+    iA = InverseMap(A; solver=(y, A, x)->ldiv!(y, lu(A), x))
+    @test ishermitian(A) == ishermitian(iA) == false
+    @test issymmetric(A) == issymmetric(iA) == false
+    @test isposdef(A) == isposdef(iA) == false
+    iA = InverseMap(lu(A))
+    @test transpose(A) \ x ≈ transpose(iA) * x
+    @test adjoint(A) \ x ≈ adjoint(iA) * x
+    ## Complex Hermitian
+    A = ComplexF64[3 2im; -2im 5]; x = rand(ComplexF64, 2)
+    iA = InverseMap(A; solver=(y, A, x)->ldiv!(y, cholesky(A), x))
+    @test ishermitian(A) == ishermitian(iA) == true
+    @test issymmetric(A) == issymmetric(iA) == false
+    @test transpose(A) \ x ≈ transpose(iA) * x
+    iA = InverseMap(cholesky(A))
+    @test isposdef(A) == isposdef(iA) == true
+    @test adjoint(A) \ x ≈ adjoint(iA) * x
+    ## Complex non-Hermitian
+    A = ComplexF64[3 2im; 3im 5]; x = rand(ComplexF64, 2)
+    iA = InverseMap(A; solver=(y, A, x)->ldiv!(y, lu(A), x))
+    @test ishermitian(A) == ishermitian(iA) == false
+    @test issymmetric(A) == issymmetric(iA) == false
+    @test isposdef(A) == isposdef(iA) == false
+    iA = InverseMap(lu(A))
+    @test transpose(A) \ x ≈ transpose(iA) * x
+    @test adjoint(A) \ x ≈ adjoint(iA) * x
+
+    # Example from https://www.dealii.org/current/doxygen/deal.II/step_20.html#SolvingusingtheSchurcomplement
+    M = sparse(2.0I, 10, 10) + sprand(10, 10, 0.1); M = M'M
+    @show M.colptr
+    @show M.rowval
+    @show M.nzval
+    iM = InverseMap(M; solver=cgz!)
+    B = sparse(5.0I, 10, 5)
+    F = rand(10)
+    G = rand(5)
+    ## Solve using Schur complement
+    G′ = B' * iM * F - G
+    S = B' * iM * B
+    P = IterativeSolvers.cg(S, G′)
+    U = IterativeSolvers.cg(M, F - B * P)
+    ## Solve using standard method and compare
+    @test [M B; B' 0I] \ [F; G] ≈ [U; P] atol=1e-8
+end

test/runtests.jl

@@ -37,3 +37,5 @@ include("nontradaxes.jl")
 include("embeddedmap.jl")
 
 include("getindex.jl")
+
+include("inversemap.jl")