Re: [Help-gsl] puzzled with result from gsl_multifit_linear...

[Patrick Alken]

> Ok, I've used the gsl_linalg_SV_decomp function to obtain the singular
> value decomposition, and indeed I have a very small (1e-33) value in my
> vector 'S'.
> 
> I'm still feeling my way around GSL: what do you suggest I do next to
> find the correct fitting plane? Can I use the remaining values in 'S' to
> form my solution? Or do I need to set up another linear fit with a
> reduced set of variables?

Well if you find a zero value in position 'i' in S, then the solution
is given by the column i of the matrix V as output by the svd
routine.

Some more info:
http://en.wikipedia.org/wiki/Singular_value_decomposition#Range.2C_null_space_and_rank

So if you find a tiny singular value, just set the coefficient vector
equal to the corresponding column of V. Otherwise use the normal
gsl_multifit_linear.

In the meantime we can start thinking if its possible to modify
gsl_multifit_linear to handle this case.

> 
> BTW if it were possible to fix gsl_multifit_linear for this case (or
> alternatively to provide gsl_multifit_mul) that would be great!
> 
> Cheers
> JP
> 
> > 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
> >>> address@hidden
> >>> http://lists.gnu.org/mailman/listinfo/help-gsl
> >>>       
> >> _______________________________________________
> >> Help-gsl mailing list
> >> address@hidden
> >> http://lists.gnu.org/mailman/listinfo/help-gsl
> >>
> >>     
> >
> >
> > _______________________________________________
> > Help-gsl mailing list
> > address@hidden
> > http://lists.gnu.org/mailman/listinfo/help-gsl
> >   
>