13
$\begingroup$

I currently trying to cheaply compute a good rank estimate for a matrix $A$. Therefore I compute a columnt pivoting QR decompostion using 

[Q,R,E]=qr(A)

in Matlab. I estimate the rank of $A$ using 

tol = size(A,n)*eps*norm(A,'fro'); 
r = sum(abs(diag(R))>tol)

This works fine and a plot over all diagonal entries of R looks like: plot(sort(abs(diag(R)),1,'descend'),'r*')

If port the whole algorithm to C/Fortran I replace [Q,R,E]=qr(A) using DGEQP3 from LAPACK, which also computes a column pivoting QR decomposition. But if I use the same estimate for the rank I mostly get something wrong. The same plot for the $R$ produced from DGEQP3 looks like

The input matrix is exactly the same for both experiments.

My question now is on which LAPACK function does the column pivoting QR decomposition from Matlab rely on?

Thanks for any help, Grisu

Edit: DGEQPF gives the same wrong result.

Edit2:

Edit3: - Using GDB I found out, that Matlab 2010b calls DGEQP3: #3 0xaa46ce2f in dgeqp3_ () from /usr/ubuntu10.04/matlabr2010b/bin/glnx86/../../bin/glnx86/../../bin/glnx86/mllapack.so Why do I get the wrong result using LAPACK(3.4.0 inlcude the fixes mentioned in Working Note 176)?

Improve this question
$\endgroup$
5
  • $\begingroup$ Can you provoke the same behaviour with a smaller matrix that you might be able to share here? $\endgroup$
    – GertVdE
    Mar 27, 2012 at 14:48
  • $\begingroup$ Does $A$ have any special structure? I want to say that MATLAB uses UMFPACK for sparse linear algebra but I'm not sure what the underlying dense linear algebra libraries are. $\endgroup$
    – Aron Ahmadia
    Mar 27, 2012 at 14:52
  • $\begingroup$ Is $A$ sparse? You can set ">> spparms('spumoni',1)" and see that something called "SuiteSparseQR" gets used by matlab in that case. $\endgroup$
    – dranxo
    Mar 27, 2012 at 16:11
  • $\begingroup$ Grisu - I would love to look at your matrix. However, the link www-e.uni-magdeburg.de/makoehle/A.mtx.gz is no longer active (at the current time, anyway). Do you have a current link to the matrix? Thanks, Les Foster $\endgroup$
    – user1723
    Jul 6, 2012 at 7:04
  • $\begingroup$ @LeslieFoster - welcome to scicomp! $\endgroup$
    – Aron Ahmadia
    Jul 6, 2012 at 8:41

2 Answers 2

Reset to default
8
$\begingroup$

There are two issues at hand here:

  • Is $A$ dense or sparse?
  • Do you have the same software stack as MATLAB's internal libraries?

Dense or sparse?

MATLAB no longer explicitly mentions the LAPACK routines it calls to obtain a QR factorization if $A$ is dense. If the information in the documentation for MATLAB R2008b also holds for later releases, then MATLAB calls DGEQP3 from LAPACK when you call [Q,R,E] = qr(A). If $A$ is sparse, then MATLAB calls SuiteSparseQR, out of Tim Davis' group, which is bundled with UMFPACK in the SuiteSparse library.

Do you have the same software stack as MATLAB's internal libraries?

Probably not, which may be one reason you're getting different results.

I ran into this issue when unit testing a library I was writing that used QR factorizations. I used MATLAB to prototype my work and got different results than using LAPACK or NumPy. As far as I can tell, because MathWorks doesn't make this information easy to find, MATLAB uses a verson of LAPACK no earlier than version 3.1.1, and Intel's MKL BLAS library (for Windows, Intel Mac, and Linux) version 9.1 or higher (see here). I couldn't find anything about the version of SuiteSparse MATLAB uses. By digging around online or looking at the library files for your system, you may be able to glean additional information. You could try changing the libraries that MATLAB links to in order to be able to compare with the same libraries across software packages; Eric Chu provides a nice write-up that shows at least how you can replace MATLAB's BLAS library with your own (of course, you do this at your own risk). He suggests that you can do the same thing with LAPACK as well. It may even be possible to replace the version of SuiteSparse that MATLAB uses with your own version.

