This is the mail archive of the
gsl-discuss@sources.redhat.com
mailing list for the GSL project.
least square problem solution
- From: Edwin Robert Tisdale <E dot Robert dot Tisdale at jpl dot nasa dot gov>
- To: gsl-discuss at sources dot redhat dot com
- Date: Tue, 22 Jan 2002 15:04:03 -0800
- Subject: least square problem solution
- Organization: Jet Propulsion Laboratory
- Reply-to: 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.