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 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 \]

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
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 LDLT factorisation for indefinite matrices etc. See the manual section on factorisation 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:

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
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 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 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 mxn 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 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 available1.

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:

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 differentiation2, 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.

You can also specify which algorithm to use. See the full list and recommendations in the documentation.

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:

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
Back to top

Footnotes

  1. This generally means you have a CUDA or OpenCL capable GPU with the appropriate drivers installed↩︎

  2. Specifically forward mode automatic differentiation. See the relevant section for more detail.↩︎