Reference

A = Matrix(op)

Materialize an operator as a dense array using op.ncol products.

Hermitian(op, uplo=:U)

Symmetric(op, uplo=:U)

DiagonalAndrei(d)

Construct a linear operator that represents a diagonal quasi-Newton approximation as described in

Andrei, N. A diagonal quasi-Newton updating method for unconstrained optimization. https://doi.org/10.1007/s11075-018-0562-7

The approximation satisfies the weak secant equation and is not guaranteed to be positive definite.

Arguments

• d::AbstractVector: initial diagonal approximation.

DiagonalAndrei(d)

Construct a linear operator that represents a diagonal quasi-Newton approximation as described in

Andrei, N. A diagonal quasi-Newton updating method for unconstrained optimization. https://doi.org/10.1007/s11075-018-0562-7

The approximation satisfies the weak secant equation and is not guaranteed to be positive definite.

Arguments

• d::AbstractVector: initial diagonal approximation.

DiagonalBFGS(d)

A diagonal approximation of the BFGS update inspired by Marnissi, Y., Chouzenoux, E., Benazza-Benyahia, A., & Pesquet, J. C. (2020). Majorize–minimize adapted Metropolis–Hastings algorithm. https://ieeexplore.ieee.org/abstract/document/9050537.

Arguments

• d::AbstractVector: initial diagonal approximation.

DiagonalBFGS(d)

A diagonal approximation of the BFGS update inspired by Marnissi, Y., Chouzenoux, E., Benazza-Benyahia, A., & Pesquet, J. C. (2020). Majorize–minimize adapted Metropolis–Hastings algorithm. https://ieeexplore.ieee.org/abstract/document/9050537.

Arguments

• d::AbstractVector: initial diagonal approximation.

DiagonalPSB(d)

Construct a linear operator that represents a diagonal PSB quasi-Newton approximation as described in

M. Zhu, J. L. Nazareth and H. Wolkowicz The Quasi-Cauchy Relation and Diagonal Updating. SIAM Journal on Optimization, vol. 9, number 4, pp. 1192-1204, 1999. https://doi.org/10.1137/S1052623498331793.

The approximation satisfies the weak secant equation and is not guaranteed to be positive definite.

Arguments

• d::AbstractVector: initial diagonal approximation.

DiagonalPSB(d)

Construct a linear operator that represents a diagonal PSB quasi-Newton approximation as described in

M. Zhu, J. L. Nazareth and H. Wolkowicz The Quasi-Cauchy Relation and Diagonal Updating. SIAM Journal on Optimization, vol. 9, number 4, pp. 1192-1204, 1999. https://doi.org/10.1137/S1052623498331793.

The approximation satisfies the weak secant equation and is not guaranteed to be positive definite.

Arguments

• d::AbstractVector: initial diagonal approximation.

A data type to hold information relative to LBFGS operators.

A type for limited-memory BFGS approximations.

LBFGSOperator(T, n; [mem=5, scaling=true])
LBFGSOperator(n; [mem=5, scaling=true])

Construct a limited-memory BFGS approximation in forward form. If the type T is omitted, then Float64 is used.

A data type to hold information relative to LSR1 operators.

A type for limited-memory SR1 approximations.

