Chapter 2

We create two vectors for data about the population of China. The first has the years of census data and the other has the population, in millions of people.

year = [1982, 2000, 2010, 2015]; 
pop = [1008.18, 1262.64, 1337.82, 1374.62];

It’s convenient to measure time in years since 1980. We use .- to subtract a scalar from every element of a vector. We will also use a floating-point value in the subtraction, so the result is also in double precision.

t = year .- 1980.0
y = pop;

Now we have four data points $(t_1,y_1),\dots,(t_4,y_4)$, so $n=4$ and we seek an interpolating cubic polynomial. We construct the associated Vandermonde matrix:

V = [ t[i]^j for i in 1:4, j in 0:3 ]

To solve for the vector of polynomial coefficients, we use a backslash to solve the linear system:

c = V \ y

The algorithms used by the backslash operator are the main topic of this chapter. As a check on the solution, we can compute the residual.

y - V * c

Using floating-point arithmetic, it is not realistic to expect exact equality of quantities; a relative difference comparable to $\macheps$ is all we can look for.

By our definitions, the elements of c are coefficients in ascending-degree order for the interpolating polynomial. We can use the polynomial to estimate the population of China in 2005:

using Polynomials
p = Polynomial(c)    # construct a polynomial
p(2005-1980)         # include the 1980 time shift

The official population value for 2005 was 1303.72, so our result is rather good.

We can visualize the interpolation process. First, we plot the data as points.

using Plots
scatter(t, y;
    label="actual",  legend=:topleft,
    xlabel="years since 1980",  ylabel="population (millions)", 
    title="Population of China")

# Choose 500 times in the interval [0,35].
tt = range(0, 35, 500)

# Evaluate the polynomial at all the vector components.
yy = p.(tt)

foreach(println, yy[1:4])

962.2387878787875

963.9282039963299
965.6118068288089
967.2896105457506

Now we use plot! to add to the current plot, rather than replacing it.

plot!(tt, yy, label="interpolant")

In Julia, vectors and matrices are one-dimensional and two-dimensional arrays, respectively. Square brackets are used to enclose elements of a matrix or vector. Use spaces for horizontal concatenation, and semicolons or new lines to indicate vertical concatenation.

A = [ 1 2 3 4 5; 50 40 30 20 10
    π sqrt(2) exp(1) (1+sqrt(5))/2 log(3) ]

m, n = size(A)

A vector is not quite the same thing as a matrix: it has only one dimension, not two. Separate its elements by commas or semicolons:

x = [ 3, 3, 0, 1, 0 ]
size(x)

For some purposes, however, an $n$-vector in Julia is treated like having a column shape. Note the difference if we use spaces instead of commas inside the brackets:

y = [ 3 3 0 1 0 ]
size(y)

This $1\times 5$ matrix is not equivalent to a vector.

Concatenated elements within brackets may be matrices or vectors for a block representation, as long as all the block sizes are compatible.

[ x  x ]

[ x; x ]

The zeros and ones functions construct matrices with entries all zero or one, respectively.

B = [ zeros(3, 2) ones(3, 1) ]

A single quote ' after a matrix returns its adjoint. For real matrices, this is the transpose; for complex-valued matrices, the elements are also conjugated.

A'

If x is simply a vector, then its transpose has a row shape.

x'

There are many convenient shorthand ways of building vectors and matrices other than entering all of their entries directly or in a loop. To get a range with evenly spaced entries between two endpoints, you have two options. One is to use a colon :.

y = 1:4              # start:stop

z = 0:3:12           # start:step:stop

(Ranges are not strictly considered vectors, but they behave identically in most circumstances.) Instead of specifying the step size, you can give the number of points in the range if you use range.

s = range(-1, 1, 5)

Accessing an element is done by giving one (for a vector) or two (for a matrix) index values within square brackets.

a = A[2, end-1]

x[2]

The indices can be vectors or ranges, in which case a block of the matrix is accessed.

