The normal equations - Fundamentals of Numerical Computation

We now solve the general linear least-squares problem in Definition 3.1.1. That is, given $\mathbf{A}\in\mathbb{R}^{m \times n}$ and $\mathbf{b}\in\mathbb{R}^m$ , with $m>n$ , find the $\mathbf{x}\in\mathbb{R}^n$ that minimizes $\| \mathbf{b} - \mathbf{A}\mathbf{x} \|_2$ .

There is a concise explicit solution. In the following proof we make use of the elementary algebraic fact that for two vectors $\mathbf{u}$ and $\mathbf{v}$ ,

(\mathbf{u}+\mathbf{v})^T(\mathbf{u}+\mathbf{v}) = \mathbf{u}^T\mathbf{u} + \mathbf{u}^T\mathbf{v} + \mathbf{v}^T\mathbf{u} + \mathbf{v}^T\mathbf{v} = \mathbf{u}^T\mathbf{u} + 2\mathbf{v}^T\mathbf{u} + \mathbf{v}^T\mathbf{v}.

(3.2.1)

Proof

Let $\mathbf{y}\in \mathbb{R}^n$ be any vector. Then $\mathbf{A}(\mathbf{x}+\mathbf{y})-\mathbf{b}=\mathbf{A}\mathbf{x}-\mathbf{b}+\mathbf{A}\mathbf{y}$ , and

\begin{split} \| \mathbf{A}(\mathbf{x}+\mathbf{y})-\mathbf{b} \|_2^2 &= [(\mathbf{A}\mathbf{x}-\mathbf{b})+(\mathbf{A}\mathbf{y})]^T[(\mathbf{A}\mathbf{x}-\mathbf{b})+(\mathbf{A}\mathbf{y})]\\ &= (\mathbf{A}\mathbf{x}-\mathbf{b})^T(\mathbf{A}\mathbf{x}-\mathbf{b}) + 2(\mathbf{A}\mathbf{y})^T(\mathbf{A}\mathbf{x}-\mathbf{b}) + (\mathbf{A}\mathbf{y})^T(\mathbf{A}\mathbf{y})\\ &= \| \mathbf{A}\mathbf{x}-\mathbf{b} \|_2^2 + 2\mathbf{y}^T \mathbf{A}^T(\mathbf{A}\mathbf{x}-\mathbf{b}) + \| \mathbf{A}\mathbf{y} \|_2^2 \\ &= \| \mathbf{A}\mathbf{x}-\mathbf{b} \|_2^2 + \| \mathbf{A}\mathbf{y} \|_2^2 \\ & \ge \| \mathbf{A}\mathbf{x}-\mathbf{b} \|_2^2. \end{split}

(3.2.2)

The normal equations have a geometric interpretation, as shown in Figure 3.2.1. The vector in the range (column space) of $\mathbf{A}$ that lies closest to $\mathbf{b}$ makes the vector difference $\mathbf{A}\mathbf{x}-\mathbf{b}$ perpendicular to the range. Thus for any $\mathbf{z}$ , we must have $(\mathbf{A} \mathbf{z})^T(\mathbf{A}\mathbf{x}-\mathbf{b})=0$ , which is satisfied if $\mathbf{A}^T(\mathbf{A}\mathbf{x}-\mathbf{b})=\boldsymbol{0}$ .

Geometry of the normal equations. The smallest residual is orthogonal to the range of the matrix \mathbf{A}. — Figure 3.2.1:Geometry of the normal equations. The smallest residual is orthogonal to the range of the matrix $\mathbf{A}$ .

3.2.1Pseudoinverse and definiteness¶

If we associate the left-hand side of the normal equations as $(\mathbf{A}^T\mathbf{A})\,\mathbf{x}$ , we recognize (3.2.3) as a square $n\times n$ linear system to solve for $\mathbf{x}$ .

Mathematically, the overdetermined least-squares problem $\mathbf{A}\mathbf{x}\approx \mathbf{b}$ has the solution $\mathbf{x}=\mathbf{A}^+\mathbf{b}$ .

Computationally we can generalize the observation about Julia from Chapter 2: backslash is equivalent mathematically to left-multiplication by the inverse (in the square case) or pseudoinverse (in the rectangular case) of a matrix. One can also compute the pseudoinverse directly using pinv, but as with matrix inverses, this is rarely necessary in practice.

The matrix $\mathbf{A}^T\mathbf{A}$ appearing in the pseudoinverse has some important properties.

Proof

The first part is left as Exercise 3.2.3. For the second part, suppose that $\mathbf{A}^T\mathbf{A}\mathbf{z}=\boldsymbol{0}$ . Note that $\mathbf{A}^T\mathbf{A}$ is singular if and only if $\mathbf{z}$ may be nonzero. Left-multiplying by $\mathbf{z}^T$ , we find that

0 = \mathbf{z}^T\mathbf{A}^T\mathbf{A}\mathbf{z}=(\mathbf{A}\mathbf{z})^T(\mathbf{A}\mathbf{z}) = \| \mathbf{A}\mathbf{z} \|_2^2,

(3.2.5)

which is equivalent to $\mathbf{A}\mathbf{z}=\boldsymbol{0}$ . Then $\mathbf{z}$ may be nonzero if and only if the columns of $\mathbf{A}$ are linearly dependent.

Finally, we can repeat the manipulations above to show that for any nonzero $n$ -vector $\mathbf{v}$ , $\mathbf{v}^T(\mathbf{A}^T\mathbf{A})\mathbf{v}=\| \mathbf{A}\mathbf{v} \|_2^2\ge 0$ , and equality is not possible thanks to the second part of the theorem.

3.2.2Implementation¶

