Linear Algebra

William F. Moss¹

1  Least Squares Linear Models

A linear model of the form

As an example, suppose a straight line is to be used to fit m ³ 2 data points. We set n = 2 and f₁(t) = 1 and f₂(t) = t. Then f has the form f(t;c) = c₁ + c₂t. Generally, this line will be an interpolant only in case m = 2. A straight line fit is often used in statistics where this line is called the regression line.

If the linear model is constructed so that it changes much more slowly than an interpolant, then curve fitting can be used as a noise filter. This approach is illustrated in the first problem of this project. Here, samples taken from a straight line are corrupted with noise generated randomly, and then a straight line is fit to the data. The original line is very nearly recovered.

In the mathematical discussion below, we will use the following matrix notation which was introduced in the Orthogonal Methods course outline. Let

The interpolation problem is to find c so that the interpolation equations

The curve fitting problem is to find c so that the graph of f passes "close to" the data points. The most commonly used method for defining "close to" is called the least squares method. In the least squares method, c is determined by minimizing the sum of the squares of the vertical deviations between the graph of f and the data points. The problem is to find c_ls so that

Figure 1 provides geometric insight into the problem of solving (3) and (8) for the case m = 3 and n = 2. Here y Î IR^{3 ×1} and A = [a₁,a₂] is 3 ×2 with two linearly independent columns, a₁ and a₂. Now Ran(A) is the span{a₁,a₂}, that is, the set of all linear combinations of the columns of A. It is a two-dimensional subspace of IR^{3 ×1} and is visualized in Figure 1 by the plane which passes through the origin and contains a₁,a₂. Also, for any c Î IR^{2 ×1}, A c = c₁ a₁+c₂ a₂ Î Ran(A). Now if y is not in the plane Ran(A), as is the case in Figure 1, the linear system (3) and the interpolation problem will not have a solution. On the other hand, the least squares problem always has a solution c_ls. The vector A c_ls is the vector in the plane Ran(A), which is closest to y. Also, the vector y - Ac_ls will be orthogonal to the the plane Ran(A). Because the columns of A are linearly independent, or equivalently because rank(A) = 2, c_ls is the only solution to the least squares problem. Now imagine a case where a₁ and a₂ lie along the same line. In this case, Ran(A) is a one-dimensional subspace of IR^{3 ×1} and can be visualized by a line passing through the origin. Now there are infinitely many ways to describe a vector in Ran(A) as a linear combination of a₁ and a₂, and so there are infinitely many solutions to the least squares problem. We would prefer to have a unique linear model which fits the data. This will happen if the basis functions used in the linear model are chosen so that the columns of A are linearly independent.

The m > 3 case works the same ways as the m = 3 case. For c, A c Î Ran(A), a subspace of IR^{m ×1}. In the least squares problem, we find the vector c_ls so that A c_ls Î Ran(A) is as close to the vector y as possible. Because we have assumed that n £ m, rank(A) £ n. If the basis functions have been well chosen, the columns of A will be linearly independent and rank(A) = n. In this case, (8) is said to be a full rank least squares problem and has a unique solution. If rank(A) < n, (8) is said to be rank deficient and has infinitely many solutions.

In Orthogonal Methods course outline, we have see that the singular value decomposition (SVD) of A can be used to solve the least squares problem. We note that in the m = n case, it can also be used to solve (3), although Gauss elimination with partial pivoting is more frequently used. In the Orthogonal Methods course outline a formula is given for c_ls. In the rank deficient case, this formula gives the least squares solution which has the smallest 2-norm. There are two other methods that are commonly used to solve the full rank least squares problem: QR factorization with column pivoting followed by the solution of a triangular system, and solution of the normal equations. However, the SVD approach is regarded as the most computationally reliable method.

As explained in the instructions for this project, you are to create an M-file with filename P1loginid.m containing MATLAB commands and statements to solve the following problems. Note that you can cut MATLAB code from this PDF document and paste it into the MATLAB editor.

Problem 1 In this problem you will fit noisy data using a line. We generate the data by taking 20 uniformly spaced samples from the line y = 1 + 2t for t Î [0,1] and adding a small amount of random noise to each sample. Add the following lines to P1loginid.m.

% Project 2. Last Name First Name. GroupXX

echo on

% Problem 1
% create m > 2 noisy data points lying close to a line
m = 20;
t = linspace(0,1,m)';
epsilon = .2; % Maximum noise amplitude
y = 1 + 2*t + epsilon*rand(m,1);
% Use the model c(1) + c(2)*t
% Construct the design matrix A.
A = [ones(m,1),t];
% Use the SVD to solve the least squares problem.
c = pinv(A)*y
if rank(A) == 2
   disp('The design matrix is full rank')
else
   disp('The design matrix is rank deficient.')
end
% Since A is full rank, we have found the unique
% solution to the full rank least squares problem.
% When the design matrix is rank deficient, the
% least squares problem has infinitely many
% solutions. The above computation would produce
% the solution that has the smallest 2-norm.
%
% Plot
tt = linspace(0,1,100)';
yy = c(1) + c(2)*tt;
plot(t,y,'o',tt,yy)

Problem 2 The water level in the North Sea is mainly determined by the so-called M₂-tide, whose period is about 12 hours. A reasonable but simple model for predicting the level at any time t is a function of the form

Suppose we have made the measurements:

Add the following lines to your M-file to solve this linear least squares problem having m = 6 data points and n = 3 basis functions.

% Problem 2

m = 6;
t = [0;2;4;6;8;10]; % hours
y = [1;1.6;1.4;0.6;0.2;0.8]; % meters
% design matrix
A = [ones(m,1),sin(pi*t/6),cos(pi*t/6)];
c = pinv(A)*y % least squares solution
% Plot the data and the model
% plot vectors
tt = linspace(0,10,100)';
yy = [ones(100,1),sin(pi*tt/6),cos(pi*tt/6)]*c;
figure % create a new plot window
plot(t,y,'o',tt,yy)

Problem 3 The interpolation equations (2) can be satisfied if and only if the linear system (3) has a solution. If n < m, there will be more equations than unknowns and this overdetermined system will not generally have a solution. This is clear in Figure 1 where m = 3 and n = 2. Show that this is the case for the example in Problem 2 by applying the rank test from the Linear Systems and Matrices Outline. You must compare the rank of A to the rank of [A,y].

Problem 4 Add lines to your M-file to find a straight line fit to the following data and plot a graph of the straight line model along with the data on the same set of axes. Label your graph. Note that here we have (x,y) data instead of (t,y) data.

Problem 5 Balancing a Chemical Reaction. ² Balance the chemical reaction

Problem 6 The following problem arises later in the course when we model heat conduction in a rod and a vibrating string. Find a scalar l (eigenvalue) and a nontrivial function X(x) (eigenfunction) so that

First, we discretize the interval [0,1] with mesh points x_i = ih,  i = 0,¼,n+1, and h = 1/(n+1). Next, we use the approximation

Add the following lines to your M-file to solve this problem. Experiment with the value of n.

% Problem 6

n = 100;
h = 1/(n+1); e = ones(n,1)/h^2;
A = spdiags([-e 2*e -e],-1:1,n,n);
z = eig(A); format long
y = sort(z)/pi^2;
y(1)

Footnotes:

¹ Department of Mathematical Sciences, Clemson University, Clemson, SC 29634-1907, U.S.A. (bmoss@math.clemson.edu). Copyright 1995,1996 William F. Moss. All rights reserved.

² This problem appears in FirstLeaves: A Tutorial Introduction, B. W. Char et. al., Springer 1992, page 197 (ISBN 0-387-97621-3).