A[1:2, end-2:end]    # first two rows, last three columns

If a dimension has only the index : (a colon), then it refers to all the entries in that dimension of the matrix.

A[:, 1:2:end]        # all of the odd columns

The matrix and vector senses of addition, subtraction, scalar multiplication, multiplication, and power are all handled by the usual symbols.

B = diagm( [-1, 0, -5] )   # create a diagonal matrix

@show size(A), size(B);
BA = B * A     # matrix product

(size(A), size(B)) = ((3, 5), (3, 3))

A * B causes an error here, because the dimensions aren’t compatible.

A * B    # throws an error

DimensionMismatch: matrix A has axes (Base.OneTo(3),Base.OneTo(5)), matrix B has axes (Base.OneTo(3),Base.OneTo(3))

Stacktrace:
 [1] _generic_matmatmul!(C::Matrix{Float64}, A::Matrix{Float64}, B::Matrix{Int64}, _add::LinearAlgebra.MulAddMul{true, true, Bool, Bool})
   @ LinearAlgebra ~/.julia/juliaup/julia-1.11.2+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/matmul.jl:883
 [2] generic_matmatmul!
   @ ~/.julia/juliaup/julia-1.11.2+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/matmul.jl:868 [inlined]
 [3] _mul!
   @ ~/.julia/juliaup/julia-1.11.2+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/matmul.jl:287 [inlined]
 [4] mul!
   @ ~/.julia/juliaup/julia-1.11.2+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/matmul.jl:285 [inlined]
 [5] mul!
   @ ~/.julia/juliaup/julia-1.11.2+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/matmul.jl:253 [inlined]
 [6] *(A::Matrix{Float64}, B::Matrix{Int64})
   @ LinearAlgebra ~/.julia/juliaup/julia-1.11.2+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/matmul.jl:114
 [7] top-level scope
   @ In[29]:1

A square matrix raised to an integer power is the same as repeated matrix multiplication.

B^3    # same as B*B*B

Sometimes one instead wants to treat a matrix or vector as a mere array and simply apply a single operation to each element of it. For multiplication, division, and power, the corresponding operators start with a dot.

C = -A;

Because both matrices are $3\times 5$, A * C would be an error here, but elementwise operations are fine.

elementwise = A .* C

The two operands of a dot operator have to have the same size—unless one is a scalar, in which case it is expanded or broadcast to be the same size as the other operand.

x_to_two = x .^ 2

two_to_x = 2 .^ x

Most of the mathematical functions, such as cos, sin, log, exp, and sqrt, expect scalars as operands. However, you can broadcast any function, including ones that you have defined, across a vector or array by using a special dot syntax.

show(cos.(π * x))    # broadcast to a function

[-1.0, -1.0, 1.0, -1.0, 1.0]

Rather than dotting multiple individual functions, you can use @. before an expression to broadcast everything within it.

show(@. cospi( (x + 1)^3) )    # broadcast an entire expression

[1.0, 1.0, -1.0, 1.0, -1.0]

For a square matrix $\mathbf{A}$, the syntax A \ b is mathematically equivalent to $\mathbf{A}^{-1} \mathbf{b}$.

A = [1 0 -1; 2 2 1; -1 -3 0]

b = [1, 2, 3]

x = A \ b

One way to check the answer is to compute a quantity known as the residual. It is (ideally) close to machine precision (relative to the elements in the data).

residual = b - A * x

If the matrix $\mathbf{A}$ is singular, you may get an error.

A = [0 1; 0 0]
b = [1, -1]
x = A \ b    # throws an error

SingularException(2)

