Least Squares for Practitioners

Abstract

In experimental science and engineering, least squares are ubiquitous in analysis and digital data processing applications. Minimizing sums of squares of some quantities can be interpreted in very different ways and confusion can arise in practice, especially concerning the optimality and reliability of the results. Interpretations of least squares in terms of norms and likelihoods need to be considered to provide guidelines for general users. Assuming minimal prerequisites, the following expository discussion is intended to elaborate on some of the mathematical characteristics of the least-squares methodology and some closely related questions in the analysis of the results, model identification, and reliability for practical applications. Examples of simple applications are included to illustrate some of the advantages, disadvantages, and limitations of least squares in practice. Concluding remarks summarize the situation and provide some indications of practical areas of current research and development.

1. Introduction

Least squares go back to Gauss and Legendre in the late 1790s. The first important publication on the topic was authored by Legendre in 1806 with the title “New Methods for Determination of a Comet’s Orbit” and had a supplement entitled “On the method of least squares”. Gauss’s first publication on least squares appeared in 1809 at the end of his Theoria motus. He mentioned there in passing that Legendre had presented the method in his work of 1806, but that he himself had already discovered it in 1795. Gauss’s correspondence and the papers found after his death proved that he was certainly the first to make the discovery, but since Legendre was first to publish it, priority rights belong to the latter. Obviously, as both of them reached the result independently of each other, both deserve the honour [1].

In these astronomical applications, “least squares” was a method of obtaining the best possible average value for a measured magnitude, given several observations of the magnitude, when the measurements are found to be unavoidably different due to (random) errors. Of course, long before the theory of errors ever saw the light of day, common sense had chosen the arithmetic average value as the most probable value, hence the dichotomy in numerous situations: independently of the nature of the errors involved, a least-squares procedure gives the arithmetic mean, or more generally some weighted average value, which may or may not be the most likely value in the probabilistic sense (see, e.g., [2] for more general discussions). These two very different interpretations of least squares, now technically often referred to as the Best Approximation Estimate (BAE) in terms of quadratic norms, and the Maximum Likelihood Estimate (MLE) in terms of distributions, respectively, are widely used in functional analysis, inference analysis, and all kinds of application areas. Notice that in general, BAEs in terms of arbitrary norms can be very different from MLEs in terms of distributions. However, it is well known that in general, BAEs for p norms and MLEs for exponential distributions coincide with most interesting implications in practice (see, e.g., [3]). The following discussions will concentrate on the fundamental characteristics of the least-squares methodology and related implementation aspects.

Least-squares parameter estimation can be applied to underdetermined just as to overdetermined linear problems. In fact, underdetermined prediction problems are generally more common than overdetermined filtering and adjustment problems. With observations and unknown parameters of unequal weights modeled using some empirical or theoretical covariance (or correlation) functions, more sophisticated estimation methods such as Kriging and least-squares collocation are employed. Correspondingly, Radial Basis Functions (RBFs) and related strategies are used for interpolations of spatially scattered data and other approximations as BAEs.

The preceding implicit assumption of model linearity is essential for practical reasons. In practice, just about any engineering problem nonlinear in terms of its unknown parameters can be linearized as follows:(a)Using a Taylor expansion about an appropriate point or parameter value, the linear term can be used to approximate the model in the neighborhood of the expansion point.(b)Using differentiation in terms of the unknown parameters, the total derivative of the function is linear in terms of the differentials corresponding to the unknown parameters, and hence differential corrections to the unknown parameters can be evaluated as a linear problem.

These two strategies for nonlinear least squares then imply iterative procedures for differential correction estimates to the unknown parameters. Convergence of such iterative procedures is usually ensured for well-chosen parametrization, and the general situation has been discussed by Pope [4] and others in the estimation literature. Obviously, in general, there are numerous types of nonlinear complex problems that cannot be treated with such simplistic strategies but for the purposes of the following discussions, linearity in terms of the unknown parameters will be assumed from here on.

From an application perspective, one exceptional class of separable nonlinear least-squares problems deserves mention in this context. These are problems for which the mathematical model function is a linear combination of nonlinear functions. Specifically, one can assume that there are two sets of unknown parameters where one set is dependent on the other and can be explicitly eliminated. The method of variable projections has proven very appropriate for such nonlinear problems in several application areas [5]. General nonlinear least-squares estimation is still the object of current research (see, e.g., [6] for further discussions and references).

The intrinsic linearity of least-squares computations implies that these can be done simultaneously or in a stepwise manner to obtain exactly the same estimation results. In other words, at the limit, one unknown parameter can be estimated at a time, or one observation or measurement can be processed at a time. This characteristic of linear computations is most useful in least-squares procedures and has led to numerous formulations such as summation of normals, and sequential adjustments. Furthermore, the quadratic computations can be avoided in critical numerically sensitive situations using orthogonal methods such as Givens rotations, Householder’s reflections, and others.

2. Least Squares and Alternatives

