06_Solving_Linear_Systems.qmd

---
title: "Solving Linear Systems"
---

For solving systems with smaller, dense matrices, you could use Julia's LinearAlgebra standard library. You can also use the SparseArrays standard library for larger sparse systems. In addition, there are several iterative methods available, such as `IterativeSolver.jl` for solving very large, sparse systems. Rather than learn how several packages work, however, you can use the SciML [`LinearSolve.jl`](https://docs.sciml.ai/LinearSolve/stable/) package that provides a uniform front-end for all the main linear system packages. You can change the solver algorithm by changing a single line of code, which is very convenient when trying to find the most optimal approach for your system.

## Built-in Methods

We shall start by looking at the methods provided by Julia directly. While it may be more convenient to use `LinearSolve.jl` for larger problems, using the Julia base methods is still very useful for quick calculations.

#### Left Division

For dividing scalar values, we use the `/` (*right division*) operator. There is however also a *left division* operator, `\`. You can solve a linear system directly as follows:

$$Ax = b$$ $$x = A \backslash b $$

``` julia
A = [1 0 3;
    -1 1 0;
     2 1 1]
# 3×3 Matrix{Int64}:
#   1  0  3
#  -1  1  0
#   2  1  1

x = [1, 2, 1]
# 3-element Vector{Int64}:
#  1
#  2
#  1

b = A*x
# 3-element Vector{Int64}:
#  4
#  1
#  5

x̂=A\b # Note that the results are not integers
# 3-element Vector{Float64}:
#  1.0
#  2.0
#  1.0

x ≈ x̂ # Use approximate check with floats
# true
```

::: callout-tip
It may help to read `A\b` as "A divides into b"
:::

This is a fairly straightforward way of solving a linear system. For dense systems, it is also quite efficient, if your matrix is well-behaved.

Internally, there is quite a bit going on. Julia analyses the matrix, `A`. It then selects an appropriate factorisation methods and uses that to solve the linear system.

Upper and lower triangular and diagonal systems are solved directly via forward or backward substitution. For non-triangular, square matrices, LU factorisation is used.

For rectangular matrices, a minimum norm, least squares solution is calculated using pivoted QR factorisation.

Similar approaches are used for sparse matrices, including use of LDL^T^ factorisation for indefinite matrices etc. See the manual section on [factorisation](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.factorize) for more detail.

If you are going to perform more than one calculation with the same matrix, it will be more efficient to store the factorisation result:

``` julia
using LinearAlgebra # to import factorize()

A = [1 0 3;
    -1 1 0;
     2 1 1]
# 3×3 Matrix{Int64}:
#   1  0  3
#  -1  1  0
#   2  1  1

A = factorize(A) # replace A with its factorised form
# LU{Float64, Matrix{Float64}, Vector{Int64}}
# L factor:
# 3×3 Matrix{Float64}:
#   1.0   0.0       0.0
#  -0.5   1.0       0.0
#   0.5  -0.333333  1.0
# U factor:
# 3×3 Matrix{Float64}:
#  2.0  1.0  1.0
#  0.0  1.5  0.5
#  0.0  0.0  2.66667

A.U # To access the U part
# 3×3 Matrix{Float64}:
#  2.0  1.0  1.0
#  0.0  1.5  0.5
#  0.0  0.0  2.66667

A.L # To access the L part
# 3×3 Matrix{Float64}:
#   1.0   0.0       0.0
#  -0.5   1.0       0.0
#   0.5  -0.333333  1.0

b = [4, 1, 5] # Our first b vector
# 3-element Vector{Int64}:
#  4
#  1
#  5

A\b # Solve using specialised methods for LU factorised matrices
# 3-element Vector{Float64}:
#  1.0
#  2.0
#  1.0

b = [1, 2, 7] # A new b vector
# 3-element Vector{Int64}:
#  1
#  2
#  7

A\b # Solve without needing to repeat the factorisation
# 3-element Vector{Float64}:
#   1.75
#   3.75
#  -0.25000000000000006
```

::: callout-important
You could of course solve a linear system by multiplying with the inverse of the matrix, and this is fairly straightforward to do, with `inv(A)`. This is however, **almost never a good choice**, especially with large, sparse matrices, where the inverse is generally dense and may not fit into your computer's memory.

Factorisation is always a better choice, and iterative methods often the best choice for large, sparse matrices.
:::

#### Specialised Matrix Types

Julia offers several specialised types for matrices with optimised algorithms. These include:

-   Symmetric
-   Hermitian
-   UpperTriangular
-   LowerTriangular
-   Tridiagonal
-   Bidiagonal
-   Diagonal
-   etc.

See the [manual](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#Special-matrices) for the full list.

#### Sparse Matrices

Models of physical systems very often result in large, sparse matrices. Instead of using a vast amount of memory to store mostly zeros, Julia has a built-in sparse array type. This allows you to save only the non-zero entries and also uses specialised linear algebra methods for sparse matrices. To access this, import the standard library `SparseArrays`.

Sparse matrices in Julia are stored in the [*Compressed Sparse Column*](https://en.wikipedia.org/wiki/Sparse_matrix#Compressed_sparse_column_.28CSC_or_CCS.29) format.

Some useful functions for working with sparse matrices include:

-   `nnz()` returns the number of entries in the array. This could include values that have been set to zero.
    -   `count(!iszero, x)` will check each entry and count only non-zero values
-   `dropzeros()` will drop entries that have been set to zero from the list of stored entries. It returns a new sparse array
    -   `dropzeros!()` is similar to `dropzeros()`, but modifies the array directly.
-   `spzeros(m,n)` is similar to `zeros()` for dense arrays. It creates a `m`x`n` sparse array with all elements set to zero.
-   `sparse(I, J, V)` creates a sparse array from thee vectors. `I` is a vector of the row indices, `J` of column indices and `V` holds the values of the entries.
    -   `sparsevec(I, V)` is the vector equivalent of `sparse()`
    -   `findnz()` is the inverse of `sparse()` and `sparsevec()`
-   `issparse()` checks is an array is sparse
-   `blockdiag()` concatenates matrices into a sparse block diagonal matrix.

See the [manual](https://docs.julialang.org/en/v1/stdlib/SparseArrays/#Sparse-Arrays) for more

## LinearSolve

The `LinearSolve.jl` package will use the built-in linear algebra methods when appropriate, but also has many other, specialised methods available from several packages. These include iterative methods that are more efficient for large, sparse matrices. It can also make use of your machines's GPU is one is available[^1].

[^1]: This generally means you have a CUDA or OpenCL capable GPU with the appropriate drivers installed

To install `LinearSolve.jl`, simply use the package manager.

To illustrate how to use the package, we shall use the same toy examples from before:

``` julia
using LinearSolve

