# LinearMaps.jl

*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.*

## Installation

`LinearMaps.jl`

is a registered package and can be installed via

`pkg> add LinearMaps`

in package mode, which can be entered by typing `]`

in the Julia REPL.

## Examples

The examples show how to define `LinearMap`

s. For more details refer to the section on Types and their constructors. A linear map is defined by the way it acts on vectors and by its dimensions, potentially also by its adjoint action on vectors. Some properties of the linear map can be defined explicitly or implicitly.

Let

```
A = LinearMap(rand(10, 10))
B = LinearMap(cumsum, reverse∘cumsum∘reverse, 10)
```

be a matrix- and a function-based linear map, respectively. For `A`

all properties are inferred from the input matrix; for `B`

the action, the adjoint action and its dimensions (for a single number, the dimension is taken as square) are defined explicitly. Given the above the following code just works, indistinguishably from the case when `A`

and `B`

are both `AbstractMatrix`

-typed objects, yet giving lazy, `LinearMap`

-typed objects whose action on vectors corresponds to the otherwise expected matrix-vector multiplication.

```
3.0A + 2B
A + I
A*B'
[A B; B A]
kron(A, B)
```

The `LinearMap`

type and corresponding methods combine well with the following packages:

- ArnoldiMethods.jl
- Arpack.jl: iterative eigensolver
`eigs`

and SVD`svds`

; - IterativeSolvers.jl: iterative solvers, eigensolvers, and SVD;
- KrylovKit.jl: Krylov-based algorithms for linear problems, singular value and eigenvalue problems
- TSVD.jl: truncated SVD
`tsvd`

.

```
using LinearMaps
import Arpack, IterativeSolvers, KrylovKit, TSVD, ArnoldiMethod
# Example 1, 1-dimensional Laplacian with periodic boundary conditions
function leftdiff!(y::AbstractVector, x::AbstractVector) # left difference assuming periodic boundary conditions
N = length(x)
axes(y) == axes(x) || throw(DimensionMismatch())
@inbounds for i in eachindex(x, y)
y[i] = x[i] - x[mod1(i-1, N)]
end
return y
end
function mrightdiff!(y::AbstractVector, x::AbstractVector) # minus right difference
N = length(x)
axes(y) == axes(x) || throw(DimensionMismatch())
@inbounds for i in eachindex(x, y)
y[i] = x[i] - x[mod1(i+1, N)]
end
return y
end
D = LinearMap(leftdiff!, mrightdiff!, 100; ismutating=true) # by default has eltype(D) = Float64
Arpack.eigs(D'D; nev=3, which=:SR) # note that D'D is recognized as symmetric => real eigenfact
Arpack.svds(D; nsv=3)
ArnoldiMethod.partialschur(D'D; nev=3, which=ArnoldiMethod.SR())
KrylovKit.eigsolve(D'D, 100, 3, :SR)
Σ, L = IterativeSolvers.svdl(D; nsv=3)
TSVD.tsvd(D, 3)
# Example 2, 3 smallest eigenvalues of 1-dimensional Laplacian
A = LinearMap(100; issymmetric=true, ismutating=true) do C, B
C[1] = -2B[1] + B[2]
for i in 2:length(B)-1
C[i] = B[i-1] - 2B[i] + B[i+1]
end
C[end] = B[end-1] - 2B[end]
return C
end
Arpack.eigs(-A; nev=3, which=:SR)
ArnoldiMethod.partialschur(-A; nev=3, which=ArnoldiMethod.SR())
KrylovKit.eigsolve(-A, size(A, 1), 3, :SR)
# Example 3, 2-dimensional Laplacian
Δ = kronsum(A, A)
Arpack.eigs(Δ; nev=3, which=:LR)
ArnoldiMethod.partialeigen(ArnoldiMethod.partialschur(Δ; nev=3, which=ArnoldiMethod.LR())[1])
KrylovKit.eigsolve(Δ, size(Δ, 1), 3, :LR)
```

In the last line above we leverage 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).

```
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

```
# 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 only require an efficient evaluation of the matrix-vector product, where the concept of a matrix can be formalized / generalized to a linear map (or linear operator in the special case of a square matrix).

The LinearMaps package provides the following functionality:

A

`LinearMap`

type that shares with the`AbstractMatrix`

type that it responds to the functions`size`

,`eltype`

,`isreal`

,`issymmetric`

,`ishermitian`

and`isposdef`

,`transpose`

and`adjoint`

and multiplication with a vector using both`*`

or the in-place version`mul!`

. Linear algebra functions that use duck-typing for their arguments can handle`LinearMap`

objects similar to`AbstractMatrix`

objects, provided that they can be written using the above methods. Unlike`AbstractMatrix`

types,`LinearMap`

objects cannot be indexed, neither using`getindex`

or`setindex!`

.A single function

`LinearMap`

that acts as a general purpose constructor (though it is only an abstract type) and allows to construct linear map objects from functions, or to wrap objects of type`AbstractMatrix`

or`LinearMap`

. The latter functionality is useful to (re)define the properties (`isreal`

,`issymmetric`

,`ishermitian`

,`isposdef`

) of the existing matrix or linear map.A framework for combining objects of type

`LinearMap`

and of type`AbstractMatrix`

using linear combinations, transposition, composition, concatenation and Kronecker product/sums, where the linear map resulting from these operations is never explicitly evaluated but only its matrix-vector product is defined (i.e. lazy evaluation). The matrix-vector product is written to minimize memory allocation by using a minimal number of temporary vectors. There is full support for the in-place version`mul!`

, which should be preferred for higher efficiency in critical algorithms. In addition, it tries to recognize the properties of combinations of linear maps. In particular, compositions such as`A'*A`

for arbitrary`A`

or even`A'*B*C*B'*A`

with arbitrary`A`

and`B`

and positive definite`C`

are recognized as being positive definite and hermitian. In case a certain property of the resulting`LinearMap`

object is not correctly inferred, the`LinearMap`

method can be called to redefine the properties.