Matrix: MathWorks/QRpivot

Description: Limitation of basic solution to x=A\b using qr(A); needs min 2-norm

 (bipartite graph drawing)

• Matrix group: MathWorks
• download as a MATLAB mat-file, file size: 13 KB. Use UFget(1894) or UFget('MathWorks/QRpivot') in MATLAB.

 Matrix properties number of rows 660 number of columns 749 nonzeros 3,808 structural full rank? yes structural rank 660 # of blocks from dmperm 9 # strongly connected comp. 1 explicit zero entries 0 nonzero pattern symmetry 0% numeric value symmetry 0% type real structure rectangular Cholesky candidate? no positive definite? no

 author P. Quillen editor T. Davis date 2008 kind counter-example problem 2D/3D problem? no

 Additional fields size and type b sparse 660-by-1

Notes:

```Counter-example problem from The MathWorks, Pat Quillen

This matrix was obtained from a MATLAB user.  It illustrates the
limitations inherent in computing a basic solution to an under-
determined system without the use of column pivoting.

With column pivoting (which can only be done in MATLAB with full
matrices) the problem is solved properly.

When finding the min 2-norm solution (ignoring fill-in):

[Q,R] = qr (A') ;
x = Q*(R'\b) ;

a good solution is found.  To reduce fill-in:

p = colamd (A') ;
[Q,R] = qr (A (p,:)') ;
x = Q*(R'\b(p)) ;

which also finds a good solution.

However, x=A\b computes a basic solution, using this algorithm:

q = colamd (A) ;
[Q,R] = qr (A (:,q)) ;
x = R\(Q'*b) ;
x (q) = q ;

which finds an error with norm(A*x-b) of 1e-9 in MATLAB 7.6.

With random permutations, and determining the cond(R1) of the leading
trianglar part (R is "squeezed" and the columns can be partitioned into
[R1 R2] where R1 is square and upper triangular) leads to the following
results.

Note that the error is high when condest(R1) is high.  Note in
particular the last trial.

So this clinches the question.  MATLAB's QR, and my new sparse QR, both
use a rank-detection method (by Heath) that does not do column pivoting,
and which is known to fail for some problems - for which Grimes & Lewis'
method will likely succeed.

The advantage of my QR is that I now always return R as upper
trapezoidal, so if the user is concerned, he/she can easily check
condest(R(:,1:m)) if m < n.

err 7.71e-07 condest R1 2.18e+12
err 1.25e-09 condest R1 9.82e+08
err 2.47e-09 condest R1 2.46e+11
err 4.00e-09 condest R1 4.03e+09
err 9.88e-10 condest R1 4.73e+09
err 2.25e-08 condest R1 5.34e+09
err 2.00e-08 condest R1 1.04e+09
err 1.09e-09 condest R1 6.83e+08
err 6.18e-08 condest R1 8.13e+10
err 3.13e-10 condest R1 4.23e+09
err 6.64e-10 condest R1 2.46e+10
err 5.76e-09 condest R1 4.31e+09
err 7.61e-07 condest R1 5.08e+10
err 2.27e-09 condest R1 4.94e+09
err 3.99e-10 condest R1 2.80e+09
err 1.37e-09 condest R1 3.13e+09
err 6.93e-05 condest R1 1.84e+14
err 1.35e-08 condest R1 7.18e+09
err 1.09e-08 condest R1 1.79e+11
err 1.81e-09 condest R1 2.99e+08
err 1.55e-01 condest R1 2.45e+18

In summary, this is a "feature" not a "bug".  If you want a reliable
solution to an underdetermined system, find the min 2norm solution
via a QR factorization of A'.
```

 Ordering statistics: result nnz(V) for QR, upper bound nnz(L) for LU, with COLAMD 15,330 nnz(R) for QR, upper bound nnz(U) for LU, with COLAMD 15,836

 SVD-based statistics: norm(A) 9.33232 min(svd(A)) 0.0320586 cond(A) 291.102 rank(A) 660 sprank(A)-rank(A) 0 null space dimension 0 full numerical rank? yes

 singular values (MAT file): click here SVD method used: s = svd (full (A)) ; status: ok