Consider a system of 𝑀 linear algebraic nonhomogeneous equations with 𝑁 unknowns 𝑥₁, 𝑥₂, 𝑥₃,…, 𝑥_𝑁, where 𝑀 is not necessarily equal to 𝑁, 𝑎₁₁𝑥₁+𝑎₁₂𝑥₂+𝑎₁₃𝑥₃+⋯+𝑎_1𝑁𝑥_𝑁=𝑓₁,𝑎₂₁𝑥₁+𝑎₂₂𝑥₂+𝑎₂₃𝑥₃+⋯+𝑎_2𝑁𝑥_𝑁=𝑓₂,𝑎₃₁𝑥₁+𝑎₃₂𝑥₂+𝑎₃₃𝑥₃+⋯+𝑎_3𝑁𝑥_𝑁=𝑓₃,⋮𝑎_𝑀1𝑥₁+𝑎_𝑀2𝑥₂+𝑎_𝑀3𝑥₃+⋯+𝑎_𝑀𝑁𝑥_𝑁=𝑓_𝑀 (2.1) with corresponding matrix representation 𝐀𝐱=𝐟, (2.2) and assuming 𝐟≠𝟎 for simplicity. When 𝑀 = 𝑁 without any rank deficiencies in the matrix 𝐀, that is, with 𝐀 nonsingular, then the unique solution is simply 𝐱=𝐀⁻¹𝐟, (2.3) which can be evaluated in practice by Gaussian elimination or any other solution method for simultaneous linear equations. Notice that such solution methods are usually more efficient than the direct matrix inversion method, which is a consideration in numerous application contexts.

When the system is overdetermined with 𝑀>𝑁, that is, more equations than unknowns, then one could use the first 𝑁 equations, or some other selection of 𝑁 equations, and assuming no rank deficiency, proceed as in the previous case of 𝑀 = 𝑁. However, this is not appropriate for most applications as all the observations should somehow contribute to some “optimal” solution. Hence, rewriting the given system of equations with an error term e, that is, 𝐀𝐱=𝐟+𝐞, to emphasize that there may not exist one 𝐱 value that would satisfy 𝐀𝐱=𝐟 exactly, one obvious strategy is to minimize some norm of 𝐞, that is, some acceptable measure of the “length” of the vector e. In practical terms, this norm of e can simply be its Euclidean length, that is, ‖𝐞‖₂=(𝑒²₁+𝑒²₂+𝑒²₃+⋯+𝑒²_𝑀)^1/2=(|𝑒₁|²+|𝑒₂|²+|𝑒₃|²+⋯+|𝑒_𝑀|²)^1/2, (2.4) but more generally, using 𝑝 norms, denoted by 𝐿_𝑝, ‖𝐞‖_𝑝=(^𝑀_𝑖=1|𝑒_𝑖|^𝑝)^1/𝑝 (2.5) for 𝑝=1,2,…,∞. The solution 𝐱 for a specified value of 𝑝, if one exists, is called an 𝐿_𝑝 BAE of 𝐱. For 𝑝=2, the 𝐿₂ estimate is the familiar least-squares estimate of 𝐱, which is going to be discussed below. When p = 1, the 𝐿₁ estimate is a least-magnitude estimate of 𝐱, a generalization of the median, and is well known in robust estimation. When 𝑝=∞, the 𝐿_∞ estimate is a least-maximum or min-max estimate of 𝐱. For other values of 𝑝, some BAEs are possible but not often used in practice, except perhaps for 1 <𝑝< 2 in multifacility location-allocation problems (e.g., [7]). Notice that for 𝑝≠2, the BAE does not necessarily exist and when it does exist, it is not necessarily unique, which can greatly complicate matters in applications.

When 𝑝=2, the least-squares estimate always exists for a finite set of linear equations, assuming linearly independent columns, and its unique value is easily obtained using basic calculus: ‖𝐞‖²₂=𝐞^𝑇𝐞=𝑒²₁+𝑒²₂+𝑒²₃+⋯+𝑒²_𝑀 (2.6) using matrix notation and to minimize 𝐞^𝑇e = (𝐀𝐱−𝐟)^𝑇(𝐀𝐱−𝐟), 𝜕𝜕𝐱(𝐞^𝑇𝜕𝐞)=𝜕𝐱((𝐀𝐱−𝐟)^𝑇(𝐀𝐱−𝐟))=𝟎, (2.7) which gives the familiar normal equations 𝐀^𝑇𝐀𝐱=𝐀^𝑇𝐟, (2.8) where the square matrix 𝐀^𝑇A is easily seen to be symmetric and positive definite. The least-squares estimate is then written as 𝐱=(𝐀^𝑇𝐀)⁻¹𝐀^𝑇𝐟 (2.9) which is easily verified to correspond to a minimum as 𝜕²𝜕𝐱²(𝐞^𝑇𝜕𝐞)=²𝜕𝐱²((𝐀𝐱−𝐟)^𝑇(𝐀𝐱−𝐟))=2𝐀^𝑇𝐀>𝟎. (2.10) The previous matrix inequality simply means that, as the matrix 𝐀^𝑇A is symmetric and positive definite, it has positive real eigenvalues and hence the situation corresponds to a minimum of 𝐞^𝑇e, as desired.

For the corresponding underdetermined system Ax = f, f≠0, assuming linearly independent rows, there are obviously infinitely many solutions in general. For an optimal solution 𝐱 with minimum quadratic norm, the easiest approach is to use unknown correlates x = 𝐀^𝑇y which imply by substitution 𝐀𝐱=𝐀𝐀^𝑇𝐲=𝐟 (2.11) and assuming no rank deficiency as before, AA^T is nonsingular and hence 𝐲=(𝐀𝐀^𝑇)⁻¹𝐟, (2.12) which gives by substitution 𝐱=𝐀^𝑇(𝐀𝐀^𝑇)⁻¹𝐟. (2.13) To see the appropriateness of this estimate 𝐱, consider 𝐱−𝐱=𝐱−𝐀^𝑇(𝐀𝐀^𝑇)⁻¹𝐟=𝐱−𝐀^𝑇(𝐀𝐀^𝑇)⁻¹𝐀𝐱=[𝐈−𝐀^𝑇(𝐀𝐀^𝑇)⁻¹𝐀]𝐱=𝐍_𝐀𝐱, (2.14) with 𝐍_𝐀 usually called the nullspace projector corresponding to A. More explicitly,(i)for a vector z with Az = 0, 𝐍_𝐀𝐳=𝐳, and conversely;(ii)𝐀𝐍_𝐀𝐳=𝟎 for all vector z.

