Types and methods

LinearMap(A::LinearMap; kwargs...)::WrappedMap
LinearMap(A::AbstractVecOrMatOrQ; kwargs...)::WrappedMap
LinearMap(J::UniformScaling, M::Int)::UniformScalingMap
LinearMap{T=Float64}(f, [fc,], M::Int, N::Int = M; kwargs...)::FunctionMap
LinearMap(A::MapOrVecOrMat, dims::Dims{2}, index::NTuple{2, AbstractVector{Int}})::EmbeddedMap
LinearMap(A::MapOrVecOrMat, dims::Dims{2}; offset::Dims{2})::EmbeddedMap

Construct a linear map object, either

1. from an existing LinearMap or AbstractVecOrMat/AbstractQ A, with the purpose of

redefining its properties via the keyword arguments kwargs, see below;

1. a UniformScaling object J with specified (square) dimension M;
2. from a function or callable object f;
3. from an existing LinearMap or AbstractVecOrMat/AbstractQ A, embedded in a larger zero map.

In the case of item 3, one also needs to specify the size of the equivalent matrix representation (M, N), i.e., for functions f acting on length N vectors and producing length M vectors (with default value N=M). Preferably, also the eltype T of the corresponding matrix representation needs to be specified, i.e., whether the action of f on a vector will be similar to, e.g., multiplying by numbers of type T. If not specified, the default value T=Float64 will be assumed. Optionally, a corresponding function fc can be specified that implements the adjoint (or transpose in the real case) of f.

The keyword arguments and their default values are:

• issymmetric::Bool = false : whether A or f act as a symmetric matrix
• ishermitian::Bool = issymmetric & T<:Real : whether A or f act as a Hermitian matrix
• isposdef::Bool = false : whether A or f act as a positive definite matrix.

For existing linear maps or matrices A, the default values will be taken by calling internal functions _issymmetric, _ishermitian and _isposdef on the existing object A. These in turn dispatch to (overloads of) LinearAlgebra's issymmetric, ishermitian, and isposdef methods whenever these checks are expected to be computationally cheap or even known at compile time as for certain structured matrices, but return false for generic AbstractMatrix types.

For the function-based constructor, there is one more keyword argument:

• ismutating::Bool : flags whether the function acts as a mutating matrix multiplication f(y,x) where the result vector y is the first argument (in case of true), or as a normal matrix multiplication that is called as y=f(x) (in case of false). The default value is guessed by looking at the number of arguments of the first occurrence of f in the method table.

For the EmbeddedMap constructors, dims specifies the total dimensions of the map. The index argument specifies two collections of indices inds1 and inds2, such that for the big zero map L (thought of as a matrix), one has L[inds1,inds2] == A. In other words, inds1 specifies the output indices, inds2 specifies the input indices. Alternatively, A may be shifted by offset, such that (thinking in terms of matrices again) L[offset[1] .+ axes(A, 1), offset[2] .+ axes(A, 2)] == A.

FunctionMap{T,iip}(f, [fc,], M, N = M; kwargs...)

Construct a FunctionMap object from a function or callable object f that represents a linear map of size (M, N), where N can be omitted for square operators of size (M,M). Furthermore, the eltype T of the corresponding matrix representation needs to be specified, i.e., whether the action of f on a vector will be similar to, e.g., multiplying by numbers of type T. Optionally, a second function fc can be specified that implements the adjoint (or transpose in the real case) of f.

Accepted keyword arguments and their default values are as in the LinearMap constructor.

Examples

