Issue 86 - mpmath - lu_solve should detect numerically singular matrices

Issue 86:	lu_solve should detect numerically singular matrices	‹ Prev 8 of 75 Next ›
	2 people starred this issue and may be notified of changes.	Back to list

Status:	Accepted
Owner:	Vinzent.Steinberg
Type-Defect
Priority-Medium

Reported by Vinzent.Steinberg, Oct 15, 2008

>>> 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

Comment 1 by Vinzent.Steinberg, Oct 15, 2008

Or qr_solve is maybe to strict.

Comment 2 by Vinzent.Steinberg, Oct 15, 2008

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

Comment 3 by jorn.baayen, Mar 24, 2009

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).

linalg.py.patch
1.3 KB Download

Comment 4 by fredrik.johansson, Mar 25, 2009

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

Comment 5 by Vinzent.Steinberg, Mar 25, 2009

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.

Comment 6 by Vinzent.Steinberg, Jun 11, 2009

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

Comment 7 by Vinzent.Steinberg, Jun 11, 2009

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