Dense under-determined system with parametric solutions

[ejspeiro]

Dear all,

I am interesting in using the LAPACK to solve for under-determined Vandermonde systems, such as (for example):



I have many related questions:

• How can I have access to the particular solution depicted in the image (the one for which the parameter is 0)?
• How can I have access to the one element in the null-space, i.e. (in the image) [-1, 5, -10, 10, -5, 1]^T

Thanks in advanced!

[Julien Langou]

DGELS might be the routine you are looking for. For underdetermined system of equations, there are several ways to determine a solution. I am not certain DGELS is what you want. J.

[ejspeiro]

Dear Julien,

Thanks for your reply. I have come across this routine, and I will give it a try, yet your reply started positive, and then:

I am not certain DGELS is what you want

Do you mind expanding on why aren't you certain?

Thanks!

[Julien Langou]

Hi,

So you have a 5-by-6 matrix A. The matrix A has full rank 5.
You have a right-hand side b of length 5. You wish to solve Ax =b.
In your case, we have

Code: Select all

a =[ -1/2 1/2 3/2 5/2 7/2 9/2 ]; A = [ a.^0; a.^1; a.^2; a.^3; a.^4 ];
b = [ 0; 1; 0; 0; 0];

In this case, the set of solution is of the form ( x1 + Null( A ) ). So you are looking for
(1) one particular solution. In your case, you chose x1 = [-11/12; 17/24; 3/8; -5/24; 1/24; 0]. Fine.
(2) and a nonzero vector in the null space. In your case you chose u = [-1; 5; -10; 10; -5; 1];
so that all solution writes x1 + k * u where k is any scalar. All good.

So now how LAPACK can help you solve your problem.

First to find a particular solution you can use LAPACK DGELS. Now LAPACK DGELS will give the particular solution with minimum norm. In general it has no zero in it. In your case, the answer returned by LAPACK is:

Code: Select all

x2 =
        -0.922949735449755
         0.739748677248692
         0.312169312169325
         -0.14550264550266
        0.0102513227513299
        0.0062830687830675


