LinearMaps
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.
What's new
Simplified interface. Only the names AbstractLinearMap
and LinearMap
are exported. AbstractLinearMap
is the root type that newly defined linear maps should be subtypes of. LinearMap
acts as a general purpose constructor for constructing linear maps from matrices, functions or altering the properties of existing AbstractLinearMap
objects, even though there is no actual type called LinearMap
.
Installation
Install with the package manager, i.e. Pkg.add("LinearMaps")
.
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:
- An
AbstractLinearMap
type that shares with theAbstractMatrix
type that it responds to the functionssize
,eltype
,isreal
,issym
,ishermitian
andisposdef
,transpose
andctranspose
and multiplication with a vector using both*
or the in-place versionA_mul_B!
. Depending on the subtype, alsoAt_mul_B
,At_mul_B!
,Ac_mul_B
andAc_mul_B!
are supported. Linear algebra functions that uses duck-typing for its arguments can handleAbstractLinearMap
objects similar toAbstractMatrix
objects, provided that they can be written using the above methods. UnlikeAbstractMatrix
types,AbstractLinearMap
objects cannot be indexed, neither usinggetindex
orsetindex!
. - A single method
LinearMap
that allows to constructAbstractLinearMap
objects from objects of typeFunction
,AbstractMatrix
orAbstractLinearMap
. This method allows to (re)define the properties (isreal
,issym
,ishermitian
,isposdef
) of the corresponding linear map. - A framework for combining objects of type
AbstractLinearMap
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 versionA_mul_B!
, 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 resultingAbstractLinearMap
object is not correctly inferred, theLinearMap
method can be called to redefine the properties.
Methods
-
LinearMap
General purpose method to construct AbstractLinearMap objects of specific types, as described in the Types section below
LinearMap(A::AbstractMatrix;isreal::Bool,issym::Bool,ishermitian::Bool,isposdef::Bool) LinearMap(A::AbstractLinearMap;isreal::Bool,issym::Bool,ishermitian::Bool,isposdef::Bool)
Create a
WrappedMap
object that will respond to the methodsisreal
,issym
,ishermitian
,isposdef
with the values provided by the keyword arguments. The default values correspond to the result of calling these methods on the argumentA
. This allows to use anAbstractMatrix
within theAbstractLinearMap
framework and to redefine the properties of an existingAbstractLinearMap
.LinearMap(f::Function,M::Int,N::Int=M;ismutating::Bool,issym::Bool,ishermitian::Bool,isposdef::Bool,ftranspose,fctranspose) LinearMap(f::Function,eltype::Type,M::Int,N::Int=M;ismutating::Bool,issym::Bool,ishermitian::Bool,isposdef::Bool,ftranspose,fctranspose)
Create
FunctionMap
object that wraps a function describing the action of the linear map on a vector. The corresponding properties of the linear map can also be specified. Here,f
represents the function implementing the action of the linear map on a vector, either as returning the result (i.e.f(src::AbstractVector) -> dest::AbstractVector
) whenismutating=false
(default) or as a mutating function that accepts a vector for the destination (i.e.f(dest::AbstractVector,src::AbstractVector) -> dest
).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
. Using the second calling convention, theeltype
of the resulting linear map can explicitly be specified. The keyword arguments and their default values are:-
ismutating [=false]
:false
if the functionf
(and if providedftranspose
and orfctranspose
) accepts a single vector argument corresponding to the input, andtrue
if they accept 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 as nonmutating matrix vector multiplication. -
isreal [=true]
(only in the first calling convention): iftrue
, it will createA::FunctionMap{Float64}
, otherwiseA::FunctionMap{Complex128}
. If the matrix representation of the function could be represented using a differenteltype
, then the second calling scheme is recommended. -
issym [=false]
: whether the function represents the multiplication with a symmetric matrix. Iftrue
, this will automatically enableA'*x
andA.'*x
. -
ishermitian [=false]
: whether the function represents the multiplication with a hermitian matrix. Iftrue
, this will automatically enableA'*x
andA.'*x
. -
isposdef [=false]
: whether the function represents the multiplication with a positive definite matrix. -
ftranspose [=nothing]
: an optional argument that can be used to pass a function that implements the multiplication with the transposed matrix -
fctranspose [=nothing]
: an optional argument that can be used to pass a function that implements the multiplication with the hermitian conjugated matrix
-
-
Base.full(linearmap)
Creates a full matrix representation of the linearmap object, by multiplying it with the successive basis vectors.
All matrix multiplication methods and the corresponding mutating versions.
Types
None of the types below need to be constructed directly; they arise from performing operations between AbstractLinearMap
objects or by calling the LinearMap
method described above.
-
AbstractLinearMap
Abstract supertype
-
FunctionMap
Type for wrapping an arbitrary function that is supposed to implement the matrix vector product as an
AbstractLinearMap
. -
WrappedMap
Type for wrapping an
AbstractMatrix
orAbstractLinearMap
and to possible redefine the propertiesisreal
,issym
,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 byAbstractLinearMap
objects anyway. While the cost of matrix addition is comparible 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. -
IdentityMap
Type for representing the identity map of a certain size
M=N
, obtained simply asIdentityMap{T}(M)
,IdentityMap(T,M)=IdentityMap(T,M,N)=IdentityMap(T,(M,N))
or evenIdentityMap(M)=IdentityMap(M,N)=IdentityMap((M,N))
. IfT
is not specified,Bool
is assumed, since operations betweenBool
and any otherNumber
will always be converted to the type of the otherNumber
. IfM!=N
, an error is returned. AnIdentityMap
of the correct size and element type will automatically be created ifLinearMap
objects are combined withI
, Julia's built in identity (UniformScaling
). -
LinearCombination
,CompositeMap
,TransposeMap
andCTransposeMap
Used to add and multiply
LinearMap
objects, don't need to be constructed explicitly.
Examples
The LinearMap
object combines well with the iterative eigensolver eigs
, which is the Julia wrapper for Arpack.
using LinearMaps
function leftdiff!(y::Vector,x::Vector) # left difference assuming periodic boundary conditions
length(y)==length(x) || throw(DimensionMismatch())
N=length(x)
@inbounds for i=1:N
y[i]=x[i]-x[mod1(i-1,N)]
end
return y
end
function mrightdiff!(y::Vector,x::Vector) # minus right difference
length(y)==length(x) || throw(DimensionMismatch())
N=length(x)
@inbounds for i=1:N
y[i]=x[i]-x[mod1(i+1,N)]
end
return y
end
D=LinearMap(leftdiff!,100;ismutating=true,fctranspose=mrightdiff!)
eigs(D'*D;nev=3,which=:SR)