Therefore, for any vector z, 𝐀(𝐱+𝐍_𝐀𝐳)=𝐀𝐱+𝐀𝐍_𝐀𝐳=𝐟 (2.15) but for a minimum quadratic norm estimate 𝐱, only z≡0 is acceptable.

The previous result can readily be generalized to the quadratic norm with weight matrices for correlated observations or measurements of different quality as follows. Let f denote an 𝑀 nonzero vector and A an 𝑀×𝑁 matrix with linearly independent columns. Then there is a unique 𝑁 vector 𝐱 which minimizes {(𝐟−𝐀𝐱)^𝑇𝐏(𝐟−𝐀𝐱)}^1/2 over all 𝐱, for some appropriate weight matrix P. Furthermore, 𝐱 = (𝐀^𝑇𝐏𝐀)⁻¹𝐀^𝑇𝐏𝐟.

More generally, the unknowns themselves may have different relevance or other characteristics requiring weight matrices such as for regularization. Letting f denote an 𝑀 nonzero vector and A an 𝑀×𝑁 matrix with linearly independent columns, then there is a unique 𝑁 vector 𝐱 which minimizes {(𝐟−𝐀𝐱)^𝑇𝐏(𝐟−𝐀𝐱)+𝐱^𝑇𝐐𝐱}^1/2 over all 𝐱, for some appropriate weight matrices P and Q. Furthermore, 𝐱 = (𝐀^𝑇𝐏𝐀+𝐐)⁻¹𝐀^𝑇𝐏𝐟 and when P and Q are nonsingular, 𝐱 = Q^-1𝐀^𝑇(AQ^-1𝐀^𝑇+𝐏⁻¹)^-1f, by algebraic duality.

The proof of this last statement is a straightforward generalization of the previous situation involving the weighted observational errors and prior information about the unknown parameters. Using the Matrix Inversion Lemma (also called the Schur Identity), one can write 𝐱=(𝐀^𝑇𝐏𝐀+𝐐)⁻¹𝐀^𝑇𝐏𝐟=𝐐⁻¹𝐀^𝑇(𝐀𝐐⁻¹𝐀^𝑇+𝐏⁻¹)⁻¹𝐟, (2.16) when the inverses of the weight matrices exist. Notice that the first RHS expression is in terms of the weight matrices while the second is in terms of their inverses, and that the preceding least-squares estimates are obtained with P = I and Q = 0 for a simple overdetermined system, and with Q^-1 = I and P^-1 = 0 for a simple underdetermined system.

There is also an extensive theory dealing with the cases of rank deficiency in the matrix A implying a singular matrix A^TA or AA^T in the above expressions. In such cases, special precautions are required to reduce the number of parameters to be estimated or constrain their estimation. In general, the singular value decomposition discussed in the next section provides the best general strategy for linear problems with rank deficiencies.

In geometrical terms, the least-squares approach corresponds to an approximation using a normal (i.e., orthogonal) projection, and as a BAE, it does not necessarily involve any statistical information. In other words, as shown explicitly above, in all cases of underdetermined, determined, and overdetermined situations, even with weights associated with the observations and parameters, the least-squares solution is simply a weighted average of the observations for each unknown parameter. This is the reason for considering the least-squares approach basically as a mathematical approximation procedure which turns out to be most appropriate for statistical applications (see, e.g., [2] for more general discussions).

However, when interpreting the measurements or observations with finite first and second moments as a Gaussian sample, the average and hence the least-squares estimate becomes the unbiased minimum-variance estimate or the MLE. This is because any sample with finite first and second moments may be identified with a Gaussian sample as the Gaussian or normal distribution is fully specified by the first two moments. This statistical interpretation of least-squares estimates is really useful for error analysis and reliability considerations as in addition to Gaussian implications for the first moment, the second moment information behaves as a Chi-Square (𝜒²) distribution. Gaussian statistics are widely used in the analysis of least-squares estimates largely because of the well-developed theory and wide-ranging practical experience.

The previously introduced weight matrices P and Q are usually interpreted in the statistical sense as inversely proportional to the covariance matrices of the measurements and unknown parameters, respectively. Using unit proportionality factors, these are explictly 𝐏={𝐄[(𝐞−𝐞₀)(𝐞−𝐞₀)^𝑇]}⁻¹,𝐐={𝐄[(𝐱−𝐱₀)(𝐱−𝐱₀)^𝑇]}⁻¹, (2.17) in which the zero subscript corresponds to the mean or expected value. Notice that 𝐞₀ = E[e] = 0 for unbiasedness while 𝐱₀ = E[𝐱] is not necessarily zero in general applications. Using the well-known covariance propagation law; that is, for any linear transformation z = Rx, one has for the corresponding second moment𝐄[(𝐳−𝐳_𝑜)(𝐳−𝐳_𝑜)^𝑇]=𝐑𝐄[(𝐱−𝐱₀)(𝐱−𝐱₀)^𝑇]𝐑^𝑇 (2.18) then the variance of the estimated parameters is readily obtained as 𝐄[(𝐱−𝐱₀)(𝐱−𝐱₀)^𝑇]=(𝐀^𝑇𝐏𝐀+𝐐)⁻¹=𝐐⁻¹−𝐐⁻¹𝐀^𝑇(𝐀𝐐⁻¹𝐀^𝑇+𝐏⁻¹)⁻¹𝐀𝐐⁻¹, (2.19) which again shows a formulation in terms of the weights P and Q, and a dual formulation in terms of their inverses. In practical applications, these two equivalent formulations can be exploited to minimize the computational efforts either in terms of weight (or information) matrices or in terms of covariance matrices. Notice that adding a diagonal term to an arbitrary matrix can be interpreted in several different ways to numerically stabilize the matrix inversion such as, for example, in ridge estimation [8], Tikhonov regularization [9], and variations thereof.

