Add InverseMap for lazy inverse of a linear map.

Author: fredrikekre

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,48 @@
+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 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)
+
+function LinearMaps._unsafe_mul!(y::AbstractVecOrMat, imap::InverseMap, x::AbstractVector)
+    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,59 @@
+using Test, LinearMaps, LinearAlgebra, SparseArrays, IterativeSolvers
+
+@testset "inversemap" begin
+    # 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=IterativeSolvers.cg!) * b
+    # 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=IterativeSolvers.cg!) * b
+    # 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
+
+
+    # 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=IterativeSolvers.cg!)
+    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]
+end

test/runtests.jl

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