## Help: Least Squares function not working as expected.

General CULA Dense (LAPACK & BLAS) support and troubleshooting. Use this forum if you are having a general problem or have encountered a bug.

### Help: Least Squares function not working as expected.

Hello.

I'm having trouble writting a least squares program to fit a 5th order polynomial to the following set of points (1,1) , (2,4), (3,9), (4,11), (5,16), (6,32).

I've made the calculations (using sympy) to build the matrices involved in the equation system a = (M^T * M)^-1 * M^T * y

Code: Select all
`M = [  6     21      91      441      2275     12201  ][                                                 ][ 21     91     441     2275     12201     67171  ][                                                 ][ 91     441    2275    12201    67171     376761 ][                                                 ][ 441   2275   12201    67171    376761   2142595 ][                                                 ][2275   12201  67171   376761   2142595   12313161][                                                 ][12201  67171  376761  2142595  12313161  71340451]y = [  73  ][      ][ 352  ][      ][ 1826 ][      ][ 9892 ][      ][55082 ][      ][312412]`

I expect to have the following vector as result of my CULA code

Code: Select all
`a = [  25  ][      ][  1579][- ----][   30 ][      ][ 467  ][ ---  ][  12  ][      ][  283 ][- --- ][   24 ][      ][  19  ][  --  ][  12  ][      ][-3/40 ]`

But instead I got the following
Code: Select all
`a =[-14.81][      ][ 30.34][      ][-21.43][      ][ 7.99 ][      ][-1.41 ][      ][ 0.09 ]`

As you can see, the expected and current results differ a lot .

This is the code I wrote to solve the problem. It throws the previous result and debuger shows the next message:

"First-chance exception at 0x768ec41f (KernelBase.dll) in mySystemSolve.exe: Microsoft C++ exception: cudaError_enum at memory location 0x0021f99c.."

Code: Select all
`//Least Squares demo with CULA #include <stdio.h>#include <stdlib.h>#include <cula_lapack.h>void checkStatus(culaStatus status){  char buf[256];  if(!status)    return;  culaGetErrorInfoString(status, culaGetErrorInfo(), buf, sizeof(buf));  printf("%s\n", buf);  culaShutdown();  exit(EXIT_FAILURE);}void printM(float* M,int rows,int cols, int ld){  int i;  int j;  for (i = 0; i < rows; ++i)  {    for (j = 0; j < cols; ++j)    {      printf("%.2f ", M[j*ld + i]);    }    printf("\n");  }}void leastSquares(){  culaStatus status;  const int ROWS = 6;  const int COLS = 6;  const int NRHS = 1;  const int LDM = 6;  const int LDY = 6;  float M[36] = {6.0, 21.0, 91.0, 441.0, 2275.0, 12201.0, 21.0, 91.0, 441.0, 2275.0, 12201.0, 67171.0, 91.0, 441.0, 2275.0, 12201.0, 67171.0, 376761.0, 441.0, 2275.0, 12201.0, 67171.0, 376761.0, 2142595.0, 2275.0, 12201.0, 67171.0, 376761.0, 2142595.0, 12313161.0, 12201.0, 67171.0, 376761.0, 2142595.0, 12313161.0, 71340451.0};  float y[6]  = {73.0, 352.0, 1826.0, 9892.0, 55082.0, 312412.0};    printf("\n\nM = \n");  printM(M, ROWS, COLS, ROWS);  printf("\n\ny = \n");  printM(y, ROWS, 1, ROWS);  status = culaSgels('N', COLS, ROWS, NRHS, M, LDM, y, LDY);  checkStatus(status);  printf("\n\na = \n");  printM(y, ROWS, 1, ROWS); }int main(int argc, char const *argv[]){  culaStatus status;  printf("Initializing CULA\n");  status = culaInitialize();  checkStatus(status);  printf("Done!\n");  leastSquares();  culaShutdown();  getchar();  return 0;}`

I will be very grateful if you can help me with this piece of software.

Thank you!
Last edited by calderov on Wed Feb 06, 2013 9:26 pm, edited 1 time in total.
calderov

Posts: 2
Joined: Wed Jan 16, 2013 12:00 pm

### Re: Help: Least Squares function not working as expected.

You shouldn't break on first-chance exceptions. You'll likely be flagging internal exceptions that are not actually errors. Both CUDA and CULA will appear as first-chance.

Since you are solving a linear system, you should use culaSgesv() rather than culaSgels().
john

Posts: 587
Joined: Thu Jul 23, 2009 2:31 pm

### Re: Help: Least Squares function not working as expected.

I'm sorry. I'm new in CULA and LAPACK concepts but I need a least squares method to find vector x in a Ax = b system where:

Code: Select all
`A = |     6    21     91     441     2275    12201 |     |    21    91    441    2275    12201    67171 |    |    91   441   2275   12201    67171   376761 |    |   441  2275  12201   67171   376761  2142595 |    |  2275 12201  67171  376761  2142595 12313161 |    | 12201 67171 376761 2142595 12313161 71340451 |b = |     73 |    |    352 |    |   1826 |    |   9892 |    |  55082 |    | 312412 |`

I need the least squares solution, I know I can use culaSgesv to solve this particular system but I'm writting a piece of software that needs to build several linear systems from a random sample, and I want to have a result even if the systems are undetermined or overdetermined.

So far this is the code I've used to try to solve the problem:

Code: Select all
`void leastSquares(){  culaStatus status;  const int ROWS = 6;  const int COLS = 6;  const int NRHS = 1;  const int LDA = ROWS;  const int LDB = ROWS;  //I'm aware of colum-major order, notice that A is symetric  float A[36] = {     6,    21,     91,     441,     2275,    12201,                      21,    91,    441,    2275,    12201,    67171,                      91,   441,   2275,   12201,    67171,   376761,                     441,  2275,  12201,   67171,   376761,  2142595,                    2275, 12201,  67171,  376761,  2142595, 12313161,                   12201, 67171, 376761, 2142595, 12313161, 71340451};  float b[6]  = {     73,                     352,                     1826,                     9892,                    55082,                   312412};  status = culaSgels('N', ROWS, COLS, NRHS, A, LDA, b, LDB);  checkStatus(status);  printf("\n\nx = \n");  printM(b, COLS, 1); }`

I get the following output:
Code: Select all
`x = | -14.81 |    |  30.34 |    | -21.43 |    |  7.99  |    | -1.41  |    |  0.09  |`

Code: Select all
`x = |       25 |    | -1579/30 |    |   467/12 |    |  -283/24 |    |    19/12 |    |    -3/40 |`

I've tried with simpler and smaller systems and it has worked fine. Is this a precision problem? I'm I doing something wrong?

calderov

Posts: 2
Joined: Wed Jan 16, 2013 12:00 pm

### Re: Help: Least Squares function not working as expected.

Yep, I can confirm this - you've overrun single precision. You'll get the right answer in double.
john