Revision history [back]
 | 1 | initial version |
My answer elaborates on using the Jordan normal form to implement a broad class of matrix functions, the matrix power included. It relies on the jordan_form() method, which is implemented in exact rings.
The mathematical preliminaries can be found in extending scalar function for matrix functions and in the excellent book Functions of Matrices by N. Higham.
Failure of a straightforward approach
(using πππππΌπππππππππππ½.πΊ,ππππππππ³πππ:πΈπΆπ·πΌβ―π·πΆβ―π·πΎ)
Let
A = matrix(QQ, [[2, -1], [1, 0]])
var('k')
A^k
we get
sage: NotImplementedError: non-integral exponents not supported
Function of a matrix using Jordan normal form
The code below implements $f(A)$ using the Jordan normal form of $A$, see [N. Higham, Functions of Matrices, Sec. 1.2] for details.
def matrix_function_Jordan(A, f):
# returns jordan matrix J and invertlbe matrix P such that A = P*J*~P
[J, P] = A.jordan_form(transformation=True);
fJ = zero_matrix(SR, J.ncols())
num_Jordan_blocks = 1+len(J.subdivisions()[0])
fJ.subdivide(J.subdivisions());
for k in range(num_Jordan_blocks):
# get Jordan block Jk
Jk = J.subdivision(k, k)
# dimension of Jordan block Jk
mk = Jk.ncols();
fJk = zero_matrix(SR, mk, mk);
# compute the first row of f(Jk)
vk = [f.derivative(x, i)(Jk[i][i])/factorial(i) for i in range(mk)]
# insert vk into each row (above the main diagonal)
for i in range(mk):
row_Jk_i = vector(SR, zero_vector(SR, i).list() + vk[0:mk-i])
fJk.set_row(i, row_Jk_i)
fJ.set_block(k, k, fJk)
fA = P*fJ*~P
return fA
Application: computing $A^k$
First, we define our matrix function:
sage: var('k'); pow_sym(x) = x^k
Let's consider the OP's example, a $2\times 2$ matrix with one Jordan block of size $2$:
sage: A = matrix(QQ, [[2, -1], [1, 0]])
This matrix is not diagonalizable.
sage: A.is_diagonalizable()
sage: False
Calling our function,
sage: matrix_function_Jordan(A, pow_sym)
gives $\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rr} k + 1 & -k \\ k & -k + 1 \end{array}\right)$.
Pitfalls
A limitation of this approach is that we need the Jordan form. Notice that it is only implemented for exact rings (the Jordan form is unstable for inexact rings). Moreover, the spectrum should belong to the same ring. Consider for instance:
sage: A = matrix(QQ, [[0,4,1],[-1,1,5],[2,0,-89]]);
sage: [J, P] = A.jordan_form(transformation=True)
gives:
RuntimeError: Some eigenvalue does not exist in Rational Field.
This is not a piece of cake for WolframAlpha either: this case givesStandard computation time exceeded...
| | 2 | No.2 Revision |
My answer elaborates on using the Jordan normal form to implement a broad class of matrix functions, the matrix power included. It relies on the jordan_form() method, which is implemented in exact rings.
The mathematical preliminaries can be found in extending scalar function for matrix functions and in the excellent book Functions of Matrices by N. Higham.
Failure of a straightforward approach
(using πππππΌπππππππππππ½.πΊ,ππππππππ³πππ:πΈπΆπ·πΌβ―π·πΆβ―π·πΎ)
Let
A = matrix(QQ, [[2, -1], [1, 0]])
var('k')
A^k
we get
sage: NotImplementedError: non-integral exponents not supported
Function of a matrix using Jordan normal form
The code below implements $f(A)$ using the Jordan normal form of $A$, see [N. Higham, Functions of Matrices, Sec. 1.2] for details.
def matrix_function_Jordan(A, f):
# returns jordan matrix J and invertlbe invertible matrix P such that A = P*J*~P
[J, P] = A.jordan_form(transformation=True);
fJ = zero_matrix(SR, J.ncols())
num_Jordan_blocks = 1+len(J.subdivisions()[0])
fJ.subdivide(J.subdivisions());
for k in range(num_Jordan_blocks):
# get Jordan block Jk
Jk = J.subdivision(k, k)
# dimension of Jordan block Jk
mk = Jk.ncols();
fJk = zero_matrix(SR, mk, mk);
# compute the first row of f(Jk)
vk = [f.derivative(x, i)(Jk[i][i])/factorial(i) for i in range(mk)]
# insert vk into each row (above the main diagonal)
for i in range(mk):
row_Jk_i = vector(SR, zero_vector(SR, i).list() + vk[0:mk-i])
fJk.set_row(i, row_Jk_i)
fJ.set_block(k, k, fJk)
fA = P*fJ*~P
return fA
Application: computing $A^k$
First, we define our matrix function:
sage: var('k'); pow_sym(x) = x^k
Let's consider the OP's example, a $2\times 2$ matrix with one Jordan block of size $2$:
sage: A = matrix(QQ, [[2, -1], [1, 0]])
This matrix is not diagonalizable.
sage: A.is_diagonalizable()
sage: False
Calling our function,
sage: matrix_function_Jordan(A, pow_sym)
gives $\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rr} k + 1 & -k \\ k & -k + 1 \end{array}\right)$.
Pitfalls
A limitation of this approach is that we need the Jordan form. Notice that it is only implemented for exact rings (the Jordan form is unstable for inexact rings). Moreover, the spectrum should belong to the same ring. Consider for instance:
sage: A = matrix(QQ, [[0,4,1],[-1,1,5],[2,0,-89]]);
sage: [J, P] = A.jordan_form(transformation=True)
gives:
RuntimeError: Some eigenvalue does not exist in Rational Field.
This is not a piece of cake for WolframAlpha either: this case givesStandard computation time exceeded...
| | 3 | No.3 Revision |
My answer elaborates on using the Jordan normal form to implement a broad class of matrix functions, the matrix power included. It relies on the jordan_form() method, which is implemented in available for exact rings.
The mathematical preliminaries can be found in extending scalar function for matrix functions and in the excellent book Functions of Matrices by N. Higham.
Failure of a straightforward approach
(using πππππΌπππππππππππ½.πΊ,ππππππππ³πππ:πΈπΆπ·πΌβ―π·πΆβ―π·πΎ)
Let
A = matrix(QQ, [[2, -1], [1, 0]])
var('k')
A^k
we get
sage: NotImplementedError: non-integral exponents not supported
Function of a matrix using Jordan normal form
The code below implements $f(A)$ using the Jordan normal form of $A$, see [N. Higham, Functions of Matrices, Sec. 1.2] for details.
def matrix_function_Jordan(A, f):
# returns jordan matrix J and invertible matrix P such that A = P*J*~P
[J, P] = A.jordan_form(transformation=True);
fJ = zero_matrix(SR, J.ncols())
num_Jordan_blocks = 1+len(J.subdivisions()[0])
fJ.subdivide(J.subdivisions());
for k in range(num_Jordan_blocks):
# get Jordan block Jk
Jk = J.subdivision(k, k)
# dimension of Jordan block Jk
mk = Jk.ncols();
fJk = zero_matrix(SR, mk, mk);
# compute the first row of f(Jk)
vk = [f.derivative(x, i)(Jk[i][i])/factorial(i) for i in range(mk)]
# insert vk into each row (above the main diagonal)
for i in range(mk):
row_Jk_i = vector(SR, zero_vector(SR, i).list() + vk[0:mk-i])
fJk.set_row(i, row_Jk_i)
fJ.set_block(k, k, fJk)
fA = P*fJ*~P
return fA
Application: computing $A^k$
First, we define our matrix function:
sage: var('k'); pow_sym(x) = x^k
Let's consider the OP's example, a $2\times 2$ matrix with one Jordan block of size $2$:
sage: A = matrix(QQ, [[2, -1], [1, 0]])
This matrix is not diagonalizable.
sage: A.is_diagonalizable()
sage: False
Calling our function,
sage: matrix_function_Jordan(A, pow_sym)
gives $\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rr} k + 1 & -k \\ k & -k + 1 \end{array}\right)$.
Pitfalls
A limitation of this approach is that we need the Jordan form. Notice that it is only implemented for exact rings (the Jordan form is unstable for inexact rings). Moreover, the spectrum should belong to the same ring. Consider for instance:
sage: A = matrix(QQ, [[0,4,1],[-1,1,5],[2,0,-89]]);
sage: [J, P] = A.jordan_form(transformation=True)
gives:
RuntimeError: Some eigenvalue does not exist in Rational Field.
This is not a piece of cake for WolframAlpha either: this case givesStandard computation time exceeded...
| | 4 | No.4 Revision |
My answer elaborates on using the Jordan normal form to implement a broad class of matrix functions, the matrix power included. It relies on the jordan_form() method, available for exact rings.
The mathematical preliminaries can be found in extending scalar function for matrix functions and in the excellent book Functions of Matrices by N. Higham.
Failure of a straightforward approach
(using πππππΌπππππππππππ½.πΊ,ππππππππ³πππ:πΈπΆπ·πΌβ―π·πΆβ―π·πΎ)
Let
A = matrix(QQ, [[2, -1], [1, 0]])
var('k')
A^k
we get
sage: NotImplementedError: non-integral exponents not supported
Function of a matrix using Jordan normal form
The code below implements $f(A)$ using the Jordan normal form of $A$, see [N. Higham, Functions of Matrices, Sec. 1.2] for details.
def matrix_function_Jordan(A, f):
# returns jordan matrix J and invertible matrix P such that A = P*J*~P
[J, P] = A.jordan_form(transformation=True);
fJ = zero_matrix(SR, J.ncols())
num_Jordan_blocks = 1+len(J.subdivisions()[0])
fJ.subdivide(J.subdivisions());
for k in range(num_Jordan_blocks):
# get Jordan block Jk
Jk = J.subdivision(k, k)
# dimension of Jordan block Jk
mk = Jk.ncols();
fJk = zero_matrix(SR, mk, mk);
# compute the first row of f(Jk)
vk = [f.derivative(x, i)(Jk[i][i])/factorial(i) for i in range(mk)]
# insert vk into each row (above the main diagonal)
for i in range(mk):
row_Jk_i = vector(SR, zero_vector(SR, i).list() + vk[0:mk-i])
fJk.set_row(i, row_Jk_i)
fJ.set_block(k, k, fJk)
fA = P*fJ*~P
return fA
Application: computing $A^k$
First, we define our matrix function:
sage: var('k'); pow_sym(x) = x^k
Let's consider the OP's example, a $2\times 2$ matrix with one Jordan block of size $2$:
sage: A = matrix(QQ, [[2, -1], [1, 0]])
This matrix is not diagonalizable.
sage: A.is_diagonalizable()
sage: False
Calling our function,
sage: matrix_function_Jordan(A, pow_sym)
gives $\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rr} k + 1 & -k \\ k & -k + 1 \end{array}\right)$.
Pitfalls
A limitation of this approach is that we need the Jordan form. Notice that it is only implemented for exact rings (the Jordan form is unstable for inexact rings). Moreover, the spectrum should belong to the same ring. Consider for instance:
sage: A = matrix(QQ, [[0,4,1],[-1,1,5],[2,0,-89]]);
sage: [J, P] = A.jordan_form(transformation=True)
gives:
RuntimeError: Some eigenvalue does not exist in Rational Field.
This is not a piece of cake for WolframAlpha either: this case giveswe get Standard computation time exceeded...
| | 5 | No.5 Revision |
My answer elaborates on using the Jordan normal form to implement a broad class of matrix functions, the matrix power included. It relies on the jordan_form() method, available for exact rings.
The mathematical preliminaries can be found in extending scalar function for matrix functions and in the excellent book Functions of Matrices by N. Higham.
Failure of a straightforward approach
(using πππππΌπππππππππππ½.πΊ,ππππππππ³πππ:πΈπΆπ·πΌβ―π·πΆβ―π·πΎ)
Let
A = matrix(QQ, [[2, -1], [1, 0]])
var('k')
A^k
we get
sage: NotImplementedError: non-integral exponents not supported
Function of a matrix using Jordan normal form
The code below implements $f(A)$ using the Jordan normal form of $A$, see [N. Higham, Functions of Matrices, Sec. 1.2] for details.
def matrix_function_Jordan(A, f):
# returns jordan matrix J and invertible matrix P such that A = P*J*~P
[J, P] = A.jordan_form(transformation=True);
A.jordan_form(transformation=True)
fJ = zero_matrix(SR, J.ncols())
num_Jordan_blocks = 1+len(J.subdivisions()[0])
fJ.subdivide(J.subdivisions());
fJ.subdivide(J.subdivisions())
for k in range(num_Jordan_blocks):
# get Jordan block Jk
Jk = J.subdivision(k, k)
# dimension of Jordan block Jk
mk = Jk.ncols();
fJk = zero_matrix(SR, mk, mk);
mk)
# compute the first row of f(Jk)
vk = [f.derivative(x, i)(Jk[i][i])/factorial(i) for i in range(mk)]
# insert vk into each row (above the main diagonal)
for i in range(mk):
row_Jk_i = vector(SR, zero_vector(SR, i).list() + vk[0:mk-i])
fJk.set_row(i, row_Jk_i)
fJ.set_block(k, k, fJk)
fA = P*fJ*~P
return fA
Application: computing $A^k$
First, we define our matrix function:
sage: var('k'); pow_sym(x) = x^k
Let's consider the OP's example, a $2\times 2$ matrix with one Jordan block of size $2$:
sage: A = matrix(QQ, [[2, -1], [1, 0]])
This matrix is not diagonalizable.
sage: A.is_diagonalizable()
sage: False
Calling our function,
sage: matrix_function_Jordan(A, pow_sym)
gives $\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rr} k + 1 & -k \\ k & -k + 1 \end{array}\right)$.
Pitfalls
A limitation of this approach is that we need the Jordan form. Notice that it is only implemented for exact rings (the Jordan form is unstable for inexact rings). Moreover, the spectrum should belong to the same ring. Consider for instance:
sage: A = matrix(QQ, [[0,4,1],[-1,1,5],[2,0,-89]]);
sage: [J, P] = A.jordan_form(transformation=True)
gives:
RuntimeError: Some eigenvalue does not exist in Rational Field.
This is not a piece of cake for WolframAlpha either: we get Standard computation time exceeded...
| | 6 | No.6 Revision |
My answer elaborates on using the Jordan normal form to implement a broad class of matrix functions, the matrix power included. It relies on the jordan_form() method, available for exact rings.
The mathematical preliminaries can be found in extending scalar function for matrix functions and in the excellent book Functions of Matrices by N. Higham.
Failure of a straightforward approach
(using πππππΌπππππππππππ½.πΊ,ππππππππ³πππ:πΈπΆπ·πΌβ―π·πΆβ―π·πΎ)
Let
sage: A = matrix(QQ, [[2, -1], [1, 0]])
sage: var('k')
sage: A^k
we get
sage: NotImplementedError: non-integral exponents not supported
Function of a matrix using Jordan normal form
The code below implements $f(A)$ using the Jordan normal form of $A$, see [N. Higham, Functions of Matrices, Sec. 1.2] for details.
def matrix_function_Jordan(A, f):
# returns jordan matrix J and invertible matrix P such that A = P*J*~P
[J, P] = A.jordan_form(transformation=True)
fJ = zero_matrix(SR, J.ncols())
num_Jordan_blocks = 1+len(J.subdivisions()[0])
fJ.subdivide(J.subdivisions())
for k in range(num_Jordan_blocks):
# get Jordan block Jk
Jk = J.subdivision(k, k)
# dimension of Jordan block Jk
mk = Jk.ncols();
fJk = zero_matrix(SR, mk, mk)
# compute the first row of f(Jk)
vk = [f.derivative(x, i)(Jk[i][i])/factorial(i) for i in range(mk)]
# insert vk into each row (above the main diagonal)
for i in range(mk):
row_Jk_i = vector(SR, zero_vector(SR, i).list() + vk[0:mk-i])
fJk.set_row(i, row_Jk_i)
fJ.set_block(k, k, fJk)
fA = P*fJ*~P
return fA
Application: computing $A^k$
First, we define our matrix function:
sage: var('k'); pow_sym(x) = x^k
Let's consider the OP's example, a $2\times 2$ matrix with one Jordan block of size $2$:
sage: A = matrix(QQ, [[2, -1], [1, 0]])
This matrix is not diagonalizable.
sage: A.is_diagonalizable()
sage: False
Calling our function,
sage: matrix_function_Jordan(A, pow_sym)
gives $\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rr} k + 1 & -k \\ k & -k + 1 \end{array}\right)$.
Pitfalls
A limitation of this approach is that we need the Jordan form. Notice that it is only implemented for exact rings (the Jordan form is unstable for inexact rings). Moreover, the spectrum should belong to the same ring. Consider for instance:
sage: A = matrix(QQ, [[0,4,1],[-1,1,5],[2,0,-89]]);
sage: [J, P] = A.jordan_form(transformation=True)
gives:
RuntimeError: Some eigenvalue does not exist in Rational Field.
This is not a piece of cake for WolframAlpha either: we get Standard computation time exceeded...