All nonzero covariance matrices and their inverses, the nonzero weight (or information) matrices, are symmetric and positive definite in least-squares estimation. The optimality of the estimates in the sense of minimum variances requires such symmetry and positive definiteness for all the nonzero weight and covariance matrices involved. In some applications, it may be required to control the dynamic range and spectral shape of the covariance of the estimation error and to that end, such methods as Covariance Shaping Least-Squares Estimation [10] can be used advantageously in practice.

For illustration purposes, consider the following situation. Given five measurements {(1,5),(2,7),(3,13),(4,8),(5,6)} with some a priori information about a desired quadratic regression model, that is, 𝑦=𝑎+𝑏𝑥+𝑐𝑥², with unknown parameters a, b, and c, the preceding least-squares formulations can be used to provide the estimates for the quadratic polynomial in the interior of the interval spanned by the observations used (see Table 1) for interpolation purposes. Notice that the mathematical model, namely, the quadratic polynomial, was assumed from the context of the experiment or exercise. In general, the model identification needs to be resolved and the situation will be briefly discussed in Section 6.

Table 1 

Estimation of quadratic polynomials for interpolation purposes.

3. Least-Squares Interpolation and Prediction

Simple linear interpolation consists in an arithmetic mean of 𝑁 quantities, that is, 1𝑢=𝑁^𝑁_𝑖=1𝑢_𝑖=^𝑁_𝑖=1𝜆_𝑖𝑢_𝑖 (3.1) with all the coefficients 𝜆_𝑖≡ 1/𝑁. In a more general context, the coefficients 𝜆_𝑖 would be estimated to optimize the interpolation or prediction. For instance writing 𝑓₀=𝑓(𝑥₀) and for the observations, 𝑓_𝑖 = 𝑓(𝑥_𝑖), i = 1,𝑁, one has the general interpolation formula 𝑓₀=^𝑁_𝑖=1𝑝_𝑖𝑓_𝑖,assuming^𝑁_𝑖=1𝑝_𝑖=1, (3.2) that is, with normalized weights 𝑝₁,…,𝑝_𝑁. In the presence of correlations between the observations, using the usual matrix notation, 𝑓₀=𝑝1⋯1¹¹⋯𝑝_1𝑁𝑝⋮⋮⋮_𝑁1⋯𝑝_𝑁𝑁𝑓₁⋮𝑓_𝑁=𝑐1⋯1¹¹⋯𝑐_1𝑁𝑐⋮⋮⋮_𝑁1⋯𝑐_𝑁𝑁⁻¹𝑓₁⋮𝑓_𝑁, (3.3) assuming the weight matrix P = [𝑝_𝑖𝑗] corresponds to the inverse of the covariance matrix C = [𝑐_𝑖𝑗].

Correspondingly, the least-squares prediction formula 𝑓₀=𝑐₀₁⋯𝑐_0𝑁𝑐₁₁⋯𝑐_1𝑁𝑐⋮⋮⋮_𝑁1⋯𝑐_𝑁𝑁⁻¹𝑓₁⋮𝑓_𝑁, (3.4) in which the quantities 𝑐_𝑖𝑗 are defined in terms of the correlations between 𝑓_𝑖 and 𝑓_𝑗, including 𝑓₀, assuming the expected value of 𝑓 to be zero, E[f(𝑥)] = 0. Notice that this expression is of the form of the least-squares solution to an underdetermined linear problem as discussed before. For prediction applications, the correlation terms are often modelled empirically using correlation functions of the separation distance 𝑑_𝑖𝑗=‖𝑥_𝑖−𝑥_𝑗‖,𝑖,𝑗=0,…,𝑁.

When E[𝑓(𝑥)] is an unknown constant 𝛼, the normal equations for the unknown coefficients 𝜆₁,…,𝜆_𝑁 and 𝛼 are written explicitly as 𝑐₁₁⋯𝑐_1𝑁1𝑐⋮⋮⋮⋮_𝑁1⋯𝑐_𝑁𝑁1𝜆1⋯10¹⋮𝜆_𝑁𝛼=𝑓₁⋮𝑓_𝑁0 (3.5) for some unknown Lagrange multiplier 𝛼. In general when E[𝑓(𝑥)] is modelled as a 𝑘th degree polynomial, these normal equations become 𝑐₁₁⋯𝑐_1𝑁1⋯𝑥^𝑘𝑐⋮⋮⋮⋮⋮⋮_𝑁1⋯𝑐_𝑁𝑁1𝑥^𝑘𝑥1⋯10⋯0⋮⋮⋮⋮⋮⋮_𝑘⋯𝑥^𝑘𝜆0⋯0¹⋮𝜆_𝑁𝛼₀⋮𝛼_𝑘=𝑓₁⋮𝑓_𝑁0⋮0 (3.6) for some unknown coefficients 𝛼₀,…,𝛼_𝑘.