Stacktrace:
 [1] generic_trimatdiv!(C::Vector{Float64}, uploc::Char, isunitc::Char, tfun::typeof(identity), A::Matrix{Int64}, B::Vector{Int64})
   @ LinearAlgebra ~/.julia/juliaup/julia-1.11.2+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/triangular.jl:1388
 [2] _ldiv!
   @ ~/.julia/juliaup/julia-1.11.2+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/triangular.jl:966 [inlined]
 [3] ldiv!
   @ ~/.julia/juliaup/julia-1.11.2+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/triangular.jl:959 [inlined]
 [4] \
   @ ~/.julia/juliaup/julia-1.11.2+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/triangular.jl:1721 [inlined]
 [5] \(A::Matrix{Int64}, B::Vector{Int64})
   @ LinearAlgebra ~/.julia/juliaup/julia-1.11.2+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/generic.jl:1130
 [6] top-level scope
   @ In[41]:3

In this case we can check that the rank of $\mathbf{A}$ is less than its number of columns, indicating singularity.

rank(A)

A linear system with a singular matrix might have no solution or infinitely many solutions, but in either case, backslash will fail. Moreover, detecting singularity is a lot like checking whether two floating-point numbers are exactly equal: because of roundoff, it could be missed. In Conditioning of linear systems we’ll find a robust way to fully describe this situation.

It’s easy to get just the lower triangular part of any matrix using the tril function.

A = rand(1.:9., 5, 5)
L = tril(A)

We’ll set up and solve a linear system with this matrix.

b = ones(5)
x = FNC.forwardsub(L,b)

It’s not clear how accurate this answer is. However, the residual should be zero or comparable to $\macheps$.

b - L * x

Next we’ll engineer a problem to which we know the exact answer. Use \alpha Tab and \beta Tab to get the Greek letters.

α = 0.3;
β = 2.2;
U = diagm( 0=>ones(5), 1=>[-1, -1, -1, -1] )
U[1, [4, 5]] = [ α - β, β ]
U

x_exact = ones(5)
b = [α, 0, 0, 0, 1]

Now we use backward substitution to solve for $\mathbf{x}$, and compare to the exact solution we know already.

x = FNC.backsub(U,b)
err = x - x_exact

Everything seems OK here. But another example, with a different value for β, is more troubling.

α = 0.3;
β = 1e12;
U = diagm( 0=>ones(5), 1=>[-1, -1, -1, -1] )
U[1, [4, 5]] = [ α - β, β ]
b = [α, 0, 0, 0, 1]

x = FNC.backsub(U,b)
err = x - x_exact

It’s not so good to get 4 digits of accuracy after starting with 16! The source of the error is not hard to track down. Solving for $x_1$ performs $(\alpha-\beta)+\beta$ in the first row. Since $|\alpha|$ is so much smaller than $|\beta|$, this a recipe for losing digits to subtractive cancellation.

We explore the outer product formula for two random triangular matrices.

L = tril( rand(1:9, 3, 3) )

U = triu( rand(1:9, 3, 3) )

Here are the three outer products in the sum in (2.4.4):

L[:, 1] * U[1, :]'

L[:, 2] * U[2, :]'

L[:, 3] * U[3, :]'

Simply because of the triangular zero structures, only the first outer product contributes to the first row and first column of the entire product.

For illustration, we work on a $4 \times 4$ matrix. We name it with a subscript in preparation for what comes.

A₁ = [
     2    0    4     3 
    -4    5   -7   -10 
     1   15    2   -4.5
    -2    0    2   -13
    ];
L = diagm(ones(4))
U = zeros(4, 4);

Now we appeal to (2.4.5). Since $L_{11}=1$, we see that the first row of $\mathbf{U}$ is just the first row of $\mathbf{A}_1$.

U[1, :] = A₁[1, :]
U

From (2.4.6), we see that we can find the first column of $\mathbf{L}$ from the first column of $\mathbf{A}_1$.

L[:, 1] = A₁[:, 1] / U[1, 1]
L

We have obtained the first term in the sum (2.4.4) for $\mathbf{L}\mathbf{U}$, and we subtract it away from $\mathbf{A}_1$.

A₂ = A₁ - L[:, 1] * U[1, :]'

