Repository: sympy/sympy
Pull request: Add CRootsOf incrementally for eigenvalues
<!-- 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 #25282

#### 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. Formerly, `log(-x)` incorrectly gave `-log(x)`.

* physics.units
  * Corrected a semantical error in the conversion between volt and statvolt which
    reported the volt as being larger than the statvolt.

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 -->
- matrices
  - Fixed some corner cases when `Matrix.eigenvals` gives wrong result when `multiplicity=False`.
<!-- END RELEASE NOTES -->


Problem statement:
incorrect multiplicities from `eigenvals`
To reproduce:
```
import numpy as np
import sympy as sp

print(sp.__version__)

dd = sd = [0] * 11 + [1]
ds = [2,0,1,0,0,0,1,0,1,0,1,0]
ss = ds.copy()
ss[8] = 2

def rotate(x, i):
    return x[i:] + x[:i]

mat = []
for i in range(12):
    mat.append(rotate(ss,i) + rotate(sd,i))
for i in range(12):
    mat.append(rotate(ds,i) + rotate(dd,i))

print(np.array(mat))

mat = sp.Matrix(mat)
eigenvalues = mat.eigenvals()
print('eigenvalues from eigenvals:')
for i in eigenvalues:
    print(eigenvalues[i], 'x', i)
print('matrix dimensions:', mat.shape)
print('sum of multiplicities:', sum([eigenvalues[i] for i in eigenvalues]))

eigenvectors = mat.eigenvects()
print('eigenvalues from eigenvects (multiplicities appear to be correct):')
for val, mult, _ in eigenvectors:
    print(mult, 'x', val)
```

