Frontmatter

If you are publishing this notebook on the web, you can set the parameters below to provide HTML metadata. This is useful for search engines and social media.

Notebook for 2023-02-08

In this lecture, we discussed forward and backward substitution for triangular matrices, and the basic Gaussian elimination algorithm without pivoting. We copy in the algorithms here and illustrate them, as well as illustrating what we would really want to do in Julia using the LinearAlgebra package rather than rolling our own.
md"""
# Notebook for 2023-02-08

In this lecture, we discussed forward and backward substitution for triangular matrices, and the basic Gaussian elimination algorithm without pivoting. We copy in the algorithms here and illustrate them, as well as illustrating what we would really want to do in Julia using the `LinearAlgebra` package rather than rolling our own.
"""
3.6 ms
using LinearAlgebra
289 μs

Triangular solves

LU factorization gives us a unit lower triangular matrix and an upper triangular matrix as factors, so let's start by figuring out how to solve systems with them.
md"""
## Triangular solves

LU factorization gives us a *unit* lower triangular matrix and an upper triangular matrix as factors, so let's start by figuring out how to solve systems with them.
"""
1.2 ms
forward_subst_unit (generic function with 1 method)
function forward_subst_unit(L, d)
y = copy(d)
n = length(d)
for j = 1:n
y[j] -= L[j,1:j-1]'*y[1:j-1]
end
y
end
715 μs
test_forward_subst (generic function with 1 method)
function test_forward_subst()
# Set up test problem: L is unit lower triangular (diagonals are all 1)
L = [1.0 0.0 0.0;
2.0 1.0 0.0;
3.0 2.0 1.0]
d = [1.0; 1.0; 3.0]
y1 = forward_subst_unit(L, d) # Solve with our routine
y2 = L\d # Solve with the Julia built-in
y1, y2 # Return the results for side-by-side viewing
end
416 μs
test_forward_subst()
1.8 ms
backward_subst (generic function with 1 method)
function backward_subst(U, d)
x = copy(d)
n = length(d)
for j = n:-1:1
x[j] = (d[j] - U[j,j+1:n]'*x[j+1:n])/U[j,j]
end
x
end
768 μs
test_backward_subst (generic function with 1 method)
function test_backward_subst()
# Set up test problem: U is upper triangular
U = [1.0 4.0 7.0;
0.0 -3.0 -6.0;
0.0 0.0 1.0]
d = [1.0; 1.0; 3.0]
y1 = backward_subst(U, d) # Solve with our routine
y2 = U\d # Solve with the Julia built-in
y1, y2 # Return the results for side-by-side viewing
end
412 μs
test_backward_subst()
43.8 μs

Gaussian elimination

Let's write a version of the LU factorization routine that overwrites $A$ with the factors.
md"""
## Gaussian elimination

Let's write a version of the LU factorization routine that overwrites $A$ with the factors.
"""
175 μs
my_lu! (generic function with 1 method)
function my_lu!(A)
n = size(A)[1]
for j = 1:n
A[j+1:n,j] /= A[j,j] # Multipliers (column of L)
A[j+1:n,j+1:n] -= A[j+1:n,j]*A[j,j+1:n]' # Schur complement update
end
A
end
850 μs

It's convenient to also have a version that does not overwrite the underlying matrix.
md"""
It's convenient to also have a version that does not overwrite the underlying matrix.
"""
132 μs
my_lu (generic function with 1 method)
function my_lu(A)
LU = my_lu!(copy(A))
L = UnitLowerTriangular(LU)
U = UpperTriangular(LU)
L, U
end
352 μs
test_my_lu (generic function with 1 method)
function test_my_lu()
A = [1.0 4.0 7.0;
2.0 5.0 8.0;
3.0 6.0 10.0]
L, U = my_lu(A)
norm(A-L*U)/norm(A) # Return relative residual error (Frobenius norm)
end
520 μs
0.0
test_my_lu()
835 μs

If we want to solve linear systems, we can do forward and backward substitution, either using the triangular views or directly on the packed form of the factorization.
md"""
If we want to solve linear systems, we can do forward and backward substitution, either using the triangular views or directly on the packed form of the factorization.
"""
120 μs
test_lu_solve1 (generic function with 1 method)
function test_lu_solve1()

# Manufacture a linear system
A = [1.0 4.0 7.0;
2.0 5.0 8.0;
3.0 6.0 10.0]
x = rand(3)
b = A*x

# Factor and solve
my_lu!(A)
x2 = backward_subst(A, forward_subst_unit(A, b))

# Return a relative error
norm(x2-x)/norm(x)
end
481 μs
1.2947551941976476e-15
test_lu_solve1()
111 μs
test_lu_solve2 (generic function with 1 method)
function test_lu_solve2()

# Manufacture a linear system
A = [1.0 4.0 7.0;
2.0 5.0 8.0;
3.0 6.0 10.0]
x = rand(3)
b = A*x

# Factor and solve (use built-in triangular solvers)
L, U = my_lu(A)
x2 = U\(L\b)

# Return a relative error
norm(x2-x)/norm(x)
end
614 μs
2.858578470698463e-15
test_lu_solve2()
55.1 μs

Finally, let's look at how to do the whole thing end-to-end with the built-in LU routines in Julia. The factorization we compute in this case includes pivoting for stability, which is why the $L$ and $U$ factors are different from those before.
md"""
Finally, let's look at how to do the whole thing end-to-end with the built-in LU routines in Julia. The factorization we compute in this case includes pivoting for stability, which is why the $L$ and $U$ factors are different from those before.
"""
130 μs
LU{Float64, Matrix{Float64}, Vector{Int64}}
L factor:
3×3 Matrix{Float64}:
 1.0       0.0  0.0
 0.333333  1.0  0.0
 0.666667  0.5  1.0
U factor:
3×3 Matrix{Float64}:
 3.0  6.0  10.0
 0.0  2.0   3.66667
 0.0  0.0  -0.5
begin
A = [1.0 4.0 7.0;
2.0 5.0 8.0;
3.0 6.0 10.0]
x = rand(3)
b = A*x

F = lu(A)
end
299 μs

The result with pivoting should look like $PA = LU$; note that Julia stores the permutation as a reordered index vector rather than as a matrix.
md"""
The result with pivoting should look like $PA = LU$; note that Julia stores the permutation as a reordered index vector rather than as a matrix.
"""
117 μs
3×3 Matrix{Float64}:
 0.0  0.0  0.0
 0.0  0.0  0.0
 0.0  0.0  0.0
A[F.p,:] - F.L*F.U
47.3 μs

We can solve the linear system either by using the factorization object or by using the pieces of the factorization (as we did before). The former approach is usually simpler.
md"""
We can solve the linear system either by using the factorization object or by using the pieces of the factorization (as we did before). The former approach is usually simpler.
"""
111 μs
x-F\b
223 μs
x-F.U\(F.L\b[F.p])
36.5 μs