Least-squares prediction problems can be classified differently depending on application and other considerations. In different contexts, deterministic and probabilistic interpretations are used and hence the inferences are different. Three methodologies are mentioned here.(a)Kriging methods with variograms or generalized covariance functions such as, Cov(𝑑) = −𝑑 or 𝑑³ for spatial distance d (see [11] for details).(b)Radial Basis Function (RBF) methods with empirical RBFs as weighing functions, such as, RBF(𝑟) = 𝑟² log 𝑟 or (𝑟²+𝑐²)^1/2 for radial distance 𝑟 and constant c (see, e.g., [12] for more details). (c)Least-Squares Collocation (LSC) methods with ordinary covariance functions C≡C_s + C_n with C_s denoting the signal part and C_n denoting the noise part of the covariance matrix C.

Furthermore, generalized covariance functions (that is with positive power spectra) include ordinary covariance functions and empirical RBFs that can often be interpreted as covariance or correlation functions. Notice that the nonsingularity of the normal equation matrices is always assumed to guarantee a solution without additional constraints.

4. Solution Methodology for Normal Equations

The normal equation matrices 𝐀^𝑇A or A𝐀^𝑇 and the like are positive definite symmetric matrices that lend themselves to 𝐋^𝑇L or L𝐋^𝑇 decompositions in terms of lower triangular matrices L. The best known decomposition algorithm for such a matrix is the Cholesky square-root algorithm which is usually applied simultaneously to the normal equation matrix C and right-hand side vector (or column matrix) F as follows: 𝑐₁₁𝑐₁₂⋯𝑐_1𝑁𝑓₁𝑐₂₁𝑐₂₂⋯𝑐_2𝑁𝑓₂𝑐⋮⋮⋮⋮⋮_𝑁1𝑐_𝑁2⋯𝑐_𝑁𝑁𝑓_𝑁𝑐Cholesky'ssquare-root−−−−−−−−−−−−−−→^′11𝑐^′₁₂⋯𝑐^′_1𝑁𝑓^′₁0𝑐^′₂₂⋯𝑐^′_2𝑁𝑓^′₂⋮⋮⋮⋮⋮00⋯𝑐^′_𝑁𝑁𝑓^′_𝑁. (4.1) The solution is then obtained as a back substitution in the resulting upper triangular system of equations. The computational effort for 𝑁 normal equations is approximately 1/3 of the effort 𝑂(𝑁³) required using inverse matrix strategies. Furthermore in the case of least-squares adjustments of blocks of stereomodels, photographs and networks of geodetic stations, the normal equation matrix is usually banded to 𝐵≪𝑁 and the Cholesky’s algorithm only requires 𝑂(𝑁𝐵²) in such cases. This is really advantageous and has been used in large geodetic network, photogrammetric block adjustments, and most other similar least-squares applications.

Furthermore, in the Cholesky square-root reduction of a normal equation matrix to an upper (or lower) triangular matrix, the numerical conditioning is usually monitored by the magnitude of the computed diagonal elements as these should remain positive for a positive definite symmetric matrix. In some applications, the procedure of monitoring the magnitude of the diagonal elements of the triangular matrix is used to decide on an optimal order for a polynomial model which may be in terms of complex variables. Other similar strategies for numerical analysis of least squares problems based on the triangular decomposition of the normal equations matrix will be mentioned below.

However as mentioned in the previous section, the least-squares prediction problems have a “normal” matrix of the general form 𝐃𝐂𝐃𝐅_𝑇𝟎𝟎 (4.2) with the matrix C symmetric and positive definite, D rectangular, and 0 denoting a zero matrix. Such a “normal” matrix is obviously nonpositive definite and hence is not really appropriate for the previous triangular matrix representation. However, the Cholesky’s square-root procedure can be applied to the first 𝑁 equations and then Givens rotations can be applied to transform the remaining 𝑀 rows into the upper triangular form. Givens rotations are applied to two row vectors at any one time to eliminate the first nonzero element of the second row vector thus transforming the system of equations into an upper triangular system for back substitution at any time. Explicit implementation details can be found in [13] and elsewhere. Graphically, the situation is as follows: 𝑐₁₁𝑐₁₂⋯𝑐_1𝑁𝑑₁₁⋯𝑑_1𝑀𝑓₁𝑐₂₁𝑐₂₂⋯𝑐_2𝑁𝑑₂₁⋯𝑑_2𝑀𝑓₂𝑐⋮⋮⋮⋮⋮⋮⋮⋮_𝑁1𝑐_𝑁2⋯𝑐_𝑁𝑁𝑑_𝑁1⋯𝑑_𝑁𝑀𝑓_𝑁𝑑₁₁𝑑₂₁⋯𝑑_𝑁1𝑑0⋯00⋮⋮⋮⋮⋮⋮⋮⋮_1𝑀𝑑_2𝑀⋯𝑑_𝑁𝑀𝑐0⋯00Cholesky'ssquare-root−−−−−−−−−−−−−−→^′11𝑐^′₁₂⋯𝑐^′_1𝑁𝑑^′₁₁⋯𝑑^′_1𝑀𝑓^′₁0𝑐^′₂₂⋯𝑐^′_2𝑁𝑑^′₂₁⋯𝑑^′_2𝑀𝑓^′₂⋮⋮⋮⋮⋮⋮⋮⋮00⋯𝑐^′_𝑁𝑁𝑑^′_𝑁1⋯𝑑^′_𝑁𝑀𝑓^′_𝑁𝑑₁₁𝑑₂₁⋯𝑑_𝑁1𝑑0⋯00⋮⋮⋮⋮⋮⋮⋮⋮_1𝑀𝑑_2𝑀⋯𝑑_𝑁𝑀𝑐0⋯00Givens'Rotations−−−−−−−−−−−−→^′′11𝑐^′′₁₂⋯𝑐^′′_1𝑁𝑑^′′₁₁⋯𝑑^′′_1𝑀𝑓₁^′′0𝑐^′′₂₂⋯𝑐^′′_2𝑁𝑑^′′₂₁⋯𝑑^′′_2𝑀𝑓₂^′′⋮⋮⋮⋮⋮⋮⋮⋮00⋯𝑐^′′_𝑁𝑁𝑑^′′_𝑁1⋯𝑑^′′_𝑁𝑀𝑓_𝑁^′′00⋯0𝑒^′′₁₁⋯𝑒^′′_1𝑀𝑔₁^′′⋮⋮⋮⋮⋮⋮⋮⋮00⋯00⋯𝑒^′′_𝑀𝑀𝑔_𝑀^′′. (4.3) Notice that Givens rotations could be applied to the full 𝑁+𝑀 equations but the preceding strategy is superior in terms of computational efficiency. Givens rotations have excellent numerical stability characteristics but often require slightly more computational efforts than the alternatives.

