Skip to content

Commit 85c65b7

Browse files
fredrikekrefredrikekre
committed
Add InverseMap for lazy inverse of a linear map.
1 parent 99bbc0e commit 85c65b7

File tree

9 files changed

+285
-1
lines changed

9 files changed

+285
-1
lines changed

‎Project.toml

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
name = "LinearMaps"
22
uuid = "7a12625a-238d-50fd-b39a-03d52299707e"
3-
version = "3.7.0"
3+
version = "3.8.0"
44

55
[deps]
66
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

‎docs/src/history.md

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,14 @@
11
# Version history
22

3+
## What's new in v3.8
4+
5+
* A new map called [`InverseMap`](@ref) is introduced. Letting an `InverseMap` act on a
6+
vector is equivalent to solving the linear system, i.e. `InverseMap(A) * b` is the same as
7+
`A \ b`. The default solver is `ldiv!`, but can be specified with the `solver` keyword
8+
argument to the constructor (see the docstring for details). Note that `A` must be
9+
compatible with the solver: `A` can, for example, be a factorization, or another
10+
`LinearMap` in combination with an iterative solver.
11+
312
## What's new in v3.7
413

514
* `mul!(M::AbstractMatrix, A::LinearMap, s::Number, a, b)` methods are provided, mimicking

‎docs/src/index.md

Lines changed: 46 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -111,6 +111,52 @@ KrylovKit.eigsolve(Δ, size(Δ, 1), 3, :LR)
111111

112112
leveraging the fact that objects of type `L <: LinearMap` are callable.
113113

