Reported by rcurtin on 10 Apr 44902838 02:43 UTC
For some reason I haven't nailed down, on i386, NoCholeskySingularityTest often fails despite the fact that the code should work fine. In my investigations, I have discovered that the reason for this is that the test will try to solve a singular matrix of the form

[a
 a a](a)

which will generally cause solve() to fail (at least on x86_64). This is done in lars.cpp at line 198. But you can still solve a singular linear system, and MATLAB will do this in the least-squares sense:

http://www.mathworks.com/help/matlab/math/systems-of-linear-equations.html

A bit of backtracing and digging reveals that what is happening is this:

• LARS passes a singular matrix of the form above to arma::solve()
• arma::solve() detects that the matrix is very small and tries to invert it by hand, but gives up when it notices the determinant is 0
• arma::solve() then calls out to arma::lapack::gesv() (that is, LAPACK dgesv()) which provides a solution to the system (alternately, atlas::clapack_gesv() may be called, depending on the system configuration)

But this is confusing to me because LAPACK dgesv() should fail:

http://www.netlib.org/lapack/explore-html/d8/d72/dgesv_8f.html (search 'singular')

dgesv() will perform an LU factorization on the input matrix and then fail if U is singular... but U should be singular.

So, there are questions I think are worth asking:

• Is solve() a good test for matrix singularity?
• Why is LAPACK dgesv() solving the matrix despite its singularity, but only on i386? (I've never seen this failure on x86_64.)

Answering these questions, or, at least the second, should be possible with a little bit of code modification in LAPACK to figure out what is going on. If the answer to the first question is no, then digging deeper into the dgesv() issue may be completely unnecessary.

This should be reproducible by calling arma::solve() with a matrix of the form above, and may end up in a bug report bubbling up to Armadillo or LAPACK.

For now, I have commented out the test.

Michael, I CC'ed you because you wrote this code, and may have a more correct perspective on what's going on.

Migrated-From: http://trac.research.cc.gatech.edu/fastlab/ticket/373