```jldoctest julia> F = FunctionMap{Int64,false}(cumsum, 2) 2×2 FunctionMap{Int64,false}(cumsum; issymmetric=false, ishermitian=false, isposdef=false)

julia> F * ones(Int64, 2) 2-element Vector{Int64}: 1 2

julia> Matrix(F) 2×2 Matrix{Int64}: 1 0 1 1

Lazy representation of a scaled map λ * A = A * λ with real or complex map A <: LinearMap{RealOrComplex} and real or complex scaling factor λ <: RealOrComplex.

+(A::LinearMap, B::LinearMap)::LinearCombination

Construct a (lazy) representation of the sum/linear combination of the two operators. Sums of LinearMap/LinearCombination objects and LinearMap/LinearCombination objects are reduced to a single LinearCombination. In sums of LinearMaps and AbstractMatrix/UniformScaling objects, the latter get promoted to LinearMaps automatically.

Examples

julia> CS = LinearMap{Int}(cumsum, 3)::LinearMaps.FunctionMap;

julia> LinearMap(ones(Int, 3, 3)) + CS + I + rand(3, 3);

*(A::LinearMap, B::LinearMap)::CompositeMap

Construct a (lazy) representation of the product of the two operators. Products of LinearMap/CompositeMap objects and LinearMap/CompositeMap objects are reduced to a single CompositeMap. In products of LinearMaps and AbstractMatrix/UniformScaling objects, the latter get promoted to LinearMaps automatically.

Examples

julia> CS = LinearMap{Int}(cumsum, 3)::LinearMaps.FunctionMap;

julia> LinearMap(ones(Int, 3, 3)) * CS * I * rand(3, 3);

transpose(A::LinearMap)

Construct a lazy representation of the transpose of A. This can be either a TransposeMap wrapper of A, or a suitably redefined instance of the same type as A. For instance, for a linear combination of linear maps $A + B$, the transpose is given by $A^⊤ + B^⊤$, i.e., another linear combination of linear maps.

adjoint(A::LinearMap)

Construct a lazy representation of the adjoint of A. This can be either a AdjointMap wrapper of A, or a suitably redefined instance of the same type as A. For instance, for a linear combination of linear maps $A + B$, the adjoint is given by $A^* + B^*$, i.e., another linear combination of linear maps.

kron(A::LinearMap, B::LinearMap)::KroneckerMap
kron(A, B, Cs...)::KroneckerMap

Construct a (lazy) representation of the Kronecker product A⊗B. One of the two factors can be an AbstractMatrix, which is then promoted to a LinearMap automatically.

To avoid fallback to the generic Base.kron in the multi-map case, there must be a LinearMap object among the first 8 arguments in usage like kron(A, B, Cs...).

For convenience, one can also use A ⊗ B or ⊗(A, B, Cs...) (typed as \otimes+TAB) to construct the KroneckerMap, even when all arguments are of AbstractMatrix type.

If A, B, C and D are linear maps of such size that one can form the matrix products A*C and B*D, then the mixed-product property (A⊗B)*(C⊗D) = (A*C)⊗(B*D) holds. Upon vector multiplication, this rule is checked for applicability.

Examples

julia> J = LinearMap(I, 2) # 2×2 identity map
2×2 LinearMaps.UniformScalingMap{Bool} with scaling factor: true

julia> E = spdiagm(-1 => trues(1)); D = E + E' - 2I;

julia> Δ = kron(D, J) + kron(J, D); # discrete 2D-Laplace operator

julia> Matrix(Δ)
4×4 Matrix{Int64}:
 -4   1   1   0
  1  -4   0   1
  1   0  -4   1
  0   1   1  -4

⊗(k::Integer)

Construct a lazy representation of the k-th Kronecker power A^⊗(k) = A ⊗ A ⊗ ... ⊗ A, where A can be an AbstractMatrix or a LinearMap.

kronsum(A, B)::KroneckerSumMap
kronsum(A, B, Cs...)::KroneckerSumMap

Construct a (lazy) representation of the Kronecker sum A⊕B = A ⊗ Ib + Ia ⊗ B of two square linear maps of type LinearMap or AbstractMatrix. Here, Ia and Ib are identity operators of the size of A and B, respectively. Arguments of type AbstractMatrix are automatically promoted to LinearMap.

For convenience, one can also use A ⊕ B or ⊕(A, B, Cs...) (typed as \oplus+TAB) to construct the KroneckerSumMap.

Examples

julia> J = LinearMap(I, 2) # 2×2 identity map
2×2 LinearMaps.UniformScalingMap{Bool} with scaling factor: true

julia> E = spdiagm(-1 => trues(1)); D = LinearMap(E + E' - 2I);

julia> Δ₁ = kron(D, J) + kron(J, D); # discrete 2D-Laplace operator, Kronecker sum

julia> Δ₂ = kronsum(D, D);

julia> Δ₃ = D^⊕(2);

julia> Matrix(Δ₁) == Matrix(Δ₂) == Matrix(Δ₃)
true

⊕(k::Integer)

Construct a lazy representation of the k-th Kronecker sum power A^⊕(k) = A ⊕ A ⊕ ... ⊕ A, where A can be a square AbstractMatrix or a LinearMap. This calls sumkronsum on the k-tuple (A, ..., A) for k ≥ 3.

Example

```jldoctest julia> Matrix([1 0; 0 1]^⊕(2)) 4×4 Matrix{Int64}: 2 0 0 0 0 2 0 0 0 0 2 0 0 0 0 2

