Edit on GitHub

Tutorial

This section of the documentation describes a few uses of LinearOperators.

Using matrices

Operators may be defined from matrices and combined using the usual operations, but the result is deferred until the operator is applied.

using LinearOperators, SparseArrays
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 Array{Float64,1}:
 0.7082889552644532
 0.6721427712594621
 0.9327904433408708
 0.9496589165078558
 0.5658463473796669

Inverse

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)
2.4757415607659417e-12

In this example, the Cholesky factor is computed only once and can be used many times transparently.

Using functions

Operators may be defined from functions. In the example below, the transposed isn't defined, but it may be inferred from the conjugate transposed. Missing operations are represented as nothing.

using FFTW
dft = LinearOperator(10, 10, false, false,
                     v -> fft(v),
                     nothing,       # will be inferred
                     w -> ifft(w))
x = rand(10)
y = dft * x
norm(dft' * y - x)  # DFT is an orthogonal operator
2.6622212897024595e-16
transpose(dft) * y
10-element Array{Complex{Float64},1}:
   0.8786013899885341 - 0.0im
   0.8570663597525118 - 0.0im
  0.19045332274089394 - 0.0im
   0.8320134987734346 - 0.0im
   0.3208171451335831 - 0.0im
   0.4785627240050898 - 0.0im
 0.040681202859401355 - 0.0im
   0.3160498443874862 - 0.0im
   0.4562318769037463 + 0.0im
   0.8602132056965561 - 0.0im

By default a linear operator defined by functions and that is neither symmetric nor hermitian will have element type Complex128. This behavior may be overridden by specifying the type explicitly, e.g.,

op = LinearOperator(Float64, 10, 10, false, false,
                    v -> [v[1] + v[2]; v[2]],
                    nothing,
                    w -> [w[1]; w[1] + w[2]])
Linear operator
  nrow: 10
  ncol: 10
  eltype: Float64
  symmetric: false
  hermitian: false
  nprod:   0
  ntprod:  0
  nctprod: 0

Make sure that the type passed to LinearOperator is correct, otherwise errors may occur.

dft = LinearOperator(Float64, 10, 10, false, false,
                     v -> fft(v),
                     nothing,
                     w -> ifft(w))
v = rand(10)
println("eltype(dft)         = $(eltype(dft))")
println("eltype(v)           = $(eltype(v))")
println("eltype(dft.prod(v)) = $(eltype(dft.prod(v)))")
# dft * v     # ERROR: expected Vector{Float64}
# Matrix(dft) # ERROR: tried to create a Matrix of Float64
eltype(dft)         = Float64
eltype(v)           = Float64
eltype(dft.prod(v)) = Complex{Float64}

Limited memory BFGS and SR1

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
4.109005201248445e-13

There is also a LSR1 operator that behaves similarly to these two.

Restriction, extension and slices

The restriction operator restricts a vector to a set of indices.

v = collect(1:5)
R = opRestriction([2;5], 5)
R * v
2-element Array{Int64,1}:
 2
 5

Notice that it corresponds to a matrix with rows of the identity given by the indices.

Matrix(R)
2×5 Array{Int64,2}:
 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 Array{Int64,1}:
 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 Array{Float64,1}:
 1.8121605051122305
 2.1451557985063983
 2.5064062796411193
opRestriction(I, 5) * opA * opExtension(J, 5) * ones(3)
3-element Array{Float64,1}:
 1.8121605051122305
 2.1451557985063983
 2.5064062796411193

A main difference with matrices, is that slices do not return vectors nor numbers.

opA[1,:] * ones(5)
1-element Array{Float64,1}:
 2.740909849361878
opA[:,1] * ones(1)
5-element Array{Float64,1}:
 0.8393629797922026 
 0.8628690936422081 
 0.42155359279030513
 0.9088682954119938 
 0.340029892097079  
opA[1,1] * ones(1)
1-element Array{Float64,1}:
 0.8393629797922026

Preallocated Operators

Operators created from matrices are very practical, however, it is often useful to reuse the memory used by the operator. For that use, we can use PreallocatedLinearOperator(A) to create an operator that reuses the memory.

m, n = 50, 30
A = rand(m, n) + im * rand(m, n)
op1 = PreallocatedLinearOperator(A)
op2 = LinearOperator(A)
v = rand(n)
al = @allocated op1 * v
println("Allocation of PreallocatedLinearOperator product = $al")
v = rand(n)
al = @allocated op2 * v
println("Allocation of LinearOperator product = $al")
Allocation of PreallocatedLinearOperator product = 34397200
Allocation of LinearOperator product = 16835493

Notice the memory overwrite:

Av = op1 * v
w = rand(n)
Aw = op1 * w
Aw === Av
true

which doesn't happen on LinearOperator.

Av = op2 * v
w = rand(n)
Aw = op2 * w
Aw === Av
false

You can also provide the memory to be used.

Mv  = Array{ComplexF64}(undef, m)
Mtu = Array{ComplexF64}(undef, n)
Maw = Array{ComplexF64}(undef, n)
op  = PreallocatedLinearOperator(Mv, Mtu, Maw, A)
v, u, w = rand(n), rand(m), rand(m)
(Mv === op * v, Mtu === transpose(op) * u, Maw === adjoint(op) * w)
(true, true, true)