least square problem solution

[Edwin Robert Tisdale <E dot Robert dot Tisdale at jpl dot nasa dot gov>]

Jose Miguel Buenaposada Biencinto wrote:

> I have the solution of a least square problem:
>
>     x = inverse(transpose(A)*A)*transpose(A)*b
>
> where x is nx1, A is mxn and b is mx1.
> I know (due to the nature of the problem I have)
> that some of the rows in A are not very significant
> and I'd like to be able to choose the best set of m rows of A
> to solve the system and to compute the error
> you make by doing so.
> Is there a direct procedure to do so?
> Or it is just a combinatorial problem (bad for me)?

Apparently, you are trying to solve
an overdetermined system of equations

    b = Ax

for x.  One way to do this is to pre-multiply both sides
of the equation by the transpose A^T

    (A^T)b = (A^T)Ax

then pre-multiply by the inverse of (A^T)A

    ((A^T)A)^{-1})(A^T)b = x

where ((A^T)A)^{-1})(A^T) is
the Moore-Penrose pseudoinverse.
This involves unnecessary computation
so we solve

    t = Sx

where t = (A^T)b and S = (A^T)A
for x instead using a Cholesky decomposion

    S = L(L^T)

so that

    t = L(L^T)x = L(u = (L^T)x)

and we use triangular solvers to solve

    t = Lu

for u then

    u = (L^T)x

for x.

This doesn't works unless

    n <= rank(A)

and S = (A^T)A is positive definite.

For an underdetermined set of equations

    b = Ax

try a Singular Value Decomposition

    A = UD(V^T)

where D >= 0 is a diagonal matrix and compute

    x = V(D^{-1})(U^T)b

if m < n.  Just substitute zero for 1/D_{jj}
for D_{jj} close to zero.

Take a look at "Numerical Recipes in C
 The Art of Scientific Computing: Second Edition",
by William H. Press, Saul A. Teukolsky,
William T. Vettering and Brian P. Flannery,
Chapter 2 Solutions of Linear Algebraic Equations,
Section 6 Singular Value Decomposition,
pages 59-70.