6
$\begingroup$

I am trying to use LAPACK geqp3 function for QR decomposition with pivoting but result is wrong. I tried almost two days but can't figure out the problem. I compared the result with Matlab and Python CvxOpt both are same and both uses the geqp3 LPACK routines as internally .

Here is input matrix and code. Any one can tell me what i am doing wrong ?

> Blockquote A =
   -0.9586   -0.8171   -0.6033    1.1895    1.3031   -0.0759
    1.9877    1.0979    0.7795   -0.7164   -0.0505   -0.2412
    0.2328    0.2210    0.9401    0.0691    0.1082   -0.8700
    0.2808    0.9539    0.0198   -0.2278   -1.1498   -0.3127
    0.1335    1.2238    0.5449    0.2712   -0.7684   -2.5430
    0.3996   -0.5725    0.7090   -1.2492    0.5508   -1.0513
   -0.3073    0.2512   -1.2444    0.1169   -1.1271    0.6893
    0.3238   -0.8187    2.7957    0.3582    1.8493   -1.5166
    0.6406   -0.4428   -0.1864   -1.0097    0.6713    0.7002


void fReadMtrx(string f, double *mtrx, int nrow, int ncol){

    ifstream data(f.c_str(), ios::in);
    if (!data) {
        cerr << "File " << f << " could not be opened." << endl;
        exit(1);
    }

    for(int i = 0; i < nrow; i++) 
        for(int j = 0; j < ncol; j++) 
            data >> mtrx[j*nrow+i]; 
}

void showMtrx(double *mtrx, int nrow, int ncol){

    for (int i = 0; i < nrow; i++) {
        for (int j = 0; j < ncol; j++) { 
            cout <<  mtrx[j*nrow+i] << "  ";
        }
        cout << endl;
    }
}


int main(int argc, char *argv[])
{
    int m = 9;
    int n = 6;

    //allocate and read in a
    double* a = new double[m*n];

    fReadMtrx("A.dat", a, m, n);

    cout << endl << "a:" << endl;
    showMtrx(a, m, n);

    int ldA = m;
    double wl;
    int lwork;
    int info;


    // QR decomposition
    lwork = -1;
    dgeqp3_(&m, &n, NULL, &ldA, NULL, NULL, &wl, &lwork, &info);

    lwork = (int)wl;
    double* work  = new double[lwork];
    int* jpvt_ptr = new int[n];
    double* tau = new double[n];

    dgeqp3_(&m, &n, a, &ldA, jpvt_ptr,tau,(double *)work, &lwork, &info);

    cout << endl << "QR:" << endl;
    showMtrx(a, m, n);


   return 0;
}

When i perform QR decomposition using geqp3 routines it gives me this result.

QR :
 2.4063   1.0778   1.6110  -1.5135  -0.0642  -0.7111 
-0.5907   2.1035  -1.1498   0.4283  -2.6023  -0.6335 
-0.0692  -0.0423   2.8526   0.7086   0.9759  -2.7723
-0.0835  -0.3745  -0.0332   1.3431   0.8472   0.5267
-0.0397  -0.5405  -0.3887  -0.3816   0.7893  -1.2573
-0.1188   0.3753  -0.1044   0.9602  -0.0370   1.1425
 0.0913  -0.1996   0.4502  -0.1306   2.1257  -0.1959
-0.0962   0.4710  -1.1228   0.0682   1.4831  -0.4424
-0.1904   0.3782   0.4133   0.4090  -1.0372   0.2011

While in python and matlab i am getting this result.

QR :
   [ 3.47196268, 1.63382425, -2.39788322, -0.26189627,  1.11652819,  -0.1965108 ],
   [-0.19128521,  2.52976731,  0.93296076,  0.49012156, -0.78217128, -2.06427021],
   [-0.23067443, -0.0120604 ,  2.26624645, -0.67846333,  0.74839188,  -0.28431637],
   [-0.00486702,  0.43558037, -0.00634783, -2.00278565,  1.22747977,  -0.01853067],
   [-0.13370508,  0.3074231 ,  0.76958444, -0.0335343 ,  1.36566733,   1.04995806],
   [-0.17396713, -0.18659867,  0.30029599, -0.6800568 ,  0.37513732,   0.3212295 ],
   [ 0.30534072,  0.38815978, -0.10102664,  0.06497129,  0.35993755,   0.72318025],
   [-0.68601098, -0.61373672,  0.14243916,  0.28777571, -0.68376182,  0.3557272 ],
   [ 0.04573902, -0.25965053, -0.15498792, -0.52485689, -0.0986465 ,  -0.38370839]
$\endgroup$
2
  • 2
    $\begingroup$ It's not clear to me from your post; what makes the result wrong? $\endgroup$ Jan 13, 2013 at 8:17
  • $\begingroup$ Also, exactly what output did your code produce? $\endgroup$ Jan 14, 2013 at 6:16

1 Answer 1

8
$\begingroup$

Interpreting the results

To those of you confused as to the difference, he is reporting the compact representation of the pivoted QR factorization, where the upper triangle gives the factor $R$ from the pivoted QR decomposition and each column of the strictly lower triangle represents a Householder reflection, say $$ H_j = I - \tau_j \left[\begin{array}{c}1\\ u_j\end{array}\right] \left[\begin{array}{cc}1 & u_j^H\end{array}\right],$$ where $u_j$ is the vector stored below the $j$'th diagonal entry of the result.

I recommend loading his input matrix into MATLAB or Octave in order to follow along at home. You will notice that the first result is essentially the same as an unpivoted QR factorization, e.g., via

X=qr(A)

but with signs of the upper triangle appropriately changed so that the diagonal entries are positive. This is a strong clue as to the problem.

Fixing the problem

If you look through the comments at the beginning of the dgeqp3 source, you will notice that if an entry of the JPVT array is nonzero, then that column is not allowed to be a pivot candidate. Notice that your code is simply allocating an integer array, which does not ensure that its contents are zero. If you explicitly zero the JPVT array before calling dgeqp3 then I suspect that your problem will go away.

$\endgroup$
2
  • 1
    $\begingroup$ You are right it was the problem of JPVT. it working fine now thanx again for suggestion $\endgroup$
    – alex
    Jan 18, 2013 at 8:15
  • $\begingroup$ Wonderful! Would you mind accepting my answer then? $\endgroup$ Jan 18, 2013 at 15:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.