GitHub

API

Modules

Types and constants

Functions and macros

Documentation

Base.:\ — Method
op::MatrixShapedOperator \ x::AbstractVector
op::MatrixShapedOperator \ X::AbstractMatrix

Eager inverse application (solve) of a square operator to a vector, or to the columns of a matrix as multiple right-hand sides.

Extended help

Like lower_cholesky and logabsdet this is structural: it follows the inverse path of the operator's representation, or throws an ArgumentError if there is none. LinearAlgebra.factorize caches a factorization for repeated solves.

source
MatrixShapedOperators.MatrixShapedOperators — Module
MatrixShapedOperators

Linear operators with matrix shape and semantics, accessed through application and operator algebra, not through element access: operators combine lazily with operators (*, +) and eagerly with arrays (op * x, column-wise op * X, solves via op \ x).

MatrixShapedOperators is designed to be tracing- and GPU-friendly. Structural and behavioral information (triangularity, symmetricity, positive-definiteness, etc.) is captured at construction time and managed through type-stable traits.

Use asoperator to wrap matrices and other compatible objects as operators.

source
MatrixShapedOperators.BlockDiagOperator — Type
struct BlockDiagOperator{T<:Number,BS<:Union{Tuple,AbstractVector}} <: MatrixShapedOperator{T}

Represents the block-diagonal composition of matrix-shaped operators. Blocks need not be square, they may be a tuple (few, heterogeneous) or a vector (many, uniform).

User code should not call BlockDiagOperator directly, but use blockdiag_operator instead.

Traits derive from the blocks via the trait algebra; lower_cholesky and logabsdet operate block-wise.

source
MatrixShapedOperators.ColGramOperator — Type
struct ColGramOperator{T<:Number,OP} <: MatrixShapedOperator{T}

Represents the column-Gram operator A' * A of a matrix or matrix-shaped operator A, storing only asoperator(A) internally.

User code should not call ColGramOperator directly, but use colgram_operator instead.

Keeps structured factors (UpperTriangular, Diagonal, ...) as such, their structure drives trait derivation and lower_cholesky.

source
MatrixShapedOperators.IdentityOperator — Type
struct IdentityOperator <: MatrixShapedOperator{Bool}

Represents the identity operator of an explicit size.

User code should not call IdentityOperator directly, but use 𝟙 instead:

id = 𝟙(size(x, 1))
id * x == x

The element type is the multiplicative neutral Bool, which promotes away; scaled identities are UniformScalingOperators.

source
MatrixShapedOperators.MatrixAsOperator — Type
struct MatrixAsOperator{T,TRS<:Tuple,M} <: MatrixShapedOperator{T}

Wraps an AbstractMatrix or LinearAlgebra.AbstractQ as a MatrixShapedOperator op.

The wrapped object is available via Base.parent(op).

User code should not call MatrixAsOperator directly, but use asoperator instead.

Extended help

Application delegates to multiplication with the wrapped matrix, so specialized matrix types (e.g. Diagonal, GPU arrays) keep their optimized code paths.

The structure traits are captured at construction as a traitset in the type parameter TRS - from the structure of the wrapped matrix by default, plus declarations via asoperator.

source
MatrixShapedOperators.MatrixShapedOperator — Type
abstract type MatrixShapedOperator{T<:Number}

Abstract supertype for linear operators with matrix shape and semantics, accessed through application and operator algebra, not element access.

Operators combine lazily with operators and eagerly with arrays, on which they act just like a matrix:

op_a * op_b    # lazy MatrixShapedProduct
op_a + op_b    # lazy MatrixShapedSum
op * x         # eager application to a vector
op * X         # eager application to the columns of a matrix
op \ x         # eager structural solve

Application is out-of-place, though identity-like operations may return the argument itself. Matrices enter the lazy algebra via asoperator; adding a matrix to an operator is an error, so implicit operators are never materialized by accident.

Implementation

Subtypes must implement

  • MatrixShapedOperators.explicit_mul_impl(op, x::AbstractVector{<:Number})
  • Base.adjoint(op) and Base.size(op)

