|
| 1 | +# LinearMaps.jl |
| 2 | + |
| 3 | +*A Julia package for defining and working with linear maps, also known as linear transformations or linear operators acting on vectors. The only requirement for a LinearMap is that it can act on a vector (by multiplication) efficiently.* |
| 4 | + |
| 5 | +## Installation |
| 6 | + |
| 7 | +`LinearMaps.jl` is a registered package and can be installed via |
| 8 | + |
| 9 | + pkg> add LinearMaps |
| 10 | + |
| 11 | +in package mode, to be entered by typing `]` in the Julia REPL. |
| 12 | + |
| 13 | +## Examples |
| 14 | + |
| 15 | +Let |
| 16 | + |
| 17 | + A = LinearMap(rand(10, 10)) |
| 18 | + B = LinearMap(cumsum, reverse∘cumsum∘reverse, 10) |
| 19 | + |
| 20 | +be a matrix- and function-based linear map, respectively. Then the following code just works, |
| 21 | +indistinguishably from the case when `A` and `B` are both `AbstractMatrix`-typed objects. |
| 22 | + |
| 23 | +``` |
| 24 | +3.0A + 2B |
| 25 | +A*B' |
| 26 | +[A B; B A] |
| 27 | +kron(A, B) |
| 28 | +``` |
| 29 | + |
| 30 | +The `LinearMap` type and corresponding methods combine well with the following packages: |
| 31 | +* [Arpack.jl](https://github.com/JuliaLinearAlgebra/Arpack.jl): iterative eigensolver |
| 32 | + `eigs` and SVD `svds`; |
| 33 | +* [IterativeSolvers.jl](https://github.com/JuliaMath/IterativeSolvers.jl): iterative |
| 34 | + solvers, eigensolvers, and SVD; |
| 35 | +* [KrylovKit.jl](https://github.com/Jutho/KrylovKit.jl): Krylov-based algorithms for linear problems, singular value and eigenvalue problems |
| 36 | +* [TSVD.jl](https://github.com/andreasnoack/TSVD.jl): truncated SVD `tsvd`. |
| 37 | + |
| 38 | +```julia |
| 39 | +using LinearMaps |
| 40 | +import Arpack, IterativeSolvers, KrylovKit, TSVD |
| 41 | + |
| 42 | +# Example 1, 1-dimensional Laplacian with periodic boundary conditions |
| 43 | +function leftdiff!(y::AbstractVector, x::AbstractVector) # left difference assuming periodic boundary conditions |
| 44 | + N = length(x) |
| 45 | + length(y) == N || throw(DimensionMismatch()) |
| 46 | + @inbounds for i=1:N |
| 47 | + y[i] = x[i] - x[mod1(i-1, N)] |
| 48 | + end |
| 49 | + return y |
| 50 | +end |
| 51 | + |
| 52 | +function mrightdiff!(y::AbstractVector, x::AbstractVector) # minus right difference |
| 53 | + N = length(x) |
| 54 | + length(y) == N || throw(DimensionMismatch()) |
| 55 | + @inbounds for i=1:N |
| 56 | + y[i] = x[i] - x[mod1(i+1, N)] |
| 57 | + end |
| 58 | + return y |
| 59 | +end |
| 60 | + |
| 61 | +D = LinearMap(leftdiff!, mrightdiff!, 100; ismutating=true) # by default has eltype(D) = Float64 |
| 62 | + |
| 63 | +Arpack.eigs(D'D; nev=3, which=:SR) # note that D'D is recognized as symmetric => real eigenfact |
| 64 | +Arpack.svds(D; nsv=3) |
| 65 | + |
| 66 | +Σ, L = IterativeSolvers.svdl(D; nsv=3) |
| 67 | + |
| 68 | +TSVD.tsvd(D, 3) |
| 69 | + |
| 70 | +# Example 2, 1-dimensional Laplacian |
| 71 | +A = LinearMap(100; issymmetric=true, ismutating=true) do C, B |
| 72 | + C[1] = -2B[1] + B[2] |
| 73 | + for i in 2:length(B)-1 |
| 74 | + C[i] = B[i-1] - 2B[i] + B[i+1] |
| 75 | + end |
| 76 | + C[end] = B[end-1] - 2B[end] |
| 77 | + return C |
| 78 | +end |
| 79 | + |
| 80 | +Arpack.eigs(-A; nev=3, which=:SR) |
| 81 | + |
| 82 | +# Example 3, 2-dimensional Laplacian |
| 83 | +Δ = kronsum(A, A) |
| 84 | + |
| 85 | +Arpack.eigs(Δ; nev=3, which=:LR) |
| 86 | +KrylovKit.eigsolve(x -> Δ*x, size(Δ, 1), 3, :LR) |
| 87 | +``` |
| 88 | + |
| 89 | +## Philosophy |
| 90 | + |
| 91 | +Several iterative linear algebra methods such as linear solvers or eigensolvers |
| 92 | +only require an efficient evaluation of the matrix-vector product, where the |
| 93 | +concept of a matrix can be formalized / generalized to a linear map (or linear |
| 94 | +operator in the special case of a square matrix). |
| 95 | + |
| 96 | +The LinearMaps package provides the following functionality: |
| 97 | + |
| 98 | +1. A `LinearMap` type that shares with the `AbstractMatrix` type that it |
| 99 | + responds to the functions `size`, `eltype`, `isreal`, `issymmetric`, |
| 100 | + `ishermitian` and `isposdef`, `transpose` and `adjoint` and multiplication |
| 101 | + with a vector using both `*` or the in-place version `mul!`. Linear algebra |
| 102 | + functions that use duck-typing for their arguments can handle `LinearMap` |
| 103 | + objects similar to `AbstractMatrix` objects, provided that they can be |
| 104 | + written using the above methods. Unlike `AbstractMatrix` types, `LinearMap` |
| 105 | + objects cannot be indexed, neither using `getindex` or `setindex!`. |
| 106 | + |
| 107 | +2. A single function `LinearMap` that acts as a general purpose |
| 108 | + constructor (though it is only an abstract type) and allows to construct |
| 109 | + linear map objects from functions, or to wrap objects of type |
| 110 | + `AbstractMatrix` or `LinearMap`. The latter functionality is useful to |
| 111 | + (re)define the properties (`isreal`, `issymmetric`, `ishermitian`, |
| 112 | + `isposdef`) of the existing matrix or linear map. |
| 113 | + |
| 114 | +3. A framework for combining objects of type `LinearMap` and of type |
| 115 | + `AbstractMatrix` using linear combinations, transposition, composition, |
| 116 | + concatenation and Kronecker product/sums, |
| 117 | + where the linear map resulting from these operations is never explicitly |
| 118 | + evaluated but only its matrix-vector product is defined (i.e. lazy |
| 119 | + evaluation). The matrix-vector product is written to minimize memory |
| 120 | + allocation by using a minimal number of temporary vectors. There is full |
| 121 | + support for the in-place version `mul!`, which should be preferred for |
| 122 | + higher efficiency in critical algorithms. In addition, it tries to recognize |
| 123 | + the properties of combinations of linear maps. In particular, compositions |
| 124 | + such as `A'*A` for arbitrary `A` or even `A'*B*C*B'*A` with arbitrary `A` |
| 125 | + and `B` and positive definite `C` are recognized as being positive definite |
| 126 | + and hermitian. In case a certain property of the resulting `LinearMap` |
| 127 | + object is not correctly inferred, the `LinearMap` method can be called to |
| 128 | + redefine the properties. |
0 commit comments