Repository: sympy/sympy
Pull request: fix(matrices): Don't use DomainMatrix.rref for RR in linsolve
The DomainMatrix.rref method as implemented by the sdm_irref function
has not been well designed for handling matrices of floats as it does
not have a pivoting strategy that can mitigate errors arising from
erroneously choosing a zero pivot (that is nonzero only because of
roundoff). This commit makes linsolve use the `ddm_irref()` for the case of
pure RR or CC matrices. The behaviour of sdm_rref for floats should be
improved. Until then it is better not to use `sdm_rref` in these
cases.

<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->

#### References to other Issues or PRs
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234" (see
https://tinyurl.com/auto-closing for more information). Also, please
write a comment on that issue linking back to this pull request once it is
open. -->

Fixes #21520

#### Brief description of what is fixed or changed


#### Other comments


#### Release Notes

<!-- Write the release notes for this release below between the BEGIN and END
statements. The basic format is a bulleted list with the name of the subpackage
and the release note for this PR. For example:

* solvers
  * Added a new solver for logarithmic equations.

* functions
  * Fixed a bug with log of integers.

or if no release note(s) should be included use:

NO ENTRY

See https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more
information on how to write release notes. The bot will check your release
notes automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->
* solvers
    * Some bugs in the new implementation of linsolve are fixed by the dense RREF implementation instead of the sparse implemention. These bugs previously led to incorrect results being returned linsolve.
<!-- END RELEASE NOTES -->


Problem statement:
linsolve fails simple system of two equations
```
import sympy
x,y = sympy.symbols('x, y')

sympy.linsolve([sympy.Eq(y, x), sympy.Eq(y, 0.0215 * x)], (x, y))
>> FiniteSet((0, 0))

sympy.linsolve([sympy.Eq(y, x), sympy.Eq(y, 0.0216 * x)], (x, y))
>> FiniteSet((-4.07992766242527e+17*y, 1.0*y))

sympy.linsolve([sympy.Eq(y, x), sympy.Eq(y, 0.0217 * x)], (x, y))
>> FiniteSet((0, 0))
```

Any thoughts on why these don't all return the same solution? Thanks!

