Linear operators

A LinearOperator is a structure representing a linear operator from a VectorSpace to an other. More precisely, a LinearOperator is comprised of the three fields domain::VectorSpace, codomain::VectorSpace and coefficients::AbsractMatrix with matching dimensions and size.

julia> A = LinearOperator(Taylor(1), Taylor(1), [1 2 ; 3 4])LinearOperator : Taylor(1) → Taylor(1) with coefficients Matrix{Int64}:
 1  2
 3  4

The three fields domain, codomain and coefficients are accessible via the respective functions of the same name.

julia> domain(A)Taylor(1)
julia> codomain(A)Taylor(1)
julia> coefficients(A)2×2 Matrix{Int64}:
 1  2
 3  4

For convenience, the methods zeros, ones, fill and fill! are available:

julia> dom, codom = Taylor(1), Taylor(2)(Taylor(1), Taylor(2))
julia> zeros(dom, codom)LinearOperator : Taylor(1) → Taylor(2) with coefficients Matrix{Float64}:
 0.0  0.0
 0.0  0.0
 0.0  0.0
julia> ones(dom, codom)LinearOperator : Taylor(1) → Taylor(2) with coefficients Matrix{Float64}:
 1.0  1.0
 1.0  1.0
 1.0  1.0
julia> fill(2, dom, codom)LinearOperator : Taylor(1) → Taylor(2) with coefficients Matrix{Int64}:
 2  2
 2  2
 2  2
julia> fill!(zeros(dom, codom), 2)LinearOperator : Taylor(1) → Taylor(2) with coefficients Matrix{Float64}:
 2.0  2.0
 2.0  2.0
 2.0  2.0

The coefficients of a LinearOperator are indexed according to the indices of the domain and codomain (as given by indices).

julia> A[0:1,0:1] # indices(domain(A)), indices(codomain(A))2×2 Matrix{Int64}:
 1  2
 3  4

When the domain and/or the codomain of a LinearOperator is a CartesianSpace, its coefficients can be thought of as a block matrix . The function block extracts a LinearOperator composing the cartesian space.

julia> B = LinearOperator(ScalarSpace() × Taylor(1)^2, ScalarSpace() × Taylor(1)^2, reshape(1:25, 5, 5))LinearOperator : 𝕂 × Taylor(1)² → 𝕂 × Taylor(1)² with coefficients Base.ReshapedArray{Int64, 2, UnitRange{Int64}, Tuple{}}:
 1   6  11  16  21
 2   7  12  17  22
 3   8  13  18  23
 4   9  14  19  24
 5  10  15  20  25
julia> B[1:5,1:5] # indices(domain(B)), indices(codomain(B))5×5 Matrix{Int64}:
 1   6  11  16  21
 2   7  12  17  22
 3   8  13  18  23
 4   9  14  19  24
 5  10  15  20  25
julia> block(B, 1, 1) # extract the linear operator associated with the domain ScalarSpace() and codomain ScalarSpace()LinearOperator : 𝕂 → 𝕂 with coefficients SubArray{Int64, 2, Base.ReshapedArray{Int64, 2, UnitRange{Int64}, Tuple{}}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, false}:
 1
julia> block(B, 2, 2) # extract the linear operator associated with the domain Taylor(1)^2 and codomain Taylor(1)^2LinearOperator : Taylor(1)² → Taylor(1)² with coefficients SubArray{Int64, 2, Base.ReshapedArray{Int64, 2, UnitRange{Int64}, Tuple{}}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, false}:
  7  12  17  22
  8  13  18  23
  9  14  19  24
 10  15  20  25
julia> block(block(B, 2, 2), 1, 1)LinearOperator : Taylor(1) → Taylor(1) with coefficients SubArray{Int64, 2, Base.ReshapedArray{Int64, 2, UnitRange{Int64}, Tuple{}}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, false}:
 7  12
 8  13
julia> block(block(B, 2, 2), 2, 2)LinearOperator : Taylor(1) → Taylor(1) with coefficients SubArray{Int64, 2, Base.ReshapedArray{Int64, 2, UnitRange{Int64}, Tuple{}}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, false}:
 19  24
 20  25

Similarly, the function eachblock returns a Generator whose iterates yield each LinearOperator composing the cartesian space.

Arithmetic

The addition and subtraction operations are implemented as the + and - functions respectively.

julia> C = LinearOperator(Taylor(1), Taylor(1), [1 2 ; 3 4])LinearOperator : Taylor(1) → Taylor(1) with coefficients Matrix{Int64}:
 1  2
 3  4
julia> D = LinearOperator(Taylor(1), Taylor(2), [1 2 ; 3 4 ; 5 6])LinearOperator : Taylor(1) → Taylor(2) with coefficients Matrix{Int64}:
 1  2
 3  4
 5  6
julia> C + DLinearOperator : Taylor(1) → Taylor(2) with coefficients Matrix{Int64}:
 2  4
 6  8
 5  6
julia> C - DLinearOperator : Taylor(1) → Taylor(2) with coefficients Matrix{Int64}:
  0   0
  0   0
 -5  -6
