Subject: Re: puzzled with result from gsl_multifit_linear... - msg#00043

Previous Message by Date: click to view message preview

Re: puzzled with result from gsl_multifit_linear...

Hi Patrick, Patrick Alken wrote: > 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. > 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? 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 >>> 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 >> >> > > > _______________________________________________ > Help-gsl mailing list > Help-gsl@xxxxxxx > http://lists.gnu.org/mailman/listinfo/help-gsl >

Next Message by Date: click to view message preview

Re: puzzled with result from gsl_multifit_linear...

Hi Patrick Patrick Alken wrote: >> 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. > OK, I managed to implement this using the gsl_linalg_SV_decomp and gsl_linalg_SV_solve functions as following. In case it's useful to anyone else. It nicely picks up the case of points being colinear as well. Cheers JP static void centroid(gsl_matrix *X, double n, double num_params, double *cent){ /* prob there is a more GSL-ish way of doing this... */ int i,j; for(j=0;j<num_params;++j)cent[j]=0; for(i=0;i<n;++i){ for(j=0;j<num_params;++j)cent[j]+=gsl_matrix_get(X,i,j); } for(j=0;j<num_params;++j)cent[j]/=n; } Plane PlaneFit::leastSquaresPlane(){ unsigned n = G.size(); if(n < 3){ throw runtime_error("Not enough points for PlaneFit"); } int num_params = 3; gsl_matrix *X = gsl_matrix_calloc(G.size(), num_params); gsl_vector *c = gsl_vector_calloc(num_params); // load data into the matrix 'X' int ii=0; for(Group::const_iterator i=G.begin(); i<G.end(); ++i){ Vector xi = (*i)->pos(); gsl_matrix_set(X,ii,0, xi.x); gsl_matrix_set(X,ii,1, xi.y); gsl_matrix_set(X,ii,2, xi.z); ++ii; } // calculate the centroid double cent[3]; centroid(X, n, num_params, cent); // offset the points by the centroid vector for(int i=0;i<n;++i){ for(int j=0;j<num_params;++j){ double v = gsl_matrix_get(X,i,j); gsl_matrix_set(X,i,j, v - cent[j]); } } // calulate the fit using SVD directly (so we catch degenerate cases) gsl_matrix *V = gsl_matrix_calloc(num_params,num_params); gsl_vector *S = gsl_vector_calloc(num_params); int res = gsl_linalg_SV_decomp_jacobi(X,V,S); // check for zero values in the vector 'S', this indicate redundant terms in the equation unsigned have_zero=0; unsigned j_zero; double tol = 1e-10; for(int j=0;j<num_params;++j){ if(gsl_vector_get(S,j) < tol){ // guaranteed greater than zero so no 'fabs' have_zero++; j_zero = j; } } if(have_zero > 1){ // perhaps this case occurs when the points are colinear throw runtime_error("PlaneFit: Points do not define a plane (MULTIPLE ZEROS IN 'S' VECTOR)"); }else if(have_zero==1){ //cerr << "SINGULAR!" << endl; for(int j=0;j<num_params;++j){ gsl_vector_set(c,j, gsl_matrix_get(V,j,j_zero)); } }else{ // the general case, use the SVD to solve X*c=0 gsl_vector *y = gsl_vector_calloc(n); for(int i=0; i<n; ++i){ gsl_vector_set(y,i,0.0); } int res = gsl_linalg_SV_solve(X, V, S, y, c); gsl_vector_free(y); } // resulting plane fit is 0 = (x - x_0)*c[0] + (y-y_0)*c[1] + (z-z_0)*c[2] /*cerr << "FIT FOUND:" << endl; for(int i=0;i<3;++i){ cerr << " c[" <<i << "]=" <<gsl_vector_get(c,i) << endl; }*/ // create plane using direction vector 'd' and location 'p' (the centroid) Vector d(gsl_vector_get(c,0),gsl_vector_get(c,1),gsl_vector_get(c,2)); Vector p(cent[0], cent[1], cent[2]); Plane P(p,d); //cerr << "CENTROID = " << p << endl; //cerr << "FIT PLANE = " << P << endl; gsl_matrix_free(X); gsl_vector_free(c); gsl_matrix_free(V); gsl_vector_free(S); return P; }

