Linear Least Squares

Suppose we are given a set of data points {(xi , fi )}, i = 1, . . . , n. These could be measurements from an experiment or obtained simply by evaluating a function at some points. You have seen that we can interpolate these points, i.e., either find a polynomial of degree ≤ (n − 1) which passes through all n points or we can use a continuous piecewise interpolant of the data which is usually a better approach. How, it might be the case that we know that these data points should lie on, for example, a line or a parabola, but due to experimental error they do not. So what we would like to do is find a line (or some other higher degree polynomial) which best represents the data. Of course, we need to make precise what we mean by a “best fit” of the data. As a concrete example suppose we have n points (x1 , f1 ), (x2 , f2 ), ··· (xn , fn )

and we expect them to lie on a straight line but due to experimental error, they don’t. We would like to draw a line and have the line be the best representation of the points. If n = 2 then the line will pass through both points and so the error is zero at each point. However, if we have more than two data points, then we can’t find a line that passes through the three points (unless they happen to be collinear) so we have to find a line which is a good approximation in some sense. Of course we need to define what we mean by a good representation. An obvious approach would be to create an error vector of length n and each component measures the difference (fi − y(xi )) where y = a1 x + a0 is the line we fit the data with. Then we can take a norm of this error vector and our goal would be to find the line which minimizes this error vector. Of course this problem is not clearly defined because we have not specified what norm to use. The linear least squares problem finds the line which minimizes this difference in the ℓ2 (Euclidean) norm. Example We want to fit a line p1 (x) = a0 + a1 x to the data points (1, 2.2), (.8, 2.4), (0, 4.25)

in a linear least squares sense. For now, we will just write the overdetermined system and determine if it has a solution. We will find the line after we investigate how to solve the linear least squares problem. Our equations are a0 + a1 ∗ 1 = 2.2 a0 + a1 ∗ .8 = 2.4 a0 + a1 ∗ 0 = 4.25 1

Writing this as a matrix problem Ax = b we have   1 1  1 0.8  1 0 a0 a1  2.1 =  2.4  4.25 

Now we know that this over-determined problem has a solution if the right hand side is in R(A) (i.e., it is a linear combination of the columns of the coefficient matrix A). Here the rank of A is clearly 2 and thus not all of I 3 . Moreover, (2.1, 2.4, 4.25)T is not in the R(A), R T T i.e., not in the span{(1, 1, 1) , (1, 0.8, 0) } and so the system doesn’t have a solution. This just means that we can’t find a line that passes through all three points.

Example

If our data had been (1, 2.1) (0.8, 2.5) (0, 4.1)

then would we have had a solution to the over-determined system? Our matrix problem Ax = b is  1 1  1 0.8  1 0  a0 a1  2.1 =  2.5  4.1 

and we notice that in this case, the right hand side is in R(A) because      2.1 1 1  2.5  = 4.1  1  − 2  0.8  4.1 1 0  and thus the system is solvable and we have the line 4.1 − 2x which passes through all three points. But, in general, we can’t solve the over-determined system so our approach is to find a vector x such that the residual r = b − Ax is as small as possible. The residual is a vector and so we take the norm. The linear least squares method uses the ℓ2 -norm. Consider the over-determined system Ax = b where A is m × n with m > n. The linear least squares problem is to find a vector x which minimizes the ℓ2 norm of the residual that is x = min b − Az n z∈I R 2

b − Az

2

for all z ∈ I n R

2

We note that minimizing the ℓ2 norm of the residual is equivalent to minimizing its square. This is often easier to work with because we avoid dealing with square roots. So we rewrite the problem as find a vector x which minimizes the square of the ℓ2 norm b − Az
2 2

for all z ∈ I n R

Example points

For our example where we want to fit a line p1 (x) = a0 + a1 x to the data (1, 2.2),   2.2 1 r =  2.4  −  1 4.25 1  (.8, 2.4),  1 0.8  0 (0, 4.25)  2.2 − z1 − z2 =  2.4 − z1 − 0.8z2  4.25 − z1 

calculate the residual vector and then use techniques from Calculus to minimize r 2 . 2 z1 z2

To minimize r 2 we take the first partials with respect to z1 and z2 and set them equal 2 to zero. We have f= r and thus ∂f = −4.4 + 2z1 + 2z2 − 4.8 + 2z1 + 1.6z2 − 8.5 + 2z1 = 17.7 + 6z1 + 3.6z2 = 0 ∂z1 ∂f = −4.4 + 2z1 + 2z2 − 3.84 + 1.6z1 + 1.28z2 = −8.24 + 3.6z1 + 3.28z2 = 0 ∂z2 So we have to solve the linear system 6 3.6 3.6 3.28 whose solution is (4.225, −2.125)T . We now want to determine 1. Does the linear least squares problem always have a solution? 2. Does the linear least squares problem always have a unique solution? 3. How can we efficiently solve the linear least squares problem? Theorem The linear least squares problem always has a solution. It is unique if A has linearly independent columns. The solution of the problem can be found by solving the normal equations AT Ay = AT b . 3 z1 z2 = −17.7 8.24
2 2

