Faster matrix exponentiation #18595

abhinav-anand-addepar · 2020-02-07T05:43:51Z

matrices
- Faster Matrix exponentiation using Cayley Hamilton Theorem

sympy-bot · 2020-02-07T05:43:59Z

✅

Hi, I am the SymPy bot (v149). 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
- Faster Matrix exponentiation using Cayley Hamilton Theorem (#18595 by @abhinav28071999)

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

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.

<!-- BEGIN RELEASE NOTES -->
* matrices
  * Faster Matrix exponentiation using Cayley Hamilton Theorem
<!-- END RELEASE NOTES -->

Update

The release notes on the wiki have been updated.

abhinav-anand-addepar · 2020-02-07T06:05:01Z

When A = Matrix([[2, 3, 4, 4, 5], [2, 1, 2, 3, 1], [4, 3, 5, 0, 2], [3, 1, 4, 0, 3], [3, 2, 3, 4, 3]])
exp = 200000
Using Cayley 2.937811851501465
Without Using Cayley 16.406556844711304
exp = 500000
Using Cayley 10.280319690704346
Without Using Cayley 69.17847156524658
exp = 1000000
Using Cayley 26.187591552734375
Without Using Cayley 208.7181031703949
[Finished in 334.3s]

abhinav-anand-addepar · 2020-02-07T06:08:30Z

I think the condition check for n can be removed. I will have to do the timing analysis for different rows

abhinav-anand-addepar · 2020-02-07T06:22:52Z

When A = Matrix([[2, 1, 3, 0, 0, 2, 3], [2, 3, 1, 3, 3, 2, 3], [3, 3, 2, 3, 1, 2, 2], [2, 0, 0, 0, 3, 0, 1], [2, 2, 3, 0, 0, 3, 3], [0, 2, 2, 3, 1, 3, 1], [3, 3, 3, 3, 0, 0, 1]])
exp = 200000
Using Cayley 6.235580921173096
Without Using Cayley 42.61778020858765
exp = 500000
Using Cayley 20.59529948234558
Without Using Cayley 174.21074295043945
exp = 1000000
Using Cayley 49.71791124343872
Without Using Cayley 519.055638551712
[Finished in 813.1s]

sympy/matrices/common.py

sylee957 · 2020-02-07T06:50:15Z

Generally, I don't want to see another polynomial long division reimplemented here
Why not use https://docs.sympy.org/latest/modules/polys/reference.html#sympy.polys.polytools.rem

Is it slow?

abhinav-anand-addepar · 2020-02-07T06:57:37Z

Generally, I don't want to see another polynomial long division reimplemented here
Why not use https://docs.sympy.org/latest/modules/polys/reference.html#sympy.polys.polytools.rem

Is it slow?

I am not using polynomial long division

sylee957 · 2020-02-07T07:09:09Z

But I mean what about contributing to the polynomial for more general use

I've already found some bug with a matrix with complex entries
Matrix([[2+I, 3, 4, 4, 5], [2, 1, 2, 3, 1], [4, 3, 5, 0, 2], [3, 1, 4, 0, 3], [3, 2, 3, 4, 3]])
And I doubt that you can use cayley exponential with a matrix of symbolic coefficients unless it's a special case

sylee957 · 2020-02-07T07:26:52Z

And I ask you what kind of method did you use for computing the polynomial remainer

abhinav-anand-addepar · 2020-02-07T07:47:01Z

And I ask you what kind of method did you use for computing the polynomial remainer

There is no polynomial remainder here. I am using this algorithm.

abhinav-anand-addepar · 2020-02-07T07:56:05Z

But I mean what about contributing to the polynomial for more general use

I've already found some bug with a matrix with complex entries
Matrix([[2+I, 3, 4, 4, 5], [2, 1, 2, 3, 1], [4, 3, 5, 0, 2], [3, 1, 4, 0, 3], [3, 2, 3, 4, 3]])
And I doubt that you can use cayley exponential with a matrix of symbolic coefficients unless it's a special case

It's working fine.
The code i used is:

Can you tell what exponent you used.

sylee957 · 2020-02-07T07:59:23Z

from sympy import *

A = Matrix([[2+I, 3, 4, 4, 5], [2, 1, 2, 3, 1], [4, 3, 5, 0, 2], [3, 1, 4, 0, 3], [3, 2, 3, 4, 3]])
A._pow_cayley(100)

I'm calling from the private function

In my second thought, I suppose it can be used for symbolic or complex coefficients, but it should be fixed.

sylee957 · 2020-02-07T11:12:12Z

And I still think that what you're doing here is computing the remainder of x**n with the characteristic polynomial of the matrix.
So if that is the case, I'd expect that this can be made to use less symbolic substitutions, but more direct coefficient manipulation, if this is the reason it's causing the problems.

abhinav-anand-addepar · 2020-02-07T11:39:19Z

And I still think that what you're doing here is computing the remainder of x**n with the characteristic polynomial of the matrix.
So if that is the case, I'd expect that this can be made to use less symbolic substitutions, but more direct coefficient manipulation, if this is the reason it's causing the problems.

the problem was that i initialised the polynomial without telling the variable. I used poly(x**4). I should have used poly(x**4,x). It will cause problem when there is another symbol in the charpoly or A. Now it should be working fine.
But there is still a problem , that is that i used {x, c1 , c2......} symbols in the function so the matrix should not contain any of these symbols. Is there any symbol in sympy that can be initialised only within a function and not by the user?

abhinav-anand-addepar · 2020-02-07T11:51:33Z

And I still think that what you're doing here is computing the remainder of x**n with the characteristic polynomial of the matrix.
So if that is the case, I'd expect that this can be made to use less symbolic substitutions, but more direct coefficient manipulation, if this is the reason it's causing the problems.

Yes i think you are right, let me check which is faster.

abhinav-anand-addepar · 2020-02-07T12:00:26Z

@sylee957 i checked. polynomial division is too slow to be used here. where exp > 1000, i found that there was significant drop in performance when we use polynomial division.

sylee957 · 2020-02-07T12:10:14Z

In my opinion, I think that the idea from here can be generalized to improve the polynomial remainder with a monic polynomial.

abhinav-anand-addepar · 2020-02-07T12:12:23Z

In my opinion, I think that the idea from here can be generalized to improve the polynomial remainder with a monic polynomial.

this is a special case. here x**n can be divided by a multivariate polynomial.

sylee957 · 2020-02-07T12:18:51Z

But there is still a problem , that is that i used {x, c1 , c2......} symbols in the function so the matrix should not contain any of these symbols. Is there any symbol in sympy that can be initialised only within a function and not by the user?

You should use dummy if you are going to use symbols in library code. But generally, I don't think that it's a right direction unless the problem is too difficult.

sympy/matrices/common.py

sympy/matrices/tests/test_commonmatrix.py

sympy/matrices/common.py

oscarbenjamin · 2020-02-07T16:57:18Z

sympy/matrices/common.py

@@ -2377,7 +2454,9 @@ def __pow__(self, exp):

        return self.pow(exp)

-    def pow(self, exp, jordan=None):
+
+    def pow(self, exp, jordan=None, cayley=False):


Having method='jordan' is better than two boolean arguments

The pow method here is new since the last release of sympy so the signature can be changed.

Yes, I’d also want to use a single keyword to control every methods

jksuom · 2020-02-07T18:24:31Z

This looks like solving the linear recursion defined by the characteristic function. The function linrec was designed for numerical sequences but it seems to work with more general objects that have the relevant arithmetic methods defined. For matrices, it is necessary to define addition with zero. I have tested with this quick fix:

--- a/sympy/matrices/common.py
+++ b/sympy/matrices/common.py
@@ -2234,6 +2234,8 @@ def __abs__(self):
     @call_highest_priority('__radd__')
     def __add__(self, other):
         """Return self + other, raising ShapeError if shapes don't match."""
+        if not other:
+            return self
         other = _matrixify(other)
         # matrix-like objects can have shapes.  This is
         # our first sanity check.

Then

>>> from sympy import Matrix, eye
>>> from sympy.discrete.recurrences import linrec
>>> A = Matrix([[1,2],[4,5]])
>>> p = A.charpoly()
>>> coeffs = (-p).all_coeffs()[1:]
>>> init = [eye(2), A]
>>> linrec(coeffs, init, 20)
Matrix([
[3428578511656497,  4683525345792216],
[9367050691584432, 12795629203240929]])
>>> _ == A**20
True

jksuom · 2020-02-08T08:51:03Z

sympy/discrete/recurrences.py

@@ -6,7 +6,7 @@
 from sympy.core import S, sympify
 from sympy.core.compatibility import as_int, iterable

-def linrec(coeffs, init, n):
+def linrec(coeffs, init, n, final_expr=False):


How about having a separate linrec_coeffs(coeffs, n) that could also be used by linrec?

Yes, I would also agree on @jksuom 's idea to have a separate function than adding up a keyword.
It is not a good design because setting the keyword True would make the variable init obscure.

can you look at the error?

codecov · 2020-02-09T08:12:23Z

Codecov Report

Merging #18595 into master will decrease coverage by 0.005%.
The diff coverage is 51.063%.

@@             Coverage Diff              @@
##            master   #18595       +/-   ##
============================================
- Coverage   75.585%   75.58%   -0.006%     
============================================
  Files          644      644               
  Lines       167452   167573      +121     
  Branches     39462    39502       +40     
============================================
+ Hits        126570   126652       +82     
- Misses       35364    35391       +27     
- Partials      5518     5530       +12

sympy/discrete/recurrences.py

jksuom · 2020-02-09T11:29:00Z

sympy/discrete/recurrences.py

+
+    return b[n] if n < k else sum(u*v for u, v in zip(_linrec_coeffs(c, n), b))
+
+def _linrec_coeffs(c, n):


Perhaps this could be explicitly public (no _) though not necessarily imported in __init__.py. The docstring could refer to linrec.

sympy/discrete/recurrences.py

sympy/matrices/common.py

oscarbenjamin · 2020-02-09T14:32:59Z

sympy/matrices/common.py

-            be used if possible and False means it should not be used unless
-            it is the only way to calculate the power.
+
+        method : jordon, cayley


It should also be possible to specify recursion

Perhaps recursion should be specified as "multiply" though since I think that makes more sense to users.

Actually since dotprodsimp is only used for recursion it should probably be:

method : multiply, mulsimp, jordon, cayley

Then mulsimp corresponds to using recursion with dotprodsimp and the dotprodsimp parameter can be removed (it was added since last release).

I have made the changes. Is is good now?

sympy/matrices/common.py

sympy/discrete/recurrences.py

sympy/matrices/common.py

sylee957 · 2020-02-10T03:18:46Z

sympy/matrices/common.py

-    def pow(self, exp, jordan=None, dotprodsimp=None):
+
+    def pow(self, exp, method=None):
+


Suggested change

sylee957 · 2020-02-10T03:19:05Z

sympy/matrices/common.py

-            If left as None then Jordan form exponentiation will be used under
-            certain conditions, True specifies that jordan_pow should always
-            be used if possible and False means it should not be used unless
-            it is the only way to calculate the power.



Suggested change

abhinav-anand-addepar · 2020-02-10T03:59:42Z

I don't know why the tests are failing, could anyone take a look

abhinav-anand-addepar · 2020-02-10T05:17:43Z

got it

oscarbenjamin · 2020-02-10T13:59:28Z

sympy/matrices/common.py

                        raise

-            if _get_intermediate_simp_bool(True, dotprodsimp):
+            if method == 'cayley':
+                return a._eval_pow_by_cayley(exp)


These things should outside the enclosing if

Done. Is this okay?

jksuom · 2020-02-13T07:17:46Z

sympy/matrices/common.py

-                return a._eval_pow_by_recursion_dotprodsimp(exp)
+        elif method == 'cayley':
+            if not exp.is_Number or exp % 1 != 0:
+                raise ValueError("cayley method is only valid for integer powers")


nonnegative integer powers?

when the power is negative, then the matrix is inverted and exp *= -1 is done in previous step.
So, even method = multiply works for negative powers.

abhinav-anand-addepar requested review from oscarbenjamin and sylee957 February 7, 2020 05:56

abhinav-anand-addepar requested review from asmeurer, jksuom and smichr February 7, 2020 06:26

Mohitbalwani26 reviewed Feb 7, 2020

View reviewed changes

sympy/matrices/common.py Outdated Show resolved Hide resolved

sylee957 reviewed Feb 7, 2020

View reviewed changes

sympy/matrices/common.py Outdated Show resolved Hide resolved

smichr reviewed Feb 7, 2020

View reviewed changes

sympy/matrices/tests/test_commonmatrix.py Outdated Show resolved Hide resolved

smichr reviewed Feb 7, 2020

View reviewed changes

sympy/matrices/common.py Outdated Show resolved Hide resolved

smichr reviewed Feb 7, 2020

View reviewed changes

sympy/matrices/common.py Outdated Show resolved Hide resolved

oscarbenjamin reviewed Feb 7, 2020

View reviewed changes

abhinav-anand-addepar closed this Feb 8, 2020

abhinav-anand-addepar reopened this Feb 8, 2020

jksuom reviewed Feb 8, 2020

View reviewed changes

abhinav-anand-addepar requested review from sylee957 and jksuom February 9, 2020 09:56

jksuom reviewed Feb 9, 2020

View reviewed changes

sympy/discrete/recurrences.py Outdated Show resolved Hide resolved

jksuom reviewed Feb 9, 2020

View reviewed changes

sympy/discrete/recurrences.py Outdated Show resolved Hide resolved

oscarbenjamin reviewed Feb 9, 2020

View reviewed changes

sympy/matrices/common.py Outdated Show resolved Hide resolved

oscarbenjamin reviewed Feb 9, 2020

View reviewed changes

sympy/matrices/common.py Show resolved Hide resolved