[mlpack-svn] [MLPACK] #373: LARSTest/NoCholeskySingularityTest often fails for i386

[MLPACK Trac]

#373: LARSTest/NoCholeskySingularityTest often fails for i386
---------------------+------------------------------------------------------
 Reporter:  rcurtin  |        Owner:                                          
     Type:  defect   |       Status:  new                                     
 Priority:  minor    |    Milestone:  mlpack 1.1.0                            
Component:  mlpack   |     Keywords:  lapack, armadillo, lars, solve, singular
 Blocking:           |   Blocked By:                                          
---------------------+------------------------------------------------------
 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.

-- 
Ticket URL: <http://trac.research.cc.gatech.edu/fastlab/ticket/373>
MLPACK <www.fast-lab.org>
MLPACK is an intuitive, fast, and scalable C++ machine learning library developed at Georgia Tech.