Ask Your Question
2

Bug in general power of a matrix

asked 2018-03-19 13:45:23 +0100

Frank Zenter gravatar image

updated 2023-01-10 02:35:20 +0100

tmonteil gravatar image

The code to question 35658 gives a wrong answer i.e. for the matrix

A=matrix([[4,1,2],[0,2,-4],[0,1,6]])

Where can I find an improvement?

For $k\in {\bf N}$ the $k$-th power of a matrix $A\in {\bf R}^{n\times n}$ satisfies by definition the recursion $A^{k+1}=A*A^k$ and the initial condition $A^0=I$. Substitution of $k=0$ and $k=1$ should therefore give the identity matrix $I$ resp. the matrix $A$. The following live code does (currently) not give the expected answer.

Btw: One of the M-programs gives the correct answer.

edit retag flag offensive close merge delete

Comments

I bet this method will work if you do A=matrix(QQ,[[4,1,2],[0,2,-4],[0,1,6]]).

kcrisman gravatar imagekcrisman ( 2018-03-19 18:25:59 +0100 )edit

Also, if you read the other ticket carefully, you will see this is already in Sage now: see here for a live example. Just define as completely normal and take the power.

kcrisman gravatar imagekcrisman ( 2018-03-19 18:29:11 +0100 )edit
1

I finally got the point of the question. (Too many hidden links in a simple question. Please always give code here, not elsewhere with a to be copied html address...)

Bug:

sage: k = var( 'k' )
sage: A = matrix( QQbar, 3, 3, [[4,1,0],[0,4,0],[0,0,4]] )
sage: A
[4 1 0]
[0 4 0]
[0 0 4]
sage: A.jordan_form()
[4 1|0]
[0 4|0]
[---+-]
[0 0|4]
sage: A^k
[        4^k 4^(k - 1)*k           0]
[          0         4^k           0]
[          0           0           0]

The Jordan blocks are not correctly parsed, the $0$ in the bottom right corner...

dan_fulea gravatar imagedan_fulea ( 2018-03-23 01:01:48 +0100 )edit

Thanks for this simpler example, which put me on the track to pin down the bug.

The problem was with FJ.set_block(k, k, ...) in the code of _matrix_power_symbolic.

See my answer, with a replacement of k, k by the appropriate index where the block needs to go.

slelievre gravatar imageslelievre ( 2018-03-23 11:09:57 +0100 )edit

1 Answer

Sort by » oldest newest most voted
4

answered 2018-03-23 11:02:24 +0100

slelievre gravatar image

updated 2018-04-03 09:13:17 +0100

Thank you for pointing out this problem.

The computation of symbolic powers of a matrix was recently introduced in Sage, based on an idea proposed on Ask Sage:

  • compute the Jordan form, and the transformation matrix
  • compute powers of each Jordan block
  • form a block-diagonal matrix with these powers
  • change back to the old basis using the transformation matrix

The idea is good. Unfortunately there was a glitch in the implementation: the powers of the Jordan blocks were being inserted at the wrong place.

To illustrate the bug in an enlightening way, consider the matrix:

sage: n = SR.var('n')
sage: A = matrix(QQ, 3, [[2, 1, 0], [0, 2, 0], [0, 0, 3]])
sage: A
[2 1 0]
[0 2 0]
[0 0 3]

The current implementation in Sage gives:

sage: B = A^n; B
[        2^n 2^(n - 1)*n           0]
[          0         3^n           0]
[          0           0           0]

Block k is inserted starting at position (k, k) instead of its correct position.

The revised implementation, suggested below, gives:

sage: B = _matrix_power_symbolic(A, n); B
[        2^n 2^(n - 1)*n           0]
[          0         2^n           0]
[          0           0         3^n]

and can deal with the matrix in the question here.

sage: A = matrix([[4, 1, 2], [0, 2, -4], [0, 1, 6]])
sage: A
[ 4  1  2]
[ 0  2 -4]
[ 0  1  6]
sage: B = _matrix_power_symbolic(A, n); B
[                 4^n          4^(n - 1)*n        2*4^(n - 1)*n]
[                   0 -2*4^(n - 1)*n + 4^n       -4*4^(n - 1)*n]
[                   0          4^(n - 1)*n  2*4^(n - 1)*n + 4^n]

Here is a revised implementation with a fix for the glitch, along with some efficiency improvements.

