Skip to content

JuliaSmoothOptimizers/LinearOperators.jl

Repository files navigation

A Julia Linear Operator Package

Documentation Linux/macOS/Windows/FreeBSD Coverage DOI
docs-stable docs-dev build-gh build-cirrus codecov doi

How to Cite

If you use LinearOperators.jl in your work, please cite using the format given in CITATION.cff.

Philosophy

Operators behave like matrices (with some exceptions - see below) 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.6 and up.

How to Install

pkg> add LinearOperators
pkg> test LinearOperators

How to use

Check the tutorial.

Operators Available

Operator Description
LinearOperator Base class. Useful to define operators from functions
TimedLinearOperator Linear operator instrumented with timers from TimerOutputs
BlockDiagonalOperator Block-diagonal linear operator
opEye Identity operator
opOnes All ones operator
opZeros All zeros operator
opDiagonal Square (equivalent to diagm()) or rectangular diagonal operator
opInverse Equivalent to \
opCholesky More efficient than opInverse for symmetric positive definite matrices
opHouseholder Apply a Householder transformation I-2hh'
opHermitian Represent a symmetric/hermitian operator based on the diagonal and strict lower triangle
opRestriction Represent a selection of "rows" when composed on the left with an existing operator
opExtension Represent a selection of "columns" when composed on the right with an existing operator
LBFGSOperator Limited-memory BFGS approximation in operator form (damped or not)
InverseLBFGSOperator Inverse of a limited-memory BFGS approximation in operator form (damped or not)
LSR1Operator Limited-memory SR1 approximation in operator form

Utility Functions

Function Description
check_ctranspose Cheap check that A' is correctly implemented
check_hermitian Cheap check that A = A'
check_positive_definite Cheap check that an operator is positive (semi-)definite
diag Extract the diagonal of an operator
Matrix Convert an abstract operator to a dense array
hermitian Determine 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
show Display basic information about an operator
size Return the size of a linear operator
symmetric Determine whether the operator is symmetric
normest Estimate the 2-norm
solve_shifted_system! Solves linear system $(B + \sigma I) x = b$, where $B$ is a forward L-BFGS operator and $\sigma \geq 0$.

Other Operations on Operators

Operators can be transposed (transpose(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 below).

Differences

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

A = rand(5,5)
opA = LinearOperator(A)
A[:,1] * 3 # Vector
opA[:,1] * 3 # LinearOperator
A[:,1] * [3] # ERROR
opA[:,1] * [3] # Vector

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

opA[1,1] # LinearOperator
A[1,1] # Number

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

full(opA[:,1]) # nx1 matrix

Other Operators

  • LimitedLDLFactorizations features a limited-memory LDLT 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.

Bug reports and discussions

If you think you found a bug, feel free to open an issue. Focused suggestions and requests can also be opened as issues. Before opening a pull request, start an issue or a discussion on the topic, please.

If you want to ask a question not suited for a bug report, feel free to start a discussion here. This forum is for general discussion about this repository and the JuliaSmoothOptimizers organization, so questions about any of our packages are welcome.