#22523 closed enhancement (fixed)
Symbolic power of a matrix
Reported by:  mforets  Owned by:  

Priority:  major  Milestone:  sage8.0 
Component:  linear algebra  Keywords:  days84 
Cc:  tmonteil  Merged in:  
Authors:  Marcelo Forets  Reviewers:  Travis Scrimshaw 
Report Upstream:  N/A  Work issues:  
Branch:  3c82b43 (Commits)  Commit:  
Dependencies:  Stopgaps: 
Description (last modified by )
Compute the kth power of a matrix, where k is a symbolic variable.
See this question.
Followup at #25082.
Change History (23)
comment:1 Changed 3 years ago by
 Branch set to u/mforets/symbolic_power_of_a_matrix
comment:2 Changed 3 years ago by
 Commit set to ce9ac0b084b0a1b2f4c00b3a8d23232da2a1acdb
comment:3 followup: ↓ 4 Changed 3 years ago by
Either _matrix_power_symbolic
is a function and then you should call _matrix_power_symbolic(self, n)
, or _matrix_power_symbolic
is a method (of a generic matrix class), and then you should call self._matrix_power_symbolic(n)
. Here, you seem to mix both points by calling self._matrix_power_symbolic(self, n)
(so that self
is used for the two arguments of _matrix_power_symbolic
).
comment:4 in reply to: ↑ 3 Changed 3 years ago by
Replying to tmonteil:
Either
_matrix_power_symbolic
is a function and then you should call_matrix_power_symbolic(self, n)
, or_matrix_power_symbolic
is a method (of a generic matrix class), and then you should callself._matrix_power_symbolic(n)
. Here, you seem to mix both points by callingself._matrix_power_symbolic(self, n)
(so thatself
is used for the first two arguments of_matrix_power_symbolic
).
comment:5 Changed 3 years ago by
 Commit changed from ce9ac0b084b0a1b2f4c00b3a8d23232da2a1acdb to 215cd88ee7d32aeec709d8b1a81f7b53f897742e
Branch pushed to git repo; I updated commit sha1. New commits:
215cd88  Add symbolic matrix power.

comment:6 Changed 3 years ago by
 Status changed from new to needs_review
comment:7 Changed 3 years ago by
 Commit changed from 215cd88ee7d32aeec709d8b1a81f7b53f897742e to 3d69f78e9165026d9e98b8b2815e39ca26997b2a
Branch pushed to git repo; I updated commit sha1. New commits:
3d69f78  amend call to toplevel var

comment:8 Changed 3 years ago by
 Commit changed from 3d69f78e9165026d9e98b8b2815e39ca26997b2a to f6de852ce1f5eb5817ff2336ec1558b2425c9356
Branch pushed to git repo; I updated commit sha1. New commits:
f6de852  get rid of var and get rid of taking the derivative

comment:9 Changed 3 years ago by
also related, this other ask.sagemath question
comment:10 Changed 3 years ago by
the patchbot is complaining about:
sage t long src/doc/en/reference/references/index.rst Error: TAB character found at line 797 [0 tests, 0.00 s]  sage t long src/doc/en/reference/references/index.rst # Tab character found 
i'll try to fix this.
comment:11 Changed 3 years ago by
 Commit changed from f6de852ce1f5eb5817ff2336ec1558b2425c9356 to 883bb4a6fd679c6f0f193db108a6189b8cddf4ad
Branch pushed to git repo; I updated commit sha1. New commits:
883bb4a  remove tab character in a reference

comment:12 Changed 3 years ago by
Some quick comments:
 Use
Jk[i,i]
instead ofJk[i][i]
as the latter creates a temporaryvector
object. However, I would doJk_ii = Jk[i,i] if hasattr(Jk_ii, 'radical_expression'): Jk_ii = Jk_ii.radical_expression()
 Use