Now $\mathbf{A}_2 = \boldsymbol{\ell}_2\mathbf{u}_2^T + \boldsymbol{\ell}_3\mathbf{u}_3^T + \boldsymbol{\ell}_4\mathbf{u}_4^T.$ If we ignore the first row and first column of the matrices in this equation, then in what remains we are in the same situation as at the start. Specifically, only $\boldsymbol{\ell}_2\mathbf{u}_2^T$ has any effect on the second row and column, so we can deduce them now.

U[2, :] = A₂[2, :]
L[:, 2] = A₂[:, 2] / U[2, 2]
L

If we subtract off the latest outer product, we have a matrix that is zero in the first two rows and columns.

A₃ = A₂ - L[:, 2] * U[2, :]'

Now we can deal with the lower right $2\times 2$ submatrix of the remainder in a similar fashion.

U[3, :] = A₃[3, :]
L[:, 3] = A₃[:, 3] / U[3, 3]
A₄ = A₃ - L[:, 3] * U[3, :]'

Finally, we pick up the last unknown in the factors.

U[4, 4] = A₄[4, 4];

We now have all of $\mathbf{L}$,

L

and all of $\mathbf{U}$,

U

We can verify that we have a correct factorization of the original matrix by computing the backward error:

A₁ - L * U

IIn floating point, we cannot expect the difference to be exactly zero as we found in this toy example. Instead, we would be satisfied to see that each element of the difference above is comparable in size to machine precision.

Here are the data for a linear system $\mathbf{A}\mathbf{x}=\mathbf{b}$.

A = [2 0 4 3; -4 5 -7 -10; 1 15 2 -4.5; -2 0 2 -13];
b = [4,9,9,4];

We apply Function 2.4.1 and then do two triangular solves.

L, U = FNC.lufact(A)
z = FNC.forwardsub(L, b)
x = FNC.backsub(U, z)

A check on the residual assures us that we found the solution.

b - A*x

Here is a straightforward implementation of matrix-vector multiplication.

n = 6
A = randn(n, n)
x = rand(n)
y = zeros(n)
for i in 1:n
    for j in 1:n
        y[i] += A[i, j] * x[j]    # 1 multiply, 1 add
    end
end

Each of the loops implies a summation of flops. The total flop count for this algorithm is

Since the matrix $\mathbf{A}$ has $n^2$ elements, all of which have to be involved in the product, it seems unlikely that we could get a flop count that is smaller than $O(n^2)$ in general.

Let’s run an experiment with the built-in matrix-vector multiplication. Note that Julia is unusual in that loops have a variable scope separate from its enclosing code. Thus, for n in n below means that inside the loop, the name n will take on each one of the values that were previously assigned to the vector n.

n = 1000:1000:5000
t = []
for n in n
    A = randn(n, n)  
    x = randn(n)
    time = @elapsed for j in 1:80; A * x; end
    push!(t, time)
end

The reason for doing multiple repetitions at each value of $n$ in the loop above is to avoid having times so short that the resolution of the timer is significant.

@pt :header = ["size", "time (sec.)"] [n t]

Looking at the timings just for $n=2000$ and $n=4000$, they have ratio

@show t[n.==4000] ./ t[n.==2000];

t[n .== 4000] ./ t[n .== 2000] = [1.4117809950868896]

If the run time is dominated by flops, then we expect this ratio to be

Let’s repeat the experiment of the previous figure for more, and larger, values of $n$.

randn(5,5)*randn(5);  # throwaway to force compilation

n = 400:200:6000
t = []
for n in n
    A = randn(n, n)  
    x = randn(n)
    time = @elapsed for j in 1:50; A * x; end
    push!(t, time)
end

Plotting the time as a function of $n$ on log-log scales is equivalent to plotting the logs of the variables.

scatter(n, t, label="data", legend=false,
    xaxis=(:log10, L"n"), yaxis=(:log10, "elapsed time (sec)"),
    title="Timing of matrix-vector multiplications")