Output:
```
1.12
[[2 0 1 0 0 0 1 0 2 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1]
 [0 1 0 0 0 1 0 2 0 1 0 2 0 0 0 0 0 0 0 0 0 0 1 0]
 [1 0 0 0 1 0 2 0 1 0 2 0 0 0 0 0 0 0 0 0 0 1 0 0]
 [0 0 0 1 0 2 0 1 0 2 0 1 0 0 0 0 0 0 0 0 1 0 0 0]
 [0 0 1 0 2 0 1 0 2 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0]
 [0 1 0 2 0 1 0 2 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0]
 [1 0 2 0 1 0 2 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0]
 [0 2 0 1 0 2 0 1 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0]
 [2 0 1 0 2 0 1 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0]
 [0 1 0 2 0 1 0 0 0 1 0 2 0 0 1 0 0 0 0 0 0 0 0 0]
 [1 0 2 0 1 0 0 0 1 0 2 0 0 1 0 0 0 0 0 0 0 0 0 0]
 [0 2 0 1 0 0 0 1 0 2 0 1 1 0 0 0 0 0 0 0 0 0 0 0]
 [2 0 1 0 0 0 1 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1]
 [0 1 0 0 0 1 0 1 0 1 0 2 0 0 0 0 0 0 0 0 0 0 1 0]
 [1 0 0 0 1 0 1 0 1 0 2 0 0 0 0 0 0 0 0 0 0 1 0 0]
 [0 0 0 1 0 1 0 1 0 2 0 1 0 0 0 0 0 0 0 0 1 0 0 0]
 [0 0 1 0 1 0 1 0 2 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0]
 [0 1 0 1 0 1 0 2 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0]
 [1 0 1 0 1 0 2 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0]
 [0 1 0 1 0 2 0 1 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0]
 [1 0 1 0 2 0 1 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0]
 [0 1 0 2 0 1 0 0 0 1 0 1 0 0 1 0 0 0 0 0 0 0 0 0]
 [1 0 2 0 1 0 0 0 1 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0]
 [0 2 0 1 0 0 0 1 0 1 0 1 1 0 0 0 0 0 0 0 0 0 0 0]]
eigenvalues from eigenvals:
1 x CRootOf(lambda**8 - 10*lambda**6 + 15*lambda**4 - 10*lambda**2 + 1, 0)
1 x CRootOf(lambda**4 - 5*lambda**2 + 1, 0)
1 x -1
1 x CRootOf(lambda**4 - 5*lambda**2 + 1, 1)
1 x CRootOf(lambda**8 - 10*lambda**6 + 15*lambda**4 - 10*lambda**2 + 1, 1)
1 x CRootOf(lambda**2 - 6*lambda - 1, 0)
1 x CRootOf(lambda**2 - 8*lambda + 1, 0)
1 x CRootOf(lambda**8 - 10*lambda**6 + 15*lambda**4 - 10*lambda**2 + 1, 2)
1 x CRootOf(lambda**4 - 5*lambda**2 + 1, 2)
1 x 1
1 x CRootOf(lambda**4 - 5*lambda**2 + 1, 3)
1 x CRootOf(lambda**8 - 10*lambda**6 + 15*lambda**4 - 10*lambda**2 + 1, 3)
1 x CRootOf(lambda**2 - 6*lambda - 1, 1)
1 x CRootOf(lambda**2 - 8*lambda + 1, 1)
1 x CRootOf(lambda**8 - 10*lambda**6 + 15*lambda**4 - 10*lambda**2 + 1, 4)
1 x CRootOf(lambda**8 - 10*lambda**6 + 15*lambda**4 - 10*lambda**2 + 1, 5)
1 x CRootOf(lambda**8 - 10*lambda**6 + 15*lambda**4 - 10*lambda**2 + 1, 6)
1 x CRootOf(lambda**8 - 10*lambda**6 + 15*lambda**4 - 10*lambda**2 + 1, 7)
matrix dimensions: (24, 24)
sum of multiplicities: 18
eigenvalues from eigenvects (multiplicities appear to be correct):
2 x -1
2 x 1
2 x -sqrt(5/2 - sqrt(21)/2)
2 x sqrt(5/2 - sqrt(21)/2)
1 x 3 - sqrt(10)
1 x 3 + sqrt(10)
1 x 4 - sqrt(15)
1 x sqrt(15) + 4
2 x -sqrt(sqrt(21)/2 + 5/2)
2 x sqrt(sqrt(21)/2 + 5/2)
1 x -cos(atan(sqrt(-33 + 20*sqrt(3))/(5 - 2*sqrt(3)))/2) - I*sin(atan(sqrt(-33 + 20*sqrt(3))/(5 - 2*sqrt(3)))/2)
1 x -cos(atan(sqrt(-33 + 20*sqrt(3))/(5 - 2*sqrt(3)))/2) + I*sin(atan(sqrt(-33 + 20*sqrt(3))/(5 - 2*sqrt(3)))/2)
1 x cos(atan(sqrt(-33 + 20*sqrt(3))/(5 - 2*sqrt(3)))/2) - I*sin(atan(sqrt(-33 + 20*sqrt(3))/(5 - 2*sqrt(3)))/2)
1 x cos(atan(sqrt(-33 + 20*sqrt(3))/(5 - 2*sqrt(3)))/2) + I*sin(atan(sqrt(-33 + 20*sqrt(3))/(5 - 2*sqrt(3)))/2)
1 x -sqrt(sqrt(3) + 5/2 + sqrt(33 + 20*sqrt(3))/2)
1 x sqrt(sqrt(3) + 5/2 + sqrt(33 + 20*sqrt(3))/2)
1 x -sqrt(-sqrt(33 + 20*sqrt(3))/2 + sqrt(3) + 5/2)
1 x sqrt(-sqrt(33 + 20*sqrt(3))/2 + sqrt(3) + 5/2)
```
The multiplicities from `eigenvals` function only sum up to 18, while the matrix is of dimension 24. The output from `eigenvects` appears to be correct (and it's impressive that it can solve the 8th degree equation symbolically).

Hints from issue discussion:
I think that some roots are ignored from `roots`. I'm trying to investigate whether the matrix code is responsible for something
but `sum(roots(mat.charpoly()).values())` already behaves suspicious. With `multiple=True`, the result seems correct
Eigenvects are computable but quite slow though.

Commit message: Use counter

Patch to review:
diff --git a/sympy/matrices/eigen.py b/sympy/matrices/eigen.py
index ea39e83a940c..cc421077a3e8 100644
--- a/sympy/matrices/eigen.py
+++ b/sympy/matrices/eigen.py
@@ -258,7 +258,7 @@ def _eigenvals_dict(
             f = charpoly.as_expr()
             x = charpoly.gen
             try:
-                eigs = {CRootOf(f, x, idx): 1 for idx in range(degree)}
+                eigs = Counter(CRootOf(f, x, idx) for idx in range(degree))
             except NotImplementedError:
                 if error_when_incomplete:
                     raise MatrixError(eigenvals_error_message)
diff --git a/sympy/matrices/tests/test_eigen.py b/sympy/matrices/tests/test_eigen.py
index 75ec80bcdd55..e952d225f74e 100644
--- a/sympy/matrices/tests/test_eigen.py
+++ b/sympy/matrices/tests/test_eigen.py
@@ -705,3 +705,21 @@ def test_issue_20752():
     b = symbols('b', nonzero=True)
     m = Matrix([[0, 0, 0], [0, b, 0], [0, 0, b]])
     assert m.is_positive_semidefinite is None
+
+
+def test_issue_25282():
+    dd = sd = [0] * 11 + [1]
+    ds = [2, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0]
+    ss = ds.copy()
+    ss[8] = 2
+
+    def rotate(x, i):
+        return x[i:] + x[:i]
+
+    mat = []
+    for i in range(12):
+        mat.append(rotate(ss, i) + rotate(sd, i))
+    for i in range(12):
+        mat.append(rotate(ds, i) + rotate(dd, i))
+
+    assert sum(Matrix(mat).eigenvals().values()) == 24


Code review comments:
diff_hunk: @@ -258,7 +258,7 @@ def _eigenvals_dict(
             f = charpoly.as_expr()
             x = charpoly.gen
             try:
-                eigs = {CRootOf(f, x, idx): 1 for idx in range(degree)}
+                eigs = Counter(CRootOf(f, x, idx) for idx in range(degree))
path: sympy/matrices/eigen.py
text: @user1:
Maybe it could be `charpoly.all_roots(multiple=False)`