`
instead of$
for latex formatting in docstrings.  In the
INPUT:
block: ``A``  a square matrix over an exact field
Note the removal of the period too.  Make
A
andA^n
in latex.  Is there some reason why the function is not a
cdef
function?  I don't understand why you try to change the base ring to
QQbar
. Could you justify this?  This seems strange:
vector(SR, i).list()
. Why not[SR.zero()]*i
?  I would use a call to
binomial
instead of(factorial_n/(factorial(ni)*factorial(i))
because it probably requires less operations whenn
is big. Perhaps some testing is needed here.  How many of those imports are needed in the function? Having them there slows the function down.
 I'm not completely sure about creating the blank (dense) matrix and then filling it in afterwards. I feel like a better solution would be to just construct the data for the matrix and then construct the matrix from that (if this is even needed).
comment:13 Changed 3 years ago by
 Status changed from needs_review to needs_work
comment:14 Changed 3 years ago by
 Commit changed from 883bb4a6fd679c6f0f193db108a6189b8cddf4ad to 32e7cf9cee379d13642eeefdf81aa89ff73efc3b
comment:15 Changed 3 years ago by
@tscrim : thanks for the feedback, i've tried to correct the formatting issues.
some testing for the algorithm before changes (before 32e7cf9):
Test in the rationals, twodimensional:
sage: k = var('k') sage: A = matrix(QQ, 2, [1/2, 1, 2, 0]) sage: timeit('A^k') 5 loops, best of 3: 105 ms per loop
Test of the stochastic matrix in wikipedia:
sage: A = matrix([[0, 0, 1/2, 0, 1/2], [0, 0, 1, 0, 0], [1/4, 1/4, 0, 1/4, 1/4], [0, 0, 1/2, 0, 1/2], [0, 0, 0, 0, 1]]) sage: timeit('A^k') 5 loops, best of 3: 192 ms per loop
Test in the symbolic ring:
sage: n = var('n') sage: A = matrix([[pi, e],[0, 2*I]]) sage: timeit('A^n') 25 loops, best of 3: 26 ms per loop
here are some partial answers:
 Use
Jk[i,i]
instead ofJk[i][i]
as the latter creates a temporaryvector
object. However, I would doJk_ii = Jk[i,i] if hasattr(Jk_ii, 'radical_expression'): Jk_ii = Jk_ii.radical_expression()
yes, it is better that way. it didn't have an impact in the timing of my quick examples above, though (+ 10ms)
 I don't understand why you try to change the base ring to
QQbar
. Could you justify this?
this is to cover the cases where the eigenvalues do not belong to the same ring as the given matrix. in general the jordan_form
method breaks with the message RuntimeError: Some eigenvalue does not exist in Rational Field.
For example consider taking the jordan form of A = matrix(2, [1, 2, 1, 2])
.
On the other hand, the example A = matrix([[pi, e],[0, 2*I]])
tests the case when it is not possible to transform to QQbar (it would break with TypeError?: Illegal initializer for algebraic number
).
 This seems strange:
vector(SR, i).list()
. Why not[SR.zero()]*i
?
sure, i'm used to thinking in matrices/vectors in the 1st place, rather than lists :)
done in the new commit, for the record it didn't have an impact on timeit for the 3 examples above (+ 5..10ms)
 I would use a call to
binomial
instead of(factorial_n/(factorial(ni)*factorial(i))
because it probably requires less operations whenn
is big. Perhaps some testing is needed here.
right, at least with the binomial method the code is more clear.
more generally, i'm not confident that jordanization will be useful at all for 'big' n, since this operation is very costly. however if the matrix is very sparse then maybe it works, i didn't try any of that.
 How many of those imports are needed in the function? Having them there slows the function down.
well yes, i was able to get rid of importing vector, based on the suggestion of above.
comment:16 followup: ↓ 17 Changed 3 years ago by
Is there some reason why the function is not a cdef function?
i tried def
> cdef
in the function's definition, but this breaks many things (compilation). where can i learn more about using cdef in Sage library code? thanks
I'm not completely sure about creating the blank (dense) matrix and then filling it in afterwards. I feel like a better solution would be to just construct the data for the matrix and then construct the matrix from that (if this is even needed).
hmm i would normally use the same matrix to store partial computations (instead of creating a new one, say). are you suggesting that it's preferable to do differently like what, use lists?
in this particular case, there are FJ_k.set_row(i, Jk_row_i)
and FJ.set_block(k, k, FJ_k)
inside a loop. yes, the matrix FJ
is needed afterwards to do the computation P * FJ * Pinv
.
comment:17 in reply to: ↑ 16 Changed 3 years ago by
Replying to mforets:
Is there some reason why the function is not a cdef function?
i tried
def
>cdef
in the function's definition, but this breaks many things (compilation). where can i learn more about using cdef in Sage library code? thanks
Hmm...interesting. I will try this and see what I can do. To learn about cdef
, you should probably go to Cython's docs. The error messages on these failures can be helpful sometimes to fix things.
I'm not completely sure about creating the blank (dense) matrix and then filling it in afterwards. I feel like a better solution would be to just construct the data for the matrix and then construct the matrix from that (if this is even needed).
hmm i would normally use the same matrix to store partial computations (instead of creating a new one, say). are you suggesting that it's preferable to do differently like what, use lists?
Ideally you would use a single list, but a list of lists also works. At least for dense matrices. If the matrix was expected to be sparse, then a dict is good. The issue is more of that element creation is a relatively heavy operation compared to creating lists. Although I would directly modify the matrix FJ
.
in this particular case, there are
FJ_k.set_row(i, Jk_row_i)
andFJ.set_block(k, k, FJ_k)
inside a loop. yes, the matrixFJ
is needed afterwards to do the computationP * FJ * Pinv
.
However, FJ_k
is not. I would just like to modify FJ
since the intermediate matrix FJ_k
doesn't really need to be there.
comment:18 Changed 3 years ago by
 Commit changed from 32e7cf9cee379d13642eeefdf81aa89ff73efc3b to 3c82b4370366942e55885cdf88b504d4caaa6529
Branch pushed to git repo; I updated commit sha1. New commits:
3c82b43  directly modify FJ

comment:19 Changed 3 years ago by
 Milestone changed from sage7.6 to sage8.0
 Status changed from needs_work to needs_review
comment:20 Changed 3 years ago by
for the 2x2 case there is a simple closed formula by K. S. Williams. since it does not require the jordan form, it can also be used with an inexact ring, and maybe it's faster (i don't plan to try it, just a comment.)
comment:21 Changed 3 years ago by
 Reviewers set to Travis Scrimshaw
 Status changed from needs_review to positive_review
Let's get this into Sage. Positive review.
comment:22 Changed 3 years ago by
 Branch changed from u/mforets/symbolic_power_of_a_matrix to 3c82b4370366942e55885cdf88b504d4caaa6529
 Resolution set to fixed
 Status changed from positive_review to closed
comment:23 Changed 22 months ago by
 Commit 3c82b4370366942e55885cdf88b504d4caaa6529 deleted
 Description modified (diff)
Followup at #25028 after a bug was reported (and a fix found) at Ask Sage question 41622.
this patch is broken:
just adding the new
_matrix_power_symbolic
and updating the generic__pow__
doesn't seem to be correct. some hints?New commits:
added the code, but the method is not found