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.
LinearMaps.jl is a registered package and can be installed via
pkg> add LinearMaps
in package mode, to be entered by typing
] in the Julia REPL.
A = LinearMap(rand(10, 10)) B = LinearMap(cumsum, reverse∘cumsum∘reverse, 10)
be a matrix- and function-based linear map, respectively. Then the following code just works, indistinguishably from the case when
B are both
3.0A + 2B A + I A*B' [A B; B A] kron(A, B)
LinearMap type and corresponding methods combine well with the following packages:
- Arpack.jl: iterative eigensolver
- IterativeSolvers.jl: iterative solvers, eigensolvers, and SVD;
- KrylovKit.jl: Krylov-based algorithms for linear problems, singular value and eigenvalue problems
- TSVD.jl: truncated SVD
using LinearMaps import Arpack, IterativeSolvers, KrylovKit, TSVD # Example 1, 1-dimensional Laplacian with periodic boundary conditions function leftdiff!(y::AbstractVector, x::AbstractVector) # left difference assuming periodic boundary conditions N = length(x) length(y) == N || throw(DimensionMismatch()) @inbounds for i=1:N 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) length(y) == N || throw(DimensionMismatch()) @inbounds for i=1:N 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) Σ, L = IterativeSolvers.svdl(D; nsv=3) TSVD.tsvd(D, 3) # Example 2, 1-dimensional Laplacian A = LinearMap(100; issymmetric=true, ismutating=true) do C, B C = -2B + B 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) # Example 3, 2-dimensional Laplacian Δ = kronsum(A, A) Arpack.eigs(Δ; nev=3, which=:LR) KrylovKit.eigsolve(x -> Δ*x, size(Δ, 1), 3, :LR)
In Julia v1.3 and above, the last line can be simplified to
julia KrylovKit.eigsolve(Δ, size(Δ, 1), 3, :LR)`
leveraging the fact that objects of type
L <: LinearMap are callable.
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:
LinearMaptype that shares with the
AbstractMatrixtype that it responds to the functions
adjointand multiplication with a vector using both
*or the in-place version
mul!. Linear algebra functions that use duck-typing for their arguments can handle
LinearMapobjects similar to
AbstractMatrixobjects, provided that they can be written using the above methods. Unlike
LinearMapobjects cannot be indexed, neither using
A single function
LinearMapthat 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
LinearMap. The latter functionality is useful to (re)define the properties (
isposdef) of the existing matrix or linear map.
A framework for combining objects of type
LinearMapand of type
AbstractMatrixusing 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
Band positive definite
Care recognized as being positive definite and hermitian. In case a certain property of the resulting
LinearMapobject is not correctly inferred, the
LinearMapmethod can be called to redefine the properties.