Previous Message by Thread: click to view message preview

Re: puzzled with result from gsl_multifit_linear...

Hi Patrick, Patrick Alken wrote: > 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. > 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? 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 >>> 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 >> >> > > > _______________________________________________ > Help-gsl mailing list > Help-gsl@xxxxxxx > http://lists.gnu.org/mailman/listinfo/help-gsl >

Next Message by Thread: click to view message preview

Re: puzzled with result from gsl_multifit_linear...

Hi Patrick Patrick Alken wrote: >> 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. > OK, I managed to implement this using the gsl_linalg_SV_decomp and gsl_linalg_SV_solve functions as following. In case it's useful to anyone else. It nicely picks up the case of points being colinear as well. Cheers JP static void centroid(gsl_matrix *X, double n, double num_params, double *cent){ /* prob there is a more GSL-ish way of doing this... */ int i,j; for(j=0;j<num_params;++j)cent[j]=0; for(i=0;i<n;++i){ for(j=0;j<num_params;++j)cent[j]+=gsl_matrix_get(X,i,j); } for(j=0;j<num_params;++j)cent[j]/=n; } Plane PlaneFit::leastSquaresPlane(){ unsigned n = G.size(); if(n < 3){ throw runtime_error("Not enough points for PlaneFit"); } int num_params = 3; gsl_matrix *X = gsl_matrix_calloc(G.size(), num_params); gsl_vector *c = gsl_vector_calloc(num_params); // load data into the matrix 'X' int ii=0; for(Group::const_iterator i=G.begin(); i<G.end(); ++i){ Vector xi = (*i)->pos(); gsl_matrix_set(X,ii,0, xi.x); gsl_matrix_set(X,ii,1, xi.y); gsl_matrix_set(X,ii,2, xi.z); ++ii; } // calculate the centroid double cent[3]; centroid(X, n, num_params, cent); // offset the points by the centroid vector for(int i=0;i<n;++i){ for(int j=0;j<num_params;++j){ double v = gsl_matrix_get(X,i,j); gsl_matrix_set(X,i,j, v - cent[j]); } } // calulate the fit using SVD directly (so we catch degenerate cases) gsl_matrix *V = gsl_matrix_calloc(num_params,num_params); gsl_vector *S = gsl_vector_calloc(num_params); int res = gsl_linalg_SV_decomp_jacobi(X,V,S); // check for zero values in the vector 'S', this indicate redundant terms in the equation unsigned have_zero=0; unsigned j_zero; double tol = 1e-10; for(int j=0;j<num_params;++j){ if(gsl_vector_get(S,j) < tol){ // guaranteed greater than zero so no 'fabs' have_zero++; j_zero = j; } } if(have_zero > 1){ // perhaps this case occurs when the points are colinear throw runtime_error("PlaneFit: Points do not define a plane (MULTIPLE ZEROS IN 'S' VECTOR)"); }else if(have_zero==1){ //cerr << "SINGULAR!" << endl; for(int j=0;j<num_params;++j){ gsl_vector_set(c,j, gsl_matrix_get(V,j,j_zero)); } }else{ // the general case, use the SVD to solve X*c=0 gsl_vector *y = gsl_vector_calloc(n); for(int i=0; i<n; ++i){ gsl_vector_set(y,i,0.0); } int res = gsl_linalg_SV_solve(X, V, S, y, c); gsl_vector_free(y); } // resulting plane fit is 0 = (x - x_0)*c[0] + (y-y_0)*c[1] + (z-z_0)*c[2] /*cerr << "FIT FOUND:" << endl; for(int i=0;i<3;++i){ cerr << " c[" <<i << "]=" <<gsl_vector_get(c,i) << endl; }*/ // create plane using direction vector 'd' and location 'p' (the centroid) Vector d(gsl_vector_get(c,0),gsl_vector_get(c,1),gsl_vector_get(c,2)); Vector p(cent[0], cent[1], cent[2]); Plane P(p,d); //cerr << "CENTROID = " << p << endl; //cerr << "FIT PLANE = " << P << endl; gsl_matrix_free(X); gsl_vector_free(c); gsl_matrix_free(V); gsl_vector_free(S); return P; }