Given an MM by NN real matrix AA and an MM by 1 vector bb, find the NN by 1 vector xx when

or, using matrix notation,

If bb does not lie in the range space of AA (the space spanned by the columns of AA), there is no exact solution to Equation 2, therefore, an approximation problem can be posed by minimizing an equation error defined by

A generalized solution (or an optimal approximate solution) to Equation 2 is usually considered to be an xx that minimizes some norm of εε. If that problem does not have a unique solution, further conditions, such as also minimizing the norm of xx, are imposed. The l2l2 or root-mean-squared error or Euclidean norm is εT*εεT*ε and minimization sometimes has an analytical solution. Minimization of other norms such as l∞l∞ (Chebyshev) or l1l1 require iterative solutions. The general lplp norm is defined as qq where

for 1<p<∞1<p<∞ and a “pseudonorm" (not convex) for 0<p<10<p<1. These can sometimes be evaluated using IRLS (iterative reweighted least squares) algorithms [3], [5], [34], [16], [11].

If there is a non-zero solution of the homogeneous equation

then Equation 2 has infinitely many generalized solutions in the sense that any particular solution of Equation 2 plus an arbitrary scalar times any non-zero solution of Equation 5 will have the same error in Equation 3 and, therefore, is also a generalized solution. The number of families of solutions is the dimension of the null space of AA.

This is analogous to the classical solution of linear, constant coefficient differential equations where the total solution consists of a particular solution plus arbitrary constants times the solutions to the homogeneous equation. The constants are determined from the initial (or other) conditions of the solution to the differential equation.

Examination of the basic problem shows there are ten cases [19] listed in Figure 1 to be considered. These depend on the shape of the MM by NN real matrix AA, the rank rr of AA, and whether bb is in the span of the columns of AA.

Figure 1. Ten Cases for the Pseudoinverse.

Here we have:

This is a setting for frames and sparse representations.

In case 1a and 3a, bb is necessarily in the span of AA. In addition to these classifications, the possible orthogonality of the columns or rows of the matrices gives special characteristics.

Case 1: Here we see a 3 x 3 square matrix which is an example of case 1 in Figure 1 and 2.

If the matrix has rank 3, then the bb vector will necessarily be in the space spanned by the columns of AA which puts it in case 1a. This can be solved for xx by inverting AA or using some more robust method. If the matrix has rank 1 or 2, the bb may or may not lie in the spanned subspace, so the classification will be 1b or 1c and minimization of ||x||22||x||22 yields a unique solution.

Case 2: If AA is 4 x 3, then we have more equations than unknowns or the overspecified or overdetermined case.

If this matrix has the maximum rank of 3, then we have case 2a or 2b depending on whether bb is in the span of AA or not. In either case, a unique solution xx exists which can be found by Equation 15 or Equation 21. For case 2a, we have a single exact solution with no equation error, ϵ=0ϵ=0 just as case 1a. For case 2b, we have a single optimal approximate solution with the least possible equation error. If the matrix has rank 1 or 2, the classification will be 2c or 2d and minimization of ||x||22||x||22 yelds a unique solution.

Case 3: If AA is 3 x 4, then we have more unknowns than equations or the underspecified case.

If this matrix has the maximum rank of 3, then we have case 3a and bb must be in the span of AA . For this case, many exact solutions xx exist, all having zero equation error and a single one can be found with minimum solution norm ||x||||x|| using Equation 17 or Equation 22. If the matrix has rank 1 or 2, the classification will be 3b or 3c.

There are several assumptions or side conditions that could be used in order to define a useful unique solution of Equation 2. The side conditions used to define the Moore-Penrose pseudo-inverse are that the l2l2 norm squared of the equation error εε be minimized and, if there is ambiguity (several solutions with the same minimum error), the l2l2 norm squared of xx also be minimized. A useful alternative to minimizing the norm of xx is to require certain entries in xx to be zero (sparse) or fixed to some non-zero value (equality constraints).

In using sparsity in posing a signal processing problem (e.g. compressive sensing), an l1l1 norm can be used (or even an l0l0 “pseudo norm”) to obtain solutions with zero components if possible [12], [31].

In addition to using side conditions to achieve a unique solution, side conditions are sometimes part of the original problem. One interesting case requires that certain of the equations be satisfied with no error and the approximation be achieved with the remaining equations.

If the l2l2 norm is used, a unique generalized solution to Equation 2 always exists such that the norm squared of the equation error εT*εεT*ε and the norm squared of the solution xT*xxT*x are both minimized. This solution is denoted by

where A+A+ is called the Moore-Penrose inverse [1] of AA (and is also called the generalized inverse [6] and the pseudoinverse [1])

Roger Penrose [28] showed that for all AA, there exists a unique A+A+ satisfying the four conditions:

There is a large literature on this problem. Five useful books are [19], [1], [6], [9], [29]. The Moore-Penrose pseudo-inverse can be calculated in Matlab [22] by the pinv(A,tol) function which uses a singular value decomposition (SVD) to calculate it. There are a variety of other numerical methods given in the above references where each has some advantages and some disadvantages.

