Linear Solvers

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

\[ Ax = y\]

is satisfied.

A linear system can be called with^[1]:

using SimpleSolvers

A = [(0. + 1e-6) 1. 2.; 3. 4. 5.; 6. 7. 8.]
y = [1., 2., 3.]
ls = LinearSystem(A, y)
update!(ls, 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 LinearSystem with an lu solver (using LU and solve) in double precision and without pivoting:

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

3-element StaticArraysCore.SizedVector{3, Float64, Vector{Float64}} with indices SOneTo(3):
  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ˢ = LinearSystem(Aˢ, yˢ)
update!(lsˢ, Aˢ, yˢ)
y² = solve(lu, lsˢ)

3-element StaticArraysCore.SizedVector{3, Float32, Vector{Float32}} with indices SOneTo(3):
 -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 StaticArraysCore.SizedVector{3, Float32, Vector{Float32}} with indices SOneTo(3):
  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.