by Geoffroy Leconte and Dominique Orban
LinearOperators.jl is a package for matrix-like operators. Linear operators are defined by how they act on a vector, which is useful in a variety of situations where you don't want to materialize the matrix.
This section of the documentation describes a few uses of LinearOperators.
Operators may be defined from matrices and combined using the usual operations, but the result is deferred until the operator is applied.
using LinearOperators, Random, SparseArrays
Random.seed!(0)
A1 = rand(5,7)
A2 = sprand(7,3,.3)
op1 = LinearOperator(A1)
op2 = LinearOperator(A2)
op = op1 * op2 # Does not form A1 * A2
x = rand(3)
y = op * x
5-element Vector{Float64}:
0.24766385713448147
0.34471427066907173
0.6399025892338569
0.246664873173519
0.5258420478788777
Operators may be defined to represent (approximate) inverses.
using LinearAlgebra
A = rand(5,5)
A = A' * A
op = opCholesky(A) # Use, e.g., as a preconditioner
v = rand(5)
norm(A \ v - op * v) / norm(v)
1.148733996146629e-13
In this example, the Cholesky factor is computed only once and can be used many times transparently.
It is often useful to reuse the memory used by the operator. For that reason, we can use mul! on operators as if we were using matrices using preallocated vectors:
m, n = 50, 30
A = rand(m, n) + im * rand(m, n)
op = LinearOperator(A)
v = rand(n)
res = zeros(eltype(A), m)
res2 = copy(res)
mul!(res2, op, v) # compile 3-args mul!
al = @allocated mul!(res, op, v) # op * v, store result in res
println("Allocation of LinearOperator mul! product = $al")
v = rand(n)
α, β = 2.0, 3.0
mul!(res2, op, v, α, β) # compile 5-args mul!
al = @allocated mul!(res, op, v, α, β) # α * op * v + β * res, store result in res
println("Allocation of LinearOperator mul! product = $al")
Allocation of LinearOperator mul! product = 0
Allocation of LinearOperator mul! product = 0
Operators may be defined from functions. They have to be based on the 5-arguments mul! function. In the example below, the transposed isn't defined, but it may be inferred from the conjugate transposed. Missing operations are represented as nothing. You will have deal with cases where β == 0 and β != 0 separately because * will allocate an uninitialized res vector that may contain NaN values, and 0 * NaN == NaN.
using FFTW
function mulfft!(res, v, α, β)
if β == 0
res .= α .* fft(v)
else
res .= α .* fft(v) .+ β .* res
end
end
function mulifft!(res, w, α, β)
if β == 0
res .= α .* ifft(w)
else
res .= α .* ifft(w) .+ β .* res
end
end
dft = LinearOperator(ComplexF64, 10, 10, false, false,
mulfft!,
nothing, # will be inferred
mulifft!)
x = rand(10)
y = dft * x
norm(dft' * y - x) # DFT is a unitary operator
3.3612925285989963e-16
transpose(dft) * y
10-element Vector{ComplexF64}:
0.7113168325963566 - 0.0im
0.2148297961020408 - 0.0im
0.9453427110792467 - 0.0im
0.7309321298104021 - 0.0im
0.23194328498409336 - 0.0im
0.9501874396162999 - 0.0im
0.5123847829172379 - 0.0im
0.9037088931078092 - 0.0im
0.005729552514365821 + 0.0im
0.2588088046728426 - 0.0im
Another example:
function customfunc!(res, v, α, β)
if β == 0
res[1] = (v[1] + v[2]) * α
res[2] = v[2] * α
else
res[1] = (v[1] + v[2]) * α + res[1] * β
res[2] = v[2] * α + res[2] * β
end
end
function tcustomfunc!(res, w, α, β)
if β == 0
res[1] = w[1] * α
res[2] = (w[1] + w[2]) * α
else
res[1] = w[1] * α + res[1] * β
res[2] = (w[1] + w[2]) * α + res[2] * β
end
end
op = LinearOperator(Float64, 2, 2, false, false,
customfunc!,
nothing,
tcustomfunc!)
Linear operator
nrow: 2
ncol: 2
eltype: Float64
symmetric: false
hermitian: false
nprod: 0
ntprod: 0
nctprod: 0
Operators can also be defined with the 3-args mul! function:
op2 = LinearOperator(Float64, 2, 2, false, false,
(res, v) -> customfunc!(res, v, 1.0, 0.),
nothing,
(res, w) -> tcustomfunc!(res, w, 1.0, 0.))
Linear operator
nrow: 2
ncol: 2
eltype: Float64
symmetric: false
hermitian: false
nprod: 0
ntprod: 0
nctprod: 0
When using the 5-args mul! with the above operator, some vectors will be allocated (only at the first call):
res, a = zeros(2), rand(2)
mul!(res, op2, a) # compile
println("allocations 1st call = ", @allocated mul!(res, op2, a, 2.0, 3.0))
println("allocations 2nd call = ", @allocated mul!(res, op2, a, 2.0, 3.0))
allocations 1st call = 80
allocations 2nd call = 0
Make sure that the type passed to LinearOperator is correct, otherwise errors may occur.
dft = LinearOperator(Float64, 10, 10, false, false,
mulfft!,
nothing,
mulifft!)
v = rand(10)
println("eltype(dft) = $(eltype(dft))")
println("eltype(v) = $(eltype(v))")
eltype(dft) = Float64
eltype(v) = Float64
try
dft * v # ERROR: expected Vector{Float64}
catch ex
println("ex = $ex")
end
ex = InexactError(:Float64, Float64, 0.47507530806279163 - 0.3832301880068187im)
try
Matrix(dft) # ERROR: tried to create a Matrix of Float64
catch ex
println("ex = $ex")
end
# Using external modules
ex = InexactError(:Float64, Float64, 0.8090169943749475 - 0.5877852522924731im)
It is possible to use certain modules made for matrices that do not need to access specific elements of their input matrices, and only use operations implemented within LinearOperators, such as mul!, *, +, ... For example, we show the solution of a linear system using Krylov.jl:
using Krylov
A = rand(5, 5)
opA = LinearOperator(A)
opAAT = opA + opA'
b = rand(5)
(x, stats) = minres(opAAT, b)
norm(b - opAAT * x)
2.3624810325472306e-13
Two other useful operators are the Limited-Memory BFGS in forward and inverse form.
B = LBFGSOperator(20)
H = InverseLBFGSOperator(20)
r = 0.0
for i = 1:100
global r
s = rand(20)
y = rand(20)
push!(B, s, y)
push!(H, s, y)
r += norm(B * H * s - s)
end
r
3.919077465054229e-13
There is also a LSR1 operator that behaves similarly to these two.
The restriction operator restricts a vector to a set of indices.
v = collect(1:5)
R = opRestriction([2;5], 5)
R * v
2-element Vector{Int64}:
2
5
Notice that it corresponds to a matrix with rows of the identity given by the indices.
Matrix(R)
2×5 Matrix{Int64}:
0 1 0 0 0
0 0 0 0 1
The extension operator is the transpose of the restriction. It extends a vector with zeros.
v = collect(1:2)
E = opExtension([2;5], 5)
E * v
5-element Vector{Int64}:
0
1
0
0
2
With these operators, we define the slices of an operator op.
A = rand(5,5)
opA = LinearOperator(A)
I = [1;3;5]
J = 2:4
A[I,J] * ones(3)
3-element Vector{Float64}:
1.043315409692812
1.126646200633167
1.2315932943898407
opRestriction(I, 5) * opA * opExtension(J, 5) * ones(3)
3-element Vector{Float64}:
1.043315409692812
1.126646200633167
1.2315932943898407
A main difference with matrices, is that slices do not return vectors nor numbers.
opA[1,:] * ones(5)
1-element Vector{Float64}:
1.9634567072673046
opA[:,1] * ones(1)
5-element Vector{Float64}:
0.02203447865143171
0.18781833339603016
0.1692361439290926
0.20276666600590854
0.24969824326802326
opA[1,1] * ones(1)
1-element Vector{Float64}:
0.02203447865143171