and declare their structural properties via the trait constructors (see MatrixShapedOperators.StructureTrait). Symmetry is typically only claimed for real element types, where transpose equals adjoint. The is* queries default to the traits and refine them by cheap value checks where those exist. batched_mul_impl may specialize batched application, the fallback applies columns individually (see BatchedMulStyle).

*, + and \ are defined centrally and dispatch to lazy_mul_impl, lazy_add_impl, explicit_mul_impl resp. ldiv_impl after argument checking. transpose, scalar scaling, -, /, mul! and AbstractMatrix/Matrix materialization are generic.

source
MatrixShapedOperators.MatrixShapedProduct — Type
struct MatrixShapedProduct{T<:Number,OPS<:Union{Tuple,AbstractVector}} <: MatrixShapedOperator{T}

Represents the product factors[1] * factors[2] * ... of matrix-shaped operators with matching adjacent sizes, applied right-to-left.

User code should not call MatrixShapedProduct directly, but use * instead:

op_a * op_b isa MatrixShapedProduct

Extended help

Factors may be a tuple (few, heterogeneous) or a vector (many, uniform), MatrixShaped values are normalized via asoperator.

Symmetry/definiteness traits are not preserved by products; use RowGramOperator for A * A'. Triangularity, unitarity and full row rank compose (see the trait algebra).

source
MatrixShapedOperators.MatrixShapedSum — Type
struct MatrixShapedSum{T<:Number,OPS<:Union{Tuple,AbstractVector}} <: MatrixShapedOperator{T}

Represents the sum terms[1] + terms[2] + ... of equal-size matrix-shaped operators.

User code should not call MatrixShapedSum directly, but use + instead:

op_a + op_b + 2 * I isa MatrixShapedSum

Extended help

Adding a UniformScaling adds a UniformScalingOperator term. Terms may be a tuple (few, heterogeneous) or a vector (many, uniform), MatrixShaped values are normalized via asoperator.

The sum is symmetric/hermitian/positive semi-definite if all terms are, and positive definite if additionally at least one term is.

source
MatrixShapedOperators.MulFuncOperator — Type
struct MulFuncOperator{T<:Number,TRS<:Tuple,F,G} <: MatrixShapedOperator{T}

Linear operator defined by a multiplication function ovp and an adjoint multiplication function vop.

User code should not call MulFuncOperator directly, but use mulfunc_operator instead.

Extended help

The structure traits declared at construction are stored as a traitset in the type parameter TRS.

op * X hands the whole matrix to ovp if BatchedMulStyle(ovp) is BatchedMul(), and calls ovp on each column separately otherwise (the default); likewise for op' * X and vop. Scalar scaling drops the definiteness declaration. Compatible with program tracing if ovp and vop are.

source
MatrixShapedOperators.RowGramOperator — Type
struct RowGramOperator{T<:Number,OP} <: MatrixShapedOperator{T}

Represents the row-Gram operator A * A' of a matrix or matrix-shaped operator A, storing only asoperator(A) internally.

User code should not call RowGramOperator directly, but use rowgram_operator instead.

Keeps structured factors (LowerTriangular, Diagonal, ...) as such, their structure drives trait derivation and lower_cholesky.

source
MatrixShapedOperators.ScaledOperator — Type
struct ScaledOperator{T<:Number,S<:Number,O} <: MatrixShapedOperator{T}

Represents s * op for a scalar s and a matrix-shaped operator op.

User code should not call ScaledOperator directly, but use scalar multiplication:

sop = s * op
sop * x == s .* (op * x)

Scaling drops definiteness, unitarity and rank structure traits, and hermitianity for complex scalars; is* queries refine by the value of s. Nested scalings collapse into a single ScaledOperator; logabsdet, inverse application and (for real non-negative scalings) Cholesky and Gram factors derive from the wrapped operator.

source
MatrixShapedOperators.UniformScalingOperator — Type
const UniformScalingOperator{T,S} = ScaledOperator{T,S,IdentityOperator}

A scaled identity λ * 𝟙(n): a type alias for dispatch on uniform scalings.

User code should not construct UniformScalingOperator directly, but use λ * 𝟙(n), or add a sizeless λ * 𝟙 resp. λ * I to an operator.