5. Singular Value Decomposition

For indepth analysis of least-squares results and other related applications, the Singular Value Decomposition (SVD) approach is essential in practice and a brief overview follows. Considering the preceding (rectangular 𝑀×𝑁) matrix A, its SVD gives 𝐀=𝐔Λ𝐕^𝑇, (5.1) where U is the matrix of (unit) eigenvectors of A𝐀^𝑇, V is the matrix of (unit) eigenvectors of 𝐀^𝑇A, and Λ is an 𝑀×𝑁 matrix with diagonal elements (called the singular values) equal to the square roots of the nonzero eigenvalues of 𝐀^𝑇A or A𝐀^𝑇. By substitution, one readily obtains 𝐀𝐀^𝑇=(𝐔Λ𝐕^𝑇)(𝐔Λ𝐕^𝑇)^𝑇=𝐔Λ𝐕^𝑇𝐕Λ^𝑇𝐔^𝑇=𝐔(ΛΛ^𝑇)𝐔^𝑇,𝐀^𝑇𝐀=(𝐔Λ𝐕^𝑇)^𝑇(𝐔Λ𝐕^𝑇)=𝐕Λ^𝑇𝐔^𝑇𝐔Λ𝐕^𝑇=𝐕(Λ^𝑇Λ)𝐕^𝑇, (5.2) where ΛΛ^𝑇 and Λ^𝑇Λ denote the diagonal matrices of the squares of the singular values of A and dimensions 𝑀×𝑀 and 𝑁×𝑁, respectively. The last step in both derivations follows from the orthogonality of the (unit) eigenvectors in U and V, respectively. Their inverses are, respectively, (𝐀𝐀^𝑇)⁻¹=𝐔(ΛΛ^𝑇)⁻¹𝐔^𝑇,(𝐀^𝑇𝐀)⁻¹=𝐕(Λ^𝑇Λ)⁻¹𝐕^𝑇, (5.3) as matrix inversion does not change the eigenvectors in the SVD of a symmetric matrix.

The previous least-squares solution for the overdetermined system A𝐱 = f is simply 𝐱=(𝐀^𝑇𝐀)^−𝟏𝐀^𝑇𝐟=𝐕(Λ^𝑇Λ)^−𝟏𝐕^𝑇𝐕Λ^𝑇𝐔^𝑇𝐟=𝐕(Λ^𝑇Λ)^−𝟏Λ^𝑇𝐔^𝑇𝐟=𝐕Λ⁻𝐔^𝑇𝐟 (5.4) with the special notation Λ⁻ = (Λ^𝑇Λ)⁻¹Λ^𝑇 and for the underdetermined case 𝐱=𝐀^𝑇(𝐀𝐀^𝑇)⁻¹𝐟=𝐕Λ^𝑇𝐔^𝑇𝐔(ΛΛ^𝑇)⁻¹𝐔^𝑇𝐟=𝐕Λ^𝑇(ΛΛ^𝑇)⁻¹𝐔^𝑇𝐟=𝐕Λ⁺𝐔^𝑇𝐟 (5.5) with the special notation Λ⁺ = Λ^𝑇(ΛΛ^𝑇)⁻¹. These Λ⁻ and Λ⁺ are usually called generalized inverses of Λ.

From a computational perspective, for an overdetermined system with 𝑀≫𝑁, it may be more efficient to first perform a QR factorization of A with Q as an 𝑀×𝑁 matrix with orthogonal columns and an upper triangular matrix R of order 𝑁, and then compute the SVD of R, since if A = QR and R = UΛ𝐕^𝑇, then the SVD of A is given by A = (QU) Λ𝐕^𝑇. Similarly, for an underdetermined system with 𝑀≪𝑁, it may be more efficient to first perform an LQ factorization of A with a lower triangular matrix L of order 𝑀 and an 𝑀×𝑁 matrix Q with orthogonal rows, and then compute the SVD of L, since if A = LQ and L = UΛ𝐕^𝑇, then the SVD of A is given by A = UΛ(𝐐^𝑇𝐕)^𝑇. The SVD approach is often used in spectral analysis and computing a minimum norm solution for (possibly) rank-deficient linear least-squares and related problems. More discussion of the computational aspects can be found, for example, in [14].

In practice, the SVD of a matrix has been described as “one of the most elegant algorithms in numerical algebra for exposing quantitative information about the structure of a system of linear equations” [15]. In current data assimilation and prediction research using spatiotemporal processes such as in global change and other environmental applications, the SVD approach has become very important for at least three problem areas.

