Sage: Ticket #22523: Symbolic power of a matrix
https://trac.sagemath.org/ticket/22523
<p>
Compute the k-th power of a matrix, where k is a symbolic variable.
</p>
<p>
See <a class="ext-link" href="https://ask.sagemath.org/question/35658/general-power-of-a-matrix/#35666"><span class="icon"></span>this question</a>.
</p>
<p>
Follow-up at <a class="closed ticket" href="https://trac.sagemath.org/ticket/25082" title="defect: Fix symbolic power of matrix (closed: fixed)">#25082</a>.
</p>
en-usSagehttps://trac.sagemath.org/chrome/site/logo_sagemath_trac.png
https://trac.sagemath.org/ticket/22523
Trac 1.1.6mforetsFri, 10 Mar 2017 18:55:07 GMTbranch set
https://trac.sagemath.org/ticket/22523#comment:1
https://trac.sagemath.org/ticket/22523#comment:1
<ul>
<li><strong>branch</strong>
set to <em>u/mforets/symbolic_power_of_a_matrix</em>
</li>
</ul>
TicketmforetsFri, 10 Mar 2017 19:06:41 GMTcommit set
https://trac.sagemath.org/ticket/22523#comment:2
https://trac.sagemath.org/ticket/22523#comment:2
<ul>
<li><strong>commit</strong>
set to <em>ce9ac0b084b0a1b2f4c00b3a8d23232da2a1acdb</em>
</li>
</ul>
<p>
this patch is broken:
</p>
<pre class="wiki">sage: A = random_matrix(RDF, 2)
sage: A
[-0.9428173549065391 -0.301642050267529]
[ 0.6212822063471177 -0.9690284523509545]
sage: A^x
...
AttributeError: 'sage.matrix.matrix_real_double_dense.Matrix_real_double_dense' object has no attribute '_matrix_power_symbolic'
</pre><p>
just adding the new <code>_matrix_power_symbolic</code> and updating the generic <code> __pow__</code> doesn't seem to be correct. some hints?
</p>
<hr />
<p>
New commits:
</p>
<table class="wiki">
<tr><td><a class="ext-link" href="https://git.sagemath.org/sage.git/commit?id=ce9ac0b084b0a1b2f4c00b3a8d23232da2a1acdb"><span class="icon"></span>ce9ac0b</a></td><td><code>added the code, but the method is not found</code>
</td></tr></table>
TickettmonteilFri, 10 Mar 2017 19:12:25 GMT
https://trac.sagemath.org/ticket/22523#comment:3
https://trac.sagemath.org/ticket/22523#comment:3
<p>
Either <code>_matrix_power_symbolic</code> is a function and then you should call <code>_matrix_power_symbolic(self, n)</code>, or <code>_matrix_power_symbolic</code> is a method (of a generic matrix class), and then you should call <code>self._matrix_power_symbolic(n)</code>. Here, you seem to mix both points by calling <code>self._matrix_power_symbolic(self, n)</code> (so that <code>self</code> is used for the two arguments of <code>_matrix_power_symbolic</code>).
</p>
TickettmonteilFri, 10 Mar 2017 19:19:11 GMT
https://trac.sagemath.org/ticket/22523#comment:4
https://trac.sagemath.org/ticket/22523#comment:4
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/22523#comment:3" title="Comment 3">tmonteil</a>:
</p>
<blockquote class="citation">
<p>
Either <code>_matrix_power_symbolic</code> is a function and then you should call <code>_matrix_power_symbolic(self, n)</code>, or <code>_matrix_power_symbolic</code> is a method (of a generic matrix class), and then you should call <code>self._matrix_power_symbolic(n)</code>. Here, you seem to mix both points by calling <code>self._matrix_power_symbolic(self, n)</code> (so that <code>self</code> is used for the first two arguments of <code>_matrix_power_symbolic</code>).
</p>
</blockquote>
TicketgitSat, 11 Mar 2017 13:30:17 GMTcommit changed
https://trac.sagemath.org/ticket/22523#comment:5
https://trac.sagemath.org/ticket/22523#comment:5
<ul>
<li><strong>commit</strong>
changed from <em>ce9ac0b084b0a1b2f4c00b3a8d23232da2a1acdb</em> to <em>215cd88ee7d32aeec709d8b1a81f7b53f897742e</em>
</li>
</ul>
<p>
Branch pushed to git repo; I updated commit sha1. New commits:
</p>
<table class="wiki">
<tr><td><a class="ext-link" href="https://git.sagemath.org/sage.git/commit/?id=215cd88ee7d32aeec709d8b1a81f7b53f897742e"><span class="icon"></span>215cd88</a></td><td><code>Add symbolic matrix power.</code>
</td></tr></table>
TicketmforetsSat, 11 Mar 2017 13:30:50 GMTstatus changed
https://trac.sagemath.org/ticket/22523#comment:6
https://trac.sagemath.org/ticket/22523#comment:6
<ul>
<li><strong>status</strong>
changed from <em>new</em> to <em>needs_review</em>
</li>
</ul>
TicketgitSat, 11 Mar 2017 19:59:01 GMTcommit changed
https://trac.sagemath.org/ticket/22523#comment:7
https://trac.sagemath.org/ticket/22523#comment:7
<ul>
<li><strong>commit</strong>
changed from <em>215cd88ee7d32aeec709d8b1a81f7b53f897742e</em> to <em>3d69f78e9165026d9e98b8b2815e39ca26997b2a</em>
</li>
</ul>
<p>
Branch pushed to git repo; I updated commit sha1. New commits:
</p>
<table class="wiki">
<tr><td><a class="ext-link" href="https://git.sagemath.org/sage.git/commit/?id=3d69f78e9165026d9e98b8b2815e39ca26997b2a"><span class="icon"></span>3d69f78</a></td><td><code>amend call to top-level var</code>
</td></tr></table>
TicketgitSun, 12 Mar 2017 15:29:38 GMTcommit changed
https://trac.sagemath.org/ticket/22523#comment:8
https://trac.sagemath.org/ticket/22523#comment:8
<ul>
<li><strong>commit</strong>
changed from <em>3d69f78e9165026d9e98b8b2815e39ca26997b2a</em> to <em>f6de852ce1f5eb5817ff2336ec1558b2425c9356</em>
</li>
</ul>
<p>
Branch pushed to git repo; I updated commit sha1. New commits:
</p>
<table class="wiki">
<tr><td><a class="ext-link" href="https://git.sagemath.org/sage.git/commit/?id=f6de852ce1f5eb5817ff2336ec1558b2425c9356"><span class="icon"></span>f6de852</a></td><td><code>get rid of var and get rid of taking the derivative</code>
</td></tr></table>
TicketmforetsMon, 10 Apr 2017 07:47:35 GMT
https://trac.sagemath.org/ticket/22523#comment:9
https://trac.sagemath.org/ticket/22523#comment:9
<p>
also related, this other <a class="ext-link" href="https://ask.sagemath.org/question/8280/an-nvarn/"><span class="icon"></span>ask.sagemath question</a>
</p>
TicketmforetsMon, 10 Apr 2017 07:51:03 GMT
https://trac.sagemath.org/ticket/22523#comment:10
https://trac.sagemath.org/ticket/22523#comment:10
<p>
the patchbot is complaining about:
</p>
<pre class="wiki">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
----------------------------------------------------------------------
</pre><p>
i'll try to fix this.
</p>
TicketgitTue, 11 Apr 2017 08:48:37 GMTcommit changed
https://trac.sagemath.org/ticket/22523#comment:11
https://trac.sagemath.org/ticket/22523#comment:11
<ul>
<li><strong>commit</strong>
changed from <em>f6de852ce1f5eb5817ff2336ec1558b2425c9356</em> to <em>883bb4a6fd679c6f0f193db108a6189b8cddf4ad</em>
</li>
</ul>
<p>
Branch pushed to git repo; I updated commit sha1. New commits:
</p>
<table class="wiki">
<tr><td><a class="ext-link" href="https://git.sagemath.org/sage.git/commit/?id=883bb4a6fd679c6f0f193db108a6189b8cddf4ad"><span class="icon"></span>883bb4a</a></td><td><code>remove tab character in a reference</code>
</td></tr></table>
TickettscrimTue, 11 Apr 2017 16:02:48 GMT
https://trac.sagemath.org/ticket/22523#comment:12
https://trac.sagemath.org/ticket/22523#comment:12
<p>
Some quick comments:
</p>
<ul><li>Use <code>Jk[i,i]</code> instead of <code>Jk[i][i]</code> as the latter creates a temporary <code>vector</code> object. However, I would do
<div class="wiki-code"><div class="code"><pre> Jk_ii <span class="o">=</span> Jk<span class="p">[</span>i<span class="p">,</span>i<span class="p">]</span>
<span class="k">if</span> <span class="nb">hasattr</span><span class="p">(</span>Jk_ii<span class="p">,</span> <span class="s">'radical_expression'</span><span class="p">):</span>
Jk_ii <span class="o">=</span> Jk_ii<span class="o">.</span>radical_expression<span class="p">()</span>
</pre></div></div></li><li>Use <code>`</code> instead of <code>$</code> for latex formatting in docstrings.
</li><li>In the <code>INPUT:</code> block: <code>- ``A`` -- a square matrix over an exact field</code> Note the removal of the period too.
</li><li>Make <code>A</code> and <code>A^n</code> in latex.
</li><li>Is there some reason why the function is not a <code>cdef</code> function?
</li><li>I don't understand why you try to change the base ring to <code>QQbar</code>. Could you justify this?
</li><li>This seems strange: <code>vector(SR, i).list()</code>. Why not <code>[SR.zero()]*i</code>?
</li><li>I would use a call to <code>binomial</code> instead of <code>(factorial_n/(factorial(n-i)*factorial(i))</code> because it probably requires less operations when <code>n</code> is big. Perhaps some testing is needed here.
</li><li>How many of those imports are needed in the function? Having them there slows the function down.
</li><li>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).
</li></ul>
TicketmforetsTue, 11 Apr 2017 18:15:07 GMTstatus changed
https://trac.sagemath.org/ticket/22523#comment:13
https://trac.sagemath.org/ticket/22523#comment:13
<ul>
<li><strong>status</strong>
changed from <em>needs_review</em> to <em>needs_work</em>
</li>
</ul>
TicketgitWed, 12 Apr 2017 19:15:48 GMTcommit changed
https://trac.sagemath.org/ticket/22523#comment:14
https://trac.sagemath.org/ticket/22523#comment:14
<ul>
<li><strong>commit</strong>
changed from <em>883bb4a6fd679c6f0f193db108a6189b8cddf4ad</em> to <em>32e7cf9cee379d13642eeefdf81aa89ff73efc3b</em>
</li>
</ul>
<p>
Branch pushed to git repo; I updated commit sha1. New commits:
</p>
<table class="wiki">
<tr><td><a class="ext-link" href="https://git.sagemath.org/sage.git/commit/?id=c473b60796d78a49cc572901dd7192d9ecea5cf5"><span class="icon"></span>c473b60</a></td><td><code>fix docstring formatting</code>
</td></tr><tr><td><a class="ext-link" href="https://git.sagemath.org/sage.git/commit/?id=32e7cf9cee379d13642eeefdf81aa89ff73efc3b"><span class="icon"></span>32e7cf9</a></td><td><code>updates in _matrix_power_symbolic</code>
</td></tr></table>
TicketmforetsWed, 12 Apr 2017 19:24:27 GMT
https://trac.sagemath.org/ticket/22523#comment:15
https://trac.sagemath.org/ticket/22523#comment:15
<p>
@tscrim : thanks for the feedback, i've tried to correct the formatting issues.
</p>
<p>
some testing for the algorithm before changes (before 32e7cf9):
</p>
<p>
Test in the rationals, two-dimensional:
</p>
<pre class="wiki"> 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
</pre><p>
Test of the <a class="ext-link" href="https://en.wikipedia.org/wiki/Stochastic_matrix"><span class="icon"></span>stochastic matrix</a> in wikipedia:
</p>
<pre class="wiki"> 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
</pre><p>
Test in the symbolic ring:
</p>
<pre class="wiki"> 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
</pre><p>
here are some partial answers:
</p>
<blockquote class="citation">
<ul><li>Use <code>Jk[i,i]</code> instead of <code>Jk[i][i]</code> as the latter creates a temporary <code>vector</code> object. However, I would do
<div class="wiki-code"><div class="code"><pre> Jk_ii <span class="o">=</span> Jk<span class="p">[</span>i<span class="p">,</span>i<span class="p">]</span>
<span class="k">if</span> <span class="nb">hasattr</span><span class="p">(</span>Jk_ii<span class="p">,</span> <span class="s">'radical_expression'</span><span class="p">):</span>
Jk_ii <span class="o">=</span> Jk_ii<span class="o">.</span>radical_expression<span class="p">()</span>
</pre></div></div></li></ul></blockquote>
<p>
yes, it is better that way. it didn't have an impact in the timing of my quick examples above, though (+- 10ms)
</p>
<blockquote class="citation">
<ul><li>I don't understand why you try to change the base ring to <code>QQbar</code>. Could you justify this?
</li></ul></blockquote>
<p>
this is to cover the cases where the eigenvalues do not belong to the same ring as the given matrix. in general the <code>jordan_form</code> method breaks with the message <code>RuntimeError: Some eigenvalue does not exist in Rational Field.</code>
</p>
<p>
For example consider taking the jordan form of <code>A = matrix(2, [1, 2, 1, -2])</code>.
</p>
<p>
On the other hand, the example <code>A = matrix([[pi, e],[0, -2*I]])</code> tests the case when it is not possible to transform to QQbar (it would break with <code></code><a class="missing wiki">TypeError?</a>: Illegal initializer for algebraic number<code></code>).
</p>
<blockquote class="citation">
<ul><li>This seems strange: <code>vector(SR, i).list()</code>. Why not <code>[SR.zero()]*i</code>?
</li></ul></blockquote>
<p>
sure, i'm used to thinking in matrices/vectors in the 1st place, rather than lists :)
</p>
<p>
done in the new commit, for the record it didn't have an impact on timeit for the 3 examples above (+- 5..10ms)
</p>
<blockquote class="citation">
<ul><li>I would use a call to <code>binomial</code> instead of <code>(factorial_n/(factorial(n-i)*factorial(i))</code> because it probably requires less operations when <code>n</code> is big. Perhaps some testing is needed here.
</li></ul></blockquote>
<p>
right, at least with the binomial method the code is more clear.
</p>
<p>
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.
</p>
<blockquote class="citation">
<ul><li>How many of those imports are needed in the function? Having them there slows the function down.
</li></ul></blockquote>
<p>
well yes, i was able to get rid of importing vector, based on the suggestion of above.
</p>
TicketmforetsThu, 13 Apr 2017 13:11:55 GMT
https://trac.sagemath.org/ticket/22523#comment:16
https://trac.sagemath.org/ticket/22523#comment:16
<blockquote class="citation">
<p>
Is there some reason why the function is not a cdef function?
</p>
</blockquote>
<p>
i tried <code>def</code> -> <code>cdef</code> in the function's definition, but this breaks many things (compilation). where can i learn more about using cdef in Sage library code? thanks
</p>
<blockquote class="citation">
<p>
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).
</p>
</blockquote>
<p>
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?
</p>
<p>
in this particular case, there are <code>FJ_k.set_row(i, Jk_row_i)</code> and <code>FJ.set_block(k, k, FJ_k)</code> inside a loop. yes, the matrix <code>FJ</code> is needed afterwards to do the computation <code>P * FJ * Pinv</code>.
</p>
TickettscrimSat, 15 Apr 2017 04:12:51 GMT
https://trac.sagemath.org/ticket/22523#comment:17
https://trac.sagemath.org/ticket/22523#comment:17
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/22523#comment:16" title="Comment 16">mforets</a>:
</p>
<blockquote class="citation">
<blockquote class="citation">
<p>
Is there some reason why the function is not a cdef function?
</p>
</blockquote>
<p>
i tried <code>def</code> -> <code>cdef</code> in the function's definition, but this breaks many things (compilation). where can i learn more about using cdef in Sage library code? thanks
</p>
</blockquote>
<p>
Hmm...interesting. I will try this and see what I can do. To learn about <code>cdef</code>, you should probably go to Cython's docs. The error messages on these failures can be helpful sometimes to fix things.
</p>
<blockquote class="citation">
<blockquote class="citation">
<p>
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).
</p>
</blockquote>
<p>
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?
</p>
</blockquote>
<p>
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 <code>FJ</code>.
</p>
<blockquote class="citation">
<p>
in this particular case, there are <code>FJ_k.set_row(i, Jk_row_i)</code> and <code>FJ.set_block(k, k, FJ_k)</code> inside a loop. yes, the matrix <code>FJ</code> is needed afterwards to do the computation <code>P * FJ * Pinv</code>.
</p>
</blockquote>
<p>
However, <code>FJ_k</code> is not. I would just like to modify <code>FJ</code> since the intermediate matrix <code>FJ_k</code> doesn't really need to be there.
</p>
TicketgitMon, 17 Apr 2017 22:31:17 GMTcommit changed
https://trac.sagemath.org/ticket/22523#comment:18
https://trac.sagemath.org/ticket/22523#comment:18
<ul>
<li><strong>commit</strong>
changed from <em>32e7cf9cee379d13642eeefdf81aa89ff73efc3b</em> to <em>3c82b4370366942e55885cdf88b504d4caaa6529</em>
</li>
</ul>
<p>
Branch pushed to git repo; I updated commit sha1. New commits:
</p>
<table class="wiki">
<tr><td><a class="ext-link" href="https://git.sagemath.org/sage.git/commit/?id=3c82b4370366942e55885cdf88b504d4caaa6529"><span class="icon"></span>3c82b43</a></td><td><code>directly modify FJ</code>
</td></tr></table>
TicketmforetsMon, 17 Apr 2017 22:38:24 GMTstatus, milestone changed
https://trac.sagemath.org/ticket/22523#comment:19
https://trac.sagemath.org/ticket/22523#comment:19
<ul>
<li><strong>status</strong>
changed from <em>needs_work</em> to <em>needs_review</em>
</li>
<li><strong>milestone</strong>
changed from <em>sage-7.6</em> to <em>sage-8.0</em>
</li>
</ul>
TicketmforetsMon, 17 Apr 2017 22:48:24 GMT
https://trac.sagemath.org/ticket/22523#comment:20
https://trac.sagemath.org/ticket/22523#comment:20
<p>
for the 2x2 case there is a simple closed formula by <a class="ext-link" href="http://people.math.carleton.ca/~williams/papers/pdf/175.pdf"><span class="icon"></span>K. S. Williams</a>. 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.)
</p>
TickettscrimMon, 17 Apr 2017 23:18:13 GMTstatus changed; reviewer set
https://trac.sagemath.org/ticket/22523#comment:21
https://trac.sagemath.org/ticket/22523#comment:21
<ul>
<li><strong>status</strong>
changed from <em>needs_review</em> to <em>positive_review</em>
</li>
<li><strong>reviewer</strong>
set to <em>Travis Scrimshaw</em>
</li>
</ul>
<p>
Let's get this into Sage. Positive review.
</p>
TicketvbraunSun, 23 Apr 2017 12:58:09 GMTstatus, branch changed; resolution set
https://trac.sagemath.org/ticket/22523#comment:22
https://trac.sagemath.org/ticket/22523#comment:22
<ul>
<li><strong>status</strong>
changed from <em>positive_review</em> to <em>closed</em>
</li>
<li><strong>resolution</strong>
set to <em>fixed</em>
</li>
<li><strong>branch</strong>
changed from <em>u/mforets/symbolic_power_of_a_matrix</em> to <em>3c82b4370366942e55885cdf88b504d4caaa6529</em>
</li>
</ul>
TicketslelievreMon, 02 Apr 2018 18:08:44 GMTdescription changed; commit deleted
https://trac.sagemath.org/ticket/22523#comment:23
https://trac.sagemath.org/ticket/22523#comment:23
<ul>
<li><strong>commit</strong>
<em>3c82b4370366942e55885cdf88b504d4caaa6529</em> deleted
</li>
<li><strong>description</strong>
modified (<a href="/ticket/22523?action=diff&version=23">diff</a>)
</li>
</ul>
<p>
Follow-up at <a class="closed ticket" href="https://trac.sagemath.org/ticket/25028" title="defect: pari('f(x)=f(x)')(0) crashes Sage (closed: worksforme)">#25028</a> after a bug was reported (and a fix found) at
<a class="ext-link" href="https://ask.sagemath.org/question/41622/bug-in-general-power-of-a-matrix"><span class="icon"></span>Ask Sage question 41622</a>.
</p>
Ticket