LEAST-SQUARES AND CHI-SQUARE FOR THE BUDDING AFICIONADO: ART AND PRACTICE

\begin{matrix}c \\ ◯ \end{matrix}Carl Heiles March 21, 2007

0. USING SINGULAR VALUE DECOMPOSITION (SVD)

Occasionally, a normal-equation matrix [α] = X^{T} · X is degenerate, or at least sufficiently ill-posed that inverting it using standard matrix inversion doesn’t work. In its invert function, IDL even provides a keyword called status to check on this (although I find that it is not perfectly reliable; the best indicator of reliability is to check that the matrix product [α^{−1}] · [α] = I, the unitary matrix). In these cases, Singular Value Decomposition (SVD) comes to the rescue.

First, we reiterate the least-squares problem. Least squares begins with equations of condition (equation ??), which are expressed in matrix form as

X · a = y  (0.1)

In our treatments above we premultiply both sides by X^{T}, on the left generating the curvature matrix [α] to obtain the normal equations (equation ?? or its equivalent derived from equations ??)

[α]·a = X^{T} · y ,  (0.2)

for which the solution is, of course,

a = ([α]^{−1}·X^{T}) ·y  (0.3)

We need to find the inverse matrix [α]^{−1}. Above in this document we did this using simple matrix inversion, which doesn’t work if [α] is degenerate.

SVD provides a bombproof and interestingly informative way do least squares. Perhaps surprising, with SVD you don’t form normal equations. Rather, you solve for the coefficients directly. In essence, SVD provides the combination ([α]^{−1}·X^{T}) without taking any inverses. By itself, this doesn’t prevent blowup degenerate cases; however, SVD provides a straightforward way to eliminate the blowup and get reasonable solutions.

For a discussion of the details of SVD, see NR §2.6; for SVD applied to least squares, see NR §15.4; if you are rusty on matrix algebra, look at NR §11.0. Below, we provide a brief description of SVD. Implementing SVD in our least-squares solutions is trivially easy, and we provide the IDL prescription below in §0.3 and §0.4. Be sure to look at §0.2.1!!

0.1. Phenomenological description of SVD

The cornerstone of SVD is that our (or any) M × N (M rows, N columns) matrix X, where M ≥ N, can be expressed as a product of three matrices:

X = U·[w]·V^{T} ,  (0.4)

where

1. U is M × N, [w] is N × N and diagonal, and V is N × N; and

2. the columns of U and V are unit vectors that are orthonormal. Because V is square, its rows are also orthonormal so that V · V^{T} = 1. Recall that, for square orthonormal vectors, the transpose equals the inverse so V^{T} = V^{−1}.

So what?

Consider a single row m of the matrix X. This row has a single value of x_{m}, with an associated value of y_{m}. For this row, we can consider the entry of each column n of X to be a basis function f_{n}(x_{m}), evaluated at x_{m}. The N columns contain N basis functions. And, of course, we have N coefficients in the vector a. These N basis functions and associated coefficients represent the entire set of M measurements y_{m} taken at positions x_{m}.

Now consider a single column n of the matrix X. This column consists of the entire set of M values of f_{n}(x_{m}) for a single value of n. We can regard the M values of x_{m} to be a vector of length M. Similarly, we can regard each column n of X to be an M-length vector consisting of the elements f_{n}(x_{m}).

Finally, consider the set of N columns of the matrix X. Each column n is an M-element vector with elements f_{n}(x_{m}). Now, there is no reason for these N column vectors of X to have any special property, such as being orthogonal. In fact, consider a typical least-squares polynomial fit; this basis set is certainly not orthogonal, as we illustrated in our simple numerical example of §??. What SVD does is to replace the set of N original nonorthogonal vectors X by an orthogonal basis set that consists of the N rows of V—in fact, they are not only orthogonal, but orthonormal.

0.2. Using SVD for Least Squares

Some of the original nonorthogonal column vectors of X might not only be nonorthogonal, they might be degenerate! This happens if the particular set of the M values x_{m} make some of the original nonorthogonal vectors linear combinations of others. This leads to the unfortunate situation of the curvature matrix α not having an inverse. In such cases, SVD comes to the rescue.

With SVD you don’t form normal equations, so you don’t explicitly calculate [α]^{−1}. We apply SVD to X, as in equation 0.4:

X = U·[w]·V^{T} ,  (0.5)

NR §15.4 shows that the covariance matrix is given by

[α^{−1}] = V· [\frac{1}{w^{2}} ] · V^{T}  (0.6)

Given this, the fact that V^{T} · V = I, and the general matrix identity (P · Q)^{T} = Q^{T} · P^{T} it’s straightforward to show that