squarekron(A₁::MapOrMatrix, A₂::MapOrMatrix, A₃::MapOrMatrix, Aᵢ::MapOrMatrix...)::CompositeMap

Construct a (lazy) representation of the Kronecker product ⨂ᵢ₌₁ⁿ Aᵢ of at least 3 square Kronecker factors. In contrast to kron, this function assumes that all Kronecker factors are square, and makes use of the following identity<sup>[1]</sup>:

\[\bigotimes_{i=1}^n A_i = \prod_{i=1}^n I_1 \otimes \ldots \otimes I_{i-1} \otimes A_i \otimes I_{i+1} \otimes \ldots \otimes I_n\]

where $I_k$ is an identity matrix of the size of $A_k$. By associativity, the Kronecker product of the identity operators may be combined to larger identity operators $I_{1:i-1}$ and $I_{i+1:n}$, which yields

\[\bigotimes_{i=1}^n A_i = \prod_{i=1}^n I_{1:i-1} \otimes A_i \otimes I_{i+1:n}\]

i.e., a CompositeMap where each factor is a Kronecker product consisting of three maps: outer UniformScalingMaps and the respective Kronecker factor. This representation is expected to yield significantly faster multiplication (and reduce memory allocation) compared to kron, but benchmarking intended use cases is highly recommended.

sumkronsum(A, B)::LinearCombination
sumkronsum(A, B, Cs...)::LinearCombination

Construct a (lazy) representation of the Kronecker sum A⊕B of two or more square objects of type LinearMap or AbstractMatrix. This function makes use of the following representation of Kronecker sums<sup>[1]</sup>:

\[\bigoplus_{i=1}^n A_i = \sum_{i=1}^n I_1 \otimes \ldots \otimes I_{i-1} \otimes A_i \otimes I_{i+1} \otimes \ldots \otimes I_n\]

where $I_k$ is the identity operator of the size of $A_k$. By associativity, the Kronecker product of the identity operators may be combined to larger identity operators $I_{1:i-1}$ and $I_{i+1:n}$, which yields

\[\bigoplus_{i=1}^n A_i = \sum_{i=1}^n I_{1:i-1} \otimes A_i \otimes I_{i+1:n},\]

i.e., a LinearCombination where each summand is a Kronecker product consisting of three maps: outer UniformScalingMaps and the respective Kronecker factor. This representation is expected to yield significantly faster multiplication (and reduce memory allocation) compared to kronsum, especially for 3 or more Kronecker summands, but benchmarking intended use cases is highly recommended.

Examples

julia> J = LinearMap(I, 2) # 2×2 identity map
2×2 LinearMaps.UniformScalingMap{Bool} with scaling factor: true