You can see that while the full story is complicated, the graph is trending to a straight line of positive slope. For comparison, we can plot a line that represents $O(n^2)$ growth exactly. (All such lines have slope equal to 2.)

plot!(n, t[end] * (n/n[end]).^2, l=:dash,
    label=L"O(n^2)", legend=:topleft)

We’ll test the conclusion of $O(n^3)$ flops experimentally, using the built-in lu function instead of the purely instructive lufact.

lu(randn(3, 3));   # throwaway to force compilation

n = 400:400:4000
t = []
for n in n
    A = randn(n, n)  
    time = @elapsed for j in 1:12; lu(A); end
    push!(t, time)
end

We plot the timings on a log-log graph and compare it to $O(n^3)$. The result could vary significantly from machine to machine, but in theory the data should start to parallel the line as $n\to\infty$.

scatter(n, t, label="data", legend=:topleft,
    xaxis=(:log10, L"n"), yaxis=(:log10, "elapsed time"))
plot!(n, t[end ]* (n/n[end]).^3, l=:dash, label=L"O(n^3)")

Here is a previously encountered matrix that factors well.

A = [2 0 4 3 ; -4 5 -7 -10 ; 1 15 2 -4.5 ; -2 0 2 -13];
L, U = FNC.lufact(A)
L

If we swap the second and fourth rows of $\mathbf{A}$, the result is still nonsingular. However, the factorization now fails.

A[[2, 4], :] = A[[4, 2], :]  
L, U = FNC.lufact(A)
L

The presence of NaN in the result indicates that some impossible operation was required. The source of the problem is easy to locate. We can find the first outer product in the factorization just fine:

U[1, :] = A[1, :]
L[:, 1] = A[:, 1] / U[1, 1]
A -= L[:, 1] * U[1, :]'

The next step is U[2, :] = A[2, :], which is also OK. But then we are supposed to divide by U[2, 2], which is zero. The algorithm cannot continue.

Here is the trouble-making matrix from Demo 2.6.1.

A₁ = [2 0 4 3 ; -2 0 2 -13; 1 15 2 -4.5 ; -4 5 -7 -10]

We now find the largest candidate pivot in the first column. We don’t care about sign, so we take absolute values before finding the max.

i = argmax( abs.(A₁[:, 1]) )

This is the row of the matrix that we extract to put into $\mathbf{U}$. That guarantees that the division used to find $\boldsymbol{\ell}_1$ will be valid.

L, U = zeros(4,4),zeros(4,4)
U[1, :] = A₁[i, :]
L[:, 1] = A₁[:, 1] / U[1, 1]
A₂ = A₁ - L[:, 1] * U[1, :]'

Observe that $\mathbf{A}_2$ has a new zero row and zero column, but the zero row is the fourth rather than the first. However, we forge on by using the largest possible pivot in column 2 for the next outer product.

@show i = argmax( abs.(A₂[:, 2]) ) 
U[2, :] = A₂[i, :]
L[:, 2] = A₂[:, 2] / U[2, 2]
A₃ = A₂ - L[:, 2] * U[2, :]'

i = argmax(abs.(A₂[:, 2])) = 3

Now we have zeroed out the third row as well as the second column. We can finish out the procedure.

@show i = argmax( abs.(A₃[:, 3]) ) 
U[3, :] = A₃[i, :]
L[:, 3] = A₃[:, 3] / U[3, 3]
A₄ = A₃ - L[:, 3] * U[3, :]'

i = argmax(abs.(A₃[:, 3])) = 2

@show i = argmax( abs.(A₄[:, 4]) ) 
U[4, :] = A₄[i, :]
L[:, 4] = A₄[:, 4] / U[4, 4];

i = argmax(abs.(A₄[:, 4])) = 1

We do have a factorization of the original matrix:

A₁ - L * U

And $\mathbf{U}$ has the required structure:

U

