Add InverseMap for lazy inverse of a linear map.

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,52 @@ KrylovKit.eigsolve(Δ, size(Δ, 1), 3, :LR)
 
 leveraging the fact that objects of type `L <: LinearMap` are callable.
 
+### Inverse map with conjugate gradient
+
+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 the linear map `S` will act on a vector `v`. Since `S` is a composed
+linear map, `S * v` is roughly equivalent to
+
+```julia
+# Apply first linear map B to v
+tmp1 = B * v
+# 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,59 @@
+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)
+
+# Two separate methods to deal with method ambiguities
+function _unsafe_mul!(y::AbstractVector, imap::InverseMap, x::AbstractVector)
+    imap.ldiv!(y, imap.A, x)
+    return y
+end
+function _unsafe_mul!(y::AbstractMatrix, imap::InverseMap, x::AbstractMatrix)
+    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,203 @@
+using Test, LinearMaps, LinearAlgebra, SparseArrays, IterativeSolvers
+
+@testset "inversemap" begin
+    # First argument to cg!/gmres! 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, 0), A, b)
+    gmresz! = (x, A, b) -> IterativeSolvers.gmres!(fill!(x, 0), 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
+    ## 3- and 5-arg mul!
+    iA = InverseMap(lu(A))
+    y = zeros(size(A, 1))
+    mul!(y, iA, b)
+    @test A \ b ≈ y
+    Y = zeros(size(A))
+    mul!(Y, iA, B)
+    @test A \ B ≈ Y
+    y = rand(size(A, 1))
+    yc = copy(y)
+    α, β = 1.2, 3.4
+    mul!(y, iA, b, α, β)
+    @test A \ b * α + yc * β ≈ y
+    Y = rand(size(A)...)
+    Yc = copy(Y)
+    mul!(Y, iA, B, α, β)
+    @test A \ B * α + Yc * β ≈ Y
+
+    # 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
+    ## 3- and 5-arg mul!
+    iA = InverseMap(lu(A))
+    y = zeros(size(A, 1))
+    mul!(y, iA, b)
+    @test A \ b ≈ y
+    y = rand(size(A, 1))
+    yc = copy(y)
+    α, β = 1.2, 3.4
+    mul!(y, iA, b, α, β)
+    @test A \ b * α + yc * β ≈ y
+
+    # Combine with another LinearMap
+    A = LinearMap(cumsum, 10, 10)
+    iA = InverseMap(A; solver=gmresz!)
+    y = zeros(size(A, 1))
+    mul!(y, iA, b)
+    @test IterativeSolvers.gmres(A, b) ≈ iA * b ≈ y
+    y = rand(size(A, 1))
+    yc = copy(y)
+    α, β = 1.2, 3.4
+    mul!(y, iA, b, α, β)
+    @test IterativeSolvers.gmres(A, b * α) + β * yc ≈ iA * b * α + yc * β ≈ y
+
+    # Interface testing: note that not all combinations of factorization and
+    # is(symmetric|hermitian|posdef) and transpose/adjoint are supported by LinearAlgebra,
+    # so we test this for a custom type just to make sure the call is forwarded correctly
+    # and then run some tests for supported combinations.
+    struct TestMap{T} <: LinearMap{T}
+        A::Matrix{T}
+    end
+    Base.size(tm::TestMap) = size(tm.A)
+    Base.transpose(tm::TestMap) = transpose(tm.A)
+    Base.adjoint(tm::TestMap) = adjoint(tm.A)
+    LinearAlgebra.issymmetric(tm::TestMap) = issymmetric(tm.A)
+    LinearAlgebra.ishermitian(tm::TestMap) = ishermitian(tm.A)
+    LinearAlgebra.isposdef(tm::TestMap) = isposdef(tm.A)
+    A = [5.0 2.0; 2.0 4.0]
+    itm = InverseMap(TestMap(A))
+    @test size(itm) == size(A)
+    @test transpose(itm).A === transpose(A)
+    @test adjoint(itm).A === adjoint(A)
+    @test issymmetric(itm)
+    @test ishermitian(itm)
+    @test isposdef(itm)
+    ## Real symmetric (and Hermitian)
+    A = Float64[3 2; 2 5]; x = rand(2)
+    ### Wrapping a matrix and factorize in solver
+    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 isposdef(A) == isposdef(iA) == true
+    @test A \ x ≈ iA * x
+    if VERSION >= v"1.8.0-"
+        @test transpose(A) \ x ≈ transpose(iA) * x
+        @test adjoint(A) \ x ≈ adjoint(iA) * x
+    end
+    ### Wrapping a factorization
+    iA = InverseMap(cholesky(A))
+    # @test ishermitian(A) == ishermitian(iA) == true
+    # @test issymmetric(A) == issymmetric(iA) == true
+    @test isposdef(A) == isposdef(iA) == true
+    @test A \ x ≈ iA * x
+    # @test transpose(A) \ x ≈ transpose(iA) * x
+    if VERSION >= v"1.7.0"
+        @test adjoint(A) \ x ≈ adjoint(iA) * x
+    end
+    ## Real non-symmetric
+    A = Float64[3 2; -2 5]; x = rand(2)
+    ### Wrapping a matrix and factorize in solver
+    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
+    @test A \ x ≈ iA * x
+    @test transpose(A) \ x ≈ transpose(iA) * x
+    @test adjoint(A) \ x ≈ adjoint(iA) * x
+    ### Wrapping a factorization
+    iA = InverseMap(lu(A))
+    # @test ishermitian(A) == ishermitian(iA) == true
+    # @test issymmetric(A) == issymmetric(iA) == true
+    # @test isposdef(A) == isposdef(iA) == true
+    @test A \ x ≈ iA * x
+    @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)
+    ### Wrapping a matrix and factorize in solver
+    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 isposdef(A) == isposdef(iA) == true
+    @test A \ x ≈ iA * x
+    if VERSION >= v"1.8.0-"
+        @test transpose(A) \ x ≈ transpose(iA) * x
+        @test adjoint(A) \ x ≈ adjoint(iA) * x
+    end
+    ### Wrapping a factorization
+    iA = InverseMap(cholesky(A))
+    # @test ishermitian(A) == ishermitian(iA) == true
+    # @test issymmetric(A) == issymmetric(iA) == true
+    @test isposdef(A) == isposdef(iA) == true
+    @test A \ x ≈ iA * x
+    # @test transpose(A) \ x ≈ transpose(iA) * x
+    if VERSION >= v"1.7.0"
+        @test adjoint(A) \ x ≈ adjoint(iA) * x
+    end
+    ## Complex non-Hermitian
+    A = ComplexF64[3 2im; 3im 5]; x = rand(ComplexF64, 2)
+    ### Wrapping a matrix and factorize in solver
+    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
+    @test A \ x ≈ iA * x
+    @test transpose(A) \ x ≈ transpose(iA) * x
+    @test adjoint(A) \ x ≈ adjoint(iA) * x
+    ### Wrapping a factorization
+    iA = InverseMap(lu(A))
+    # @test ishermitian(A) == ishermitian(iA) == true
+    # @test issymmetric(A) == issymmetric(iA) == true
+    # @test isposdef(A) == isposdef(iA) == true
+    @test A \ x ≈ iA * x
+    @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
+    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")