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

t = year .- 1980.0
y = pop;

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

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.

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.

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,length=500)   
# Evaluate the polynomial at all the vector components.
yy = p.(tt)
foreach(println,yy[1:4])

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

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)

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

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

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

A * B    # throws an error

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

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

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

show(@. cos(π*(x+1)^3))    # broadcast an entire expression

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

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.

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

α = 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) )

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.

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.

pretty_table([n t], header=(["size", "time"], ["", "(sec)"]))

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

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)

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]

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, :]'

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, :]'

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

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)

# 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) )

# 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)

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);

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))")

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();

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);

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

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

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

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

If we use an ordinary or dense matrix, then there’s no way to exploit a banded structure such as tridiagonality.

n = 10000
A = diagm(0=>1:n, 1=>n-1:-1:1, -1=>ones(n-1))
lu(rand(3, 3))  # force compilation
@time lu(A);

If instead we construct a proper sparse matrix, though, the speedup can be dramatic.

A = spdiagm(0=>1:n, 1=>n-1:-1:1, -1=>ones(n-1))
lu(A);    # compile for sparse case
@time lu(A);

You can also see above that far less memory was used in the sparse case.

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

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'

cholesky(B)    # throws an error

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)