A Julia Linear Operator Package

Operators behave like matrices (with exceptions) but are defined by their effect when applied to a vector. They can be transposed, conjugated, or combined with other operators cheaply. The costly operation is deferred until multiplied with a vector.

Compatibility

Julia 1.3 and up.

How to Install

pkg> add LinearOperators
pkg> test LinearOperators

Operators Available

OperatorDescription
LinearOperatorBase class. Useful to define operators from functions
TimedLinearOperatorLinear operator instrumented with timers from TimerOutputs
BlockDiagonalOperatorBlock-diagonal linear operator
opEyeIdentity operator
opOnesAll ones operator
opZerosAll zeros operator
opDiagonalSquare (equivalent to diagm()) or rectangular diagonal operator
opInverseEquivalent to \
opCholeskyMore efficient than opInverse for symmetric positive definite matrices
opLDLSimilar to opCholesky, for general sparse symmetric matrices
opHouseholderApply a Householder transformation I-2hh'
opHermitianRepresent a symmetric/hermitian operator based on the diagonal and strict lower triangle
opRestrictionRepresent a selection of "rows" when composed on the left with an existing operator
opExtensionRepresent a selection of "columns" when composed on the right with an existing operator
LBFGSOperatorLimited-memory BFGS approximation in operator form (damped or not)
InverseLBFGSOperatorInverse of a limited-memory BFGS approximation in operator form (damped or not)
LSR1OperatorLimited-memory SR1 approximation in operator form
kronKronecker tensor product in linear operator form

Utility Functions

FunctionDescription
check_ctransposeCheap check that A' is correctly implemented
check_hermitianCheap check that A = A'
check_positive_definiteCheap check that an operator is positive (semi-)definite
diagExtract the diagonal of an operator
MatrixConvert an abstract operator to a dense array
hermitianDetermine whether the operator is Hermitian
push!For L-BFGS or L-SR1 operators, push a new pair {s,y}
reset!For L-BFGS or L-SR1 operators, reset the data
shapeReturn the size of a linear operator
showDisplay basic information about an operator
sizeReturn the size of a linear operator
symmetricDetermine whether the operator is symmetric
normestEstimate the 2-norm

Other Operations on Operators

Operators can be transposed (A.'), conjugated (conj(A)) and conjugate-transposed (A'). Operators can be sliced (A[:,3], A[2:4,1:5], A[1,1]), but unlike matrices, slices always return operators (see differences).

Differences

Unlike matrices, an operator never reduces to a vector or a number.

A = rand(5,5)
opA = LinearOperator(A)
A[:,1] * 3 # Vector
5-element Vector{Float64}:
 1.4471654905502307
 2.0720315298108285
 0.09493729757887093
 0.6397410183094988
 2.502579300790608
opA[:,1] * 3 # LinearOperator
Linear operator
  nrow: 5
  ncol: 1
  eltype: Float64
  symmetric: false
  hermitian: false
  nprod:   0
  ntprod:  0
  nctprod: 0

# A[:,1] * [3] # ERROR
opA[:,1] * [3] # Vector
5-element Vector{Float64}:
 1.4471654905502307
 2.0720315298108285
 0.09493729757887093
 0.6397410183094988
 2.502579300790608

This is also true for A[i,:], which returns vectors on Julia 0.6, and for the scalar A[i,j]. Similarly, opA[1,1] is an operator of size (1,1):"

(opA[1,1] * [3])[1] - A[1,1] * 3
0.0

In the same spirit, the operator Matrix always returns a matrix.

Matrix(opA[:,1])
5×1 Matrix{Float64}:
 0.4823884968500769
 0.6906771766036095
 0.03164576585962364
 0.21324700610316627
 0.834193100263536

Other Operators

  • LLDL features a limited-memory LDLᵀ factorization operator that may be used as preconditioner in iterative methods
  • MUMPS.jl features a full distributed-memory factorization operator that may be used to represent the preconditioner in, e.g., constraint-preconditioned Krylov methods.

Testing

julia> Pkg.test("LinearOperators")