However, the triangularity of $\mathbf{L}$ has been broken.

L

Here again is the matrix from Demo 2.6.2.

A = [2 0 4 3 ; -2 0 2 -13; 1 15 2 -4.5 ; -4 5 -7 -10]

As the factorization proceeded, the pivots were selected from rows 4, 3, 2, and finally 1. If we were to put the rows of $\mathbf{A}$ into that order, then the algorithm would run exactly like the plain LU factorization from LU factorization.

B = A[[4, 3, 2, 1], :]
L, U = FNC.lufact(B);

We obtain the same $\mathbf{U}$ as before:

U

And $\mathbf{L}$ has the same rows as before, but arranged into triangular order:

L

The third output of plufact is the permutation vector we need to apply to $\mathbf{A}$.

A = rand(1:20, 4, 4)
L, U, p = FNC.plufact(A)
A[p,:] - L * U   # should be ≈ 0

Given a vector $\mathbf{b}$, we solve $\mathbf{A}\mathbf{x}=\mathbf{b}$ by first permuting the entries of $\mathbf{b}$ and then proceeding as before.

b = rand(4)
z = FNC.forwardsub(L,b[p])
x = FNC.backsub(U,z)

A residual check is successful:

b - A*x

With the syntax A \ b, the matrix A is PLU-factored, followed by two triangular solves.

A = randn(500, 500)   # 500x500 with normal random entries
A \ rand(500)          # force compilation
@elapsed for k=1:50; A \ rand(500); end

In Efficiency of matrix computations we showed that the factorization is by far the most costly part of the solution process. A factorization object allows us to do that costly step only once per unique matrix.

factored = lu(A)     # store factorization result
factored \ rand(500)   # force compilation
@elapsed for k=1:50; factored \ rand(500); end

We construct a linear system for this matrix with $\epsilon=10^{-12}$ and exact solution $[1,1]$:

ϵ = 1e-12
A = [-ϵ 1; 1 -1]
b = A * [1, 1]

We can factor the matrix without pivoting and solve for $\mathbf{x}$.

L, U = FNC.lufact(A)
x = FNC.backsub( U, FNC.forwardsub(L, b) )

Note that we have obtained only about 5 accurate digits for $x_1$. We could make the result even more inaccurate by making ε even smaller:

ϵ = 1e-20; A = [-ϵ 1; 1 -1]
b = A * [1, 1]
L, U = FNC.lufact(A)
x = FNC.backsub( U, FNC.forwardsub(L, b) )

This effect is not due to ill conditioning of the problem—a solution with PLU factorization works perfectly:

A \ b

In Julia the LinearAlgebra package has a norm function for vector norms.

x = [2, -3, 1, -1]
twonorm = norm(x)         # or norm(x,2)

infnorm = norm(x, Inf)

onenorm = norm(x, 1)

There is also a normalize function that divides a vector by its norm, making it a unit vector.

normalize(x, Inf)

A = [ 2 0; 1 -1 ]

In Julia, one uses norm for vector norms and for the Frobenius norm of a matrix, which is like stacking the matrix into a single vector before taking the 2-norm.

Fronorm = norm(A)

Most of the time we want to use opnorm, which is an induced matrix norm. The default is the 2-norm.

twonorm = opnorm(A)

You can get the 1-norm as well.

onenorm = opnorm(A, 1)

According to (2.7.15), the matrix 1-norm is equivalent to the maximum of the sums down the columns (in absolute value).

# Sum down the rows (1st matrix dimension):
maximum( sum(abs.(A), dims=1) )

Similarly, we can get the ∞-norm and check our formula for it.

infnorm = opnorm(A, Inf)

 # Sum across columns (2nd matrix dimension):
maximum( sum(abs.(A), dims=2) )

Next we illustrate a geometric interpretation of the 2-norm. First, we will sample a lot of vectors on the unit circle in $\mathbb{R}^2$.

