lu_solve should detect numerically singular matrices

[GoogleCodeExporter]

>>> A = matrix([[5.6, 1.2], [7./15, .1]])
>>> det(A)
mpf('7.3617327511765746e-18')
>>> cond(A)
mpf('5.6037531825289503e+18')
>>> lu_solve(A,[1,2])
matrix(
[['-3.12426446020118e+17'],
 ['1.45799008142722e+18']])
>>> qr_solve(A,[1,2])
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/vinzent/src/mpmath/mpmath/settings.py", line 135, in g
    return f(*args, **kwargs)
  File "mpmath/linalg.py", line 272, in qr_solve
    H, p, x, r = householder(extend(A, b))
  File "mpmath/linalg.py", line 202, in householder
    raise ValueError('matrix is numerically singular')
ValueError: matrix is numerically singular

Original issue reported on code.google.com by Vinzent.Steinberg@gmail.com on 15 Oct 2008 at 7:38

[GoogleCodeExporter]

Or qr_solve is maybe to strict.

Original comment by Vinzent.Steinberg@gmail.com on 15 Oct 2008 at 7:48

[GoogleCodeExporter]

This should be fixed too.

>>> lu_solve(zeros(2),[0,0])
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/vinzent/src/mpmath/mpmath/settings.py", line 135, in g
    return f(*args, **kwargs)
  File "mpmath/linalg.py", line 110, in lu_solve
    A, p = LU_decomp(A)
  File "mpmath/linalg.py", line 37, in LU_decomp
    * absmin(A[k,j])
  File "/home/vinzent/src/mpmath/mpmath/mptypes.py", line 198, in __rdiv__
    return make_mpf(mpf_rdiv_int(t, s._mpf_, prec, rounding))
  File "mpmath/libmpf.py", line 802, in mpf_rdiv_int
    return mpf_div(from_int(n), t, prec, rnd)
  File "mpmath/libmpf.py", line 761, in mpf_div
    raise ZeroDivisionError
ZeroDivisionError

Original comment by Vinzent.Steinberg@gmail.com on 15 Oct 2008 at 7:50

[GoogleCodeExporter]

The function LU_decomp has several issues:

 - 'tol' can be 0;
 - 'fsum([absmin(A[k,l]) for l in xrange(j, n)])' can be 0;
 - A[n-1,n-1] remains unchecked even though it is assumed non-zero in U_solve.

The attached patch fixes these (in a somewhat crude way).

Original comment by jorn.baa...@gmail.com on 24 Mar 2009 at 10:35

Attachments:

• linalg.py.patch

[GoogleCodeExporter]

Adding these checks seems fine to me. Vinzent, can you review the patch (and 
commit
if you have time)?

Original comment by fredrik....@gmail.com on 25 Mar 2009 at 9:19

[GoogleCodeExporter]

Thank you a lot, I committed your patch and some tests.

Should singular matrices raise a ZeroDivisionError or a ValueError?

In fact, raising a ZeroDivisionError was a workaround for det(), because 
LU_decomp()
could not recognize singular matrices and raised such an error. This is fixed 
now.

However, trying to inverse a singular matrix is analog to scalar zero division 
(it's
like trying to inverse 0), so maybe ZeroDivisionError is more accurate.

Original comment by Vinzent.Steinberg@gmail.com on 25 Mar 2009 at 3:07

[GoogleCodeExporter]

Jorn Baayen, I credited you in mpmath svn, please let me know if I did not 
spell your
name correctly.

Original comment by Vinzent.Steinberg@gmail.com on 11 Jun 2009 at 8:47

[GoogleCodeExporter]

Also, please tell me your email address (if you want it in mpmath's README).

Original comment by Vinzent.Steinberg@gmail.com on 11 Jun 2009 at 8:50