Re: puzzled with result from gsl_multifit_linear...

Hmm, this is an interesting problem. As John noted d = 0 in your
example so you need to use the form ax + by + cd + d = 0, which means
the RHS of your least squares fit is 0.

This means you are trying to solve the least squares problem

X c = 0

which means the solution is the null space of the matrix X. The
null-space of X corresponds to a singular value of 0 in the SVD
decomposition (the algorithm GSL uses to solve least squares).
Unfortunately, according to the GSL manual gsl_multifit_linear:

"Any components which have zero singular value (to machine precision)
are discarded from the fit."

This means GSL will return the solution c = 0, which is wrong since
there is a well defined non-zero solution to X c = 0 for your example.
(the solution to your problem is: 0x + y - z + 0 = 0, or y = z), so
c = [0 1 -1 0] (the null space of the matrix X).

One solution to your problem is to not use the function
gsl_multifit_linear, and instead call the gsl_linalg_svd_decomp routine
yourself, search for singular values in the matrix S and when you find
them the corresponding vector in V is a null space vector.

Perhaps the routine gsl_multifit_linear should be modified to
provide a correct solution in the case of zero singular values, but
maybe Brian will have more to say about this. Right now the routine
gsl_multifit_linear computes the SVD of X = U S V^t and then simply
says c = V S^{-1} U^t y to find c, but of course if there is a zero
singular value then S^{-1} does not exist, so the result returned
by gsl is wrong.

Off the top of my head I think gsl_multifit_linear can be modified
to find the correct solution in the case of a zero RHS, simply by
searching the singular values and if it finds one that is zero, compute
the null-space of the matrix X instead of S^{-1}.

On Tue, Nov 20, 2007 at 01:07:03PM +1100, John Gehman wrote:
> G'day John,
> 
> I haven't looked at the examples all that closely, but just quickly,  
> I believe a general equation for a plane is ax+by+cz+d=0, i.e. in  
> your equation you're effectually fixing d at -1, while I think it  
> should be d=0 (zero) in your example?
> 
> Cheers,
> john
> 
> ---------------------------------------------------------
> 
> John Gehman
> Research Fellow
> School of Chemistry
> University of Melbourne
> Australia
> 
> On 20/11/2007, at 12:53 PM, John Pye wrote:
> 
> >Hi all
> >
> >I have a question about the use of the gsl_multifit_linear routine  
> >that
> >perhaps is a question more about geometry/algebra than coding, but I'm
> >not sure.
> >
> >I want to construct a routine that fits a plane (in three dimensions)
> >through a set of data points (x,y,z). I have set up  
> >gsl_multifit_linear
> >to fit the plane equation a*x + b*y + c*z = 1to my data, and for most
> >cases that seems to work OK. However there are a few degenerate cases
> >that don't work, and I'm trying to work out what I should do. Is  
> >there a
> >better equation that describes a plane?
> >
> >The following example shows what happens when I try to fit a plane to
> >the points (0,0,0), (1,0,0), (1,1,1), (0,1,1). These points  
> >represent an
> >inclined rectangle with one side on the positive x-axis.
> >
> >I would expect the fit to these points to be 0*x + INF*y - INF*z =  
> >1, or
> >else for it to give an error. But instead, gsl_multifit_linear  
> >gives me
> >the result 0.666*x + 0.333*y + 0.333*z = 1, which is wrong.
> >
> >I am obviously making a logic error here, but I'm a bit stuck on  
> >what it
> >is, and what the correct general way of working around it would be?  
> >Can
> >anyone here offer any suggestions?
> >
> >Cheers
> >JP
> >
> >
> >// for the fitting of a plane through a group of points, testfit.cpp
> >#include <gsl/gsl_vector.h>
> >#include <gsl/gsl_matrix.h>
> >#include <gsl/gsl_multifit.h>
> >#include <iostream>
> >
> >using namespace std;
> >
> >int main(void){
> >
> >    int n = 4;
> >    int num_params = 3;
> >
> >    gsl_multifit_linear_workspace *work = gsl_multifit_linear_alloc(n,
> >num_params);
> >
> >    gsl_vector *y = gsl_vector_calloc(n);
> >    for(unsigned i=0; i < n; ++i){
> >        gsl_vector_set(y,i,1.0);
> >    }
> >
> >    gsl_matrix *X = gsl_matrix_calloc(n, num_params);
> >    gsl_vector *c = gsl_vector_calloc(num_params);
> >    gsl_matrix *cov = gsl_matrix_calloc(num_params, num_params);
> >
> >#define SET(I,A,B,C) \
> >    gsl_matrix_set(X,I,0,A); gsl_matrix_set(X,I,1,B);
> >gsl_matrix_set(X,I,2,C)
> >
> >    SET(0, 0,0,0);
> >    SET(1, 1,0,0);
> >    SET(2, 1,1,1);
> >    SET(3, 0,1,1);
> >
> >    double chisq;
> >    int res = gsl_multifit_linear(X,y,c,cov,&chisq,work);
> >
> >    cerr << "FIT FOUND:" << endl;
> >    for(int i=0;i<3;++i){
> >        cerr << "  c[" <<i << "]=" <<gsl_vector_get(c,i) << endl;
> >    }
> >
> >}
> >
> >
> >_______________________________________________
> >Help-gsl mailing list
> >Help-gsl@xxxxxxx
> >http://lists.gnu.org/mailman/listinfo/help-gsl
> 
> _______________________________________________
> Help-gsl mailing list
> Help-gsl@xxxxxxx
> http://lists.gnu.org/mailman/listinfo/help-gsl
>