# Construct 601 unit column vectors.
θ = 2π * (0:1/600:1)   # type \theta then Tab
x = [ fun(t) for fun in [cos, sin], t in θ ];

To create an array of plots, start with a plot that has a layout argument, then do subsequent plot! calls with a subplot argument.

plot(aspect_ratio=1, layout=(1, 2),
    xlabel=L"x_1",  ylabel=L"x_2")
plot!(x[1, :], x[2, :], subplot=1, title="Unit circle")

The linear function $\mathbf{f}(\mathbf{x}) = \mathbf{A}\mathbf{x}$ defines a mapping from $\mathbb{R}^2$ to $\mathbb{R}^2$. We can apply A to every column of x by using a single matrix multiplication.

Ax = A * x;

The image of the transformed vectors is an ellipse.

plot!(Ax[1, :], Ax[2, :], 
    subplot=2, title="Image under x → Ax")

That ellipse just touches the circle of radius $\|\mathbf{A}\|_2$.

plot!(twonorm*x[1, :], twonorm*x[2, :], subplot=2, l=:dash)

Julia has a function cond to compute matrix condition numbers. By default, the 2-norm is used. As an example, the family of Hilbert matrices is famously badly conditioned. Here is the $6\times 6$ case.

A = [ 1 / (i + j) for i in 1:6, j in 1:6 ]
κ = cond(A)

Because $\kappa\approx 10^8$, it’s possible to lose nearly 8 digits of accuracy in the process of passing from $\mathbf{A}$ and $\mathbf{b}$ to $\mathbf{x}$. That fact is independent of the algorithm; it’s inevitable once the data are expressed in finite precision.

Let’s engineer a linear system problem to observe the effect of a perturbation. We will make sure we know the exact answer.

x = 1:6
b = A * x

Now we perturb the system matrix and vector randomly by 10^-10 in norm.

# type \Delta then Tab to get Δ
ΔA = randn(size(A));  ΔA = 1e-10 * (ΔA / opnorm(ΔA));
Δb = randn(size(b));  Δb = 1e-10 * normalize(Δb);

We solve the perturbed problem using pivoted LU and see how the solution was changed.

new_x = ((A + ΔA) \ (b + Δb))
Δx = new_x - x

Here is the relative error in the solution.

@show relative_error = norm(Δx) / norm(x);

relative_error = norm(Δx) / norm(x) = 2.624224063874282e-5

And here are upper bounds predicted using the condition number of the original matrix.

println("Upper bound due to b: $(κ * norm(Δb) / norm(b))")
println("Upper bound due to A: $(κ * opnorm(ΔA) / opnorm(A))")

Upper bound due to b: 0.0006723667714371329
Upper bound due to A: 0.004566989873939067

Even if we didn’t make any manual perturbations to the data, machine roundoff does so at the relative level of $\macheps$.

Δx = A\b - x
@show relative_error = norm(Δx) / norm(x);
@show rounding_bound = κ * eps();

relative_error = norm(Δx) / norm(x) = 7.822650774976615e-10
rounding_bound = κ * eps() = 1.134607141116935e-8

Larger Hilbert matrices are even more poorly conditioned:

A = [ 1 / (i + j) for i=1:14, j=1:14 ];
κ = cond(A)

Note that κ exceeds $1/\macheps$. In principle we therefore may end up with an answer that has relative error greater than 100%.

rounding_bound = κ*eps()

Let’s put that prediction to the test.

x = 1:14
b = A * x  
Δx = A\b - x
@show relative_error = norm(Δx) / norm(x);

relative_error = norm(Δx) / norm(x) = 4.469466154206132

As anticipated, the solution has zero accurate digits in the 2-norm.

Here is a small tridiagonal matrix. Note that there is one fewer element on the super- and subdiagonals than on the main diagonal.

A = diagm( -1 => [4, 3, 2, 1, 0], 
    0 => [2, 2, 0, 2, 1, 2], 
    1 => fill(-1, 5) )