Hints from issue discussion:
It seems that in rref the pivot is not fully cancelled due to a rounding error so e.g. we have something like:
```python
In [1]: M = Matrix([[1.0, 1.0], [3.1, 1.0]])

In [2]: M
Out[2]: 
⎡1.0  1.0⎤
⎢        ⎥
⎣3.1  1.0⎦
```
Then one step of row reduction gives:
```python
In [3]: M = Matrix([[1.0, 1.0], [1e-16, -2.1]])

In [4]: M
Out[4]: 
⎡  1.0    1.0 ⎤
⎢             ⎥
⎣1.0e-16  -2.1⎦
```
With exact arithmetic the 1e-16 would have been 0 but the rounding error makes it not zero and then it throws off subsequent steps. I think that the solution is to make it exactly zero:
```diff
diff --git a/sympy/polys/matrices/sdm.py b/sympy/polys/matrices/sdm.py
index cfa624185a..647eb6af3d 100644
--- a/sympy/polys/matrices/sdm.py
+++ b/sympy/polys/matrices/sdm.py
@@ -904,6 +904,8 @@ def sdm_irref(A):
             Ajnz = set(Aj)
             for k in Ajnz - Ainz:
                 Ai[k] = - Aij * Aj[k]
+            Ai.pop(j)
+            Ainz.remove(j)
             for k in Ajnz & Ainz:
                 Aik = Ai[k] - Aij * Aj[k]
                 if Aik:
```
That gives:
```python
In [1]: import sympy

In [2]: sympy.linsolve([sympy.Eq(y, x), sympy.Eq(y, 0.0215 * x)], (x, y))
Out[2]: {(0, 0)}

In [3]: sympy.linsolve([sympy.Eq(y, x), sympy.Eq(y, 0.0216 * x)], (x, y))
Out[3]: {(0, 0)}

In [4]: sympy.linsolve([sympy.Eq(y, x), sympy.Eq(y, 0.0217 * x)], (x, y))
Out[4]: {(0, 0)}
```
@tyler-herzer-volumetric can you try out that diff?
Haven't been able to replicate the issue after making that change. Thanks a lot @oscarbenjamin!
This example still fails with the diff:
```python
In [1]: linsolve([0.4*x + 0.3*y + 0.2, 0.4*x + 0.3*y + 0.3], [x, y])
Out[1]: {(1.35107988821115e+15, -1.8014398509482e+15)}
```
In this case although the pivot is set to zero actually the matrix is singular but row reduction leads to something like:
```python
In [3]: Matrix([[1, 1], [0, 1e-17]])
Out[3]: 
⎡1     1   ⎤
⎢          ⎥
⎣0  1.0e-17⎦
```
That's a trickier case. It seems that numpy can pick up on it e.g.:
```python
In [52]: M = np.array([[0.4, 0.3], [0.4, 0.3]])

In [53]: b = np.array([0.2, 0.3])

In [54]: np.linalg.solve(M, b)
---------------------------------------------------------------------------
LinAlgError: Singular matrix 
```
A slight modification or rounding error leads to the same large result though:
```python
In [55]: M = np.array([[0.4, 0.3], [0.4, 0.3-3e-17]])

In [56]: np.linalg.solve(M, b)
Out[56]: array([ 1.35107989e+15, -1.80143985e+15])
```
I'm not sure it's possible to arrange the floating point calculation so that cases like this are picked up as being singular without introducing some kind of heuristic threshold for the determinant. This fails in numpy with even fairly simple examples:
```python
In [83]: b
Out[83]: array([0.2, 0.3, 0.5])

In [84]: M
Out[84]: 
array([[0.1, 0.2, 0.3],
       [0.4, 0.5, 0.6],
       [0.7, 0.8, 0.9]])

In [85]: np.linalg.solve(M, b)
Out[85]: array([-4.50359963e+14,  9.00719925e+14, -4.50359963e+14])

In [86]: np.linalg.det(M)
Out[86]: 6.661338147750926e-18
```
Maybe the not full-rank case for float matrices isn't so important since it can't be done reliably with floats. I guess that the diff shown above is good enough then since it fixes the calculation in the full rank case.
There is another problematic case. This should have a unique solution but a parametric solution is returned instead:
```python
In [21]: eqs = [0.8*x + 0.8*z + 0.2, 0.9*x + 0.7*y + 0.2*z + 0.9, 0.7*x + 0.2*y + 0.2*z + 0.5]

In [22]: linsolve(eqs, [x, y, z])
Out[22]: {(-0.32258064516129⋅z - 0.548387096774194, 1.22033022161007e+16⋅z - 5.37526407137769e+15, 1.0⋅z)}
```
That seems to be another bug in the sparse rref routine somehow:
```python
In [34]: M = Matrix([
    ...: [0.8,   0, 0.8, -0.2],
    ...: [0.9, 0.7, 0.2, -0.9],
    ...: [0.7, 0.2, 0.2, -0.5]])

In [35]: from sympy.polys.matrices import DomainMatrix

In [36]: dM = DomainMatrix.from_Matrix(M)

In [37]: M.rref()
Out[37]: 
⎛⎡1  0  0  -0.690476190476191⎤           ⎞
⎜⎢                           ⎥           ⎟
⎜⎢0  1  0  -0.523809523809524⎥, (0, 1, 2)⎟
⎜⎢                           ⎥           ⎟
⎝⎣0  0  1  0.440476190476191 ⎦           ⎠

In [38]: dM.rref()[0].to_Matrix()
Out[38]: 
⎡1.0  5.55111512312578e-17    0.32258064516129      -0.548387096774194  ⎤
⎢                                                                       ⎥
⎢0.0          1.0           -1.22033022161007e+16  -5.37526407137769e+15⎥
⎢                                                                       ⎥
⎣0.0          0.0                    0.0                    0.0         ⎦

In [39]: dM.to_dense().rref()[0].to_Matrix()
Out[39]: 
⎡1.0  0.0  0.0  -0.69047619047619 ⎤
⎢                                 ⎥
⎢0.0  1.0  0.0  -0.523809523809524⎥
⎢                                 ⎥
⎣0.0  0.0  1.0   0.44047619047619 ⎦
```
The last one was a similar problem to do with cancelling above the pivot:
```diff
diff --git a/sympy/polys/matrices/sdm.py b/sympy/polys/matrices/sdm.py
index cfa624185a..7c4ad43660 100644
--- a/sympy/polys/matrices/sdm.py
+++ b/sympy/polys/matrices/sdm.py
@@ -904,6 +904,8 @@ def sdm_irref(A):
             Ajnz = set(Aj)
             for k in Ajnz - Ainz:
                 Ai[k] = - Aij * Aj[k]
+            Ai.pop(j)
+            Ainz.remove(j)
             for k in Ajnz & Ainz:
                 Aik = Ai[k] - Aij * Aj[k]
                 if Aik:
@@ -938,6 +940,8 @@ def sdm_irref(A):
             for l in Ainz - Aknz:
                 Ak[l] = - Akj * Ai[l]
                 nonzero_columns[l].add(k)
+            Ak.pop(j)
+            Aknz.remove(j)
             for l in Ainz & Aknz:
                 Akl = Ak[l] - Akj * Ai[l]
                 if Akl:
diff --git a/sympy/polys/matrices/tests/test_linsolve.py b/sympy/polys/matrices/tests/test_linsolve.py
index eda4cdbdf3..6b79842fa7 100644
--- a/sympy/polys/matrices/tests/test_linsolve.py
+++ b/sympy/polys/matrices/tests/test_linsolve.py
@@ -7,7 +7,7 @@
 from sympy.testing.pytest import raises
 
 from sympy import S, Eq, I
-from sympy.abc import x, y
+from sympy.abc import x, y, z
 
 from sympy.polys.matrices.linsolve import _linsolve
 from sympy.polys.solvers import PolyNonlinearError
@@ -23,6 +23,14 @@ def test__linsolve():
     raises(PolyNonlinearError, lambda: _linsolve([x*(1 + x)], [x]))
 
 
+def test__linsolve_float():
+    assert _linsolve([Eq(y, x), Eq(y, 0.0216 * x)], (x, y)) == {x:0, y:0}
+
+    eqs = [0.8*x + 0.8*z + 0.2, 0.9*x + 0.7*y + 0.2*z + 0.9, 0.7*x + 0.2*y + 0.2*z + 0.5]
+    sol = {x:-0.69047619047619047, y:-0.52380952380952395, z:0.44047619047619047}
+    assert _linsolve(eqs, [x,y,z]) == sol
+
+
 def test__linsolve_deprecated():
     assert _linsolve([Eq(x**2, x**2+y)], [x, y]) == {x:x, y:S.Zero}
     assert _linsolve([(x+y)**2-x**2], [x]) == {x:-y/2}
```
Another problem has emerged though:
```python
In [1]: eqs = [0.9*x + 0.3*y + 0.4*z + 0.6, 0.6*x + 0.9*y + 0.1*z + 0.7, 0.4*x + 0.6*y + 0.9*z + 0.5]

In [2]: linsolve(eqs, [x, y, z])
Out[2]: {(-0.5, -0.4375, -0.0400000000000001)}

In [3]: solve(eqs, [x, y, z])
Out[3]: {x: -0.502857142857143, y: -0.438095238095238, z: -0.04}
```
> Another problem has emerged though:
> 
> ```python
> In [1]: eqs = [0.9*x + 0.3*y + 0.4*z + 0.6, 0.6*x + 0.9*y + 0.1*z + 0.7, 0.4*x + 0.6*y + 0.9*z + 0.5]
> ```

