Skip to article frontmatterSkip to article content

Chapter 3

Functions

Solution of least squares by the normal equations
lsnormal.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
"""
    lsnormal(A, b)

Solve a linear least-squares problem by the normal equations.
Returns the minimizer of ||b-Ax||.
"""
function lsnormal(A, b)
    N = A' * A
    z = A' * b
    R = cholesky(N).U
    w = forwardsub(R', z)                   # solve R'z=c
    x = backsub(R, w)                       # solve Rx=z
    return x
end
Solution of least squares by QR factorization
lsqrfact.jl
1
2
3
4
5
6
7
8
9
10
11
12
"""
    lsqrfact(A, b)

Solve a linear least-squares problem by QR factorization. Returns
the minimizer of ||b-Ax||.
"""
function lsqrfact(A, b)
    Q, R = qr(A)
    c = Q' * b
    x = backsub(R, c)
    return x
end
QR factorization by Householder reflections
qrfact.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
"""
    qrfact(A)

QR factorization by Householder reflections. Returns Q and R.
"""
function qrfact(A)
    m, n = size(A)
    Qt = diagm(ones(m))
    R = float(copy(A))
    for k in 1:n
        z = R[k:m, k]
        w = [-sign(z[1]) * norm(z) - z[1]; -z[2:end]]
        nrmw = norm(w)
        if nrmw < eps()
            continue    # already in place; skip this iteration
        end
        v = w / nrmw
        # Apply the reflection to each relevant column of R and Q
        for j in k:n
            R[k:m, j] -= v * (2 * (v' * R[k:m, j]))
        end
        for j in 1:m
            Qt[k:m, j] -= v * (2 * (v' * Qt[k:m, j]))
        end
    end
    return Qt', triu(R)
end

Examples

3.1 Fitting functions to data

Example 3.1.1

Here are 5-year averages of the worldwide temperature anomaly as compared to the 1951–1980 average (source: NASA).

year = 1955:5:2000
temp = [ -0.0480, -0.0180, -0.0360, -0.0120, -0.0040,
       0.1180, 0.2100, 0.3320, 0.3340, 0.4560 ]
    
scatter(year, temp, label="data",
    xlabel="year", ylabel="anomaly (degrees C)", 
    legend=:bottomright)
Loading...

A polynomial interpolant can be used to fit the data. Here we build one using a Vandermonde matrix. First, though, we express time as decades since 1950, as it improves the condition number of the matrix.

t = @. (year - 1950) / 10
n = length(t)
V = [ t[i]^j for i in 1:n, j in 0:n-1 ]
c = V \ temp
10-element Vector{Float64}: -14.114000001832462 76.36173810552113 -165.45597224550528 191.96056669514388 -133.27347224319684 58.015577787494486 -15.962888891734785 2.6948063497166928 -0.2546666667177082 0.010311111113288083

The coefficients in vector c are used to create a polynomial. Then we create a function that evaluates the polynomial after changing the time variable as we did for the Vandermonde matrix.

using Polynomials, Plots
p = Polynomial(c)
f = yr -> p((yr - 1950) / 10)
plot!(f, 1955, 2000, label="interpolant")
Loading...

As you can see, the interpolant does represent the data, in a sense. However it’s a crazy-looking curve for the application. Trying too hard to reproduce all the data exactly is known as overfitting.

Example 3.1.2

Here are the 5-year temperature averages again.

year = 1955:5:2000
t = @. (year - 1950) / 10
temp = [ -0.0480, -0.0180, -0.0360, -0.0120, -0.0040,
          0.1180, 0.2100, 0.3320, 0.3340, 0.4560 ]
10-element Vector{Float64}: -0.048 -0.018 -0.036 -0.012 -0.004 0.118 0.21 0.332 0.334 0.456

The standard best-fit line results from using a linear polynomial that meets the least-squares criterion.

V = [ t.^0 t ]    # Vandermonde-ish matrix
@show size(V)
c = V \ temp
p = Polynomial(c)
Loading...
f = yr -> p((yr - 1955) / 10)
scatter(year, temp, label="data",
    xlabel="year", ylabel="anomaly (degrees C)", leg=:bottomright)
plot!(f, 1955, 2000, label="linear fit")
Loading...

If we use a global cubic polynomial, the points are fit more closely.

V = [ t[i]^j for i in 1:length(t), j in 0:3 ]   
@show size(V);
size(V) = (10, 4)

Now we solve the new least-squares problem to redefine the fitting polynomial.

p = Polynomial( V \ temp )
plot!(f, 1955, 2000, label="cubic fit")
Loading...

If we were to continue increasing the degree of the polynomial, the residual at the data points would get smaller, but overfitting would increase.

Example 3.1.3
a = [1/k^2 for k=1:100] 
s = cumsum(a)        # cumulative summation
p = @. sqrt(6*s)

scatter(1:100, p;
    title="Sequence convergence",
    xlabel=L"k",  ylabel=L"p_k")
Loading...

This graph suggests that maybe pkπp_k\to \pi, but it’s far from clear how close the sequence gets. It’s more informative to plot the sequence of errors, ϵk=πpk\epsilon_k= |\pi-p_k|. By plotting the error sequence on a log-log scale, we can see a nearly linear relationship.

ϵ = @. abs(π - p)    # error sequence
scatter(1:100, ϵ;
    title="Convergence of errors",
    xaxis=(:log10,L"k"),  yaxis=(:log10,"error"))
Loading...

The straight line on the log-log scale suggests a power-law relationship where ϵkakb\epsilon_k\approx a k^b, or logϵkb(logk)+loga\log \epsilon_k \approx b (\log k) + \log a.

k = 1:100
V = [ k.^0 log.(k) ]     # fitting matrix
c = V \ log.(ϵ)          # coefficients of linear fit
2-element Vector{Float64}: -0.18237524972829994 -0.9674103233127929

In terms of the parameters aa and bb used above, we have

a, b = exp(c[1]), c[2];
@show b;
b = -0.9674103233127929

It’s tempting to conjecture that the slope b1b\to -1 asymptotically. Here is how the numerical fit compares to the original convergence curve.

plot!(k, a * k.^b, l=:dash, label="power-law fit")
Loading...

3.2 The normal equations

Example 3.2.1

Because the functions sin2(t)\sin^2(t), cos2(t)\cos^2(t), and 1 are linearly dependent, we should find that the following matrix is somewhat ill-conditioned.

t = range(0, 3, 400)
f = [ x -> sin(x)^2, x -> cos((1 + 1e-7) * x)^2, x -> 1. ]
A = [ f(t) for t in t, f in f ]
@show κ = cond(A);
κ = cond(A) = 1.8253225426741675e7

Now we set up an artificial linear least-squares problem with a known exact solution that actually makes the residual zero.

x = [1., 2, 1]
b = A * x;

Using backslash to find the least-squares solution, we get a relative error that is well below κ times machine epsilon.

x_BS = A \ b
@show observed_error = norm(x_BS - x) / norm(x);
@show error_bound = κ * eps();
observed_error = norm(x_BS - x) / norm(x) = 1.0163949045357309e-10
error_bound = κ * eps() = 4.053030228488391e-9

If we formulate and solve via the normal equations, we get a much larger relative error. With κ21014\kappa^2\approx 10^{14}, we may not be left with more than about 2 accurate digits.

N = A' * A
x_NE = N \ (A'*b)
@show observed_err = norm(x_NE - x) / norm(x);
@show digits = -log10(observed_err);
observed_err = norm(x_NE - x) / norm(x) = 0.021745909192780664
digits = -(log10(observed_err)) = 1.6626224298403076

3.3 The QR factorization

Example 3.3.1

Julia provides access to both the thin and full forms of the QR factorization.

A = rand(1.:9., 6, 4)
@show m,n = size(A);
(m, n) = size(A) = (6, 4)

Here is a standard call:

Q,R = qr(A);
Q
6×6 LinearAlgebra.QRCompactWYQ{Float64, Matrix{Float64}, Matrix{Float64}}
R
4×4 Matrix{Float64}: -12.1244 -13.7739 -11.0521 -8.74273 0.0 -8.44268 -8.5007 -4.09564 0.0 0.0 7.58871 2.57729 0.0 0.0 0.0 5.66991

If you look carefully, you see that we seemingly got a full Q\mathbf{Q} but a thin R\mathbf{R}. However, the Q\mathbf{Q} above is not a standard matrix type. If you convert it to a true matrix, then it reverts to the thin form.

Q̂ = Matrix(Q)
6×4 Matrix{Float64}: -0.742307 0.381927 0.00560904 0.363319 -0.329914 -0.0539855 0.645015 -0.664534 -0.164957 -0.559999 0.318433 0.254603 -0.494872 -0.140201 -0.614228 -0.408772 -0.247436 -0.188546 0.21908 0.440904 -0.0824786 -0.694559 -0.239277 0.00898257

We can test that Q\mathbf{Q} is an orthogonal matrix:

opnorm(Q' * Q - I)
6.690337603028481e-16

The thin Q^\hat{\mathbf{Q}} cannot be an orthogonal matrix, because it is not square, but it is still ONC:

Q̂' * Q̂ - I
4×4 Matrix{Float64}: -2.22045e-16 5.97028e-17 1.35881e-17 1.00192e-16 5.97028e-17 -4.44089e-16 2.94959e-16 -1.17144e-16 1.35881e-17 2.94959e-16 -1.11022e-16 -4.092e-17 1.00192e-16 -1.17144e-16 -4.092e-17 6.66134e-16
Example 3.3.2

We’ll repeat the experiment of Demo 3.2.1, which exposed instability in the normal equations.

t = range(0, 3, 400)
f = [ x -> sin(x)^2, x -> cos((1 + 1e-7) * x)^2, x -> 1. ]
A = [ f(t) for t in t, f in f ]
x = [1., 2, 1]
b = A * x;

The error in the solution by Function 3.3.2 is similar to the bound predicted by the condition number.

observed_error = norm(FNC.lsqrfact(A, b) - x) / norm(x);
@show observed_error;
@show error_bound = cond(A) * eps();
observed_error = 4.665273501889628e-9

error_bound = cond(A) * eps() = 4.053030228488391e-9

3.4 Computing QR factorizations

Example 3.4.1

We will use Householder reflections to produce a QR factorization of a random matrix.

A = rand(float(1:9), 6, 4)
m,n = size(A)
(6, 4)

Our first step is to introduce zeros below the diagonal in column 1 by using (3.4.4) and (3.4.1).

z = A[:, 1];
v = normalize(z - norm(z) * [1; zeros(m-1)])
P₁ = I - 2v * v'   # reflector
6×6 Matrix{Float64}: 0.181818 0.545455 0.181818 0.545455 0.363636 0.454545 0.545455 0.636364 -0.121212 -0.363636 -0.242424 -0.30303 0.181818 -0.121212 0.959596 -0.121212 -0.0808081 -0.10101 0.545455 -0.363636 -0.121212 0.636364 -0.242424 -0.30303 0.363636 -0.242424 -0.0808081 -0.242424 0.838384 -0.20202 0.454545 -0.30303 -0.10101 -0.30303 -0.20202 0.747475

We check that this reflector introduces zeros as it should:

P₁ * z
6-element Vector{Float64}: 11.000000000000002 1.1102230246251565e-16 -5.412337245047638e-16 -3.3306690738754696e-16 -8.604228440844963e-16 -1.3322676295501878e-15

Now we replace A\mathbf{A} by PA\mathbf{P}\mathbf{A}.

A = P₁ * A
6×4 Matrix{Float64}: 11.0 14.9091 12.2727 11.0 1.11022e-16 1.39394 -0.515152 0.333333 -5.41234e-16 1.46465 3.16162 -0.555556 -3.33067e-16 0.393939 0.484848 3.33333 -8.60423e-16 5.92929 0.323232 3.88889 -1.33227e-15 4.16162 2.40404 -2.88889

We are set to put zeros into column 2. We must not use row 1 in any way, lest it destroy the zeros we just introduced. So we leave it out of the next reflector.

z = A[2:m, 2]
v = normalize(z - norm(z) * [1; zeros(m-2)])
P₂ = I - 2v * v'
5×5 Matrix{Float64}: 0.185089 0.194477 0.0523076 0.787297 0.552583 0.194477 0.953588 -0.0124831 -0.187887 -0.131873 0.0523076 -0.0124831 0.996642 -0.0505351 -0.0354693 0.787297 -0.187887 -0.0505351 0.239382 -0.533858 0.552583 -0.131873 -0.0354693 -0.533858 0.625299

We now apply this reflector to rows 2 and below only.

A[2:m, :] = P₂ * A[2:m, :]
A
6×4 Matrix{Float64}: 11.0 14.9091 12.2727 11.0 -1.51573e-15 7.5312 2.12779 1.59337 -1.53013e-16 2.46165e-16 2.53088 -0.856261 -2.28649e-16 -2.54927e-17 0.315203 3.25245 7.11202e-16 3.20968e-16 -2.23015 2.67155 -2.29185e-16 4.41626e-16 0.611891 -3.74331

We need to iterate the process for the last two columns.

for j in 3:n
    z = A[j:m, j]
    v = normalize(z - norm(z) * [1; zeros(m-j)])
    P = I - 2v * v'
    A[j:m, :] = P * A[j:m, :]
end

We have now reduced the original to an upper triangular matrix using four orthogonal Householder reflections:

R = triu(A)
6×4 Matrix{Float64}: 11.0 14.9091 12.2727 11.0 0.0 7.5312 2.12779 1.59337 0.0 0.0 3.44277 -2.72756 0.0 0.0 0.0 5.00216 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0