def _matrix_power_symbolic(A, n):
    r"""
    Return the symbolic matrix power `A^n`

    This function implements the computation of `A^n` for symbolic `n`,
    relying on the Jordan normal form of `A`, available for exact rings
    as :meth:`jordan_form`. See [Hig2008]_, §1.2, for further details.

    INPUT:

    - ``A`` -- a square matrix over an exact field

    - ``n`` -- the symbolic exponent

    OUTPUT:

    The matrix `A^n` (with symbolic entries).

    EXAMPLES:

    Powers of two by two matrix::

        sage: n = SR.var('n')
        sage: A = matrix(QQ, [[2, -1], [1,  0]])
        sage: B = A^n
        sage: B
        [ n + 1     -n]
        [     n -n + 1]
        sage: all(A^k == B.subs({n: k}) for k in range(8))
        True

    Powers of a three by three matrix in Jordan form::

        sage: n = SR.var('n')
        sage: A = matrix(QQ, 3, [[2, 1, 0], [0, 2, 0], [0, 0, 3]])
        sage: A
        [2 1 0]
        [0 2 0]
        [0 0 3]
        sage: B = A^n; B
        [        2^n 2^(n - 1)*n           0]
        [          0         2^n           0]
        [          0           0         3^n]
        sage: all(A^k == B.subs({n: k}) for k in range(8))
        True

    Powers of a three by three matrix not in Jordan form::

        sage: A = matrix([[4, 1, 2], [0, 2, -4], [0, 1, 6]])
        sage: A
        [ 4  1  2]
        [ 0  2 -4]
        [ 0  1  6]
        sage: B = A^n
        sage: B
        [                 4^n          4^(n - 1)*n        2*4^(n - 1)*n]
        [                   0 -2*4^(n - 1)*n + 4^n       -4*4^(n - 1)*n]
        [                   0          4^(n - 1)*n  2*4^(n - 1)*n + 4^n]
        sage: [B.subs({n: k}) for k in range(4)]
        [
        [1 0 0]  [ 4  1  2]  [ 16   8  16]  [  64   48   96]
        [0 1 0]  [ 0  2 -4]  [  0   0 -32]  [   0  -32 -192]
        [0 0 1], [ 0  1  6], [  0   8  32], [   0   48  160]
        ]
        sage: all(A^k == B.subs({n: k}) for k in range(8))
        True

    TESTS:

    Testing exponentiation in the symbolic ring::

        sage: n = var('n')
        sage: A = matrix([[pi, e],[0, -2*I]])
        sage: (A^n).list()
        [pi^n,
         -(-2*I)^n/(pi*e^(-1) + 2*I*e^(-1)) + pi^n/(pi*e^(-1) + 2*I*e^(-1)),
         0,
         (-2*I)^n]

    If the base ring is inexact, the Jordan normal form is not available::

        sage: A = matrix(RDF, [[2, -1], [1,  0]])
        sage: A^n
        Traceback (most recent call last):
        ...
        ValueError: Jordan normal form not implemented over inexact rings.

    Testing exponentiation in the integer ring::

        sage: A = matrix(ZZ, [[1,-1],[-1,1]])
        sage: A^(2*n+1)
        [ 1/2*2^(2*n + 1) -1/2*2^(2*n + 1)]
        [-1/2*2^(2*n + 1)  1/2*2^(2*n + 1)]

    Check if :trac:`23215` is fixed::

        sage: a, b, k = var('a, b, k')
        sage: (matrix(2, [a, b, -b, a])^k).list()
        [1/2*(a + I*b)^k + 1/2*(a - I*b)^k,
         -1/2*I*(a + I*b)^k + 1/2*I*(a - I*b)^k,
         1/2*I*(a + I*b)^k - 1/2*I*(a - I*b)^k,
         1/2*(a + I*b)^k + 1/2*(a - I*b)^k]
    """
    from sage.rings.qqbar import AlgebraicNumber
    from sage.matrix.constructor import matrix
    from sage.functions.other import binomial
    from sage.symbolic.ring import SR
    from sage.rings.qqbar import QQbar

    got_SR = A.base_ring() == SR

    # Change to QQbar if possible
    try:
        A = A.change_ring(QQbar)
    except (TypeError, NotImplementedError):
        pass

    # Get Jordan matrix J and invertible matrix P such that A = P*J*~P
    # From that, we will compute M = J^n, and obtain A^n = P*J^n*~P
    J, P = A.jordan_form(transformation=True)

    # Where each Jordan block starts, and number of blocks
    block_start = [0] + J.subdivisions()[0]
    num_blocks = len(block_start)

    # Prepare matrix M to store `J^n`, computed by Jordan block
    M = matrix(SR, J.ncols())
    M.subdivide(J.subdivisions())

    for k in range(num_blocks):

        # Jordan block Jk, its dimension nk, the eigenvalue m
        Jk = J.subdivision(k, k)
        nk = Jk.ncols()
        mk = Jk[0,0]

        # First row of block Mk; its entries are of the form
        # D^i(f) / i! with f = x^n and D = differentiation wrt x
        if hasattr(mk, 'radical_expression'):
            mk = mk.radical_expression()
        vk = [(binomial(n, i) * mk**(n-i)).simplify_full()
              for i in range(nk)]

        # Form block Mk and insert it in M
        Mk = matrix(SR, [[SR.zero()]*i + vk[:nk-i] for i in range(nk)])
        M.set_block(block_start[k], block_start[k], Mk)

    # Change entries of P and P^-1 into symbolic expressions
    if not got_SR:
        Pinv = (~P).apply_map(AlgebraicNumber.radical_expression)
        P = P.apply_map(AlgebraicNumber.radical_expression)
    else:
        Pinv = ~P

    return P * M * Pinv

I will open a ticket on Sage's issue tracker to get the fix in Sage.

Edit: The issue is now tracked at Sage Trac ticket #25082.

edit flag offensive delete link more

Comments

1

Thanks for finding and fixing the bug. Feel free to CC me when you create the ticket.

mforets gravatar imagemforets ( 2018-03-27 20:48:55 +0100 )edit

Yes, great detective work. Sorry for confusing the issue earlier since I thought it was just that an error was raised.

kcrisman gravatar imagekcrisman ( 2018-03-28 03:56:44 +0100 )edit

Which ticket is that ?

tmonteil gravatar imagetmonteil ( 2018-04-01 02:37:54 +0100 )edit
1

Sorry it took me a while. The issue is now tracked at Sage Trac ticket #25082.

slelievre gravatar imageslelievre ( 2018-04-02 20:04:15 +0100 )edit

Fixed now !

tmonteil gravatar imagetmonteil ( 2023-01-10 02:35:03 +0100 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2018-03-19 13:45:23 +0100

Seen: 565 times

Last updated: Jun 06 '18