Section 8.1.3. Implementation of linear least squares

The linear least-squares problem can be viewed geometrically as the problem of finding the point in a p-dimensional subspace, defined as the set of points that can be reached by a linear combination of the columns of A, closest to a given point, y, in an n-dimensional observation space. Since this is equivalent to finding the orthogonal projection of point y into that subspace, it is not surprising that an orthogonal decomposition of A helps to solve the problem. For convenience in this discussion, let us remove the weight matrix from the problem by defining the standardized design matrix by where U is the upper triangular Cholesky factor of W.

Consider the least-squares problem with the QR factorization of Z, as given in Subsection 8.1.1.1. For y′ = U(y − b), (8.1.2.5) becomes which reduces to The second term in (8.1.3.3) is independent of x, and is therefore the sum of squared residuals. The first term vanishes if which, because R is upper triangular, is easily solved for x. The QR decomposition of Z therefore leads naturally to the following algorithm for solving the linear least-squares problem:

Let us now consider the relationship of the QR procedure for solving the linear least-squares problem to the classical method based on the normal equations. The normal equations can be derived by differentiating (8.1.3.2) and equating the result to a null vector. This yields The algorithm is therefore to compute the cross-product matrix, , and the right-hand side, d = Z^Ty′, and to solve the resulting system of equations, Bx = d. This is usually accomplished by computing the Cholesky decomposition of B, that is , where C is upper triangular, and then solving the two triangular systems and . Because , equation (8.1.3.5) becomes or It is clear that R is the Cholesky factor of , although it is formed in a different way. This procedure requires of order operations to form the product and operations for the Cholesky decomposition. In some situations, the extra time to compute the QR factorization is justified because of greater stability, as will be discussed below. Most other quantities of statistical interest can be computed directly from the QR factorization.

The condition number of Z, which is defined (Subsection 8.1.1.1) as the square root of the ratio of the largest to the smallest eigenvalue of , is an indicator of the effect a small change in an element of Z will have on the elements of and of . A large value of the condition number means that small errors in computing an element of Z, owing possibly to truncation or roundoff in the computer, can introduce large errors into the elements of the inverse matrix. Also, when the condition number is large, the standard uncertainties of some estimated parameters will be large. A large condition number, as defined in this way, can result from either scaling or correlation or some combination of these. To illustrate this, consider the matrices and where represents machine precision, which can be defined as the smallest number in machine representation that, when added to 1, gives a result different from 1. By the conventional definition, both of these matrices have a condition number for Z of . Because numbers of order can be perfectly well represented, however, the first one can be inverted without loss of precision, whereas an inverse for the second would be totally meaningless. It is good practice, therefore, to factor the design matrix, Z, into the form where S is a p × p diagonal matrix whose elements define some kind of `natural' unit appropriate to the parameter represented in each column of Z. The ideal natural unit would be the standard uncertainty of that parameter, but this is not available until after the calculation has been completed. If correlation is not too severe, suitable values for the elements of S, of the same order of magnitude as those derived from the standard uncertainty, are the column Euclidean norms, that is where denotes the jth column of Z. This scaling causes all diagonal elements of Z^TZ to be equal to one, and errors in the elements of Z will have roughly equal effects.

Ill conditioning that results from correlation, as in the second example above, is more difficult to deal with. It is an indication that some linear combination of parameters, some eigenvector of the normal equations matrix, is poorly determined by the available data. Use of the QR factorization of Z to compute the Cholesky factor of Z^TZ may be advantageous, in spite of the additional computation time, because better numerical stability is obtained in marginal situations. As a practical matter, however, it is important to recognize that an ill conditioned matrix is a symptom of a flaw in the model or in the experimental design (or both). Use can be made of the fact that, although determining the entire set of eigenvalues and eigenvectors of a large matrix is computationally an inherently difficult problem, a relatively simple algorithm, known as a condition estimator (Anderson et al., 1992), can produce a good approximation to the eigenvector that corresponds to the smallest eigenvalue of a nearly singular matrix. This information can be used in either or both of two ways. First, without any fundamental modification to the model or the experiment, a simple, linear transformation of the parameters so that the problem eigenvector is one of the independent parameters, followed by rescaling, can resolve the numerical difficulties in computing the estimates. A common example is the situation where a phase transition results in the doubling of a unit cell, with pairs of atoms almost but not quite related by a lattice translation. A transformation that makes the estimated parameters the sums and differences of corresponding parameters in related pairs of atoms can make a dramatic improvement in the condition number. Alternatively, the problem eigenvector can be set to some value determined from theory or from some other experiment (see Section 8.3.1 ), or additional data can be collected that are selected to make that combination of parameters determinate.