Our previous discussions of least squares always led us to matrices of the form which then needed to be inverted. Golub's method of using Householder transformations works directly with the matrix and has the advantage that it is considerably more accurate than methods with invert .It seems that about twice as much precision is required to invert than is needed to deal directly with .Another reason for learning about Golub's method is that the calculation is organized in a completely different way; therefore, it will often turn out to have other advantages or disadvantages which differ from one application to the next.
A reflection transformation is a matrix of the form where is an arbitrary vector. Obviously is symmetric, that is, .It also turns out that the reflection transformation is its own inverse, that is, .To see this, we verify by substitution that .
A matrix transformation is said to be unitary if .When a matrix is unitary it means that the vector has the same length as the vector .These lengths are and which are the same. Reflection transformations are unitary because .They have a simple physical interpretation. Consider an orthogonal coordinate system in which one of the coordinate axes is aligned along the vector. Reflection transformation reverses the sign of this coordinate axis vector (since ) but it leaves unchanged all the other coordinate axis vectors. Thus it is obvious geometrically that reflection transformations preserve lengths and that applying the transformation twice returns any original vector to itself. Now, we seek a special reflection transformation called the Householder transformation which converts a matrix of the form on the left side to the form on the right where a is an arbitrary element
Having determined the required transformation, we will know how to convert any matrix to an upper triangular form like
by a succession of Householder transforms. Golub recognized the value of this technique in solving overdetermined sets of simultaneous equations. He noted that when the error vector is transformed by a unitary matrix the problem of minimizing the length of by variation of reduces to exactly the same problem as minimizing the length of with respect to variation of .Thus a succession of Householder transforms could be found to reduce to the form
Now we turn to the task of finding the special reflection transformation, called the Householder transformation, which accomplishes (51). Observe that the left-hand operator below is a reflection transformation for any numerical choice of s.
Alternatively, if (54) is to be valid, then s must take a particular value such that
This will be true only for s given by
Now let us see why the left-hand operator in (54) can achieve (51). Choice of the a vector as the third column in the matrix of (51) introduces the desired zeros on the right-hand side. Finally, it is necessary also to observe that this choice of does not destroy any of the zeros which already existed on the left-hand side in (51).
A subroutine for least squares fitting programmed by Don C. Riley.
(This program does not do the square matrix case.
It is necessary that M > N.)
SUBROUTINE GOLUB (A,X,B,M,N)
C A(M,N) ; B(M) GIVEN WITH M>N SOLVES FOR X(N) SUCH THAT
C || B - AX || = MINIMUM
C METHOD OF G.GOLUB, NUMERISCHE MATHEMATIK 7,206-216 (1965)
IMPLICIT DOUBLE PRECISION (D)
C.........PERFORM N ORTHOGONAL TRANSFORMATIONS TO A(.,.) TO
C.........UPPER TRIANGULARIZE THE MATRIX
DO 3010 K=1,N
DO 1010 I=K,M
DO 1020 I=KPLUS,M
C.........I - UU' IS A SYMMETRIC, ORTHOGONAL MATRIX WHICH WHEN APPLIED
C......... TO A(.,.) WILL ANNIHILATE THE ELEMENTS BELOW THE DIAGONAL K
DO 2030 J=K,N
c.........APPLY THE ORTHOGONAL TRANSFORMATION
DO 2010 I=K,M
DO 2020 I=K,M
DO 2040 I=K,M
DO 2050 I=K,M
C.........BACK SUBSTITUTE TO RECURSIVELY YIELD X(.)
DO 4020 I=1,LIM
IROW=N-I ! error in first edition, was =N-1
DO 4010 J=1,I
Householder transformations can also be used in problems with constraints. In the set
one may desire to satisfy the top block exactly and the bottom block only in the least-squares sense. Define as a succession of Householder transforms on ;for example, .Then substitute into (58). Householder transforms used as postmultipliers on the matrix of (58) can be chosen to introduce zeros in the top two rows of (58), for example
Now we could use premultiplying Householder transforms on (59) to bring it to the form
Since the top two equations of (59) or of (60) are to be satisfied exactly, then y1 and y2 are uniquely determined. They cannot be adjusted to help attain minimum error in the bottom three equations. Likewise the top two equations place no restraint on y3 and y4, so they may be adjusted to produce minimum error in the bottom three equations. No amount of adjustment in y3 and y4 can change the amount of error in the last equation, so we can ignore the last equation in the determination of y3 and y4. The third and fourth equations can be satisfied with zero error by suitable choice of y3 and y4. This must be the minimum-squared-error answer. Given we can go back and get with .