# Ticket #2932(closed defect: fixed)

Opened 5 years ago

## [with patch, positive review] matrix.is_invertible() has inconsisten behavior over CDF

Reported by: Owned by: dfdeshom broune major sage-3.0.4 linear algebra editor_malb

### Description

From Alex Ghitza:

```sage: M = matrix(CDF, 2, 2, [(-1 - 2*I, 5 - 6*I), (-2 - 4*I, 10 - 12*I)])
sage: M.is_invertible()
True
sage: M.determinant()
5.3290705182e-15 + 1.7763568394e-15*I
sage: M.inverse()
[ 1.01330991616e+15 - 2.58956978574e+15*I -5.06654958079e+14 +
1.29478489287e+15*I]
[ 5.62949953421e+14 + 5.62949953421e+14*I -2.81474976711e+14 -
2.81474976711e+14*I]

So because of roundoff errors, Sage thinks that we have an invertible
matrix.  But the code for echelon_form knows that it's not invertible:

sage: M.echelon_form()
[                                    1.0                             1.4
+ 3.2*I]
[-2.22044604925e-16 - 4.4408920985e-16*I
~       0]
sage: M.rank()
1
```

## Change History

### comment:1 Changed 5 years ago by jason

Using numpy gives the following:

```sage: M = matrix(CDF, 2, 2, [(-1 - 2*I, 5 - 6*I), (-2 - 4*I, 10 - 12*I)])
sage: n=M.numpy()
sage: import numpy
sage: numpy.linalg.det(n)
0.0
sage: numpy.linalg.inv(n)
---------------------------------------------------------------------------
<class 'numpy.linalg.linalg.LinAlgError'> Traceback (most recent call last)

/home/grout/sage/devel/sage-main/sage/misc/<ipython console> in <module>()

/home/grout/sage/local/lib/python2.5/site-packages/numpy/linalg/linalg.py in inv(a)
244 def inv(a):
245     a, wrap = _makearray(a)
--> 246     return wrap(solve(a, identity(a.shape[0], dtype=a.dtype)))
247
248 # Cholesky decomposition

/home/grout/sage/local/lib/python2.5/site-packages/numpy/linalg/linalg.py in solve(a, b)
187     results = lapack_routine(n_eq, n_rhs, a, n_eq, pivots, b, n_eq, 0)
188     if results['info'] > 0:
--> 189         raise LinAlgError, 'Singular matrix'
190     if one_eq:
191         return b.ravel().astype(result_t)

<class 'numpy.linalg.linalg.LinAlgError'>: Singular matrix
```

### comment:2 Changed 5 years ago by broune

• Owner changed from was to broune
• Priority changed from critical to major
• Summary changed from matrix.is_invertible() has inconsisten behavior over CDF to [with patch, needs review] matrix.is_invertible() has inconsisten behavior over CDF

I propose that this code is working as intended - it cannot be expected that floating point calculations always provide the same results as would be obtained using precise calculations. It cannot either be assumed that two calculations that are equivalent using precise numbers are equivalent when using approximations. The attached patch adds docstring explanations of this to real and complex number classes and fields.

Other than that, I get this using Sage version 3.0:

```sage: M = matrix(CDF, 2, 2, [(-1 - 2*I, 5 - 6*I), (-2 - 4*I, 10 - 12*I)])
sage: M.is_invertible()
False
```

### comment:3 Changed 5 years ago by broune

• Status changed from new to assigned

### comment:4 Changed 5 years ago by jason

I am curious what algorithms are used for is_invertible and inverse(). Apparently determinant() uses GSL; how does GSL compare with numpy generally for numeric algorithms? Numpy uses the standard blas libraries an awful lot; what does GSL do?

If is_invertible and inverse() do not use standard numerical linear algebra functions (e.g., use ATLAS), then why not? If it's an issue of NotImplemented?, then let's open a ticket. I would probably get to it by sometime this summer, since we have a numerical linear algebra workshop where I'm at.

### comment:5 Changed 5 years ago by jason

From IRC (with editing out of other unrelated conversation):

```[20:15] <jason--> re: 2932; what algorithms are used for numerical linear algebra over CDF for inverse() and is_invertible()?
[20:15] <jason--> I don't find those functions in matrix_complex_double_dense.pyx
[20:16] <jason--> are they the generic algorithms in matrix*.pyx?
[20:16] <wstein-3053> If they aren't in matrix_complex_double_dense then they are the generic ones.
[20:17] <jason--> thanks; also, does anyone have any feelings for GSL compared to numpy for numerical linear algebra?
[20:17] <wstein-3053> numpy is usually better.
[20:18] <jason--> okay; do you mind if we switch the determinant function to numpy for matrices over CDF?
[20:19] <jason--> okay; ticket coming right up.
```

### comment:7 Changed 5 years ago by craigcitro

• Keywords editor_malb added; editor_ncalexan removed

### comment:8 Changed 5 years ago by malb

What is the current status of this ticket? I just to over being the editor for this ticket so please speak up :-)

### comment:9 Changed 5 years ago by jason

I just got back into town (been gone for a month). I'd still like to make numpy the default backend for CDF and RDF and I'll probably work on that this week.

### comment:10 Changed 5 years ago by jason

• Summary changed from [with patch, needs review] matrix.is_invertible() has inconsisten behavior over CDF to [with patch, positive review] matrix.is_invertible() has inconsisten behavior over CDF

Since the patch is just documentation telling the user not to expect exact results, I like it. +1 from me. Even when I convert things to numpy, it would still be good to have such a disclaimer, I think.

### comment:11 Changed 5 years ago by jason

See ticket #3498 for switching the backend to numpy.

### comment:12 Changed 5 years ago by mabshoff

• Status changed from assigned to closed
• Resolution set to fixed
• Milestone changed from sage-3.1.1 to sage-3.0.4

Merged in Sage 3.0.4.alpha1

Note: See TracTickets for help on using tickets.