The definition of the pseudoinverse involves taking the inverse of a matrix, so it is not advisable to use the pseudoinverse computationally. Instead, we use the definition of the normal equations to set up a linear system, which we already know how to solve. In summary, the steps for solving the linear least squares problem $\mathbf{A}\mathbf{x}\approx\mathbf{b}$ are:

Steps 1 and 3 of Algorithm 3.2.1 dominate the flop count.

In the last step we can exploit the fact, proved in Theorem 3.2.2, that $\mathbf{N}$ is symmetric and positive definite, and use Cholesky factorization as in Exploiting matrix structure. This detail is included in Function 3.2.2.

Algorithm 3.2.2 (lsnormal)

Julia

MATLAB

Python

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 the normal equations

lsnormal.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
function x = lsnormal(A, b)
% LSNORMAL   Solve linear least squares by normal equations.
% Input: 
%   A     coefficient matrix (m by n, m > n)
%   b     right-hand side (m by 1)
% Output:
%   x     minimizer of || b - Ax ||

    N = A' * A;
    z = A' * b;
    R = chol(N);
    w = forwardsub(R', z);    % solve R'z=c
    x = backsub(R, w);        % solve Rx=z
end

Solution of least squares by the normal equations

lsnormal.py

1
2
3
4
5
6
7
8
9
10
11
12
13
def lsnormal(A, b):
    """
    lsnormal(A, b)
    
    Solve a linear least squares problem by the normal equations. Returns the
    minimizer of ||b-Ax||.
    """
    N = A.T @ A
    z = A.T @ b
    R = scipy.linalg.cholesky(N)
    w = forwardsub(R.T, z)                   # solve R'z=c
    x = backsub(R, w)                        # solve Rx=z
    return x

3.2.3Conditioning and stability¶

We have already used A\b as the native way to solve the linear least-squares problem $\mathbf{A}\mathbf{x}\approx\mathbf{b}$ in Julia. The algorithm employed by the backslash does not proceed through the normal equations, because of instability.

The conditioning of the linear least-squares problem relates changes in the solution $\mathbf{x}$ to those in the data, $\mathbf{A}$ and $\mathbf{b}$ . A full accounting of the condition number is too messy to present here, but we can get the main idea. We start by generalizing our previous definition of the matrix condition number.

Provided that the minimum residual norm $\|\mathbf{b}-\mathbf{A}\mathbf{x}\|$ is relatively small, the conditioning of the linear least-squares problem is close to $\kappa(\mathbf{A})$ .

As an algorithm, the normal equations begin by computing the data for the $n\times n$ system $(\mathbf{A}^T\mathbf{A})\mathbf{x} = \mathbf{A}^T \mathbf{b}$ . When these equations are solved, perturbations to the data can be amplified by a factor $\kappa(\mathbf{A}^T\mathbf{A})$ .

The following can be proved using results in Chapter 7.

This squaring of the condition number in the normal equations is the cause of instability. If $\kappa(\mathbf{A})$ is large, the squaring of it can destabilize the normal equations: while the solution of the least-squares problem is sensitive, finding it via the normal equations makes it doubly so.

Example 3.2.1 (Instability in the normal equations)

Julia

MATLAB

Python

Example 3.2.1

Because the functions $\sin^2(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 $\kappa$ 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 $\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

Example 3.2.1

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

t = linspace(0, 3, 400)';
A = [ sin(t).^2, cos((1+1e-7)*t).^2, t.^0 ];
kappa = cond(A)

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 $\kappa$ times machine epsilon.

x_BS = A \ b;
observed_err = norm(x_BS - x) / norm(x)
max_err = kappa * eps

If we formulate and solve via the normal equations, we get a much larger relative error. With $\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);
observed_err = norm(x_NE - x) / norm(x)
digits = -log10(observed_err)

Example 3.2.1

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

from numpy.linalg import cond
t = linspace(0, 3, 400)
A = array([ [sin(t)**2, cos((1+1e-7)*t)**2, 1] for t in t ])
kappa = cond(A)
print(f"cond(A) is {kappa:.3e}")

cond(A) is 1.825e+07

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

x = array([1, 2, 1])
b = A @ x

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

from scipy.linalg import lstsq
x_BS = lstsq(A, b)[0]
print(f"observed error: {norm(x_BS - x) / norm(x):.3e}")
print(f"conditioning bound: {kappa * finfo(float).eps:.3e}")

observed error: 2.019e-10
conditioning bound: 4.053e-09

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

from scipy import linalg
N = A.T @ A
x_NE = linalg.solve(N, A.T @ b)
relative_err = norm(x_NE - x) / norm(x)
print(f"observed error: {relative_err:.3e}")
print(f"accurate digits: {-log10(relative_err):.2f}")

observed error: 5.991e-02
accurate digits: 1.22

3.2.4Exercises¶

Exercise 3.2.8

⌨ Let $t_1,\ldots,t_m$ be $m$ equally spaced points in $[0,2\pi]$ . In this exercise, use $m=500$ .

(a) Let $\mathbf{A}_\epsilon$ be the matrix in (3.1.2) that corresponds to fitting data with the function $c_1 + c_2 \sin(t) + c_3 \cos\bigl((1+\epsilon)t\bigr)$ . Using the identity (3.2.7), make a table of the condition numbers of $\mathbf{A}_\epsilon$ for $\epsilon = 10^0, 10^{-1}, \ldots, 10^{-8}$ .

(b) Repeat part (a) using the fitting function $c_1 + c_2 \sin^2(t) + c_3 \cos^2\bigl( (1+\epsilon) t\bigr).$

(c) Why does it make sense that $\kappa\bigl(\mathbf{A}_\epsilon\bigr)\to \infty$ as $\epsilon\to 0$ in part (b) but not in part (a)?

Preface

Fitting functions to data

Preface

The QR factorization