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 LAPACK's getrf 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

Solving a Rank-Deficient System

Every method above answers a singular matrix with a SingularException, which is what a caller wants when a singular matrix means something has gone wrong. It is not what a caller wants when the deficiency is a property of the problem — a residual with a continuous symmetry, an over-parametrised ansatz, a constraint that is implied by the others. There the system is still consistent: the right-hand side lies in the range of the matrix, so solutions exist, and there is a whole affine family of them.

PivotedQR and SVDSolver return the minimum-norm member of that family instead of refusing:

Ad = [1. 2.; 2. 4.]       # rank 1
bd = [1., 2.]             # consistent: bd = Ad * [0.2, 0.4]
solve(PivotedQR(), Ad, bd)

2-element Vector{Float64}:
 0.20000000000000004
 0.40000000000000013

An LU cannot do anything with this matrix at all:

try
    solve(LapackLU(), Ad, bd)
catch e
    e
end

LinearAlgebra.SingularException(2)

Every rank-revealing method determines a numerical rank as it factorizes, at a relative tolerance that is theirs to carry — SimpleSolvers.rank_tolerance documents the default and why it is looser than LinearAlgebra.rank's. LinearAlgebra.rank reads it back off the factorization, and SVDSolver additionally hands back the spectrum it was read from:

using LinearAlgebra: rank

ls = LinearSolver(SVDSolver(), Ad)
factorize!(ls, Ad)
rank(ls), SimpleSolvers.singular_values(ls)

(1, [5.0, 0.0])

That last pair is the distinction no tolerance can make on its own. A spectrum with a gap is rank deficient and the directions below it are genuinely absent; a spectrum that decays smoothly is ill-conditioned, has no correct rank, and truncating it is a modelling decision rather than a numerical one.

Which of the two to use follows the ratio of solves to factorizations rather than the size of the matrix — a complete orthogonal factorization has the cheaper factorization and the more expensive solve, a singular value decomposition the reverse — and the measured table is in LapackPivotedQR's docstring.

Each comes in two kernels, exactly as LU and LapackLU do:

The Lapack pair is faster and restricted to Float32, Float64, ComplexF32 and ComplexF64. The other two are written in plain Julia and take any floating-point element type, which is what makes a rank-deficient system in Float16 or BigFloat solvable at all:

A16 = Float16[1 2; 2 4]          # rank 1, and LAPACK has no Float16
solve(SVDSolver(), A16, Float16[1, 2])

2-element Vector{Float16}:
 0.2001
 0.4001

They are also the allocation-free pair. LAPACK's Julia wrappers allocate their own factors and workspace on every call — ormrz alone asks for 98496 bytes per solve — where the pure-Julia caches are written once at construction and only written into thereafter.

At Float16 the rank tolerance deserves a second look before it is trusted: the default sqrt(eps(Float16)) is 0.031, and SimpleSolvers.rank_tolerance sets out how much room there is around it. All four are restricted to square matrices.

Choosing a Method

Nine methods. For the five that require a non-singular matrix 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. Nor is any of the four rank-revealing methods — see the warning above.

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.