We can extract the elements on any diagonal using the diag function. The main or central diagonal is numbered zero, above and to the right of that is positive, and below and to the left is negative.

@show diag_main = diag(A);
@show diag_minusone = diag(A, -1);

diag_main = diag(A) = [2, 2, 0, 2, 1, 2]

diag_minusone = diag(A, -1) = [4, 3, 2, 1, 0]

The lower and upper bandwidths of $\mathbf{A}$ are repeated in the factors from the unpivoted LU factorization.

L, U = FNC.lufact(A)
L

U

We begin with a symmetric $\mathbf{A}$.

A₁ = [  2     4     4     2
        4     5     8    -5
        4     8     6     2
        2    -5     2   -26 ];

We won’t use pivoting, so the pivot element is at position (1,1). This will become the first element on the diagonal of $\mathbf{D}$. Then we divide by that pivot to get the first column of $\mathbf{L}$.

L = diagm(ones(4))
d = zeros(4)
d[1] = A₁[1, 1]
L[:, 1] = A₁[:, 1] / d[1]
A₂ = A₁ - d[1] * L[:, 1] * L[:, 1]'

We are now set up the same way for the submatrix in rows and columns 2–4.

d[2] = A₂[2, 2]
L[:, 2] = A₂[:, 2] / d[2]
A₃ = A₂ - d[2] * L[:, 2] * L[:, 2]'

We continue working our way down the diagonal.

d[3] = A₃[3, 3]
L[:, 3] = A₃[:, 3] / d[3]
A₄ = A₃ - d[3] * L[:, 3] * L[:, 3]'
d[4] = A₄[4, 4]
@show d;
L

d = [2.0, -3.0, -2.0, 1.0]

We have arrived at the desired factorization, which we can validate:

opnorm(A₁ - (L * diagm(d) * L'))

A randomly chosen matrix is extremely unlikely to be symmetric. However, there is a simple way to symmetrize one.

A = rand(1.0:9.0, 4, 4)
B = A + A'

Similarly, a random symmetric matrix is unlikely to be positive definite. The Cholesky algorithm always detects a non-PD matrix by quitting with an error.

cholesky(B)    # throws an error

PosDefException: matrix is not positive definite; Factorization failed.

Stacktrace:
  [1] checkpositivedefinite
    @ ~/.julia/juliaup/julia-1.11.2+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/factorization.jl:68 [inlined]
  [2] #cholesky!#163
    @ ~/.julia/juliaup/julia-1.11.2+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/cholesky.jl:269 [inlined]
  [3] cholesky!
    @ ~/.julia/juliaup/julia-1.11.2+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/cholesky.jl:267 [inlined]
  [4] #cholesky!#164
    @ ~/.julia/juliaup/julia-1.11.2+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/cholesky.jl:301 [inlined]
  [5] cholesky! (repeats 2 times)
    @ ~/.julia/juliaup/julia-1.11.2+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/cholesky.jl:295 [inlined]
  [6] _cholesky
    @ ~/.julia/juliaup/julia-1.11.2+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/cholesky.jl:411 [inlined]
  [7] cholesky(A::Matrix{Float64}, ::NoPivot; check::Bool)
    @ LinearAlgebra ~/.julia/juliaup/julia-1.11.2+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/cholesky.jl:401
  [8] cholesky
    @ ~/.julia/juliaup/julia-1.11.2+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/cholesky.jl:401 [inlined]
  [9] cholesky(A::Matrix{Float64})
    @ LinearAlgebra ~/.julia/juliaup/julia-1.11.2+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/cholesky.jl:401
 [10] top-level scope
    @ In[139]:1

It’s not hard to manufacture an SPD matrix to try out the Cholesky factorization.

B = A' * A
cf = cholesky(B)

What’s returned is a factorization object. Another step is required to extract the factor as a matrix:

R = cf.U

Here we validate the factorization:

opnorm(R' * R - B) / opnorm(B)