In this case the problem is that something like `1e-17` is chosen (incorrectly) as a pivot. It doesn't look easy to resolve this because the `sdm_rref` routine doesn't have an easy way of incorporating pivoting and is really designed for exact domains.

Probably float matrices should be handled by mpmath or at least a separate routine.

Commit message: fix(matrices): Don't use DomainMatrix.rref for RR in linsolve

The DomainMatrix.rref method as implemented by the sdm_irref function
has not been well designed for handling matrices of floats as it does
not have a pivoting strategy that can mitigate errors arising from
erroneously choosing a zero pivot (that is nonzero only because of
roundoff). This commit makes linsolve use Matrix.rref() for the case of
pure float matrices. The behaviour of sdm_rref for floats should be
improved. Until then it is better not to use DomainMatrix in these
cases.

Patch to review:
diff --git a/sympy/polys/matrices/linsolve.py b/sympy/polys/matrices/linsolve.py
index 62e6d90aae36..387ea9c79cc6 100644
--- a/sympy/polys/matrices/linsolve.py
+++ b/sympy/polys/matrices/linsolve.py
@@ -35,6 +35,7 @@
 from sympy.polys.constructor import construct_domain
 from sympy.polys.solvers import PolyNonlinearError
 
+from .domainmatrix import DomainMatrix
 from .sdm import (
     SDM,
     sdm_irref,
@@ -73,6 +74,13 @@ def _linsolve(eqs, syms):
     Aaug = sympy_dict_to_dm(eqsdict, rhs, syms)
     K = Aaug.domain
 
+    # sdm_irref has issues with float matrices. This uses the Matrix.rref()
+    # method. When DomainMatrix.rref() can handle float matrices reasonably
+    # this should be removed...
+    if K.is_RealField:
+        Aaug_Matrix = DomainMatrix.from_rep(Aaug).to_Matrix().rref()[0]
+        Aaug = DomainMatrix.from_Matrix(Aaug_Matrix).to_field().to_sparse().rep
+
     # Compute reduced-row echelon form (RREF)
     Arref, pivots, nzcols = sdm_irref(Aaug)
 
diff --git a/sympy/polys/matrices/sdm.py b/sympy/polys/matrices/sdm.py
index cfa624185a5d..7c4ad43660cc 100644
--- a/sympy/polys/matrices/sdm.py
+++ b/sympy/polys/matrices/sdm.py
@@ -904,6 +904,8 @@ def sdm_irref(A):
             Ajnz = set(Aj)
             for k in Ajnz - Ainz:
                 Ai[k] = - Aij * Aj[k]
+            Ai.pop(j)
+            Ainz.remove(j)
             for k in Ajnz & Ainz:
                 Aik = Ai[k] - Aij * Aj[k]
                 if Aik:
@@ -938,6 +940,8 @@ def sdm_irref(A):
             for l in Ainz - Aknz:
                 Ak[l] = - Akj * Ai[l]
                 nonzero_columns[l].add(k)
+            Ak.pop(j)
+            Aknz.remove(j)
             for l in Ainz & Aknz:
                 Akl = Ak[l] - Akj * Ai[l]
                 if Akl:
diff --git a/sympy/polys/matrices/tests/test_linsolve.py b/sympy/polys/matrices/tests/test_linsolve.py
index eda4cdbdf3fa..354801a83816 100644
--- a/sympy/polys/matrices/tests/test_linsolve.py
+++ b/sympy/polys/matrices/tests/test_linsolve.py
@@ -7,7 +7,7 @@
 from sympy.testing.pytest import raises
 
 from sympy import S, Eq, I
-from sympy.abc import x, y
+from sympy.abc import x, y, z
 
 from sympy.polys.matrices.linsolve import _linsolve
 from sympy.polys.solvers import PolyNonlinearError
@@ -23,6 +23,33 @@ def test__linsolve():
     raises(PolyNonlinearError, lambda: _linsolve([x*(1 + x)], [x]))
 
 
+def test__linsolve_float():
+    eqs = [
+        y - x,
+        y - 0.0216 * x
+    ]
+    sol = {x:0, y:0}
+    assert _linsolve(eqs, (x, y)) == sol
+
+    eqs = [
+        0.8*x +         0.8*z + 0.2,
+        0.9*x + 0.7*y + 0.2*z + 0.9,
+        0.7*x + 0.2*y + 0.2*z + 0.5
+    ]
+    #sol = {x:-0.69047619047619047, y:-0.52380952380952395, z:0.44047619047619047}
+    sol = {x:-0.69047619047619069, y:-0.52380952380952361, z:0.44047619047619069}
+    assert _linsolve(eqs, [x,y,z]) == sol
+
+    eqs = [
+        0.9*x + 0.3*y + 0.4*z + 0.6,
+        0.6*x + 0.9*y + 0.1*z + 0.7,
+        0.4*x + 0.6*y + 0.9*z + 0.5
+    ]
+    #sol = {x:-0.50285714285714289, y:-0.43809523809523804, z:-0.040000000000000063}
+    sol = {x:-0.50285714285714278, y:-0.43809523809523798, z:-0.040000000000000022}
+    assert _linsolve(eqs, [x,y,z]) == sol
+
+
 def test__linsolve_deprecated():
     assert _linsolve([Eq(x**2, x**2+y)], [x, y]) == {x:x, y:S.Zero}
     assert _linsolve([(x+y)**2-x**2], [x]) == {x:-y/2}


Code review comments:
diff_hunk: @@ -73,6 +74,13 @@ def _linsolve(eqs, syms):
     Aaug = sympy_dict_to_dm(eqsdict, rhs, syms)
     K = Aaug.domain
 
+    # sdm_irref has issues with float matrices. This uses the Matrix.rref()
+    # method. When DomainMatrix.rref() can handle float matrices reasonably
+    # this should be removed...
+    if K.is_RealField:
path: sympy/polys/matrices/linsolve.py
text: @user1:
Complex field should treated similarly.

@author:
Yes, it should although it isn't actually used here since it always goes to EX:
```python
In [2]: construct_domain([1.0+1.0*I])
Out[2]: (EX, [EX(1.0 + 1.0*I)])

In [3]: Poly(1.0+1.0*I, x)
Out[3]: Poly(1.0 + 1.0*I, x, domain='EX')
```
I think the only way to get CC is by unifying `ZZ_I/QQ_I` with `RR`:
```python
In [4]: ZZ_I.unify(RR)
Out[4]: CC
```

@author:
I've changed it so that `construct_domain` can identify `CC`. I would have expected to see some test failures from that but there don't seem to be any which means that there is not good test coverage for Poly/linsolve etc over CC.

It should be possible to test `if K.is_Exact` here but that apparently fails for a polynomial ring:
```python
In [1]: RR.frac_field(x).is_Exact
Out[1]: False

In [2]: RR[x].is_Exact
Out[2]: True
```