Linear Solve with Caching Interface

Often, one may want to cache information that is reused between different linear solves. For example, if one is going to perform:

A \ b1
A \ b2

then it would be more efficient to LU-factorize one time and reuse the factorization:

A_lu = LA.lu!(A)
A_lu \ b1
A_lu \ b2

LinearSolve.jl's caching interface automates this process to use the most efficient means of solving and resolving linear systems. To do this with LinearSolve.jl, you simply init a cache, solve, replace b, and solve again. This looks like:

import LinearSolve as LS
import LinearAlgebra as LA

n = 4
A = rand(n, n)
b1 = rand(n);
b2 = rand(n);
prob = LS.LinearProblem(A, b1)

linsolve = LS.init(prob)
sol1 = LS.solve!(linsolve)

retcode: Success
u: 4-element Vector{Float64}:
  0.22493008877212145
 -1.4645610053561515
  1.8680183086605542
  0.42941138057983713

linsolve.b = b2
sol2 = LS.solve!(linsolve)

sol2.u

4-element Vector{Float64}:
  0.4023828510606506
 -0.5950100132572225
  1.0521611341573458
  0.32130759566461314

Then refactorization will occur when a new A is given:

A2 = rand(n, n)
linsolve.A = A2
sol3 = LS.solve!(linsolve)

sol3.u

4-element Vector{Float64}:
  1.0600771644188434
  1.7968904486444568
 -1.4707858563961487
  0.5801102236596681

The factorization occurs on the first solve, and it stores the factorization in the cache. You can retrieve this cache via sol.cache, which is the same object as the init, but updated to know not to re-solve the factorization.

The advantage of course with import LinearSolve.jl in this form is that it is efficient while being agnostic to the linear solver. One can easily swap in iterative solvers, sparse solvers, etc. and it will do all the tricks like caching the symbolic factorization if the sparsity pattern is unchanged.

Reusing the Symbolic Factorization with Sparse Matrices

For sparse factorizations such as UMFPACKFactorization and KLUFactorization, a significantly expensive part of a solve is often the symbolic factorization: the analysis of the sparsity pattern that determines the fill-in and elimination ordering. The subsequent numeric factorization, which computes the actual L and U factor values, is not cheap but must be recomputed for every change in the values of A. Thus storing the symbolic factorization can lead to some sizable gains.

A very common case (for example, the Newton steps of a nonlinear solve or the time steps of an implicit ODE/PDE integrator) is that A keeps the same sparsity pattern across solves while only its stored values change. In that case, the symbolic factorization can be reused and only the numeric factorization needs to be redone. This is controlled by the reuse_symbolic keyword on the sparse factorization algorithms, which defaults to true — so you get this behavior automatically and do not need to set it:

solver = LS.UMFPACKFactorization()   # reuse_symbolic = true by default

When a new A with an unchanged pattern is given, LinearSolve.jl reuses the cached symbolic factorization and performs only the numeric refactorization. By default it also runs a check_pattern pass to confirm the pattern really is unchanged; if your pattern is guaranteed constant, you can set check_pattern = false to skip that check for a little extra speed. Set reuse_symbolic = false only if the sparsity pattern may change between solves.

Updating the values of A

The caching interface tracks whether the factorization is stale through the assignment linsolve.A = A. Mutating the stored values of A in place does not, on its own, tell the cache that A has changed, so the next solve! will silently reuse the old factorization and return a stale (wrong) result:

linsolve.A.nzval .*= 2   # in-place edit: NOT seen by the cache
LS.solve!(linsolve)      # reuses the old factorization -> wrong answer

The supported way to signal that the values changed is to assign A back into the cache. Reassigning the same (now-mutated) matrix object is fine and still reuses the symbolic factorization when the pattern is unchanged:

import LinearSolve as LS
import LinearAlgebra as LA
import SparseArrays

n = 100
A = SparseArrays.sprand(n, n, 0.05) + 10.0 * LA.I
b = rand(n)
prob = LS.LinearProblem(A, b)

linsolve = LS.init(prob, LS.UMFPACKFactorization())
sol1 = LS.solve!(linsolve)

# Update the values of A (same sparsity pattern), then signal the change:
A.nzval .*= 2
linsolve.A = A      # triggers refactorization; symbolic part is reused
sol2 = LS.solve!(linsolve)

sol2.u

100-element Vector{Float64}:
  0.009641204824247024
  0.025945614581122104
  0.008929854148420605
  0.015362207502316571
  0.03462078799278962
  0.025749168213308166
 -0.001708566780397459
  0.006735787308315144
  0.038258658918430444
  0.019912373268050114
  ⋮
  0.040237222740472124
  0.04488068648076178
  0.01014360797619324
  0.024069562158054933
 -0.0004443506636360997
  0.006629958745332541
  0.011127402725001494
  0.02869098213810123
  0.034036938952320656

Because the pattern is unchanged, sol2 reuses the cached symbolic factorization and only redoes the (cheaper) numeric factorization.