A = Float64[1 0 3; # explicitly define as Float64
    -1 1 0;
     2 1 1]
# 3×3 Matrix{Float64}:
#   1.0  0.0  3.0
#  -1.0  1.0  0.0
#   2.0  1.0  1.0

b = Float64[4, 1, 5] # explicitly define as Float64
# 3-element Vector{Float64}:
#  4.0
#  1.0
#  5.0

prob = LinearProblem(A, b) # Define problem
# LinearProblem. In-place: true
# b: 3-element Vector{Float64}:
#  4.0
#  1.0
#  5.0

sol = solve(prob) # and solve
# u: 3-element Vector{Float64}:
#  1.0
#  2.0
#  1.0
```

Since LinearSolve has to allow for automatic differentiation[^2], which may require passing *dual variables* to the function, it will implicitly use the type of the variables as specified in `A` in `b`. This will result in an error the inputs are arrays of integers and the if the solution contains fractions. To avoid this, we explicitly define the inputs as `Float64`.

[^2]: Specifically *forward mode* automatic differentiation. See the relevant section for more detail.

You can also specify which algorithm to use. See the full list and recommendations in the [documentation](https://docs.sciml.ai/LinearSolve/stable/solvers/solvers/).

``` julia
sol = solve(prob, KrylovJL_GMRES()) #GMRES is an iterative solver, from the Krylov.jl package
# u: 3-element Vector{Float64}:
#  1.0000000000000004
#  2.0000000000000004
#  0.9999999999999998
```

In our example for using built-in methods, we stored the result of the factorisation to re-use this when solving the problem with another `b` vector. We can do the same with LinearSolve by adding one step:

``` julia
using LinearSolve

A = Float64[1 0 3;
    -1 1 0;
     2 1 1]
# 3×3 Matrix{Float64}:
#   1.0  0.0  3.0
#  -1.0  1.0  0.0
#   2.0  1.0  1.0

b1 = Float64[4, 1, 5]
# 3-element Vector{Float64}:
#  4.0
#  1.0
#  5.0

prob = LinearProblem(A, b1)
# LinearProblem. In-place: true
# b: 3-element Vector{Float64}:
#  4.0
#  1.0
#  5.0

linsolve = init(prob) # initialize a cache to store the linear system intermediates
# LinearSolve.LinearCache{Matrix{Float64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, Nothing, LU{Float64, Matrix{Float64}, Vector{Int64}}, SciMLOperators.IdentityOperator, SciMLOperators.IdentityOperator, Float64, true}([1.0 0.0 3.0; -1.0 1.0 0.0; 2.0 1.0 1.0], [4.0, 1.0, 5.0], [0.0, 0.0, 0.0], SciMLBase.NullParameters(), nothing, LU{Float64, Matrix{Float64}, Vector{Int64}}(Matrix{Float64}(undef, 0, 0), Int64[], 0), true, SciMLOperators.IdentityOperator(3), SciMLOperators.IdentityOperator(3), 1.4901161193847656e-8, 1.4901161193847656e-8, 3, false, LinearSolve.OperatorAssumptions{true}())

sol1 = solve(linsolve) # solve the problem, using the cache
# u: 3-element Vector{Float64}:
#  1.0
#  2.0
#  1.0

b2 = Float64[1, 2, 7] # another b vector
# 3-element Vector{Float64}:
#  1.0
#  2.0
#  7.0

linsolve = LinearSolve.set_b(sol1.cache, b2) # update the cache with the new b vector
# LinearSolve.LinearCache{Matrix{Float64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, Nothing, LU{Float64, Matrix{Float64}, Vector{Int64}}, SciMLOperators.IdentityOperator, SciMLOperators.IdentityOperator, Float64, true}([2.0 1.0 1.0; -0.5 1.5 0.5; 0.5 -0.3333333333333333 2.6666666666666665], [1.0, 2.0, 7.0], [1.0, 2.0, 1.0], SciMLBase.NullParameters(), nothing, LU{Float64, Matrix{Float64}, Vector{Int64}}([2.0 1.0 1.0; -0.5 1.5 0.5; 0.5 -0.3333333333333333 2.6666666666666665], [3, 2, 3], 0), false, SciMLOperators.IdentityOperator(3), SciMLOperators.IdentityOperator(3), 1.4901161193847656e-8, 1.4901161193847656e-8, 3, false, LinearSolve.OperatorAssumptions{true}())

sol2 = solve(linsolve) # and solve again
# u: 3-element Vector{Float64}:
#   1.75
#   3.75
#  -0.25000000000000006
```