Opened 2 years ago

Closed 2 years ago

Last modified 13 months ago

#22523 closed enhancement (fixed)

Symbolic power of a matrix

Reported by: mforets Owned by:
Priority: major Milestone: sage-8.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 slelievre)

Compute the k-th power of a matrix, where k is a symbolic variable.

See this question.

Follow-up at #25082.

Change History (23)

comment:1 Changed 2 years ago by mforets

  • Branch set to u/mforets/symbolic_power_of_a_matrix

comment:2 Changed 2 years ago by mforets

  • Commit set to ce9ac0b084b0a1b2f4c00b3a8d23232da2a1acdb

this patch is broken:

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'

just adding the new _matrix_power_symbolic and updating the generic __pow__ doesn't seem to be correct. some hints?


New commits:

ce9ac0badded the code, but the method is not found

comment:3 follow-up: Changed 2 years ago by 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 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 2 years ago by tmonteil

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 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 first two arguments of _matrix_power_symbolic).

comment:5 Changed 2 years ago by git

  • Commit changed from ce9ac0b084b0a1b2f4c00b3a8d23232da2a1acdb to 215cd88ee7d32aeec709d8b1a81f7b53f897742e

Branch pushed to git repo; I updated commit sha1. New commits:

215cd88Add symbolic matrix power.

comment:6 Changed 2 years ago by mforets

  • Status changed from new to needs_review

comment:7 Changed 2 years ago by git

  • Commit changed from 215cd88ee7d32aeec709d8b1a81f7b53f897742e to 3d69f78e9165026d9e98b8b2815e39ca26997b2a

Branch pushed to git repo; I updated commit sha1. New commits:

3d69f78amend call to top-level var

comment:8 Changed 2 years ago by git

  • Commit changed from 3d69f78e9165026d9e98b8b2815e39ca26997b2a to f6de852ce1f5eb5817ff2336ec1558b2425c9356

Branch pushed to git repo; I updated commit sha1. New commits:

f6de852get rid of var and get rid of taking the derivative

comment:9 Changed 2 years ago by mforets

also related, this other ask.sagemath question

comment:10 Changed 2 years ago by mforets

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 2 years ago by git

  • Commit changed from f6de852ce1f5eb5817ff2336ec1558b2425c9356 to 883bb4a6fd679c6f0f193db108a6189b8cddf4ad

Branch pushed to git repo; I updated commit sha1. New commits:

883bb4aremove tab character in a reference

comment:12 Changed 2 years ago by tscrim

Some quick comments:

  • Use Jk[i,i] instead of Jk[i][i] as the latter creates a temporary vector object. However, I would do
              Jk_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 and A^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(n-i)*factorial(i)) because it probably requires less operations when n 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 2 years ago by mforets

  • Status changed from needs_review to needs_work

comment:14 Changed 2 years ago by git

  • Commit changed from 883bb4a6fd679c6f0f193db108a6189b8cddf4ad to 32e7cf9cee379d13642eeefdf81aa89ff73efc3b

Branch pushed to git repo; I updated commit sha1. New commits:

c473b60fix docstring formatting
32e7cf9updates in _matrix_power_symbolic

comment:15 Changed 2 years ago by mforets

@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, two-dimensional:

    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 of Jk[i][i] as the latter creates a temporary vector object. However, I would do
              Jk_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(n-i)*factorial(i)) because it probably requires less operations when n 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 follow-up: Changed 2 years ago by 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

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 2 years ago by tscrim

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

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 2 years ago by git

  • Commit changed from 32e7cf9cee379d13642eeefdf81aa89ff73efc3b to 3c82b4370366942e55885cdf88b504d4caaa6529

Branch pushed to git repo; I updated commit sha1. New commits:

3c82b43directly modify FJ

comment:19 Changed 2 years ago by mforets

  • Milestone changed from sage-7.6 to sage-8.0
  • Status changed from needs_work to needs_review

comment:20 Changed 2 years ago by mforets

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 2 years ago by tscrim

  • Reviewers set to Travis Scrimshaw
  • Status changed from needs_review to positive_review

Let's get this into Sage. Positive review.

comment:22 Changed 2 years ago by vbraun

  • 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 13 months ago by slelievre

  • Commit 3c82b4370366942e55885cdf88b504d4caaa6529 deleted
  • Description modified (diff)

Follow-up at #25028 after a bug was reported (and a fix found) at Ask Sage question 41622.

Note: See TracTickets for help on using tickets.