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.

Postby calderov » Tue Feb 05, 2013 3:19 pm

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 :cry: .

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.

Postby john » Wed Feb 06, 2013 7:14 am

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
Administrator
 
Posts: 587
Joined: Thu Jul 23, 2009 2:31 pm

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

Postby calderov » Fri Feb 08, 2013 8:57 am

Thank you for answering John.

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  |


Instead of the correct one:
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?

Please help! :(
calderov
 
Posts: 2
Joined: Wed Jan 16, 2013 12:00 pm

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

Postby john » Fri Feb 08, 2013 11:22 am

Yep, I can confirm this - you've overrun single precision. You'll get the right answer in double.
john
Administrator
 
Posts: 587
Joined: Thu Jul 23, 2009 2:31 pm


Return to CULA Dense Support

Who is online

Users browsing this forum: No registered users and 1 guest

cron