# Wrong Result QR decomposition with Pivoting in LAPACK

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;

double* a = new double[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]

-
It's not clear to me from your post; what makes the result wrong? –  Geoff Oxberry Jan 13 '13 at 8:17
Also, exactly what output did your code produce? –  Brian Borchers Jan 14 '13 at 6:16

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.

-
You are right it was the problem of JPVT. it working fine now thanx again for suggestion –  alex Jan 18 '13 at 8:15
Wonderful! Would you mind accepting my answer then? –  Jack Poulson Jan 18 '13 at 15:22