Linear Solvers

Objects of type LinearSolver are used to solve LinearProblems, i.e. we want to find $x$ for given $A$ and $y$ such that

\[ Ax = y\]

is satisfied.

A linear system can be created with:

using SimpleSolvers

A = [(0. + 1e-6) 1. 2.; 3. 4. 5.; 6. 7. 8.]
y = [1., 2., 3.]
ls = LinearProblem(A, y)

Note that we here use the matrix:

\[A = \begin{pmatrix} 0 + \varepsilon & 1 & 2 \\ 3 & 4 & 5 \\ 6 & 7 & 8 \end{pmatrix}.\]

This matrix would be singular if we had $\varepsilon = 0$ because $2\cdot\begin{pmatrix} 3 \\ 4 \\ 5 \end{pmatrix} - \begin{pmatrix} 6 \\ 7 \\ 8 \end{pmatrix} = \begin{pmatrix} 0 \\ 1 \\ 2 \end{pmatrix}.$ So by choosing $\varepsilon = 10^{-6}$ the matrix is ill-conditioned.

We first solve LinearProblem with an lu solver (using LU and solve) in double precision and without pivoting:

lu = LU(; pivot = false)
y¹ = solve(lu, ls)

3-element Vector{Float64}:
  0.0
 -0.33333333333333326
  0.6666666666666666

We check the result:

A * y¹

3-element Vector{Float64}:
 1.0
 2.0
 3.0

We now do the same in single precision:

Aˢ = Float32.(A)
yˢ = Float32.(y)
lsˢ = LinearProblem(Aˢ, yˢ)
y² = solve(lu, lsˢ)

3-element Vector{Float32}:
 -0.11920929
  1.666669f-7
  0.5

and again check the result:

Aˢ * y²

3-element Vector{Float32}:
 1.0
 2.1423728
 3.2847455

As we can see the computation of the factorization returns a wrong solution. If we use pivoting however, the problem can also be solved with single precision:

lu = LU(; pivot = true)
y³ = solve(lu, lsˢ)

3-element Vector{Float32}:
  0.33333373
 -1.0000004
  1.0

Aˢ * y³

3-element Vector{Float32}:
 1.0
 1.9999998
 3.0

Solving the System with Built-In Functionality from the LinearAlgebra Package

We further try to solve the system with the inv operator from the LinearAlgebra package. First in double precision:

inv(A) * y

3-element Vector{Float64}:
  0.0
 -0.33333333395421505
  0.6666666669771075

And also in single precision

inv(Aˢ) * yˢ

3-element Vector{Float32}:
  0.25
 -0.5
  0.75

In single precision the result is completely wrong as can also be seen by computing:

inv(Aˢ) * Aˢ

3×3 Matrix{Float32}:
 1.0   0.5   1.0
 0.0   1.0  -2.0
 0.0  -0.5   1.0

If we however write:

Aˢ \ yˢ

3-element Vector{Float32}:
  0.08333349
 -0.5000001
  0.75

we again obtain a correct-looking result, as LinearAlgebra.\ uses an algorithm very similar to factorize! in SimpleSolvers.

Delegating the Factorization to LAPACK

LU is a self-contained scalar implementation. That is what makes the comparison above possible — the pivoting strategy is ours to choose — and for small systems its static-matrix cache means a factorization allocates nothing at all. It does not scale, though: the factorization is $\mathcal{O}(n^3)$ scalar operations with no blocking, so for a large dense matrix it is an order of magnitude slower than a LAPACK kernel.

LapackLU is the same interface with LinearAlgebra.lu! underneath:

solve(LapackLU(), ls)

3-element Vector{Float64}:
  1.4802958858695092e-10
 -0.3333333336293925
  0.6666666668146962

It is restricted to the element types LAPACK provides (Float32, Float64, ComplexF32 and ComplexF64) and throws an ArgumentError naming the type for anything else, so LU() remains the only dense option for e.g. BigFloat (a sparse one of those goes to SparspakLU). What it is not is a trade of allocation for speed: like LU, it allocates nothing per factorization or solve once the LinearSolver has been built. Everything else is interchangeable — factorize!, LinearAlgebra.ldiv!, solve! and solve behave the same way, and either method can be handed to a nonlinear solver as its linear_solver_method:

F(y, x, params) = y .= x .^ 3 .- 2
x = [1.5]
solve!(x, NonlinearProblem(F, zeros(1)), Newton(); linear_solver_method = LapackLU())

1-element Vector{Float64}:
 1.2599210498948732

Choosing a Method

Five methods, and the choice is made by two things: whether the matrix is sparse, and whether its element type is one LAPACK knows. SimpleSolvers.default_linear_solver_method encodes the answer, and it is what a nonlinear solver uses when no linear_solver_method is given:

RecursiveLU is never chosen automatically; see below.

Dense

Measured on an Apple M4 Max against OpenBLAS, factorize! in microseconds including the copy-in:

LU's MMatrix path stops at SimpleSolvers.N_STATIC_THRESHOLD = 10; above that it is a scalar triple loop, and the triangular solve is a further 3.5–4.5× behind getrs throughout. That is why LapackLU rather than LU is the default for the element types it covers.

RecursiveLU wins in the middle — but only against OpenBLAS. With AppleAccelerate loaded on the same machine, LapackLU factorizes a 384 × 384 in 285 µs rather than 531 and a 128 × 128 in 26.5 µs rather than 59.6, which moves the crossover down from n ≈ 200 to n ≈ 64. It also needs a package extension and a heavy dependency, and covers only Float32/Float64. Hence: opt in explicitly, after measuring on the machine that matters.

Sparse

A sparse method needs the sparsity pattern up front, so it is fixed when the LinearSolver is built — which is exactly what makes the ordering and symbolic factorization reusable across refactorizations, and where the saving comes from. A dense matrix is refused rather than converted.

Periodic banded matrices of bandwidth 2, same machine, against a dense LapackLU on the same matrix:

Two things to read off. Sparse and dense are a wash around n = 64 and sparse wins by ~7× at n = 384, so sparsity is worth exploiting only once the matrix is big enough. And SparspakLU has the faster factorization but a ~9× slower solve, which reverses the comparison in a nonlinear solve — where one factorization is followed by one or more solves. So UmfpackLU is the default.

What SparspakLU is for is element types UMFPACK cannot do at all:

So every element type outside Float64/ComplexF64 has no default at all, and that is deliberate: a sparse matrix is never densified for you. SimpleSolvers.default_linear_solver_method raises an ArgumentError naming the two things you might have meant — SparspakLU, which keeps the matrix sparse, or a dense method (LapackLU for a 32-bit float, LU otherwise), which discards the sparsity. Both are legitimate; which one is right depends on how large the matrix is and whether you can depend on the Sparspak extension, and neither is something a fallback should decide. Pass one as linear_solver_method.

An exact Rational solve goes through factorize! and LinearAlgebra.ldiv! rather than the allocating solve, whose NaN-filled solution vector those element types cannot represent.

Neither sparse method is allocation-free, and that is inside the backends rather than in the wrapper: UmfpackLU allocates ~374 kB per factorization but nothing per solve; SparspakLU allocates ~11 kB and ~10 kB respectively.

Sparse Jacobians in a nonlinear solve

To run a sparse Jacobian through a NewtonSolver or DogLegSolver, pass the pattern as jacobian_prototype together with a DF! that assembles into it:

solver = NewtonSolver(x, y; F = F!, DF! = DF!, jacobian_prototype = J0)

The prototype's storage is adopted by the Jacobian, the LinearProblem and the LinearSolver's cache, and SimpleSolvers.default_linear_solver_method then selects UmfpackLU. DF! is required: JacobianAutodiff and JacobianFiniteDifferences produce dense matrices and would write to structurally-zero positions, so that combination is refused at construction.

There is no linear_solver_method to pass for a Float64/ComplexF64 pattern — that is the one case with a default — but every other element type needs one, since a sparse matrix is never densified on your behalf.

The pattern must not change from iteration to iteration — the symbolic factorization was built for one — and DF! writing into positions outside it is an error rather than a silent reallocation.