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.
- Fully Julia v0.7/v1.0/v1.1 compatible.
- Full support of noncommutative number types such as quaternions.
- Fully Julia v0.7/v1.0 compatible.
- A
convert(SparseMatrixCSC, A::LinearMap)
function, that calls thesparse
matrix generating function.
- Fully Julia v0.7 compatible; dropped compatibility for previous versions of Julia from LinearMaps.jl v2.0.0 on.
- A 5-argument version for
mul!(y, A::LinearMap, x, α=1, β=0)
, which computesy := α * A * x + β * y
and implements the usual 3-argumentmul!(y, A, x)
for the defaultα
andβ
. - Synonymous
convert(Matrix, A::LinearMap)
andconvert(Array, A::LinearMap)
functions, that call theMatrix
constructor and return the matrix representation ofA
. - Multiplication with matrices, interpreted as a block row vector of vectors:
mul!(Y::AbstractArray, A::LinearMap, X::AbstractArray, α=1, β=0)
: appliesA
to each column ofX
and stores the result in-place in the corresponding column ofY
;- for the out-of-place multiplication, the approach is to compute
convert(Matrix, A * X)
; this is equivalent to applyingA
to each column ofX
. In generic code which handles bothA::AbstractMatrix
andA::LinearMap
, the additional call toconvert
is a noop whenA
is a matrix.
- Full compatibility with Arpack.jl's
eigs
andsvds
; previously onlyeigs
was working. For more, nicely collaborating packages see the Example section.
Install with the package manager, i.e., ]add LinearMaps
(or Pkg.add("LinearMaps")
in Julia versions below 0.7).
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 theAbstractMatrix
type that it responds to the functionssize
,eltype
,isreal
,issymmetric
,ishermitian
andisposdef
,transpose
andadjoint
and multiplication with a vector using both*
or the in-place versionmul!
. Linear algebra functions that use duck-typing for their arguments can handleLinearMap
objects similar toAbstractMatrix
objects, provided that they can be written using the above methods. UnlikeAbstractMatrix
types,LinearMap
objects cannot be indexed, neither usinggetindex
orsetindex!
. -
A single method
LinearMap
function 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 typeAbstractMatrix
orLinearMap
. 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 typeAbstractMatrix
using linear combinations, transposition and composition, 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 versionmul!
, 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 asA'*A
for arbitraryA
or evenA'*B*C*B'*A
with arbitraryA
andB
and positive definiteC
are recognized as being positive definite and hermitian. In case a certain property of the resultingLinearMap
object is not correctly inferred, theLinearMap
method can be called to redefine the properties.
-
LinearMap
General purpose method to construct
LinearMap
objects of specific types, as described in the Types section belowLinearMap{T}(A::AbstractMatrix[; issymmetric::Bool, ishermitian::Bool, isposdef::Bool]) LinearMap{T}(A::LinearMap[; issymmetric::Bool, ishermitian::Bool, isposdef::Bool])
Create a
WrappedMap
object that will respond to the methodsissymmetric
,ishermitian
,isposdef
with the values provided by the keyword arguments, and toeltype
with the valueT
. The default values correspond to the result of calling these methods on the argumentA
; in particular{T}
does not need to be specified and is set aseltype(A)
. This allows to use anAbstractMatrix
within theLinearMap
framework and to redefine the properties of an existingLinearMap
.LinearMap{T}(f, [fc = nothing], M::Int, [N::Int = M]; ismutating::Bool, issymmetric::Bool, ishermitian::Bool, isposdef::Bool])
Create a
FunctionMap
instance that wraps an object describing the action of the linear map on a vector as a function call. Here,f
can be a function or any other object for which either the callf(src::AbstractVector) -> dest::AbstractVector
(whenismutating = false
) orf(dest::AbstractVector,src::AbstractVector) -> dest
(whenismutating = true
) is supported. The value ofismutating
can be spefified, by default its value is guessed by looking at the number of arguments of the first method in the method list off
.A second function or object can optionally be provided that implements the action of the adjoint (transposed) linear map. Here, it is always assumed that this represents the adjoint/conjugate transpose, though this is of course equivalent to the normal transpose for real linear maps. Furthermore, the conjugate transpose also enables the use of
mul!(y, transpose(A), x)
using some extra conjugation calls on the input and output vector. If no second function is provided, thanmul!(y, transpose(A), x)
andmul!(y, adjoint(A), x)
cannot be used with this linear map, unless it is symmetric or hermitian.M
is the number of rows (length of the output vectors) andN
the number of columns (length of the input vectors). When the latter is not specified,N = M
.Finally, one can specify the
eltype
of the resulting linear map using the type parameterT
. If not specified, a default value ofFloat64
is assumed. Use a complex typeT
if the function represents a complex linear map.In summary, the keyword arguments and their default values are:
ismutating
:false
if the functionf
accepts a single vector argument corresponding to the input, andtrue
if it accepts two vector arguments where the first will be mutated so as to contain the result. In both cases, the resultingA::FunctionMap
will support both the mutating and non-mutating matrix vector multiplication. Default value is guessed based on the number of arguments for the first method in the method list off
; it is not possible to usef
andfc
where only one of the two is mutating and the other is not.issymmetric [=false]
: whether the function represents the multiplication with a symmetric matrix. Iftrue
, this will automatically enableA' * x
andtranspose(A) * x
.ishermitian [=T<:Real && issymmetric]
: whether the function represents the multiplication with a hermitian matrix. Iftrue
, this will automatically enableA' * x
andtranspose(A) * x
.isposdef [=false]
: whether the function represents the multiplication with a positive definite matrix.
-
Array(A::LinearMap)
,Matrix(A::LinearMap)
,convert(Matrix, A::LinearMap)
andconvert(Array, A::LinearMap)
Create a dense matrix representation of the
LinearMap
object, by multiplying it with the successive basis vectors. This is mostly for testing purposes or if you want to have the explicit matrix representation of a linear map for which you only have a function definition (e.g. to be able to use itstranspose
oradjoint
). This way, one may conveniently makeA
act on the columns of a matrixX
, instead of interpretingA * X
as a composed linear map:Matrix(A * X)
. For generic code, that is supposed to handle bothA::AbstractMatrix
andA::LinearMap
, it is recommended to useconvert(Matrix, A*X)
. -
SparseArrays.sparse(A::LinearMap)
andconvert(SparseMatrixCSC, A::LinearMap)
Create a sparse matrix representation of the
LinearMap
object, by multiplying it with the successive basis vectors. This is mostly for testing purposes or if you want to have the explicit sparse matrix representation of a linear map for which you only have a function definition (e.g. to be able to use itstranspose
oradjoint
). -
The following matrix multiplication methods.
A * x
: appliesA
tox
and returns the result;mul!(y::AbstractVector, A::LinearMap, x::AbstractVector)
: appliesA
tox
and stores the result iny
;mul!(Y::AbstractMatrix, A::LinearMap, X::AbstractMatrix)
: appliesA
to each column ofX
and stores the results in the corresponding columns ofY
;mul!(y::AbstractVector, A::LinearMap, x::AbstractVector, α::Number=1, β::Number=0)
: computesα * A * x + β * y
and stores the result iny
. Analogously forX,Y::AbstractMatrix
.
Applying the adjoint or transpose of
A
(if defined) tox
works exactly as in the usual matrix case:transpose(A) * x
andmul!(y, A', x)
, for instance.
None of the types below need to be constructed directly; they arise from performing
operations between LinearMap
objects or by calling the LinearMap
constructor
described above.
-
LinearMap
Abstract supertype
-
FunctionMap
Type for wrapping an arbitrary function that is supposed to implement the matrix vector product as a
LinearMap
. -
WrappedMap
Type for wrapping an
AbstractMatrix
orLinearMap
and to possible redefine the propertiesisreal
,issymmetric
,ishermitian
andisposdef
. AnAbstractMatrix
will automatically be converted to aWrappedMap
when it is combined with otherAbstractLinearMap
objects via linear combination or composition (multiplication). Note thatWrappedMap(mat1)*WrappedMap(mat2)
will never evaluatemat1*mat2
, since this is more costly then evaluatingmat1*(mat2*x)
and the latter is the only operation that needs to be performed byLinearMap
objects anyway. While the cost of matrix addition is comparable to matrix vector multiplication, this too is not performed explicitly since this would require new storage of the same amount as of the original matrices. -
UniformScalingMap
Type for representing a scalar multiple of the identity map (a.k.a. uniform scaling) of a certain size
M=N
, obtained simply asUniformScalingMap(λ, M)
. The typeT
of the resultingLinearMap
object is inferred from the type ofλ
. AUniformScalingMap
of the correct size will be automatically created ifLinearMap
objects are multiplied by scalars from the left or from the right, respecting the order of multiplication. -
LinearCombination
,CompositeMap
,TransposeMap
andAdjointMap
Used to add and multiply
LinearMap
objects, don't need to be constructed explicitly.
The LinearMap
object combines well with methods of the following packages:
- Arpack.jl: iterative eigensolver
eigs
and SVDsvds
; - IterativeSolvers.jl: iterative solvers, eigensolvers, and SVD;
- TSVD.jl: truncated SVD
tsvd
.
using LinearMaps
import Arpack, IterativeSolvers, TSVD
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)