LSR1Operator(T, n; [mem=5, scaling=false)
LSR1Operator(n; [mem=5, scaling=false)

Construct a limited-memory SR1 approximation in forward form. If the type T is omitted, then Float64 is used.

Base type to represent a linear operator. The usual arithmetic operations may be applied to operators to combine or otherwise alter them. They can be combined with other operators, with matrices and with scalars. Operators may be transposed and conjugate-transposed using the usual Julia syntax.

LinearOperator(M::AbstractMatrix{T}; symmetric=false, hermitian=false, S = Vector{T}) where {T}

Construct a linear operator from a dense or sparse matrix. Use the optional keyword arguments to indicate whether the operator is symmetric and/or hermitian. Change S to use LinearOperators on GPU.

LinearOperator(type::Type{T}, nrow, ncol, symmetric, hermitian, prod!,
                [tprod!=nothing, ctprod!=nothing],
                S = Vector{T}) where {T}

Construct a linear operator from functions where the type is specified as the first argument. Change S to use LinearOperators on GPU. Notice that the linear operator does not enforce the type, so using a wrong type can result in errors. For instance,

A = [im 1.0; 0.0 1.0] # Complex matrix
function mulOp!(res, M, v, α, β)
  mul!(res, M, v, α, β)
end
op = LinearOperator(Float64, 2, 2, false, false, 
                    (res, v, α, β) -> mulOp!(res, A, v, α, β), 
                    (res, u, α, β) -> mulOp!(res, transpose(A), u, α, β), 
                    (res, w, α, β) -> mulOp!(res, A', w, α, β))
Matrix(op) # InexactError

The error is caused because Matrix(op) tries to create a Float64 matrix with the contents of the complex matrix A.

Using * may generate a vector that contains NaN values. This can also happen if you use the 3-args mul! function with a preallocated vector such as Vector{Float64}(undef, n). To fix this issue you will have to deal with the cases β == 0 and β != 0 separately:

d1 = [2.0; 3.0]
function mulSquareOpDiagonal!(res, d, v, α, β::T) where T
  if β == zero(T)
    res .= α .* d .* v
  else 
    res .= α .* d .* v .+ β .* res
  end
end
op = LinearOperator(Float64, 2, 2, true, true, 
                    (res, v, α, β) -> mulSquareOpDiagonal!(res, d, v, α, β))

It is possible to create an operator with the 3-args mul!. In this case, using the 5-args mul! will generate storage vectors.

A = rand(2, 2)
op = LinearOperator(Float64, 2, 2, false, false, 
                    (res, v) -> mul!(res, A, v),
                    (res, w) -> mul!(res, A', w))

The 3-args mul! also works when applying the operator on a matrix.

LinearOperator(M::Hermitian{T}, S = Vector{T}) where {T}

Constructs a linear operator from a Hermitian matrix. If its elements are real, it is also symmetric. Change S to use LinearOperators on GPU.

LinearOperator(M::SymTridiagonal{T}, S = Vector{T}) where {T}

Constructs a linear operator from a symmetric tridiagonal matrix. If its elements are real, it is also Hermitian, otherwise complex symmetric. Change S to use LinearOperators on GPU.

LinearOperator(M::Symmetric{T}, S = Vector{T}) where {T <: Real} =

Construct a linear operator from a symmetric real square matrix M. Change S to use LinearOperators on GPU.

A data type to hold information relative to shifted operators.

ShiftedOperator(H, σ=0)

Construct a linear operator representing op = H + σI.

Implementation of a spectral gradient quasi-Newton approximation described in

Birgin, E. G., Martínez, J. M., & Raydan, M. Spectral Projected Gradient Methods: Review and Perspectives. https://doi.org/10.18637/jss.v060.i03

    SpectralGradient(σ, n)

Construct a spectral gradient Hessian approximation. The approximation is defined as σI.

Arguments

• σ::Real: initial positive multiple of the identity;
• n::Int: operator size.

TimedLinearOperator(op)

Creates a linear operator instrumented with timers from TimerOutputs.

opEye()

Identity operator.

opI = opEye()
v = rand(5)
@assert opI * v === v

opEye(T, n; S = Vector{T})
opEye(n)

Identity operator of order n and of data type T (defaults to Float64). Change S to use LinearOperators on GPU.

opEye(T, nrow, ncol; S = Vector{T})
opEye(nrow, ncol)

Rectangular identity operator of size nrowxncol and of data type T (defaults to Float64). Change S to use LinearOperators on GPU.

kron(A, B)

Kronecker tensor product of A and B in linear operator form, if either or both are linear operators. If both A and B are matrices, then Base.kron is used.

push!(op, s, y)

Push a new {s,y} pair into a L-SR1 operator.

push!(op, s, y)
push!(op, s, y, Bs)
push!(op, s, y, α, g)
push!(op, s, y, α, g, Bs)

Push a new {s,y} pair into a L-BFGS operator. The second calling sequence is used for forward updating damping, using the preallocated vector Bs. If the operator is damped, the first call will create Bs and call the second call. The third and fourth calling sequences are used in inverse LBFGS updating in conjunction with damping, where α is the most recent steplength and g the gradient used when solving d=-Hg.

show(io, op)

Display basic information about a linear operator.

m = size(op, d)

Return the size of a linear operator along dimension d.

m, n = size(op)

Return the size of a linear operator as a tuple.

diag(op)
diag!(op, d)

Extract the diagonal of a L-BFGS operator in forward mode.

diag(op)
diag!(op, d)

Extract the diagonal of a L-SR1 operator in forward mode.

ishermitian(op)

Determine whether the operator is Hermitian.

issymmetric(op)

Determine whether the operator is symmetric.

BlockDiagonalOperator(M1, M2, ..., Mn; S = promote_type(storage_type.(M1, M2, ..., Mn)))

Creates a block-diagonal linear operator:

[ M1           ]
[    M2        ]
[       ...    ]
[           Mn ]

Change S to use LinearOperators on GPU.

InverseLBFGSOperator(T, n, [mem=5; scaling=true])
InverseLBFGSOperator(n, [mem=5; scaling=true])

Construct a limited-memory BFGS approximation in inverse form. If the type T is omitted, then Float64 is used.

check_ctranspose(op)

Cheap check that the operator and its conjugate transposed are related.

check_hermitian(op)

Cheap check that the operator is Hermitian.

check_positive_definite(op; semi=false)

Cheap check that the operator is positive (semi-)definite.

has_args5(op)

Determine whether the operator can work with the 5-args mul!. If false, storage vectors will be generated at the first call of the 5-args mul!. No additional vectors are generated when using the 3-args mul!.

normest(S) estimates the matrix 2-norm of S. This function is an adaptation of Matlab's built-in NORMEST. This method allocates.

Inputs: S –- Matrix or LinearOperator type, tol –- relative error tol, default(or -1) Machine eps maxiter –- maximum iteration, default 100

Returns: e –- the estimated norm cnt –- the number of iterations used

opCholesky(M, [check=false])

Inverse of a Hermitian and positive definite matrix as a linear operator using its Cholesky factorization. The factorization is computed only once. The optional check argument will perform cheap hermicity and definiteness checks. This Operator is not in-place when using mul!.

opDiagonal(d)

Diagonal operator with the vector d on its main diagonal.

opDiagonal(nrow, ncol, d)

Rectangular diagonal operator of size nrow-by-ncol with the vector d on its main diagonal.

Z = opExtension(I, ncol)
Z = opExtension(:, ncol)

Creates a LinearOperator extending a vector of size length(I) to size ncol, where the position of the elements on the new vector are given by the indices I. The operation w = Z * v is equivalent to w = zeros(ncol); w[I] = v.

Z = opExtension(k, ncol)

Alias for opExtension([k], ncol).

opHermitian(A)

A symmetric/hermitian operator based on a matrix.

opHermitian(d, A)

A symmetric/hermitian operator based on the diagonal d and lower triangle of A.

opHouseholder(h)

Apply a Householder transformation defined by the vector h. The result is x -> (I - 2 h hᵀ) x.

opInverse(M; symm=false, herm=false)

Inverse of a matrix as a linear operator using \. Useful for triangular matrices. Note that each application of this operator applies \. This Operator is not in-place when using mul!.

opLDL(M, [check=false])

Inverse of a symmetric matrix as a linear operator using its LDLᵀ factorization if it exists. The factorization is computed only once. The optional check argument will perform a cheap hermicity check.

If M is sparse and real, then only the upper triangle should be stored in order to use LDLFactorizations.jl:

using LDLFactorizations
triu!(M)
opLDL(Symmetric(M, :U))

opOnes(T, nrow, ncol; S = Vector{T})
opOnes(nrow, ncol)

Operator of all ones of size nrow-by-ncol of data type T (defaults to Float64). Change S to use LinearOperators on GPU.

Z = opRestriction(I, ncol)
Z = opRestriction(:, ncol)

Creates a LinearOperator restricting a ncol-sized vector to indices I. The operation Z * v is equivalent to v[I]. I can be :.

Z = opRestriction(k, ncol)

Alias for opRestriction([k], ncol).

opZeros(T, nrow, ncol; S = Vector{T})
opZeros(nrow, ncol)

Zero operator of size nrow-by-ncol, of data type T (defaults to Float64). Change S to use LinearOperators on GPU.

reset!(op)

Reset the product counters of a linear operator.

reset!(op)

Resets the LBFGS data of the given operator.

reset!(op)

Resets the LSR1 data of the given operator.

reset!(op::AbstractDiagonalQuasiNewtonOperator)

Reset the diagonal data of the given operator.

reset!(data)

Resets the given LBFGS data.

reset!(data)

Reset the given LSR1 data.

solveshiftedsystem!(x, B, b, σ)

Solve linear system (B + σI) x = b, where B is a forward L-BFGS operator and σ ≥ 0.

Parameters

• x::AbstractVector{T}: preallocated vector of length n that is used to store the solution x.
• B::LBFGSOperator: forward L-BFGS operator that models a matrix of size n x n.
• b::AbstractVector{T}: right-hand side vector of length n.
• σ::T: nonnegative shift.

Returns

• x::AbstractVector{T}: solution vector x of length n.

Method

The method uses a two-loop recursion-like approach with modifications to handle the shift σ.

Example

using Random

# Problem setup
n = 100  # size of the problem
mem = 10   # L-BFGS memory size
scaling = true  # enable scaling

# Create an L-BFGS operator
B = LBFGSOperator(n, mem = mem, scaling = scaling)

# Add random {s, y} pairs to the L-BFGS operator
for _ = 1:10
    s = rand(n)   
    y = rand(n)   
    push!(B, s, y)  # Add the {s, y} pair to B
end

# Prepare vectors for the system
x = zeros(n)   # Preallocated solution vector
b = rand(n)        # Right-hand side vector
σ = 0.1            # Small shift value

# Solve the shifted system
result = solve_shifted_system!(x, B, b, σ)

# Check that the solution is close enough (residual test)
@assert norm(B * x + σ * x - b) / norm(b) < 1e-8

References

Erway, J. B., Jain, V., & Marcia, R. F. Shifted L-BFGS Systems. Optimization Methods and Software, 29(5), pp. 992-1004, 2014.