In operator products a scaled identity melts into a scalar scaling of the other factor.

source
MatrixShapedOperators.Unitarity — Type
abstract type MatrixShapedOperators.Unitarity

Structural unitarity trait: Unitarity(A) returns IsUnitary() if A is unitary (orthogonal, for real element types) for any element values of its representation, UnknownUnitarity() otherwise.

Unitary operators are square with full rank, their inverse is their adjoint and their logabsdet is zero.

source
MatrixShapedOperators.WoodburyFactorization — Type
struct WoodburyFactorization{T<:Number} <: LinearAlgebra.Factorization{T}

Cached Gram factorization op == F * F' of a positive-definite WoodburyOperator.

User code should not call WoodburyFactorization directly, but use LinearAlgebra.factorize instead. The result supports \, logabsdet and rowgram_factor without refactorizing:

op = woodbury_operator(A, B, D)
FW = factorize(op)
FW \ b
logabsdet(FW)
F = rowgram_factor(FW)  # e.g. for sampling and whitening
source
MatrixShapedOperators.WoodburyOperator — Type
struct WoodburyOperator{T<:Number,TA,TB,TD} <: MatrixShapedOperator{T}

Represents A + B * D * B'.

User code should not call WoodburyOperator directly, but use woodbury_operator instead.

Application costs one application per component, O(n k) for a structured base.

source
MatrixShapedOperators.ZeroOperator — Type
struct ZeroOperator <: MatrixShapedOperator{Bool}

Represents the zero operator of an explicit, possibly rectangular size.

User code should not call ZeroOperator directly, but use 𝟘 instead; Base.zero of any matrix-shaped operator also returns one:

z = 𝟘(3, 2)
z * x == zeros(3)  # x of length 2
zero(op) isa ZeroOperator

Adding a zero operator is a type-level no-op and multiplication with one collapses to a zero operator of the resulting shape. The element type is the additive neutral Bool, which promotes away in operator algebra.

source
MatrixShapedOperators.MatrixShaped — Type
MatrixShaped{T<:Number} = Union{AbstractMatrix{T},MatrixShapedOperator{T}}

Objects with matrix shape and algebra: matrices and matrix-shaped operators. asoperator and asmatrix normalize between the two.

Extended help

Covers the foreign AbstractMatrix, so methods of Base/LinearAlgebra functions must not dispatch on it (type piracy) - it is meant for the arguments of package-owned functions.

source
MatrixShapedOperators.𝟘 — Constant
const 𝟘 = MatrixShapedOperators.Zero()

The sizeless zero: a special algebra object (not a MatrixShapedOperator) that takes its size from context.

Zero is a singleton type; 𝟘(m, n) resp. 𝟘(n) construct the sized (square) zero:

𝟘(3, 2) isa ZeroOperator

op + 𝟘 == op
op * 𝟘 == zero(op)
source
MatrixShapedOperators.𝟙 — Constant
const 𝟙 = MatrixShapedOperators.One()

The sizeless identity: a special algebra object (not a MatrixShapedOperator) that takes its size from context.

Similar to LinearAlgebra.I but without a scaling factor. One is a singleton type.

𝟙(n) constructs the sized identity:

𝟙(size(x, 1)) * x == x

op + 𝟙 == op + 𝟙(size(op, 1))  # op must be square
op * 𝟙 == op
source
Base.:\ — Method
op::MatrixShapedOperator \ x::AbstractVector
op::MatrixShapedOperator \ X::AbstractMatrix

Eager inverse application (solve) of a square operator to a vector, or to the columns of a matrix as multiple right-hand sides.

Extended help

Like lower_cholesky and logabsdet this is structural: it follows the inverse path of the operator's representation, or throws an ArgumentError if there is none. LinearAlgebra.factorize caches a factorization for repeated solves.

source
MatrixShapedOperators.NumberDomain — Method
MatrixShapedOperators.NumberDomain(A)::Union{Type{Real},Type{Complex}}

Returns Real or Complex, the number domain of the element type of A.

Extended help