LAPACK versions change, and so do BLAS versions; they may use different algorithms from version to version, or different ordering conventions, even though $A = QR$ with $Q$ orthogonal and $R$ is upper triangular, regardless of version. These changes make reproducibility a challenge. Even the tolerance you use for your numerical rank determination is a judgment call; you appear to be using a standard tolerance.

I ended up using NumPy to prototype my results for QR factorization, because it uses the system BLAS and LAPACK libraries. NumPy and SciPy isn't a drop-in replacement for MATLAB, because the two libraries combined lack some of MATLAB's functionality, but for this particular linear algebra task, Python + NumPy + SciPy + Matplotlib should work well.

Improve this answer
$\endgroup$
6
  • $\begingroup$ Getting the same the software stack as Matlab is impossible I think. Using another environment for prototyping is also not wanted. The problem is the code works in Matlab correctly, but not in C. $\endgroup$
    – M.K. aka Grisu
    Mar 28, 2012 at 8:37
  • $\begingroup$ @Grisu: I think it'd be very difficult to get the same software stack, short of attempting to link in their libraries. What I'm confused about is how you know the result in MATLAB is correct and the result in C is wrong. Is this some sort of test matrix that has known properties? More to the point, AronAhmadia is right; beyond replicating the architecture and the software stack, you can't expect to get the same results with an unstable algorithm. I was basically told the same thing in the MATLAB forums two years ago. $\endgroup$
    – Geoff Oxberry
    Mar 28, 2012 at 17:57
  • $\begingroup$ In my opinion the QR is not unstable. I can not directly check which QR decomposition is right, but from the rank and the Q matrix I compute a projector and there I can easily verify if a get good or bad results and the one from Matlab are good. But I try linking against the Matlab libraries. $\endgroup$
    – M.K. aka Grisu
    Mar 29, 2012 at 8:19
  • $\begingroup$ @Grisu: There's a distinct difference between good results and correct results. I recently implemented a calculation incorrectly that made my results look wonderful. Nevertheless, the calculation was wrong, and the correct calculation made my results look less impressive (but thankfully, illustrates that my results are correct). Are you attempting to calculate an orthogonal projector or an oblique projector? (I ask because significant parts of my thesis are on oblique projectors.) $\endgroup$
    – Geoff Oxberry
    Mar 29, 2012 at 19:19
  • 1
    $\begingroup$ @GeoffOxberry: fwiw, on my version of MATLAB, I can call internal.matlab.language.versionPlugins.blas and internal.matlab.language.versionPlugins.lapack to get BLAS and LAPACK versions $\endgroup$
    – Amro
    Apr 26, 2013 at 0:59
8
$\begingroup$

See Leslie Foster's page on rank-revealing software. See also this LAPACK Working Note analyzing failures of rank-revealing QR xGEQP3.

You should be able to find out what routines MATLAB uses by setting breakpoints in a debugger and examining the stack. Last I looked, admittedly several years ago, MATLAB used shared libraries, in which case the symbol names cannot be stripped, so you will see the function names on the call stack (but not arguments because it definitely isn't keeping debugging information).

Improve this answer
$\endgroup$
3
  • $\begingroup$ The page on rank-revealing software did not help. The RRQR procedure described there was the first thing I use in my idea but it gives even worse results than the column pivoting idea does. $\endgroup$
    – M.K. aka Grisu
    Mar 28, 2012 at 10:17
  • 3
    $\begingroup$ @Grisu - It should have helped you. The xGEQP3 algorithm is not completely safe for revealing rank. If you want to guarantee that you get the right result, you should use the SVD or a safer QR such as xGEQPX or xGEQPY. You cannot expect an unstable algorithm to return the same result on different architectures or in different implementations (MATLAB is probably using an older LAPACK). $\endgroup$
    – Aron Ahmadia
    Mar 28, 2012 at 14:00
  • 1
    $\begingroup$ I know that the GEQP3 is not rank-revealing but, it gives more correct results than the RRQR subroutines do. Using an SVD is too expensive in my outer algorithm. I'll also talk to one of the authors of LAWN-176 and he thinks this error is not covered by the bug. I also tried DGEQPF/DGEQP3 from LAPACK 3.0.0 with the same results. $\endgroup$
    – M.K. aka Grisu
    Mar 28, 2012 at 14:26

Your Answer

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

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