> 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 > >>> [email protected] > >>> http://lists.gnu.org/mailman/listinfo/help-gsl > >>> > >> _______________________________________________ > >> Help-gsl mailing list > >> [email protected] > >> http://lists.gnu.org/mailman/listinfo/help-gsl > >> > >> > > > > > > _______________________________________________ > > Help-gsl mailing list > > [email protected] > > http://lists.gnu.org/mailman/listinfo/help-gsl > > > _______________________________________________ Help-gsl mailing list [email protected] http://lists.gnu.org/mailman/listinfo/help-gsl