= (2.2 − z1 − z2 )2 + (2.4 − z1 − .8z2 )2 + (4.25 − z1 )2

Before we prove this, recall that the matrix AT A is symmetric because (AT A)T = AT (AT )T = AT A and is positive semi-definite because xT (AT A)x = (xT AT )(Ax) = (Ax)T (Ax) = y T y ≥ 0 where y = Ax

Now y T y is just the square of the Euclidean length of y so it is only zero if y = 0. Can y ever be zero? Remember that y = Ax so if x ∈ N (A) then y = 0. When can the rectangular matrix A have something in the null space other than the zero vector? If we can take a linear combination of the columns of A (with coefficients nonzero) and get zero, i.e., if the columns of A are linearly dependent. Another way to say this is that if the columns of A are linearly independent, then AT A is positive definite; otherwise it is positive semi-definite (meaning that xT AT Ax ≥ 0). Notice in our theorem we have that the solution is unique if A has linearly independent columns. Another equivalent statement would be to require N (A) = 0. Proof First we show that the problem always has a solution. Recall that R(A) and N (AT ) are orthogonal complements in I m . This tells us that we can write any vector in R m I R as the sum of a vector in R(A) and one in N (AT ). To this end we write b = b1 + b2 where b1 ∈ R(A) and b2 ∈ R(A)⊥ = N (AT ). Now we have the residual is given by b − Ax = (b1 + b2 ) − Ax Now b1 ∈ R(A) and so the equation Ax = b1 is always solvable which says the residual is r = b2 When we take the r 2 we see that it is b2 2 ; we can never get rid of this term unless b ∈ R(A) entirely. The problem is always solvable and is the vector x such that Ax = b1 where b1 ∈ R(A). When does Ax = b1 have a unique solution? It is unique when the columns of A are linearly independent or equivalently N (A) = 0. Lastly we must show that the way to find the solution x is by solving the normal equations; note that the normal equations are a square n × n system and when A has linearly independent columns the coefficient matrix AT A is invertible with rank n. If we knew what b1 was, then we could simply solve Ax = b1 but we don’t know what the decomposition of b = b1 + b2 is, simply that it is guaranteed to exist. To demonstrate that the x which minimizes b − Ax 2 is found by solving AT Ax = AT b we first note that these normal equations can be written as AT (b − Ax) = 0 which is just AT times the residual vector so we need to show AT r = 0 to prove the result. From what we have already done we know that AT (b − Ax) = AT (b1 + b2 − Ax) = AT (b2 ) 4

Recall that b2 ∈ R(A)⊥ = N (AT ) which means that AT b2 = 0 and we have that AT (b − Ax) = 0 ⇒ AT Ax = AT b The proof relies upon the fact that R(A) and N (AT ) are orthogonal complements and that this implies we can write any vector as the sum of a vector in R(A) and its orthogonal complement. Example We return to our previous example and now determine the line which fits the data in the linear least squares sense; after we obtain the line we will compute the ℓ2 norm of the residual. We now know that the linear least squares problem has a solution and in our case it is unique because A has linearly independent columns. All we have to do is form the normal equations and solve as usual. The normal equations     1 1 2.2 1 1 1  a0 1 1 1  1 0.8  = 2.4  1 0.8 0 a1 1 0.8 0 1 0 4.25 are simplified as 3.0 1.8 1.8 1.64 a0 a1 = 8.85 4.12

which has the solution (4.225, −2.125) giving the line y(x) = 4.225−2.125x. If we calculate the residual vector we have     2.2 − y(1) 0.1  2.4 − y(0.8)  =  −0.125  4.25 − y(0) 0.025 which has an ℓ2 norm of 0.162019. We said that we only talked about the inverse of a square matrix. However, one can define a pseudo-inverse of a rectangular matrix. If A is an m×n matrix with linearly independent columns then a pseudo-inverse (or sometimes called left inverse of A ) is A† = (AT A)−1 AT which is the matrix in our solution to the normal equations x = (AT A)−1 AT b It is called the pseudo-inverse of the rectangular matrix A because (AT A)−1 AT A = (AT A)−1 (AT A) = I Note that if A is square and invertible the pseudo-inverse reduces to A−1 because (AT A)−1 AT = A−1 (AT )−1 AT = A−1 . 5

We can also find a polynomial of higher degree which fits a set of data. The following example illustrates this. Example State the linear least squares problem of finding the quadratic polynomial which fits the following data in a linear least squares sense; determine if it has a unique solution; calculate the solution and calculate the ℓ2 norm of the residual vector. (0, 0) (1, 1) (3, 2) (4, 5)