114+
### Inverse map with conjugate gradient
115+
116+
The `InverseMap` type can be used to lazily represent the inverse of an operator.
117+
When this map acts on a vector the linear system is solved. This can be used to solve a
118+
system of the form ``Sx = (C A^{-1} B) x = b`` without explicitly computing ``A^{-1}``
119+
(see for example solving a linear system using the
120+
[Schur complement](https://en.wikipedia.org/wiki/Schur_complement#Application_to_solving_linear_equations)).
121+
122+
```julia
123+
using LinearMaps, IterativeSolvers
124+
125+
A = [2.0 1.5 0.0
126+
1.5 3.0 0.0
127+
0.0 0.0 4.0]
128+
B = [2.0 0.0
129+
0.0 1.0
130+
0.0 0.0]
131+
C = B'
132+
b = [2.0, 3.0]
133+
134+
# Use IterativeSolvers.cg! to solve the system with 0 as the initial guess
135+
linsolve = (x, A, b) -> IterativeSolvers.cg!(fill!(x, 0), A, b)
136+
137+
# Construct the linear map S
138+
S = C * InverseMap(A; solver=linsolve) * B
139+
140+
# Solve the system
141+
IterativeSolvers.cg(S, b)
142+
```
143+
144+
In every CG iteration the linear map `S` will act on a vector `v`. Since `S` is a composed
145+
linear map, `S * b` is roughly equivalent to
146+
147+
```julia
148+
# Apply first linear map B to v
149+
tmp1 = B * v
150+
# Apply second linear map: solve linear system with vector tmp1 as RHS
151+
tmp2 = A \ tmp1
152+
# Apply third linear map C to tmp2
153+
result = C * tmp2
154+
```
155+
156+
i.e. inside the CG solver for solving `Sx = b` we use CG to solve another inner linear
157+
system.
158+
159+
114160
## Philosophy
115161

116162
Several iterative linear algebra methods such as linear solvers or eigensolvers

‎docs/src/types.md

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -96,6 +96,15 @@ LinearMaps.FillMap
9696

9797
Type for representing linear maps that are embedded in larger zero maps.
9898

99+
100+
### `InverseMap`
101+
102+
Type for lazy inverse of another linear map.
103+
104+
```@docs
105+
LinearMaps.InverseMap
106+
```
107+
99108
## Methods
100109

101110
### Multiplication methods

‎src/LinearMaps.jl

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@ module LinearMaps
33
export LinearMap
44
export ⊗, kronsum, ⊕
55
export FillMap
6+
export InverseMap
67

78
using LinearAlgebra
89
import LinearAlgebra: mul!
@@ -344,6 +345,7 @@ include("embeddedmap.jl") # embedded linear maps
344345
include("conversion.jl") # conversion of linear maps to matrices
345346
include("show.jl") # show methods for LinearMap objects
346347
include("getindex.jl") # getindex functionality
348+
include("inversemap.jl")
347349

348350
"""
349351
LinearMap(A::LinearMap; kwargs...)::WrappedMap

‎src/inversemap.jl

Lines changed: 54 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,54 @@
1+
struct InverseMap{T, F, S} <: LinearMap{T}
2+
A::F
3+
ldiv!::S
4+
end
5+
6+
"""
7+
InverseMap(A; solver = ldiv!)
8+
9+
Lazy inverse of `A` such that `InverseMap(A) * x` is the same as `A \\ x`. Letting an
10+
`InverseMap` act on a vector thus requires solving a linear system.
11+
12+
A solver function can be passed with the `solver` keyword argument. The solver should
13+
be of the form `f(y, A, x)` where `A` is the wrapped map, `x` the right hand side, and
14+
`y` a preallocated output vector in which the result should be stored. The default solver
15+
is `LinearAlgebra.ldiv!`.
16+
17+
Note that `A` must be compatible with the solver function. `A` can, for example`, be a
18+
factorization of a matrix, or another `LinearMap` (in combination with an iterative solver
19+
such as conjugate gradient).
20+
21+
# Examples
22+
```julia
23+
julia> using LinearMaps, LinearAlgebra
24+
25+
julia> A = rand(2, 2); b = rand(2);
26+
27+
julia> InverseMap(lu(A)) * b
28+
2-element Vector{Float64}:
29+
1.0531895201271027
30+
-0.4718540250893251
31+
32+
julia> A \\ b
33+
2-element Vector{Float64}:
34+
1.0531895201271027
35+
-0.4718540250893251
36+
```
37+
"""
38+
function InverseMap(A::F; solver::S=LinearAlgebra.ldiv!) where {F, S}
39+
T = eltype(A)
40+
InverseMap{T,F,S}(A, solver)
41+
end
42+
43+
Base.size(imap::InverseMap) = size(imap.A)
44+
Base.transpose(imap::InverseMap) = InverseMap(transpose(imap.A); solver=imap.ldiv!)
45+
Base.adjoint(imap::InverseMap) = InverseMap(adjoint(imap.A); solver=imap.ldiv!)
46+
47+
LinearAlgebra.issymmetric(imap::InverseMap) = issymmetric(imap.A)
48+
LinearAlgebra.ishermitian(imap::InverseMap) = ishermitian(imap.A)
49+
LinearAlgebra.isposdef(imap::InverseMap) = isposdef(imap.A)
50+
51+
function LinearMaps._unsafe_mul!(y::AbstractVecOrMat, imap::InverseMap, x::AbstractVector)
52+
imap.ldiv!(y, imap.A, x)
53+
return y
54+
end

‎test/Project.toml

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,7 @@
22
Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
33
BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
44
InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
5+
IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
56
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
67
Quaternions = "94ee1d12-ae83-5a48-8b1c-48b8ff168ae0"
78
SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

‎test/inversemap.jl

Lines changed: 161 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,161 @@
1+
using Test, LinearMaps, LinearAlgebra, SparseArrays, IterativeSolvers
2+
3+
@testset "inversemap" begin
4+
# First argument to cg! doubles as the initial guess so we make sure it is 0 instead
5+
# of potential garbage since we don't control the allocation of the output vector.
6+
cgz! = (x, A, b) -> IterativeSolvers.cg!(fill!(x, 0), A, b)
7+
# Dense test data
8+
A = rand(10, 10) + 5I; A = A'A
9+
B = rand(10, 10)
10+
b = rand(10)
11+
# LU
12+
@test A \ b ≈ InverseMap(lu(A)) * b
13+
# Cholesky
14+
@test A \ b ≈ InverseMap(cholesky(A)) * b
15+
# Specify solver
16+
@test A \ b ≈ InverseMap(A; solver=cgz!) * b atol=1e-8
17+
# Composition with other maps
18+
@test A \ B * b ≈ InverseMap(lu(A)) * B * b
19+
# Composition: make sure solvers called with vector B * b and not matrix B
20+
my_ldiv! = (y, A, x) -> begin
21+
@test x isa AbstractVector
22+
@test x ≈ B * b
23+
return ldiv!(y, lu(A), x)
24+
end
25+
@test A \ B * b ≈ InverseMap(A; solver=my_ldiv!) * B * b
26+
27+
# Sparse test data
28+
A = sprand(10, 10, 0.2) + 5I; A = A'A
29+
B = sprand(10, 10, 0.5)
30+
@test A \ b ≈ InverseMap(lu(A)) * b
31+
# Cholesky (CHOLMOD doesn't support inplace ldiv!)
32+
my_ldiv! = (y, A, x) -> copy!(y, A \ x)
33+
@test A \ b ≈ InverseMap(cholesky(A); solver=my_ldiv!) * b
34+
# Specify solver
35+
@test A \ b ≈ InverseMap(A; solver=cgz!) * b atol=1e-8
36+
# Composition with other maps
37+
@test A \ (B * b) ≈ InverseMap(lu(A)) * B * b
38+
# Composition: make sure solver is called with vector B * b and not matrix B
39+
my_ldiv! = (y, A, x) -> begin
40+
@test x isa AbstractVector
41+
@test x ≈ B * b
42+
return ldiv!(y, lu(A), x)
43+
end
44+
@test A \ (B * b) ≈ InverseMap(A; solver=my_ldiv!) * B * b
45+
46+
# Interface testing: note that not all combinations of factorization and
47+
# is(symmetric|hermitian|posdef) and transpose/adjoint are supported by LinearAlgebra,
48+
# so we test this for a custom type just to make sure the call is forwarded correctly
49+
# and then run some tests for supported combinations.
50+
struct TestMap{T} <: LinearMap{T}
51+
A::Matrix{T}
52+
end
53+
Base.size(tm::TestMap) = size(tm.A)
54+
Base.transpose(tm::TestMap) = transpose(tm.A)
55+
Base.adjoint(tm::TestMap) = adjoint(tm.A)
56+
LinearAlgebra.issymmetric(tm::TestMap) = issymmetric(tm.A)
57+
LinearAlgebra.ishermitian(tm::TestMap) = ishermitian(tm.A)
58+
LinearAlgebra.isposdef(tm::TestMap) = isposdef(tm.A)
59+
A = [5.0 2.0; 2.0 4.0]
60+
itm = InverseMap(TestMap(A))
61+
@test size(itm) == size(A)
62+
@test transpose(itm).A === transpose(A)
63+
@test adjoint(itm).A === adjoint(A)
64+
@test issymmetric(itm)
65+
@test ishermitian(itm)
66+
@test isposdef(itm)
67+
## Real symmetric (and Hermitian)
68+
A = Float64[3 2; 2 5]; x = rand(2)
69+
### Wrapping a matrix and factorize in solver
70+
iA = InverseMap(A; solver=(y, A, x)->ldiv!(y, cholesky(A), x))
71+
@test ishermitian(A) == ishermitian(iA) == true
72+
@test issymmetric(A) == issymmetric(iA) == true
73+
@test isposdef(A) == isposdef(iA) == true
74+
@test A \ x ≈ iA * x
75+
if VERSION >= v"1.8.0-"
76+
@test transpose(A) \ x ≈ transpose(iA) * x
77+
@test adjoint(A) \ x ≈ adjoint(iA) * x
78+
end
79+
### Wrapping a factorization
80+
iA = InverseMap(cholesky(A))
81+
# @test ishermitian(A) == ishermitian(iA) == true
82+
# @test issymmetric(A) == issymmetric(iA) == true
83+
@test isposdef(A) == isposdef(iA) == true
84+
@test A \ x ≈ iA * x
85+
# @test transpose(A) \ x ≈ transpose(iA) * x
86+
if VERSION >= v"1.7.0"
87+
@test adjoint(A) \ x ≈ adjoint(iA) * x
88+
end
89+
## Real non-symmetric
90+
A = Float64[3 2; -2 5]; x = rand(2)
91+
### Wrapping a matrix and factorize in solver
92+
iA = InverseMap(A; solver=(y, A, x)->ldiv!(y, lu(A), x))
93+
@test ishermitian(A) == ishermitian(iA) == false
94+
@test issymmetric(A) == issymmetric(iA) == false
95+
@test isposdef(A) == isposdef(iA) == false
96+
@test A \ x ≈ iA * x
97+
@test transpose(A) \ x ≈ transpose(iA) * x
98+
@test adjoint(A) \ x ≈ adjoint(iA) * x
99+
### Wrapping a factorization
100+
iA = InverseMap(lu(A))
101+
# @test ishermitian(A) == ishermitian(iA) == true
102+
# @test issymmetric(A) == issymmetric(iA) == true
103+
# @test isposdef(A) == isposdef(iA) == true
104+
@test A \ x ≈ iA * x
105+
@test transpose(A) \ x ≈ transpose(iA) * x
106+
@test adjoint(A) \ x ≈ adjoint(iA) * x
107+
## Complex Hermitian
108+
A = ComplexF64[3 2im; -2im 5]; x = rand(ComplexF64, 2)
109+
### Wrapping a matrix and factorize in solver
110+
iA = InverseMap(A; solver=(y, A, x)->ldiv!(y, cholesky(A), x))
111+
@test ishermitian(A) == ishermitian(iA) == true
112+
@test issymmetric(A) == issymmetric(iA) == false
113+
@test isposdef(A) == isposdef(iA) == true
114+
@test A \ x ≈ iA * x
115+
if VERSION >= v"1.8.0-"
116+
@test transpose(A) \ x ≈ transpose(iA) * x
117+
@test adjoint(A) \ x ≈ adjoint(iA) * x
118+
end
119+
### Wrapping a factorization
120+
iA = InverseMap(cholesky(A))
121+
# @test ishermitian(A) == ishermitian(iA) == true
122+
# @test issymmetric(A) == issymmetric(iA) == true
123+
@test isposdef(A) == isposdef(iA) == true
124+
@test A \ x ≈ iA * x
125+
# @test transpose(A) \ x ≈ transpose(iA) * x
126+
if VERSION >= v"1.7.0"
127+
@test adjoint(A) \ x ≈ adjoint(iA) * x
128+
end
129+
## Complex non-Hermitian
130+
A = ComplexF64[3 2im; 3im 5]; x = rand(ComplexF64, 2)
131+
### Wrapping a matrix and factorize in solver
132+
iA = InverseMap(A; solver=(y, A, x)->ldiv!(y, lu(A), x))
133+
@test ishermitian(A) == ishermitian(iA) == false
134+
@test issymmetric(A) == issymmetric(iA) == false
135+
@test isposdef(A) == isposdef(iA) == false
136+
@test A \ x ≈ iA * x
137+
@test transpose(A) \ x ≈ transpose(iA) * x
138+
@test adjoint(A) \ x ≈ adjoint(iA) * x
139+
### Wrapping a factorization
140+
iA = InverseMap(lu(A))
141+
# @test ishermitian(A) == ishermitian(iA) == true
142+
# @test issymmetric(A) == issymmetric(iA) == true
143+
# @test isposdef(A) == isposdef(iA) == true
144+
@test A \ x ≈ iA * x
145+
@test transpose(A) \ x ≈ transpose(iA) * x
146+
@test adjoint(A) \ x ≈ adjoint(iA) * x
147+
148+
# Example from https://www.dealii.org/current/doxygen/deal.II/step_20.html#SolvingusingtheSchurcomplement
149+
M = sparse(2.0I, 10, 10) + sprand(10, 10, 0.1); M = M'M
150+
iM = InverseMap(M; solver=cgz!)
151+
B = sparse(5.0I, 10, 5)
152+
F = rand(10)
153+
G = rand(5)
154+
## Solve using Schur complement
155+
G′ = B' * iM * F - G
156+
S = B' * iM * B
157+
P = IterativeSolvers.cg(S, G′)
158+
U = IterativeSolvers.cg(M, F - B * P)
159+
## Solve using standard method and compare
160+
@test [M B; B' 0I] \ [F; G] ≈ [U; P] atol=1e-10
161+
end

‎test/runtests.jl

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -37,3 +37,5 @@ include("nontradaxes.jl")
3737
include("embeddedmap.jl")
3838

3939
include("getindex.jl")
40+
41+
include("inversemap.jl")

0 commit comments

Comments
 (0)