gsl_multifit_linear fails in a certain case

gsl_multifit_linear fails in a certain case

[Patrick Alken]

This post is a continuation of a discussion on help-gsl. Basically
gsl_multifit_linear will fail if the following 2 conditions are met:

(1) the matrix X has no inverse (ie: it has a singular value of 0 in its
SVD decomp)

and

(2) the RHS of the fit is a zero vector

In this case, we are solving the equation X c = 0 and X has no inverse,
but gsl_multifit_linear still tries to compute the psuedo-inverse
and fudges the zero singular value to force it to make an inverse.
The result returned by the routine is c = 0 which is not the desired
solution.

A possible fix would be to check for a singular value of zero and if
one is found, check to see if the RHS vector is 0. If it is, then
don't try to compute S^{-1} but instead compute the null space of
X which is the corresponding column of the matrix V.

I'm not really sure what to do if there is more than one zero singular
value.

Re: gsl_multifit_linear fails in a certain case

[Brian Gough]

At Mon, 19 Nov 2007 23:07:32 -0700,
Patrick Alken wrote:
> 
> This post is a continuation of a discussion on help-gsl. Basically
> gsl_multifit_linear will fail if the following 2 conditions are met:
> 
> (1) the matrix X has no inverse (ie: it has a singular value of 0 in its
> SVD decomp)
> 
> and
> 
> (2) the RHS of the fit is a zero vector
> 
> In this case, we are solving the equation X c = 0 and X has no inverse,
> but gsl_multifit_linear still tries to compute the psuedo-inverse
> and fudges the zero singular value to force it to make an inverse.
> The result returned by the routine is c = 0 which is not the desired
> solution.
> 
> A possible fix would be to check for a singular value of zero and if
> one is found, check to see if the RHS vector is 0. If it is, then
> don't try to compute S^{-1} but instead compute the null space of
> X which is the corresponding column of the matrix V.
> 
> I'm not really sure what to do if there is more than one zero singular
> value.
> 

Perhaps the manual needs to be clearer about what is computed, which
is the unique least squares solution that has minimum norm solution
coefficients |c|^2 (this only comes into play if there's a null space
- otherwise there's only one solution anyway).  The idea is to
minimise any spurious components of the solution.  This is usually
what is needed when fitting data and pretty standard.

When the rhs=0 and there is one singular value there is still an
infinite number of solutions, as the overall scale is
undetermined. The least squares choice is to return the minimum norm
solution, i.e. zero. 

In the help-gsl question, which was about fitting the equation for a
plane, the appropriate method is orthogonal regression rather than
linear regression, as there's no dependent variable. 
  
-- 
Brian Gough

Re: gsl_multifit_linear fails in a certain case

[Patrick Alken]

> Perhaps the manual needs to be clearer about what is computed, which
> is the unique least squares solution that has minimum norm solution
> coefficients |c|^2 (this only comes into play if there's a null space
> - otherwise there's only one solution anyway).  The idea is to
> minimise any spurious components of the solution.  This is usually
> what is needed when fitting data and pretty standard.
> 
> When the rhs=0 and there is one singular value there is still an
> infinite number of solutions, as the overall scale is
> undetermined. The least squares choice is to return the minimum norm
> solution, i.e. zero. 
> 
> In the help-gsl question, which was about fitting the equation for a
> plane, the appropriate method is orthogonal regression rather than
> linear regression, as there's no dependent variable. 


Yes, after thinking about the problem some more I realized its a bit
more complicated than I originally thought.