In this case we seek a polynomial of the form p(x) = a0 + a1 x + a2 x2 . Our over determined system is     0 1 0 0   a0 1 1 1   1  a1 =    2 1 3 9 a2 5 1 4 16 So the linear least squares problem is to find a vector x in R3 which minimizes     0 z1 1  z2   −A 2 z3 5
2

2

for all z ∈ R3 where A is the 4 × 3 matrix given above. We see that A has linearly independent columns so its rank is 3 and thus the linear least squares problem has a unique solution. The normal equations are  1 1 0 1 0 1 1 1 1 1 3 4  1 9 16 1  0 1 3 4  0    a 1 1  0   a1 = 0 9 a2 0 16 1 1 1 3 1 9 0 1 1 4   2 16 5

leading to the square system 4  8 26     8 26 a0 8   a1  =  27  26 92 a2 92 338 99

Solving this we get a0 = 3/10, a1 = 7/30, a2 = 1/3. Our residual vector is    0.3 0 − p(0)  1 − p(1)   0.6  r=  = 0.6 2 − p(3) 0.3 5 − p(4)  6

and the square of its ℓ2 norm is r
2 2

= .32 + .62 + .62 + .32 = 0.9

Now it seems as if we are done because we know when the solution is unique and have a method for determining the solution when it is unique. What else do we need? Unfortunately, finding the normal equations works well for hand calculations but is not the preferred method for computations. Why is this? To form the normal equations we must compute AT A. This can cause problems as the following theorem tells us. Theorem Let A have linearly independent columns. Then
2 K2 (A) = K2 (AT A)

where K(A) = A A† Thus when we form AT A we are squaring the condition number of the original matrix. This is the major reason that solving the normal equations is not a preferred computational method. A more subtle problem is that the computed AT A may not be positive definite even when A has linearly independent columns so we can’t use Cholesky’s method. Can we use any of our previous results from linear algebra to help us solve the linear least squares problem? We looked at three different decompositions: LU and its variants, QR and the SVD. We use LU (or its variants) to solve the normal equations. Can we use QR or the SVD of A? In fact, we can use both. Recall that an m × n matrix with m > n and rank n has the QR decomposition A=Q R 0

where Q is an m × m orthogonal matrix and R is an n × n upper triangular matrix and 0 represents an (m − n) × n zero matrix. Now to see how we can use the QR decomposition to solve the linear least squares problem, we take QT r where r = b − Ax to get QT r = QT b − QT Ax = QT b − QT Q R 0 x

Now Q is orthogonal so QT Q = I so if we let QT b = (c, d)T , we have QT r = QT b − R 0 x= c d 7 − Rx 0 = c − Rx d

Now also recall that an orthogonal matrix preserves the ℓ2 length of any vector, i.e., Qy 2 = y 2 for Q orthogonal. Thus we have QT r and hence r
2 2 2

= r

2

= QT r

2 2

= c − Rx

2 2

+ d

2 2

So to minimize the residual we must find x which solves Rx = c and thus the minimum value of the residual is r 2 = d 2 . In conclusion, once we have a 2 2 QR decomposition of A with linearly independent columns then the solution to the linear least squares problem is the solution to the upper triangular system Rx = c where c is the first n entries of QT b and the residual is the remaining entries of QT b. Now we want to see how we can use the SVD to solve the linear least squares problem. Recall that you had a homework problem to use the SVD to compute the pseudo-inverse of an m × n matrix A so we essentially need this. Recall that the SVD of an m × n matrix A is given by A = U ΣV T where U is an m×m orthogonal matrix, V is an n ×n orthogonal matrix and Σ is an m×n diagonal matrix (i.e., Σij = 0 for all i = j). Note that this also says that U T AV = Σ. Because here m > n we write Σ as Σ= Σ 0 0 0

where Σ is a square invertible diagonal matrix. The following result gives us the solution to the linear least squares problem. Theorem given by Let A have the singular value decomposition given above. Then the vector x x=V Σ−1 0 0 0 UT b

minimizes b − Az 2 , i.e., x is the solution of the linear least squares problem. We compute our residual and use the fact that V V T = I to get r = b − Ax = b − AV V T x Now once again using the fact that an orthogonal matrix preserves the ℓ2 length of a vector, we have r
2 2

= b−AV V T x

2 2

= U T (b−AV V T x)

2 2

= U T b−(U T AV )V T x) 8

2 2

= U T b−ΣV T x)

2 2

Writing U T b = (c1 , c2 )T and V T x = (z1 , z2 )T we have r
2 2

=

c1 c2

−

Σ 0

0 0

z1 z2

2

=
2

c1 − Σz1 c2

2 2

So the residual is minimized when c1 − Σz1 = 0; note that z2 is arbitrary so we set it to zero. We have V Tx = z = Σ−1 c1 0 ⇒x=V Σ−1 c1 0 ⇒x=V Σ−1 0 0 0 UT b

because U T b = (c1 , c2 )T .

9