julia> E = spdiagm(-1 => trues(1)); D = LinearMap(E + E' - 2I);

julia> Δ₁ = kron(D, J) + kron(J, D); # discrete 2D-Laplace operator, Kronecker sum

julia> Δ₂ = sumkronsum(D, D);

julia> Δ₃ = D^⊕(2);

julia> Matrix(Δ₁) == Matrix(Δ₂) == Matrix(Δ₃)
true

hcat(As::Union{LinearMap,UniformScaling,AbstractVecOrMatOrQ}...)::BlockMap

Construct a (lazy) representation of the horizontal concatenation of the arguments. All arguments are promoted to LinearMaps automatically.

Examples

julia> CS = LinearMap{Int}(cumsum, 3)::LinearMaps.FunctionMap;

julia> L = [CS LinearMap(ones(Int, 3, 3))]::LinearMaps.BlockMap;

julia> L * ones(Int, 6)
3-element Vector{Int64}:
 4
 5
 6

vcat(As::Union{LinearMap,UniformScaling,AbstractVecOrMatOrQ}...)::BlockMap

Construct a (lazy) representation of the vertical concatenation of the arguments. All arguments are promoted to LinearMaps automatically.

Examples

julia> CS = LinearMap{Int}(cumsum, 3)::LinearMaps.FunctionMap;

julia> L = [CS; LinearMap(ones(Int, 3, 3))]::LinearMaps.BlockMap;

julia> L * ones(Int, 3)
6-element Vector{Int64}:
 1
 2
 3
 3
 3
 3

hvcat(rows::Tuple{Vararg{Int}}, As::Union{LinearMap,UniformScaling,AbstractVecOrMatOrQ}...)::BlockMap

Construct a (lazy) representation of the horizontal-vertical concatenation of the arguments. The first argument specifies the number of arguments to concatenate in each block row. All arguments are promoted to LinearMaps automatically.

Examples

julia> CS = LinearMap{Int}(cumsum, 3)::LinearMaps.FunctionMap;

julia> L = [CS CS; CS CS]::LinearMaps.BlockMap;

julia> L.rows
(2, 2)

julia> L * ones(Int, 6)
6-element Vector{Int64}:
 2
 4
 6
 2
 4
 6

cat(As::Union{LinearMap,AbstractVecOrMatOrQ}...; dims=(1,2))::BlockDiagonalMap

Construct a (lazy) representation of the diagonal concatenation of the arguments. To avoid fallback to the generic Base.cat, there must be a LinearMap object among the first 8 arguments.

blockdiag(As::Union{LinearMap,AbstractVecOrMatOrQ}...)::BlockDiagonalMap

Construct a (lazy) representation of the diagonal concatenation of the arguments. To avoid fallback to the generic blockdiag, there must be a LinearMap object among the first 8 arguments.

FillMap(λ, (m, n))::FillMap
FillMap(λ, m, n)::FillMap

Construct a (lazy) representation of an operator whose matrix representation would be an m×n-matrix filled constantly with the value λ.

InverseMap(A; solver = ldiv!)

Lazy inverse of A such that InverseMap(A) * x is the same as A \ x. Letting an InverseMap act on a vector thus requires solving a linear system.

A solver function can be passed with the solver keyword argument. The solver should be of the form f(y, A, x) where A is the wrapped map, x the right hand side, and y a preallocated output vector in which the result should be stored. The default solver is LinearAlgebra.ldiv!.

Note that A must be compatible with the solver function. A can, for example, be a factorization of a matrix, or another LinearMap (in combination with an iterative solver such as conjugate gradient).

Examples

julia> using LinearMaps, LinearAlgebra

julia> A = rand(2, 2); b = rand(2);

julia> InverseMap(lu(A)) * b
2-element Vector{Float64}:
  1.0531895201271027
 -0.4718540250893251

julia> A \ b
2-element Vector{Float64}:
  1.0531895201271027
 -0.4718540250893251

khatrirao(A::MapOrVecOrMat, B::MapOrVecOrMat) -> KhatriRaoMap

Construct a lazy representation of the Khatri-Rao (or column-wise Kronecker) product of two maps or arrays A and B. For the application to vectors, the transpose action of A on vectors needs to be defined.

facesplitting(A::AbstractMatrix, B::AbstractMatrix) -> FaceSplittingMap

Construct a lazy representation of the face-splitting (or row-wise Kronecker) product of two matrices A and B.

*(A::LinearMap, x::AbstractVector)::AbstractVector

Compute the action of the linear map A on the vector x.

Examples

julia> A=LinearMap([1.0 2.0; 3.0 4.0]); x=[1.0, 1.0];

julia> A*x
2-element Vector{Float64}:
 3.0
 7.0

julia> A(x)
2-element Vector{Float64}:
 3.0
 7.0

*(A::LinearMap, X::Union{AbstractMatrix,AbstractQ})::CompositeMap

Return the CompositeMap A*LinearMap(X), interpreting the matrix X as a linear operator, rather than a collection of column vectors. To compute the action of A on each column of X, call Matrix(A*X) or use the in-place multiplication mul!(Y, A, X[, α, β]) with an appropriately sized, preallocated matrix Y.

Examples

julia> A=LinearMap([1.0 2.0; 3.0 4.0]); X=[1.0 1.0; 1.0 1.0];

julia> A*X isa LinearMaps.CompositeMap
true

*(X::Union{AbstractMatrix,AbstractQ}, A::LinearMap)::CompositeMap

Return the CompositeMap LinearMap(X)*A, interpreting the matrix X as a linear operator. To compute the right-action of A on each row of X, call Matrix(X*A) or mul!(Y, X, A) for the in-place version.

Examples

julia> X=[1.0 1.0; 1.0 1.0]; A=LinearMap([1.0 2.0; 3.0 4.0]);

julia> X*A isa LinearMaps.CompositeMap
true

mul!(Y::AbstractVecOrMat, A::LinearMap, B::AbstractVector) -> Y
mul!(Y::AbstractMatrix, A::LinearMap, B::AbstractMatrix) -> Y

Calculates the action of the linear map A on the vector or matrix B and stores the result in Y, overwriting the existing value of Y. Note that Y must not be aliased with either A or B.

Examples

julia> A = LinearMap([1.0 2.0; 3.0 4.0]); B = ones(2); Y = similar(B); mul!(Y, A, B)
2-element Vector{Float64}:
 3.0
 7.0

julia> A = LinearMap([1.0 2.0; 3.0 4.0]); B = ones(2,2); Y = similar(B); mul!(Y, A, B)
2×2 Matrix{Float64}:
 3.0  3.0
 7.0  7.0

mul!(C::AbstractVecOrMat, A::LinearMap, B::AbstractVector, α, β) -> C
mul!(C::AbstractMatrix, A::LinearMap, B::AbstractMatrix, α, β) -> C

Combined inplace multiply-add $A B α + C β$. The result is stored in C by overwriting it. Note that C must not be aliased with either A or B.

Examples

julia> A=LinearMap([1.0 2.0; 3.0 4.0]); B=[1.0, 1.0]; C=[1.0, 3.0];

julia> mul!(C, A, B, 100.0, 10.0) === C
true

julia> C
2-element Vector{Float64}:
 310.0
 730.0

julia> A=LinearMap([1.0 2.0; 3.0 4.0]); B=[1.0 1.0; 1.0 1.0]; C=[1.0 2.0; 3.0 4.0];

julia> mul!(C, A, B, 100.0, 10.0) === C
true

julia> C
2×2 Matrix{Float64}:
 310.0  320.0
 730.0  740.0

mul!(C::AbstractMatrix, A::AbstractMatrix, B::LinearMap) -> C

Calculates the matrix representation of A*B and stores the result in C, overwriting the existing value of C. Note that C must not be aliased with either A or B. The computation C = A*B is performed via C' = B'A'.

Examples

julia> A=[1.0 1.0; 1.0 1.0]; B=LinearMap([1.0 2.0; 3.0 4.0]); C = similar(A); mul!(C, A, B);

julia> C
2×2 Matrix{Float64}:
 4.0  6.0
 4.0  6.0

mul!(C::AbstractMatrix, A::AbstractMatrix, B::LinearMap, α, β) -> C

Combined inplace multiply-add $A B α + C β$. The result is stored in C by overwriting it. Note that C must not be aliased with either A or B.

Examples

julia> A=[1.0 1.0; 1.0 1.0]; B=LinearMap([1.0 2.0; 3.0 4.0]); C = copy(A);

julia> mul!(C, A, B, 1, 1)
2×2 Matrix{Float64}:
 5.0  7.0
 5.0  7.0

mul!(Y::AbstractMatrix, A::LinearMap, b::Number) -> Y

Scales the matrix representation of the linear map A by b and stores the result in Y, overwriting the existing value of Y.

Examples

julia> A = LinearMap{Int}(cumsum, 3); b = 2; Y = Matrix{Int}(undef, (3,3));

julia> mul!(Y, A, b)
3×3 Matrix{Int64}:
 2  0  0
 2  2  0
 2  2  2

mul!(Y::AbstractMatrix, A::LinearMap, b::Number, α::Number, β::Number) -> Y

Scales the matrix representation of the linear map A by b*α, adds the result to Y*β and stores the final result in Y, overwriting the existing value of Y.

Examples

julia> A = LinearMap{Int}(cumsum, 3); b = 2; Y = ones(Int, (3,3));

julia> mul!(Y, A, b, 2, 1)
3×3 Matrix{Int64}:
 5  1  1
 5  5  1
 5  5  5

*(x::LinearAlgebra.AdjointAbsVec, A::LinearMap)::AdjointAbsVec

Compute the right-action of the linear map A on the adjoint vector x and return an adjoint vector.

Examples

julia> A=LinearMap([1.0 2.0; 3.0 4.0]); x=[1.0, 1.0]; x'A
1×2 adjoint(::Vector{Float64}) with eltype Float64:
 4.0  6.0

*(x::LinearAlgebra.TransposeAbsVec, A::LinearMap)::TransposeAbsVec

Compute the right-action of the linear map A on the transpose vector x and return a transpose vector.

Examples

julia> A=LinearMap([1.0 2.0; 3.0 4.0]); x=[1.0, 1.0]; transpose(x)*A
1×2 transpose(::Vector{Float64}) with eltype Float64:
 4.0  6.0