Judged by real(T) === T rather than T <: Real, so that tracing number types count as real. Definitive and constant-folding, dispatch on ::Type{Real}/::Type{Complex} to scope real-only behavior (e.g. structural symmetry claims, transpose == adjoint).

source
MatrixShapedOperators.asmatrix — Method
asmatrix(A::MatrixShaped)::AbstractMatrix

Returns A itself if it is an AbstractMatrix, otherwise a cheap, backend-preserving matrix representation of the operator A (e.g. the wrapped matrix of a MatrixAsOperator), which may alias A.

Idempotent dual of asoperator. Not defined for implicit operators - materialize those explicitly via AbstractMatrix(op) (the operator's own storage backend where implemented, dense CPU by default) or Matrix(op) (always dense CPU).

source
MatrixShapedOperators.asoperator — Method
asoperator(A::MatrixShaped, decls::MatrixShapedOperators.StructureTrait...)::MatrixShapedOperator

Returns A if it already is a MatrixShapedOperator, wraps an AbstractMatrix in a MatrixAsOperator.

Additional trait values declare structural knowledge beyond what the matrix type expresses, e.g. asoperator(Σ, IsPosDef()) for a dense covariance matrix known to be positive definite. Declarations are trusted, not verified; they must stay valid if the wrapped matrix is mutated.

source
MatrixShapedOperators.blockdiag_operator — Function
blockdiag_operator(blocks::MatrixShaped...)
blockdiag_operator(blocks::AbstractVector{<:MatrixShaped})

Returns a matrix-shaped operator with the given block-diagonal structure, a BlockDiagOperator in general.

The blocks are MatrixShaped values and need not be square, given either as individual arguments (best for a few blocks of different type) or as a vector (best for many blocks of equal type).

Extended help

Diagonal-like blocks (Diagonal, diagonal_operator, UniformScalingOperator, the sized identity) collapse into a single diagonal operator, sparse blocks into a single wrapped sparse matrix (SparseArrays extension), and RowGramOperator resp. ColGramOperator blocks into a single Gram operator of the block-diagonal of their factors.

source
MatrixShapedOperators.check_mulfunc_operator_support — Method
MatrixShapedOperators.check_mulfunc_operator_support(::Type{OP})::Nothing

Throws an ArgumentError if mulfunc_operator can't construct operators of type OP with the packages currently loaded, returns nothing otherwise.

Specializations of mulfunc_operator for custom operator types should be accompanied by a specialization of check_mulfunc_operator_support that returns nothing.

source
MatrixShapedOperators.colgram_operator — Method
colgram_operator(A::MatrixShaped)

Returns a matrix-shaped operator G == A' * A, the column-Gram construction of a matrix or matrix-shaped operator A.

Typically returns a ColGramOperator:

G = colgram_operator(A)
G * x == A' * (A * x)

G is hermitian (symmetric for real element types) and positive semi-definite by construction; isposdef is true only if A is structurally known to have full column rank.

colgram_factor(colgram_operator(A)) must return an equivalent of asoperator(A).

Equivalent to rowgram_operator(A').

source
MatrixShapedOperators.explicit_mul_impl — Function
MatrixShapedOperators.explicit_mul_impl(a, b)

Eager multiplication: application of an operator to a vector, or to the columns of a matrix as a batch, used by op * x/op * X - in contrast to the lazy operator-operator *.

Subtypes of MatrixShapedOperator implement explicit_mul_impl(op, x::AbstractVector{<:Number}). Matrix application dispatches on BatchedMulStyle to batched_mul_impl or to a generic column-wise fallback. Unchecked - sizes are verified centrally by op * x - so do not call *_impl methods directly.

source
MatrixShapedOperators.gettrait — Method
MatrixShapedOperators.gettrait(traits, C::Type)

Returns the trait of category C from a traitset-created trait set (a tuple or its type), unknowntrait(C) if absent.

Lookup is invariant to the order of the tuple parameters, the canonical order only serves type identity.

source
MatrixShapedOperators.lazy_mul_impl — Function
MatrixShapedOperators.lazy_mul_impl(a, b)

Lazy multiplication: combines operators into a MatrixShapedProduct by default, scalar scaling stays an operator. Base.:(*) involving matrix-shaped operators dispatches to lazy_mul_impl after argument checking, subtypes specialize lazy_mul_impl instead, e.g. for eager collapses.

source
MatrixShapedOperators.ldiv_impl — Function
MatrixShapedOperators.ldiv_impl(op, x::AbstractVecOrMat{<:Number})

Eager inverse application of a square operator to a vector or to the columns of a matrix, the implementation hook behind op \ x. Unchecked, do not call directly.

source
MatrixShapedOperators.lower_cholesky — Function
lower_cholesky(A)

Returns the lower Cholesky factor L of A, so that A == L * L'.

Unlike a general rowgram_factor, the Cholesky factor has lower-triangular structure with non-negative diagonal - positive for positive definite A, where the factor is unique.

Defined where the factor is available structurally or cheaply, implicit operators are never materialized or factorized. Structurally triangular Gram factors are canonicalized by column sign resp. phase factors. Matrix arguments yield matrix factors, operator arguments yield operator factors.

See upper_cholesky for the adjoint factor.

source
MatrixShapedOperators.mulfunc_operator — Function
mulfunc_operator(::Type{T}, sz::Dims{2}, ovp, vop, decls...)
mulfunc_operator(::Type{OP}, ::Type{T}, sz::Dims{2}, ovp, vop, decls...)

Returns a matrix-shaped operator with element type T and size sz, defined by a multiplication function ovp and an adjoint multiplication function vop. Either may be nothing, the missing direction then fails on application; materialization (Matrix target) and conversion to foreign operator types require ovp.

Typically returns a MulFuncOperator:

op = mulfunc_operator(T, sz, ovp, vop)
eltype(op) == T
size(op) == sz
op * x_r == ovp(x_r)
op' * x_l == vop(x_l)

Extended help

Structural properties, which the multiplication functions can't be inspected for, are declared via trait values decls... (see MatrixShapedOperators.StructureTrait), e.g. mulfunc_operator(T, sz, ovp, vop, IsHermitian()); non-square operators can only carry rank declarations. Declarations are trusted, not verified - MatrixShapedOperators.test_operator checks them numerically.

The form with a leading operator type OP builds that operator type instead, including foreign types like LinearMaps.LinearMap (package extensions) and Matrix, which materializes. It lets packages that generate multiplication functions, e.g. autodiff packages, target a requested operator type without depending on its package. Support can be checked via MatrixShapedOperators.check_mulfunc_operator_support.

source
MatrixShapedOperators.rowgram_factor — Method
rowgram_factor(op)

Returns a factor F with op == F * F', always as a matrix-shaped operator, with no guarantee of triangular structure, canonical form or uniqueness - the loose "matrix square root" op^(1/2) of the sampling/coloring convention in the statistics literature.

Use lower_cholesky if canonical triangular structure is required.

For a RowGramOperator this is the stored factor, returned without computation. WoodburyOperators compute their stable factorization; wrapped positive-definite matrices, uniform scalings and structured semi-definite operators (diagonal, zero, block-diagonal) fall back to lower_cholesky.

source
MatrixShapedOperators.rowgram_factor — Method
rowgram_factor(op::WoodburyOperator)

Computes a factor F with op == F * F' via the numerically stable factorization L_A * Q * blockdiag(V', I) with Q R = L_A \ B and V' V = I + R D R' (Zhang et al. 2022, appendix A), in O(n k²).

Extended help

The result is a lazy operator product with structural logabsdet. Requires a positive-definite base with a structural lower_cholesky and a positive-definite total A + B D B' (D itself may be indefinite or singular; throws PosDefException otherwise), with update components B and D that have explicit matrix representations (asmatrix).

source
MatrixShapedOperators.rowgram_operator — Method
rowgram_operator(A::MatrixShaped)

Returns a matrix-shaped operator G == A * A', the row-Gram construction of a matrix or matrix-shaped operator A.

Typically returns a RowGramOperator:

G = rowgram_operator(A)
G * x == A * (A' * x)

G is hermitian (symmetric for real element types) and positive semi-definite by construction, and positive definite if A is structurally known to have full row rank; isposdef also establishes full rank cheaply for square triangular explicit factors.

rowgram_factor(rowgram_operator(A)) must return an equivalent of asoperator(A).

Equivalent to colgram_operator(A'). Use this instead of the plain product A * A', which does not carry the Gram structure.

source
MatrixShapedOperators.test_operator — Function
MatrixShapedOperators.test_operator(
    op::MatrixShapedOperator;
    directions = :both, x, y, u, compare = isapprox, materialize = true, kwargs...
)

Tests the operator contract of op via Test.@test: linearity, adjoint consistency, batched against column-wise application and, with materialize = true (default), agreement of Matrix(op) with application and with the declared structure traits.

Meant for the test suites of operator implementations, in particular of mulfunc_operator seams. Probe arrays default to deterministic values and can be given via the keyword arguments x, y (operator domain) and u (operator codomain); compare (default isapprox) and further keyword arguments are forwarded to the comparison.

For one-directional operators directions = :forward skips the adjoint-consistency test and directions = :adjoint tests adjoint(op) that way instead; probe arrays refer to the tested direction.

!!! Note Requires the Test standard library to be loaded.

source
MatrixShapedOperators.trait_add — Method
MatrixShapedOperators.trait_add(x, y)
MatrixShapedOperators.trait_mul(x, y)
MatrixShapedOperators.trait_blockdiag(x, y)
MatrixShapedOperators.trait_adjoint(x)
MatrixShapedOperators.trait_scale(x)

Trait algebra: returns the resulting trait of one category under operator addition, multiplication (composition), block-diagonal concatenation, adjoint and scaling by a real scalar of unknown value. Conservative Unknown* fallbacks; each trait category defines its preservation rules in place.

source
MatrixShapedOperators.trait_category — Method
MatrixShapedOperators.trait_category(t)::Type

Returns the abstract trait type a trait value belongs to, e.g. PosDefNess for IsPosDef(): the last supertype before Any.

source
MatrixShapedOperators.traitset — Method
MatrixShapedOperators.traitset(decls...)::Tuple

Returns the canonical trait set for operator types that store their structural properties as a tuple (type) of traits.

The result is a tuple of the non-Unknown* trait values, deduplicated (the first non-Unknown* declaration per category wins, Unknown* values carry no information), sorted by category name. Equal trait sets have equal types. Implications are closed: unitarity and positive definiteness imply full row rank, positive semi-definiteness implies hermitian, diagonality implies symmetric.

source
MatrixShapedOperators.traitset_add — Method
MatrixShapedOperators.traitset_add(ta::Tuple, tb::Tuple)::Tuple
MatrixShapedOperators.traitset_mul(ta::Tuple, tb::Tuple)::Tuple
MatrixShapedOperators.traitset_blockdiag(ta::Tuple, tb::Tuple)::Tuple
MatrixShapedOperators.traitset_adjoint(ts::Tuple)::Tuple
MatrixShapedOperators.traitset_scale(ts::Tuple)::Tuple

Trait algebra on whole trait sets: applies trait_add etc. per category and re-canonicalizes the result.

source
MatrixShapedOperators.unknowntrait — Function
MatrixShapedOperators.unknowntrait(::Type{C})

Returns the Unknown* value of trait category C. Implemented by each trait category, including third-party ones.

source
MatrixShapedOperators.woodbury_operator — Method
woodbury_operator(A::MatrixShaped, B::MatrixShaped, D::MatrixShaped)

Returns a matrix-shaped operator op that represents A + B * D * B' for a square base A, an n × k update factor B (typically tall, k small) and a square symmetric (hermitian for complex element types) middle matrix D.

Typically returns a WoodburyOperator:

op = woodbury_operator(A, B, D)
op * x == A * x + B * (D * (B' * x))

A 𝟘 base yields a pure low-rank B * D * B', which supports application and traits but no solves or determinants - the push-through formulas require an invertible base.

The operator is symmetric/hermitian if A is, and positive (semi-)definite if A is and D is positive semi-definite. Components are MatrixShaped; keep them structured (Diagonal base, triangular D) to enable structural operations.

source