For cases 2a and 2b in Figure 1, the following NN by NN system of equations called the normal equations[1], [19] have a unique minimum squared equation error solution (minimum ϵTϵϵTϵ). Here we have the over specified case with more equations than unknowns. A derivation is outlined in "Derivations", equation Equation 28 below.

The solution to this equation is often used in least squares approximation problems. For these two cases ATAATA is non-singular and the NN by MM pseudo-inverse is simply,

A more general problem can be solved by minimizing the weighted equation error, ϵTWTWϵϵTWTWϵ where WW is a positive semi-definite diagonal matrix of the error weights. The solution to that problem [6] is

For the case 3a in Figure 1 with more unknowns than equations, AATAAT is non-singular and has a unique minimum norm solution, ||x||||x||. The NN by MM pseudoinverse is simply,

with the formula for the minimum weighted solution norm ||x||||x|| is

For these three cases, either Equation 15 or Equation 17 can be directly calculated, but not both. However, they are equal so you simply use the one with the non-singular matrix to be inverted. The equality can be shown from an equivalent definition [1] of the pseudo-inverse given in terms of a limit by

For the other 6 cases, SVD or other approaches must be used. Some properties [1], [9] are:

It is informative to consider the range and null spaces [9] of AA and A+A+

The four Penrose equations in Equation 11 are remarkable in defining a unique pseudoinverse for any A with any shape, any rank, for any of the ten cases listed in Figure 1. However, only four cases of the ten have analytical solutions (actually, all do if you use SVD).

Figure 2. Four Cases with Analytical Solutions

Fortunately, most practical cases are one of these four but even then, it is generally faster and less error prone to use special techniques on the normal equations rather than directly calculating the inverse matrix. Note the matrices to be inverted above are all rr by rr (rr is the rank) and nonsingular. In the other six cases from the ten in Figure 1, these would be singular, so alternate methods such as SVD must be used [19], [1], [6].

In addition to these four cases with “analytical” solutions, we can pose a more general problem by asking for an optimal approximation with a weighted norm [6] to emphasize or de-emphasize certain components or range of equations.

Figure 3. Three Cases with Analytical Solutions and Weights

These solutions to the weighted approxomation problem are useful in their own right but also serve as the foundation to the Iterative Reweighted Least Squares (IRLS) algorithm developed in the next chapter.

A particularly useful application of the pseudo-inverse of a matrix is to various least squared error approximations [19], [7]. A geometric view of the derivation of the normal equations can be helpful. If bb does not lie in the range space of AA, an error vector is defined as the difference between AxAx and bb. A geometric picture of this vector makes it clear that for the length of εε to be minimum, it must be orthogonal to the space spanned by the columns of AA. This means that A*ε=0A*ε=0. If both sides of Equation 2 are multiplied by A*A*, it is easy to see that the normal equations of Equation 14 result in the error being orthogonal to the columns of AA and, therefore its being minimal length. If bb does lie in the range space of AA, the solution of the normal equations gives the exact solution of Equation 2 with no error.

For cases 1b, 1c, 2c, 2d, 3a, 3b, and 3c, the homogeneous equation Equation 5 has non-zero solutions. Any vector in the space spanned by these solutions (the null space of AA) does not contribute to the equation error εε defined in Equation 3 and, therefore, can be added to any particular generalized solution of Equation 2 to give a family of solutions with the same approximation error. If the dimension of the null space of AA is dd, it is possible to find a unique generalized solution of Equation 2 with dd zero elements. The non-unique solution for these seven cases can be written in the form [6].

where yy is an arbitrary vector. The first term is the minimum norm solution given by the Moore-Penrose pseudo-inverse A+A+ and the second is a contribution in the null space of AA. For the minimum ||x||||x||, the vector y=0y=0.

To deal with measurement error and data noise, a process called “regularization" is sometimes used [15], [7], [25].

The solution of the overdetermined simultaneous equations is generally a least squared error approximation problem. A particularly interesting and useful variation on this problem adds inequality and/or equality constraints. This formulation has proven very powerful in solving the constrained least squares approximation part of FIR filter design [32]. The equality constraints can be taken into account by using Lagrange multipliers and the inequality constraints can use the Kuhn-Tucker conditions [14], [33], [20]. The iterative reweighted least squares (IRLS) algorithm described in the next chapter can be modified to give results which are an optimal constrained least p-power solution [4], [8], [5].

There is remarkable structure and subtlety in the apparently simple problem of solving simultaneous equations and considerable insight can be gained from these finite dimensional problems. These notes have emphasized the l2l2 norm but some other such as l∞l∞ and l1l1 are also interesting. The use of sparsity [31] is particularly interesting as applied in Compressive Sensing [2], [13] and in the sparse FFT [18]. There are also interesting and important applications in infinite dimensions. One of particular interest is in signal analysis using wavelet basis functions [10]. The use of weighted error and weighted norm pseudoinverses provide a base for iterative reweighted least squares (IRLS) algorithms.