DOCSEARCHbyAlgolia
Submit Docs

LinearOperators.jl

repository·main·Indexed 41 days ago

https://github.com/juliasmoothoptimizers/linearoperators.jl

A Julia package for efficient linear algebra that represents linear maps as operators. It defers expensive computations until operators are applied to vectors, enabling low-cost transformations like transposition and conjugation. It provides specialized types including BlockDiagonalOperator, LBFGSOperator, and opCholesky, and requires Julia 1.6 or higher.

Tokens
1.6K
Snippets
3
Records
8
Agent score
3%

What's inside LinearOperators.jl

  1. Understand the philosophy of LinearOperators.jl

    main

    The core philosophy of LinearOperators.jl is that operators behave like matrices but are defined by their effect when applied to a vector.

    Key characteristics:

    • Deferred Cost: The heavy computation is deferred until the operator is multiplied with a vector.
    • Cheap Transformations: Operators can be transposed, conjugated, or combined with other operators at a low computational cost.
    • Abstraction: They provide a way to represent linear maps without necessarily storing a dense matrix.

  2. How operators differ from matrices

    main

    Linear operators in this package behave similarly to matrices but have key differences in how they handle slicing and scalar multiplication. Crucially, an operator never reduces to a vector or a number through slicing or scalar operations; it always remains a LinearOperator.

    Slicing Behavior

    Unlike matrices, slicing an operator returns another operator. For example, opA[:, 1] is a LinearOperator, not a Vector.

    Scalar Multiplication

    Multiplying a slice by a scalar results in a LinearOperator. However, multiplying a slice by a single-element vector (e.g., opA[:, 1] * [3]) will result in an error, whereas in standard matrix math this would return a vector.

    Conversion

    To obtain a standard dense array from an operator, you must explicitly use the Matrix() constructor.

    using LinearOperators
A = rand(5,5)
opA = LinearOperator(A)

# In standard matrices:
A[:,1] * 3 isa Vector

# In LinearOperators:
opA[:,1] * 3 isa LinearOperator

# To get a matrix back:
Matrix(opA[:,1])

  3. How slicing and indexing differ from Matrices

    main

    A critical difference between LinearOperators.jl and standard Matrix types is that an operator never reduces to a vector or a scalar through slicing. Slices always return new operators.

    • Slicing: opA[:, 1] returns a LinearOperator (a column operator), not a Vector.
    • Scalar Indexing: opA[i, j] returns a (1,1) LinearOperator, whereas for a matrix it returns a number.
    • Conversion: To get a dense matrix from an operator slice, use the full() function.

    Example of behavior differences:

    A = rand(5,5)
opA = LinearOperator(A)

# Matrix behavior
A[:,1] * 3 # Returns a Vector

# Operator behavior
opA[:,1] * 3 # Returns a LinearOperator
opA[:,1] * [3] # Returns a Vector

# Scalar indexing
A[1,1] # Returns a Number
opA[1,1] # Returns a (1,1) LinearOperator

# Converting to dense
full(opA[:,1]) # Returns an nx1 matrix

  5. Perform common operations on operators

    main

    You can perform standard linear algebra operations on operators. Note that these operations are designed to be cheap, deferring the heavy computation until the operator is applied to a vector.

    • Transpose: A.'
    • Conjugate: conj(A)
    • Conjugate Transpose: A'
    • Slicing: A[:, 3], A[2:4, 1:5], or A[1, 1] (Note: slices return operators, not scalars or vectors).

  7. Available Linear Operators

    main

    The package provides several specialized operator types for different mathematical needs:

    OperatorDescription
    LinearOperatorBase class. Useful to define operators from functions
    TimedLinearOperatorLinear operator instrumented with timers from TimerOutputs.jl
    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
    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

  8. Utility functions for Linear Operators

    main

    Use these functions to inspect, validate, or manipulate operators:

    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
    showDisplay basic information about an operator
    sizeReturn the size of a linear operator
    symmetricDetermine whether the operator is symmetric
    normestEstimate 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$.