Problems with results for SGELS for least squares in C++

[Kadence]

I'm trying to test using sgels_ in C++ to solve least squares problems. However the results I get don't seem correct, and I don't know what I'm doing wrong. Probably something to do with FORTRAN/C++ quirks (e.g. matrices being inverted). Here's my code:

Code: Select all

#include <iostream>

using namespace std;

extern "C" void sgels_ ( ... );

int main(){
   float mA[8][2];
   float v[8];
   float mWork[16];

   mA[0][0]=100;   mA[1][0]=90;   v[0]=0.80;
   mA[0][1]=97;   mA[1][1]=90;   v[1]=0.72;
   mA[0][2]=104;   mA[1][2]=99;   v[2]=0.63;
   mA[0][3]=96;   mA[1][3]=96;   v[3]=0.55;
   mA[0][4]=99;   mA[1][5]=99;   v[4]=0.52;
   mA[0][5]=100;   mA[1][5]=97;   v[5]=0.50;
   mA[0][6]=96;   mA[1][5]=96;   v[6]=0.49;
   mA[0][7]=98;   mA[1][5]=100;   v[7]=0.45;

   char trans = 'N';
   int m=8;
   int n=2;
   int nrhs = 2;
   int lwork = 16;
   int info = 0;

   sgels_(&trans, &m, &n, &nrhs, mA, &m, v, &m, mWork, &lwork, &info);

   if(info) printf("Error - sgels returned: %d\n", info);
   for(int i=0; i<m; i++){
      printf("%f\n", v[i]);
   }
}

And I'm compiling with "g++ linear.cpp -llapack -lblas -lg". The output I get is:

Code: Select all

0.005857
-0.011565
0.026700
-0.075284
-0.124824
-0.151338
-0.135284
-0.201338

As I understand it, the first two values should be the coefficients. But using those values gives extremely wrong results. Using Excel, the actual coefficients should be B1=0.03196 and B2=-0.02684 with no intercept, and B0=1.183, B1=0.0229, B2=-0.0299 with an intercept. What is my mistake?

Some other questions:

• What should NRHS be? It says, "The number of right hand sides, i.e., the number of columns of the matrices B and X. NRHS >=0." But I don't know if that means the number of columns each has (i.e. 1), or the sum (i.e. 2). It doesn't seem to matter here though, they both give the same values.
• For LWORK, what's NB, optimum block size?
• Is there any way one can set a gels routine to return an intercept? I didn't see anything about that on the man page, or at the netlib LLS page.

[Julien Langou]

If you want to solve:

Code: Select all

   ( 100   90 ) ( x1 )   ( 0.80 )
   (  97   90 ) ( x2 )   ( 0.72 )
   ( 104   99 )          ( 0.63 )
   (  96   96 )          ( 0.55 )
   (  99   99 )        = ( 0.52 )
   ( 100   97 )          ( 0.50 )
   (  96   96 )          ( 0.49 )
   (  98  100 )          ( 0.45 )


Try something like this:

Code: Select all

#include <iostream>

extern "C" void sgels_ ( ... );

int main(){
   float mA[8*2];
   float v[8];

   mA[ 0]=100;   mA[ 8]= 90;   v[0]=0.80;
   mA[ 1]= 97;   mA[ 9]= 90;   v[1]=0.72;
   mA[ 2]=104;   mA[10]= 99;   v[2]=0.63;
   mA[ 3]= 96;   mA[11]= 96;   v[3]=0.55;
   mA[ 4]= 99;   mA[12]= 99;   v[4]=0.52;
   mA[ 5]=100;   mA[13]= 97;   v[5]=0.50;
   mA[ 6]= 96;   mA[14]= 96;   v[6]=0.49;
   mA[ 7]= 98;   mA[15]=100;   v[7]=0.45;

   char trans = 'N';
   int m=8;
   int n=2;
   int nrhs = 1;
   int lwork;
   float *mWork;
   int info = 0;

   lwork=-1;
   mWork=(float *)malloc(1*sizeof(float));
   sgels_(&trans, &m, &n, &nrhs, mA, &m, v, &m, mWork, &lwork, &info);
   lwork = (int) mWork[0];
   free(mWork);
   mWork=(float *)malloc(lwork*sizeof(float));
   sgels_(&trans, &m, &n, &nrhs, mA, &m, v, &m, mWork, &lwork, &info);
   free(mWork);

   if(info) printf("Error - sgels returned: %d\n", info);
   for(int i=0; i<m; i++){
      printf("%f\n", v[i]);
   }
}


nrhs is the number of right-hand sides, it's one in your case
you get lwork by doing a workspace query (call the routine once with
lwork=-1 and get the optimal lwork in work[0])

[Kadence]

Thank you very much, that worked :)

Is there any way to calculate an intercept using an xgels routine?

And what about standard error values for the coefficients - are those available in any way?

[Julien Langou]

Is there any way to calculate an intercept using an xgels routine?

I guess so, if you want to find the line y = a x + b, that matches
the best (in the linear least squares sense) a set of (x_i,y_i), and if
you cal b the intercept, then you can get b (and a ) by solving:

Code: Select all

( 1 x1 ) ( b )      ( y1 )
( 1 x2 ) ( a )      ( y2 )
( 1 x3 )         =  ( y3 )
( 1 x4 )            ( y4 )
( 1 x5 )            ( y5 )


And what about standard error values for the coefficients - are those available in any way?

Look at section 4 in
Marc Baboulin, Jack Dongarra, Serge Gratton, and Julien Langou.
LAWN193 -- Computing the Conditioning of the Components of a Linear
Least Squares Solution.
http://www.netlib.org/lapack/lawnspdf/lawn193.pdf
some sample code are available.

-j

[Kadence]

Thanks once again :)