julia> C + ILinearOperator(Taylor(1), Taylor(1), [1 2; 3 4]) + UniformScalingOperator{Bool}(true)
julia> C - ILinearOperator(Taylor(1), Taylor(1), [1 2; 3 4]) + UniformScalingOperator{Int64}(-1)

The product between LinearOperator is implemented as the * and ^ functions. The division between LinearOperator is implemented as the \ method.

julia> C * DLinearOperator : Taylor(1) → Taylor(1) with coefficients Matrix{Int64}:
  7  10
 15  22
julia> C ^ 3LinearOperator : Taylor(1) → Taylor(1) with coefficients Matrix{Int64}:
 37   54
 81  118
julia> C \ CLinearOperator : Taylor(1) → Taylor(1) with coefficients Matrix{Float64}:
 1.0  0.0
 0.0  1.0

The action of a LinearOperator is performed by the right product * of a LinearOperator with a Sequence; alternatively, LinearOperator defines a method on a Sequence representing *. Naturally, the resulting sequence is an element of the codomain of the LinearOperator.

Conversely, the operator \ between a LinearOperator and a Sequence corresponds to the action of the inverse of the LinearOperator; the output sequence is an element of the domain of the LinearOperator.

julia> x = Sequence(Taylor(2), [1, 1, 1])Sequence in Taylor(2) with coefficients Vector{Int64}:
 1
 1
 1
julia> C * x # C(x)Sequence in Taylor(1) with coefficients Vector{Int64}:
 3
 7
julia> D \ xSequence in Taylor(1) with coefficients Vector{Float64}:
 -1.0000000000000018
  1.0000000000000013

API

AbstractLinearOperator

LinearOperator{T<:VectorSpace,S<:VectorSpace,R<:AbstractMatrix} <: AbstractLinearOperator

julia> LinearOperator(Taylor(1), Taylor(1), [1 2 ; 3 4])
LinearOperator : Taylor(1) → Taylor(1) with coefficients Matrix{Int64}:
 1  2
 3  4

julia> LinearOperator(Taylor(2), ScalarSpace(), [1.0 0.5 0.25])
LinearOperator : Taylor(2) → 𝕂 with coefficients Matrix{Float64}:
 1.0  0.5  0.25

julia> LinearOperator([1 2 3 ; 4 5 6])
LinearOperator : 𝕂³ → 𝕂² with coefficients Matrix{Int64}:
 1  2  3
 4  5  6

block(A::LinearOperator{<:CartesianSpace,<:CartesianSpace}, i, j)

julia> A = LinearOperator(Taylor(1)^2, Taylor(1)^2, [i+j for i = 1:4, j = 1:4])
LinearOperator : Taylor(1)² → Taylor(1)² with coefficients Matrix{Int64}:
 2  3  4  5
 3  4  5  6
 4  5  6  7
 5  6  7  8

julia> block(A, 1, 1)
LinearOperator : Taylor(1) → Taylor(1) with coefficients SubArray{Int64, 2, Matrix{Int64}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, false}:
 2  3
 3  4

julia> block(A, 1, 2)
LinearOperator : Taylor(1) → Taylor(1) with coefficients SubArray{Int64, 2, Matrix{Int64}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, false}:
 4  5
 5  6

block(a::LinearOperator{<:CartesianSpace,{<:CartesianSpace})

julia> A = LinearOperator(Taylor(1)^2, Taylor(1)^2, [i+j for i = 1:4, j = 1:4])
LinearOperator : Taylor(1)² → Taylor(1)² with coefficients Matrix{Int64}:
 2  3  4  5
 3  4  5  6
 4  5  6  7
 5  6  7  8

julia> block(A)
2×2 Matrix{LinearOperator{Taylor, Taylor, SubArray{Int64, 2, Matrix{Int64}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, false}}}:
 [2 3; 3 4]  [4 5; 5 6]
 [4 5; 5 6]  [6 7; 7 8]

eachblock(A::LinearOperator{<:CartesianSpace,<:CartesianSpace})

julia> A = LinearOperator(Taylor(1)^2, Taylor(1)^2, [i+j for i = 1:4, j = 1:4])
LinearOperator : Taylor(1)² → Taylor(1)² with coefficients Matrix{Int64}:
 2  3  4  5
 3  4  5  6
 4  5  6  7
 5  6  7  8

julia> m = eachblock(A)
Base.Generator{Base.Iterators.ProductIterator{Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, RadiiPolynomial.var"#eachblock##2#eachblock##3"{LinearOperator{CartesianPower{Taylor}, CartesianPower{Taylor}, Matrix{Int64}}}}(RadiiPolynomial.var"#eachblock##2#eachblock##3"{LinearOperator{CartesianPower{Taylor}, CartesianPower{Taylor}, Matrix{Int64}}}(LinearOperator(Taylor(1)², Taylor(1)², [2 3 4 5; 3 4 5 6; 4 5 6 7; 5 6 7 8])), Base.Iterators.ProductIterator{Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}((Base.OneTo(2), Base.OneTo(2))))

julia> [v for v = m]
2×2 Matrix{LinearOperator{Taylor, Taylor, SubArray{Int64, 2, Matrix{Int64}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, false}}}:
 [2 3; 3 4]  [4 5; 5 6]
 [4 5; 5 6]  [6 7; 7 8]