Apply connected components decomposition for computing eigenvalues

[sylee957]

References to other Issues or PRs

Brief description of what is fixed or changed

I think that matrix eigenvalues computations should use connected component decomposition by default because this brings some advantage over polynomials that it can factor out some polynomial terms as much as in matrix form as possible.

However, the converse may not hold true because the companion matrix is connected as a whole, even if the polynomial can be factored.

I've also removed the usage of EigenOnlyMatrix because I think that such designs make it difficult to organize matrix methods in this way.

I've also expanded the eigenvalue computation routines to _eigenvals_dict and _eigenvals_list.

Other comments

Release Notes

• matrices
  • Matrix([]).eigenvals(multiple=True) will give an empty list instead of an empty dict.

[sympy-bot]

✅

Hi, I am the SymPy bot (v158). I'm here to help you write a release notes entry. Please read the guide on how to write release notes.

Your release notes are in good order.

Here is what the release notes will look like:

• matrices
  • Matrix([]).eigenvals(multiple=True) will give an empty list instead of an empty dict. (#19355 by @sylee957)

This will be added to https://github.com/sympy/sympy/wiki/Release-Notes-for-1.7.

Note: This comment will be updated with the latest check if you edit the pull request. You need to reload the page to see it.

Click here to see the pull request description that was parsed.

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


#### Brief description of what is fixed or changed

I think that matrix eigenvalues computations should use connected component decomposition by default because this brings some advantage over polynomials that it can factor out some polynomial terms as much as in matrix form as possible.

However, the converse may not hold true because the companion matrix is connected as a whole, even if the polynomial can be factored.

I've also removed the usage of `EigenOnlyMatrix` because I think that such designs make it difficult to organize matrix methods in this way.

I've also expanded the eigenvalue computation routines to `_eigenvals_dict` and `_eigenvals_list`.

#### Other comments


#### Release Notes

<!-- Write the release notes for this release below. 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
  - `Matrix([]).eigenvals(multiple=True)` will give an empty list instead of an empty dict.
<!-- END RELEASE NOTES -->

Update

The release notes on the wiki have been updated.

[codecov]

Codecov Report

Merging #19355 into master will increase coverage by 0.000%.
The diff coverage is 92.857%.

@@            Coverage Diff            @@
##            master    #19355   +/-   ##
=========================================
  Coverage   75.618%   75.619%           
=========================================
  Files          651       651           
  Lines       169418    169511   +93     
  Branches     39973     39998   +25     
=========================================
+ Hits        128112    128183   +71     
- Misses       35703     35717   +14     
- Partials      5603      5611    +8     

[oscarbenjamin]

This looks good. I guess the same can be applied to many other functions like det, charpoly, inv, solve, rref etc. Maybe there are a few key places to apply the optimisation so it doesn't need to be explicitly listed everywhere. For example charpoly uses det so doing it for det works for both. Then again eigenvals uses charpoly so maybe it's better to apply this in det rather than in eigenvals. The det could be returned in factored form...

[sylee957]

I have to find some good benchmark examples for this

[oscarbenjamin]

How about the one in #16207?

Also this method applies to any permutation matrix if the permutation factors. Actually there are better ways to get the eigenvalues of a permutation matrix though:
https://math.stackexchange.com/a/1970819/538434

[sylee957]

Actually, I see it is only just twice faster for millisecons examples.

class TimePermutedBlockDiagonalEigenvals:
    """Benchmark examples obtained by similarly transforming random
    block diagonal matrices with permutation matrices to make it look
    like non block diagonal."""
    def setup(self):
        self.m22 = Matrix([[26, 0, 0, 7], [0, 27, 21, 0], [0, 18, 89, 0], [13, 0, 0, 28]])
        self.m222 = Matrix([[37, 0, 0, 0, 0, 5], [0, 32, 0, 0, 33, 0], [0, 0, 78, 91, 0, 0], [0, 0, 51, 97, 0, 0], [0, 97, 0, 0, 77, 0], [37, 0, 0, 0, 0, 61]])
        self.m2222 = Matrix([[87, 0, 12, 0, 0, 0, 0, 0], [0, 35, 0, 0, 0, 0, 0, 51], [31, 0, 47, 0, 0, 0, 0, 0], [0, 0, 0, 84, 0, 41, 0, 0], [0, 0, 0, 0, 70, 0, 57, 0], [0, 0, 0, 56, 0, 30, 0, 0], [0, 0, 0, 0, 54, 0, 55, 0], [0, 61, 0, 0, 0, 0, 0, 0]])
        self.m33 = Matrix([[48, 0, 44, 0, 0, 67], [0, 16, 0, 28, 61, 0], [5, 0, 5, 0, 0, 52], [0, 28, 0, 78, 13, 0], [0, 3, 0, 52, 35, 0], [98, 0, 86, 0, 0, 70]])
        self.m333 = Matrix([[60, 0, 74, 0, 0, 0, 0, 39, 0], [0, 36, 0, 0, 14, 0, 10, 0, 0], [33, 0, 32, 0, 0, 0, 0, 46, 0], [0, 0, 0, 9, 0, 46, 0, 0, 7], [0, 51, 0, 0, 92, 0, 46, 0, 0], [0, 0, 0, 86, 0, 21, 0, 0, 16], [0, 55, 0, 0, 28, 0, 12, 0, 0], [6, 0, 1, 0, 0, 0, 0, 31, 0], [0, 0, 0, 61, 0, 59, 0, 0, 57]])

    def time_eigenvals_22(self):
        self.m22.eigenvals()

    def time_eigenvals_222(self):
        self.m222.eigenvals()

    def time_eigenvals_2222(self):
        self.m2222.eigenvals()

    def time_eigenvals_33(self):
        self.m33.eigenvals()

    def time_eigenvals_333(self):
        self.m333.eigenvals()

       before           after         ratio
     [36075578]       [c82e87d7]
     <master>         <eigen_connected_components>
         70.8±1μs       73.1±0.8μs     1.03  matrices.TimeDiagonalEigenvals.time_eigenvals
-      5.59±0.1ms      2.61±0.02ms     0.47  matrices.TimePermutedBlockDiagonalEigenvals.time_eigenvals_22
-      9.42±0.1ms      4.35±0.01ms     0.46  matrices.TimePermutedBlockDiagonalEigenvals.time_eigenvals_222
-      14.0±0.3ms       5.83±0.2ms     0.42  matrices.TimePermutedBlockDiagonalEigenvals.time_eigenvals_2222
      4.97±0.07ms       3.59±0.1ms     0.72  matrices.TimePermutedBlockDiagonalEigenvals.time_eigenvals_33
-      13.2±0.7ms       6.61±0.2ms     0.50  matrices.TimePermutedBlockDiagonalEigenvals.time_eigenvals_333

[oscarbenjamin]

The advantage will be much bigger for larger matrices that break down into much smaller matrices:

from sympy import *
from time import time
from random import shuffle

x = Symbol('x')

for m in range(5):
    n = 5*2**m
    B = randMatrix(2)
    M = BlockMatrix([[B[i, j]*eye(n) for i in range(2)] for j in range(2)]).as_explicit()
    p = list(range(2*n))
    shuffle(p)
    Ms = M[p, p]
    start = time()
    Ms.eigenvals()
    print('n = %d,  T = %.3gsecs' % (n, time() - start))

With the PR this gives:

$ python g.py 
n = 5,  T = 0.098secs
n = 10,  T = 0.0305secs
n = 20,  T = 0.0706secs
n = 40,  T = 0.144secs
n = 80,  T = 0.492secs

On master I have

$ python g.py 
n = 5,  T = 0.0771secs
n = 10,  T = 0.189secs
n = 20,  T = 2.04secs
n = 40,  T = 29.5secs
n = 80,  T = 368secs

[sylee957]

from sympy import *
from sympy.combinatorics import Permutation

x = Symbol('x')

n = 5
B = randMatrix(2)
M = BlockDiagMatrix(*[randMatrix(2) for _ in range(n)]).as_explicit()
p = Permutation.random(n).array_form
Ms = M[p, p]
%timeit Ms.eigenvals()

n = 10
B = randMatrix(2)
M = BlockDiagMatrix(*[randMatrix(2) for _ in range(n)]).as_explicit()
p = Permutation.random(n).array_form
Ms = M[p, p]
%timeit Ms.eigenvals()

n = 15
B = randMatrix(2)
M = BlockDiagMatrix(*[randMatrix(2) for _ in range(n)]).as_explicit()
p = Permutation.random(n).array_form
Ms = M[p, p]
%timeit Ms.eigenvals()

n = 20
B = randMatrix(2)
M = BlockDiagMatrix(*[randMatrix(2) for _ in range(n)]).as_explicit()
p = Permutation.random(n).array_form
Ms = M[p, p]
%timeit Ms.eigenvals()

n = 25
B = randMatrix(2)
M = BlockDiagMatrix(*[randMatrix(2) for _ in range(n)]).as_explicit()
p = Permutation.random(n).array_form
Ms = M[p, p]
%timeit Ms.eigenvals()

8.36 ms ± 236 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
27.6 ms ± 443 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
63.7 ms ± 536 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
146 ms ± 7.86 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
488 ms ± 79 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

3.98 ms ± 84.9 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
8.65 ms ± 1.05 ms per loop (mean ± std. dev. of 7 runs, 100 loops each)
13 ms ± 1.12 ms per loop (mean ± std. dev. of 7 runs, 100 loops each)
17.3 ms ± 1.15 ms per loop (mean ± std. dev. of 7 runs, 100 loops each)
21.8 ms ± 817 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

Okay, I see the performance differs more greater for larger matrices with smaller block size.
But I can't add the benchmarks for now because it is too slow.
I've slightly modified your examples because I'm aware that there are some booting time for the first run of the benchmark, so things should be done statistically.

[sylee957]

Maybe there are a few key places to apply the optimisation so it doesn't need to be explicitly listed everywhere. For example charpoly uses det so doing it for det works for both.

As long as I have found, I don't think that charpoly actually uses det at all, but has its own routine.
And right now, I don't think that it's very consistent to modify the original charpoly to compute some already factored form, so improvements for determinants or inverses can be added independently than this.

[oscarbenjamin]

I think this is fine to merge if you're done.