Linear least squares approximation

Least squares method is usually mentioned in two contexts. First, its application is widely known in regression analysis. It is used to construct models on the basis of noisy experimental data. At that, besides constructing a model, it usually performed an estimate of the error with which parameters were calculated. Some other problems could be solved as well. Secondly, the least-squares method is often applied as an approximation method without any connection to statistics. Here this technique is considered as an approximation method.

General least squares method

When approximating by using the least squares method, the function to be approximated is given by N points (x_i , y_i ). The approximation function is built as linear combination of basis functions F_j  (the number of functions M is usually less than N).

The coefficients cj are chosen in such a way as to minimize the sum of the squared deviations of the approximating function from their given values

Sometimes, if the points have weights, the items are multiplied by the corresponding w_i 

This problem is called weighted approximation by linear least squares. It should be noted that such an approximation method (linear combination of basis functions) and such an estimate function (sum of squared deviations) is not the only one which someone can choose. Basis functions can enter into g nonlinearly, and estimate function E could be changed to a maximum of deviation or any other estimate function. However, it is this approximation method and this estimate function that let us get the best c_j  in finite number of operations reducing a problem to set of linear equations solution.

Methods for finding the coefficients

Before considering the methods for finding the coefficients c_j , let's introduce the following designations:

The first way of finding optimal c_j  consists of solving a set of equations (so-called normal equations method). Taking into account problem properties, this condition is not only a necessary, but also a sufficient condition of minimum. Moreover, this set of equations is linear in c_j  and it is written as (A ^TA)c - A ^Tb = 0. The main disadvantage of this method is that the set of equations found could be ill-conditioned. The condition number of matrix A could be very big (e.g., when the polynomial basis is used), its square is even greater. That's why the normal equations method is not used in practice.

The second method is based on the fact that problem of minimizing E could be written in matrix form as finding . Thus, the problem is reduced to the problem of finding a pseudosolution of a set of linear equations. Such a solution could be solved by using the singular value decomposition (SVD). This method is more natural though it requires more complex program codes to implement it. The advantage of this method is that there is no squaring of the condition number of a basis. We solve a system whose condition number is equal to the basis condition number, which lets us get rather accurate values even if the normal equations method results in a set with a singular coefficients matrix. It should be noted that although we could use the general algorithm of the singular value decomposition, we used the more compact and faster version of this algorithm which considers the fact that we need not get SVD of matrix A but a product of the decomposition components and vector b.

There is another method based on the QR decomposition. It works a little faster than the preceding method. At first, matrix A is represented by a product of orthogonal matrix Q and upper triangular matrix R. Then, a set of equations Rx = Q ^Tb is solved. The advantage of this method is its speed, the disadvantage is that it is applicable only for non-degenerate problems (i.e. having one, and one only solution). The approximation problem is non-degenerate if the two conditions are met. First, the set of basis functions should be linear-independent. Secondly, the number of points should be not less than the number of the basis functions. As a rule, practical problems satisfy these conditions. But sometimes the algorithm on the basis of QR decomposition is not applicable because one of these conditions is not met. Moreover, the number of rows in A is usually much greater than the number of columns, so correctly implemented singular value decomposition of such a matrix is almost as fast as QR decomposition. Thus, it is more practical to use the second method (SVD).

Linear approximation

f₀  = 1 è f₁  = x are used as a set of basis functions. Linear combinations of the functions of this basis allow us to get an arbitrary line on a plane. The peculiarity of this basis is that the normal equations method is applicable for this case. The basis is small, well-conditioned, so the approximation problem is easily reduced to solving a 2x2 set of linear equations. To solve this set rotation method modification is used.

This algorithm is implemented in the BuildLinearLeastSquares subroutine.

Polynomial basis

It is a common error to use basis functions f_i  = x ⁱ for polynomial approximation. This basis is the most natural. But, although this basis is easy to work with, it is very ill-conditioned. There could be problems even if the number of basis functions is not very big (about 10). This problem was solved by choosing the Chebyshev polynomials basis instead of powers of x. Chebyshev polynomials could be represented as powers of x and vice versa linearly, so these bases are equivalent. However, the Chebyshev polynomial basis is much better conditioned. This allows us to approximate by polynomials of almost any degree.

The algorithm works as follows. The given points are scaled and reduced to [-1,1] interval, then the Chebyshev polynomials basis is built on this interval. After the least squares approximation, the user gets a set of Chebyshev polynomials coefficients. This problem is solved by the BuildChebyshevLeastSquares subroutine. These coefficients could be passed into the CalculateChebyshevLeastSquares subroutine which calculates the value of the approximating function.

If the problem requires using power basis, you could use the BuildPolynomialLeastSquares subroutine. This subroutine solves a problem by using the Chebyshev polynomials, and then the result is transformed into power basis. This enables to keep the accuracy of the calculations, because all intermediate computations are performed in a good-conditioned basis. However, the Chebyshev polynomial basis is preferable.

Sometimes it is required to solve a constrained problem. For example, we could know that the function to be approximated equals 1 in x = 0. If we have a rather big number of points in the neighborhood of x = 0, the approximating function will pass through the neighborhood of (0,1). However, it is likely not to pass through the point itself. In many cases it is not required to meet boundary conditions exactly. But there are problems which require exact passing through given points or exact equality of a derivative to the given values. At that, special demands are made only for some of x_k . In other points we just minimize the sum of the squared deviations.

In this case we have an approximation problem with the following conditions:

or

Such problems are solved by using a special algorithm implemented in the BuildChebyshevLeastSquaresConstrained subroutine. At present this subroutine can process an arbitrary number of constraints (but strictly less than M) for values of function and its first derivative. Higher derivatives constraints are not supported.

Cubic spline basis

Cubic splines are one more possible set of basis functions. A set of splines with common nodes forms a linear space, so the linear least squares approximation is applicable. Splines which are chosen as basis functions should satisfy the following conditions:

For simplicity purposes, we assume that the approximation is built for [0,M-1] with an equidistant grid. Of course, in practice the interval for which to perform the approximation can be arbitrary (but it still will be an equidistant grid). On the graph below you can find an example of such a set of basis functions.

The approximation function is built by the BuildSplineLeastSquares subroutine. It returns a coefficients array by specifying a cubic spline. This array is passed to a SplineInterpolation subroutine (see description of a cubic spline interpolation module), which calculates the value of the approximating spline in the given point.

Arbitrary basis

We've considered three sets of basis functions which can be used for approximation: lines, polynomials and cubic splines. However, these sets are not the only ones to be used for approximation. Of course, we can't foresee all possible situations and all function sets. Therefore, any program package for approximation should contain a subroutine which approximates by using an arbitrary set of basis functions. This problem is solved by using the BuildGeneralLeastSquares subroutine.

Since the set of basis functions is unknown beforehand, you have to specify the values of the basis functions in points x_i . This information is passed as FMatrix matrix which contains the value of i-th function in j-th point. Also, you have to pass arrays y è w containing ordinates of given points and their weights.

It should be noted that abscissas x_i  are not passed into the subroutine because these values are used only as basis functions arguments, and this information is already contained in FMatrix. Therefore, the approximation is independent of dimension of space. You only have to number the points, calculate the function values and pass this information into the algorithm. This makes the subroutine applicable to both one-dimensional and multidimensional problems.

Source codes