First, considering sequences of discrete data 𝐱=(𝑥₁,𝑥₂,𝑥₃,…,𝑥_𝑁) with zero mean for simplicity associated with discrete times 𝑡₁,𝑡₂,…,𝑡_𝑀, and written in matrix form as 𝑥𝐗=¹¹𝑥²₁…𝑥^𝑀₁𝑥¹₂𝑥²₂…𝑥^𝑀₂𝑥⋮⋮⋮⋮_1𝑁𝑥²_𝑁…𝑥^𝑀_𝑁 (5.6) in which the superscripts correspond to the times 𝑡₁,𝑡₂,…,𝑡_𝑀. Such a convention with columns corresponding to the data sequences at discrete times is quite common in environmental applications. Then a SVD of this data matrix X yields X = UΛ𝐕^𝑇, as discussed above. The columns of U are the Empirical Orthogonal Functions (EOFs) for the data matrix X while the columns of V are the corresponding principal components. The data transformation 𝐔^𝑇x or more generally U*𝐱 for a (complex) data vector 𝐱 is usually called a Karhunen-Loève transformation. Such resulting data sequences z = 𝐔^𝑇x are easily seen to be uncorrelated as 𝐄[𝐳𝐳^𝑇]=𝐄[𝐔^𝑇𝐱𝐱^𝑇𝐔]=𝐔^𝑇𝐄[𝐱𝐱^𝑇]𝐔=𝐔^𝑇𝐂_𝐱𝐔=𝐔^𝑇𝐔Λ𝐔^𝑇𝐔=Λ (5.7) which is most useful in practical applications. It is also important to notice that the EOFs can be described as eigenvectors of the corresponding covariance matrix of the available data. These are often called normal modes of the measured spatiotemporal process [16]. Since the (power) spectrum of a data sequence is well known to correspond to the spectrum of its (auto) covariance matrix, such normal modes have interesting interpretations in the context of dynamical systems driven by noise (e.g., [17]). Further discussions can be found in [18] and the references therein. An example of simulated application is shown in Figure 1(a) with a spatial pattern of a box followed by two cones over a sinusoidal path, with the resulting time series in Figure 1(b) of 20 occurrences of the pattern of 36 observations. Figure 1(c) shows the first modal (spatial) pattern and the corresponding first singular values which are identical to the input information except for the different scaling in amplitude and sign convention of Mathematica 7 [19]. Analogous simulations can readily be done in two dimensions with various patterns (see [20]). Such simulations show the potential of this methodology in the analysis of environmental and other geophysical time series.

(a)

(b)

(c)

(d)

(e)

Figure 1 

EOF application to template times sinusoidal trend function and extraction thereof.

Second, in numerical conditioning analysis for any linear algebraic system of equations, the singular values of the matrix give most relevant information about the propagation of numerical errors from observations to estimated unknown parameters. For instance, considering some symmetric and positive definite matrix B, consider the linear algebraic system 𝐁𝐮=𝐯, (5.8) and for some small perturbations 𝜹u and 𝜹v such that 𝐁𝜹𝐮=𝜹𝐯, (5.9) then using the spectral norm, it is well known that ‖𝜹𝐮‖‖𝐮‖≤‖𝐁‖⋅‖𝐁^−𝟏‖⋅‖𝜹𝐯‖=𝜆‖𝐯‖^max(𝐁)𝜆_min⋅(𝐁)‖𝜹𝐯‖‖𝐯‖, (5.10) in which 𝜆_max(𝐁)/𝜆_min(𝐁)≡𝜅(B) is the condition number of the matrix B in terms of its maximum and minimum eigenvalues, 𝜆_max(B) and 𝜆_min(B), respectively. This provides a powerful tool for the analysis of relative changes in the unknown parameters implied by some relative perturbations in the observations.

Third, as the nonzero eigenvalues of 𝐀^𝑇A and A𝐀^𝑇 are identical, computations, and numerical analysis in filtering and smoothing can take advantage of this fact in using the smaller matrix in the least-squares computations. In several areas of environmental research, enormous quantities of data and complex mathematical models lead to very large normal equation systems that require much computational efforts. Substantial reductions in the Kalman filtering and data assimilation computations are implied by the proper choice of covariance/information formulations but for indepth error analysis of the results, SVD-based techniques become critical (see, e.g., [21]).

6. Model Identification and Reliability Considerations

As previously mentioned, alternatives to least squares (𝐿₂), such as the least magnitude (𝐿₁) and least maximum (𝐿_∞), have advantages and disadvantages. First, the contributions of errors and especially outliers will increase in the estimation procedure from 𝐿₁ to 𝐿_∞. Therefore, in the absence of outliers, 𝐿_∞ would be best while in the presence of large errors, 𝐿₁ or robust estimation is more appropriate. However considering the linear normal equations with least squares, least squares are selected for most practical applications. Outliers are a well-known problem with least squares and much literature exists to mitigate the implications.

Furthermore, in numerous application contexts, assuming an algebraic formulation, one has to decide on degrees and orders in regression modeling such as in curve fitting and spectrum estimation. Considering a simple set of measurements or observations {(𝑡_𝑖,𝑓(𝑡_𝑖));𝑖=1,2,3,…,𝑁} to be modelled using an algebraic polynomial of the form 𝑓(𝑡)=𝑎₀+𝑎₁𝑡+𝑎₂𝑡²+⋯. If the degree of the polynomial is unknown, then the least-squares approach to estimating 𝑎₀,𝑎₁,𝑎₂,… can always achieve a perfect fit to the measurements by selecting the polynomial degree equal to 𝑁−1. Actually, the variance of the residuals will decrease with higher and higher degrees to become zero for degree 𝑁−1. Hence, the least-squares approach cannot be used to decide on the polynomial degree for such regression applications. Furthermore, the approach alone can hardly be used to decide on some other possible mathematical models such as 𝑓(𝑡)=𝑏₀+𝑏₁cos(𝑡)+𝑏₂sin(𝑡)+⋯ as additional information is necessary for model identification. Such model identification problems have been studied extensively by Akaike [22, 23] and others (see [24] for further discussions and references).