a = (V· [\frac{1}{w} ] ·U^{T}) ·y  (0.7a)

or

a = (V · [\frac{1}{w}]) · (U^{T}·y)  (0.7b)

Here in equation 0.7a, we see that (V· [\frac{1}{w}] ·U^{T}) is identical to ([α]^{−1}·X^{T}) in equation 0.3. In the rewritten equation 0.7b we see that we can regard the coefficient vector being defined by the orthogonal vectors of V. Moreover, the covariance matrix of equation 0.6 is diagonal, so that the principal axes of the χ^{2} ellipse are defined by the orthonormal column vectors of V with lengths proportional to [\frac{1}{w^{2}} ]. See NR Figure 15.6.5 and the associated discussion.

The N orthonormal vectors in V define an N-dimensional space for a with N orthogonal directions. Suppose that a particular value w_{n} is small. In equation 0.5, this means that the associated orthonormal vector v_{n}—i.e., the associated direction n—in V is not well-represented by the set of original nonorthogonal vectors in X. This, in turn, means that you can’t represent that direction in a of equation 0.7a without amplifying that orthonormal vector by a large factor; these amplification factors ∝ w_{n}^{−1}. This is an unsatisfactory situation because it means some combinations of the original m measurements are highly weighted. As w_{n} → 0 this situation becomes not only unsatisfactory, but numerically impossible.

Consider the limiting case, w_{n} = 0. In this case, the original x_{m} values do not have any projection along associated orthonormal vector v_{n} in V. (These v_{n} are called “null” orthonormal vectors.) In equation 0.7a, you can add any multiple of the null v_{n} to the solution for a and it won’t change a at all (!) because it has absolutely no effect on the fit to the data (because of the particular set of values x_{m}).

What multiple is appropriate? Common sense says that, because these V vectors have no meaning for the solution, the multiple should be zero. So in equation 0.7a, instead of trying, in vain, to include this null vector v_{n} by using a huge multiple, you toss in the towel and eliminate it altogether by replacing its corresponding w_{n}^{−1} = ∞ by w_{n}^{−1} = 0. So we have the rule: wherever w_{n} = 0 (or sufficiently small), replace w_{n}^{−1} by 0! This replacement provides the minimum length for x and thereby constitutes the least-squares solution.

0.2.1. Important Conclusion for Least Squares!!!

Suppose you have degeneracy, or near-degeneracy. This, in turn, means that the formulation of the least-squares model is faulty: some of the original basis functions represent (or nearly) represent the same physical quantity, or at least one of the functions is (nearly) a linear combination of others. Sometimes you can discover the problem and fix it by imposing additional constraints, or by finding two unknown coefficients that are nearly identical. If so, you can reformulate the model and try again.

Even if you can’t discover the root cause(s) and remove the degeneracy, the SVD solution allows you to bypass problems associated with degeneracy and provides reasonable best-fit parameter values.

0.2.2. How Small is “Small”?

When looking at w_{n}, just exactly how small is “small”?

0.2.3. Strictly Speaking…

Strictly speaking, this is governed by numerical machine accuracy. NR suggests measuring this by the ratio of the largest w_{n} to the smallest. This ratio is the condition number. Single precision 32-bit floats have accuracy ∼ 10^{−6}, so if the condition number is larger than 10^{6} you certainly have problems.

0.2.4. Practically Speaking…

Practically speaking it’s not machine accuracy that counts. Rather, it’s the accuracy of your data and the ability of those uncertain data to define the parameters. In any real problem, plot the vector w. If the values span a large range, then the data with small w_{n} might not be defining the associated vectors well enough. If so, then practically speaking, these vectors are degenerate.

How to tell what constitutes “small w_{n}”? As far as I know, and as far NR says, it’s an art.

0.3. Doing SVD in IDL

Doing SVD in IDL is easy. Ida’s la_svd procedure^{1} provides the SVD decomposition. Thus for equation 0.4, the IDL SVD decomposition is given by

la svd, X, w, U, V ;  (0.8)

the notation is identical to that in equation 0.4.

0.4. My routine lsfit svd

I’ve written an IDL routine lsfit_svd that makes the above process of dealing with the weights easy. When used without additional inputs, it returns the standard least-squares results such as the derived coefficients and the covariance matrix; it also returns [w], U, and V. It allows you to input those matrices and, also, the [\frac{1}{w}] matrix so that you can tailor the weights to your heart’s content.

1Old versions of IDL (before 5.6) had a significantly worse algorithm called svdc. Don’t use this unless you have to.