Note that Ax2 = b (as Ax1 =b) but || x2 || < || x1 ||. This is a consequence of the fact that the x2 returned by LAPACK not only is a solution, it is the solution with minimum 2-norm.
A point of reference to understand the lapack answer as opposed to say the matlab A\b answer (the matlab's answer will have some zeros in it) is my post on the forum: viewtopic.php?f=5&t=160

Second to get the null space, several ways would be OK. You can do an SVD of A but that would be expensive. The best way is to do a FULL QR factorization of A^T and get the Null space from the last column of Q. (In your case, we know that the null space is of dimension 1, so we know we are interested in the last column of Q.) Then Null(A) = Span( w ).

So now you can both (1) and (2) at once. By "at once", I mean with only one QR factorization of A.
This is pretty cool. The problem is that you do not have an LAPACK code to do this. I can provide a Matlab pseudocode below:

Code: Select all

[Q,R]=qr(A');
w=Q(:,6);
x2 = Q(1:6,1:5)*((R(1:5,1:5)')\b);


Then solutions are of the form x2 + k*w where k is any scalar. x2 is a particular solution of Ax =b, to be precise, x2 is the solution of Ax=b with minimum 2-norm. w is a nonzero vector in Null(A).

Julien.

[ejspeiro]

Dear Julien,

Thanks for such a well elaborated response to this. The theoretical implications of my problem regarding the LAPACK have now been clarified.

The issues I am facing now are of a computational nature. Specifically, the fact that I am invoking the DGELS routine from C++.

I keep getting segmentation faults because of the C++/Fortran transposition issues. Considering the actual form of the matrix I want so solve for, given on my first post, I want to double check how am I invoking DGELS.

1. Constructing the matrix

First, my matrix it's being built already transposed:

Code: Select all

double* Vandermonde(double *gen, int length_gen, int num_rows, bool transpose) {
   
   double *TT; // Output Vandermonde matrix.
   int ii;         // Iterator.
   int jj;         // Iterator.
   
   TT = NULL;
   if (length_gen < 1) {
      return NULL;
   }
   if (num_rows < 1) {
      return NULL;
   }
   TT = new double[num_rows*length_gen];
  if (TT == (double*) NULL) {
    cerr << "Could not calloc at line " << __LINE__ - 2 << endl;
      return NULL;
  }
  for (ii = 0; ii < num_rows; ii++) {
    for (jj = 0; jj < length_gen; jj++) {
         if (transpose) {
            TT[jj*num_rows + ii] = pow(gen[jj], (double) ii);
         } else {
            TT[ii*length_gen + jj] = pow(gen[jj], (double) ii);
         }
    }
  }
  return TT;
}


Such a routine, when called as:

Code: Select all

TT = Vandermonde([ -1/2 1/2 3/2 5/2 7/2 9/2 ], 6, 5, true);


Yields:

Code: Select all

Vandermonde matrix 1 out of 1

          1       -0.5       0.25     -0.125     0.0625
          1        0.5       0.25      0.125     0.0625
          1        1.5       2.25      3.375     5.0625
          1        2.5       6.25     15.625    39.0625
          1        3.5      12.25     42.875    150.062
          1        4.5      20.25     91.125    410.062


2. Invoking DGELS

This is how I do it... I first inquire the right size of the work array, by doing:

Code: Select all

      work = new double[1];
      if (work == (double*) NULL) {
         cerr << "Could not calloc at line " << __LINE__ - 2 << endl;
          return EXIT_FAILURE;
      }
      lwork = -1;
      dgels_(&trans, &TT_num_rows, &TT_num_cols, &nrhs, TT, &TT_ld, ob, &ob_ld, work, &lwork, &info);
      lwork = (int) work[0];
               cout << "lwork = " << endl << endl << setw(12)<< lwork << endl<< endl;


I then reallocate the work array, and invoke the routine doing the following:

Code: Select all

      trans = 'N';  // Not transposed, because I already build it transposed...
      info = 0;
      TT_num_rows = 5;
      TT_num_cols = 6;
      TT_ld = 5;
      ob_ld = 6;

      delete[] work;
      work = new double[lwork];
      if (work == (double*) NULL) {
         cerr << "Could not calloc at line " << __LINE__ - 2 << endl;
         return EXIT_FAILURE;
      }
      
      // We now invoke the solver again:
      dgels_(&trans, &TT_num_rows, &TT_num_cols, &nrhs, TT, &TT_ld, ob, &ob_ld, work, &lwork, &info);


BUT executing leads to:

Code: Select all

lwork =

         165

System successfully solved!

    -0.92295    0.739749    0.312169   -0.145503   0.0102513

*** glibc detected *** ./mdivergence: double free or corruption (out): 0x000000000fd915e0 ***


ALMOST THERE! But, why is the solution of dimension 5 I specified its leading dimension to be 6, and that is how your reference solution reads like.

Thanks in advanced!

[ejspeiro]

UPDATE.

Because LAPACK uses the right-hand side, to overwrite it with the solution!

Then my question becomes: how can I work this problem, where the original RHS vector is a 5 x 1 vector, BUT the solution is supposed to be a 6 x 1 vector?

Thanks.

[ejspeiro]

UPDATE:

Solved. I allocated an extra element for it, and then I made sure I was printing the right amount of elements out of it:

Code: Select all

System successfully solved!

    -0.92295    0.739749    0.312169   -0.145503   0.0102513  0.00628307


Now up to computing the null-space!

[ejspeiro]

Dear Dr. Julien,

Once again thank you for all of your help.

I am trying to code the QR factorization using LAPACK to get the null space. The MATLAB driver I created based on your pseudo-code works great. In fact, I as able to scale the solutions getting those matching my original example:

Code: Select all

a =[ -1/2 1/2 3/2 5/2 7/2 9/2 ];
A = [ a.^0; a.^1; a.^2; a.^3; a.^4 ]
A'
b = [ 0; 1; 0; 0; 0]
[Q,R] = qr(A')
w = Q(:,6);
w = w/w(6)
% w is now the kernel that I showed on my original post.
x1 = A\b;
x1 = x1 - x1(6)*w
% x1 is now the solution that I showed on my original post.


I want to duplicate this using the LAPACK. For which I am aware that I should:

1. Invoke dgeqrf_ to inquire for the value of lwork.
2. Invoke it again to attain both A and tau, containing the results of the factorization.
3. I then invoke dorgqr_ to get the matrix Q, from which I can get the required null space.

Notice that I am assuming that the same value for lwork works for both dgeqrf_ and dorgqr_. Based on all of this, my questions are:

• Could you show me the values you are getting for: A and tau, after executing dgeqrf_? It was a tremendous help when you did it with the vector x2, on previous posts within this thread
• Could you tell me if once lwork has been computed, I can use it for both dgeqrf_ and dorgqr_?

Thank you very much in advanced, and have a great day!