For example, the sample measurements {(1,5),(2,7),(3,13),(4,8),(5,6)} with an algebraic polynomial of the form 𝑦=𝑎+𝑏𝑥+𝑐𝑥²+⋯ would lead to an exact fit for degree 4. However, such high degree would likely be unacceptable because of the oscillations between the data points which would imply uncertainty in any prediction. Considering the least-squares estimates for degrees 0, 1, 2, 3, and 4, 𝑞₀𝑞(𝑥)=7.8,₁𝑞(𝑥)=6.9+0.3𝑥,₂(𝑥)=−2.6+8.44286𝑥−1.35714𝑥²,𝑞₃(𝑥)=−1.2+6.47619𝑥−0.60714𝑥²−0.08333𝑥³,𝑞₄(𝑥)=51.0−91.9167𝑥+59.2917𝑥²−14.5833𝑥³+1.20833𝑥⁴ (6.1) with corresponding error variances 𝜎²₀=9.7,𝜎²₁=9.475,𝜎²₂=3.02857,𝜎²₃=3.00357,𝜎²₄=0.0. (6.2) However, considering the (normalized) Akaike Information Criterion (AIC), defined as AIC=log𝜎²+2𝑀+3𝑁 (6.3) for approximating 𝑀 (model) parameters using 𝑁 measurements, one obtains 3.27213,3.64866,2.90809,3.2998, (6.4) corresponding to each degree 0, 1, 2, and 3, respectively, implying an optimal degree 2 for the modeling as 2.90809 is the minimum AIC value. Other examples can be found in [24] and the references therein.

In this context of least squares, the assumption of finite first and second moments with Gaussian statistics interpretation have implications in terms of expected moments for the estimated parameters. Essentially, assuming a given mathematical model, the given variances for the measurements and/or observations and a priori variances for the unknown parameters can be propagated using the variance propagation law into the estimated parameters and interpreted at some confidence level such as 95%. This is the familiar approach in geomatics with error ellipses reflecting the accuracy of measurements and/or observations and the geometrical strength of a network in positioning.

For example, given positional information at two discrete points 𝑃₁ and 𝑃₂ with a variance 𝜎², the mid point located by the arithmetic average of the coordinates of 𝑃₁ and 𝑃₂ has a predicted variance 𝜎²/2 when a linear model is known for any intermediate point. However, when the location and/or definition of the mid point is ambiguous or unknown, such as along some fuzzy line, its predicted uncertainty is likely to be much more than 𝜎²/2. In other words, the uncertainty in any estimation results is attributable to uncertainty in the assumed mathematical model and in the observational information used.

Another illustrative example is in the prediction of some quantity g, such as gravity, as a function of known data at the discrete points 𝑃₁ and 𝑃₂ with observational variance 𝜎². At the mid point between 𝑃₁ and 𝑃₂, the average 𝑔 value, that is, [𝑔(𝑃₁)+𝑔(𝑃₂)]/2, is usually an adequate prediction of the 𝑔 value there but its variance is likely to be greater than 𝜎²/2. Otherwise, why bother with measurements! It should be noticed that in nonlinear and/or non-Gaussian situations, the error propagation is much more complex, and only numerical Monte Carlo simulations offer a general strategy for uncertainty modeling (see [25] for further details and references).

7. Concluding Remarks

Least squares are ubiquitous in applied science and engineering data processing. From a mathematical perspective, a least-squares estimate is a (weighted) mean solution which may be interpreted differently depending on the application context. Furthermore, any linear finite problem, even an ill-posed one, has a unique solution in the “average sense”. Such a solution is a BAE using a minimum quadratic norm or minimum variance, with the flexibility of possible statistical interpretation as MLE for optimal and reliable predictions.

Advantages of the least-squares approach are essentially in the simple assumptions (i.e., finite first and second moments), the unique estimates from linear normal equations, with excellent computational and applicability characteristics. Disadvantages of the least-squares approach are mainly in terms of oversmoothing properties (such as in curve or surface fitting) and relative overemphasis of outlier observations or measurements.

In terms of numerical computations, the least-squares approach has excellent characteristics in terms of stability and efficiency. This is best seen using the SVD approach for any indepth analysis of least-squares results. The readjustment of the geodetic networks in the North American Datum of 1983 has demonstrated that nearly one million unknowns can be handled reliably with only 32 bit arithmetic on conventional computer platforms (see [26, 27]). A better numerical approach would be difficult to find!

Furthermore, it is also important to emphasize that least squares are not appropriate for all types of estimation problems as there are numerous application contexts where a best observation or measurement value needs to be selected among the available ones (as the most frequent value or the one with minimum error). In other application contexts when dealing with observations or measurements likely to be affected by outliers, a more robust estimate such as a median value or 𝐿₁ estimate may be preferable. No single estimation method can be considered best or optimal for all applications as data characteristics and desired estimates need to be considered.

Finally, some areas of current research and development in least squares and computational analysis include multiresolution analysis and synthesis, data regularization and fusion, EOFs of multidimensional time sequences, RBFs, and related techniques for optimal data assimilation and prediction, especially in spatiotemporal processes. With scientific data generally considered to be increasing faster than computational power, real challenges in the analysis of current observations and measurements abound and strategies have to become more sophisticated.

Acknowledgments

The author would like to acknowledge the sponsorship of the Natural Science and Engineering Research Council in the form of a Research Grant on Computational Tools for the Geosciences. Comments and suggestions from the reviewers are also acknowledged with gratitude.