fredrikj.net / blog /
Exact matrix functions in Calcium
November 10, 2020
The latest git version of Calcium includes functions for computing the exact Jordan decomposition, exponential and logarithm of complex matrices. In this post, I will show some examples using the Python wrapper.
Jordan decomposition
Any complex $n \times n$ matrix $A$ has a Jordan decomposition $A = P J P^{-1}$ where $J$ is a block diagonal matrix (called the Jordan form of $A$) with blocks of the form $$J_i = \left(\begin{smallmatrix} \lambda_i & 1 & \\ & \lambda_i & \ddots \\ & & \ddots \end{smallmatrix}\right).$$ The Jordan form $J$ is uniquely defined up to the order of the blocks, and can be used to determine whether two matrices are similar. This decomposition also allows us to evaluate any holomorphic function $f$ of $A$ using the identity $f(A) = P f(J) P^{-1}$. Indeed, $f(J)$ can be evaluated blockwise, each block $f(J_i)$ being expressible in terms of the jet of $f$ at the eigenvalue $\lambda_i$.
Here is an example of computing the Jordan form of a matrix with Calcium:
>>> A = ca_mat([[20,77,59,40], [0,-2,-3,-2], [-10,-35,-23,-15], [2,7,3,1]])
>>> A.jordan_form()
ca_mat of size 4 x 4
[1.00000*I {a where a = I [a^2+1=0]}, 0, 0, 0]
[ 0, -1.00000*I {-a where a = I [a^2+1=0]}, 0, 0]
[ 0, 0, -2, 1]
[ 0, 0, 0, -2]
Rendered in LaTeX for easier viewing:
$$J = \begin{pmatrix} +i & 0 & 0 & 0 \\ 0 & -i & 0 & 0 \\ 0 & 0 & -2 & 1 \\ 0 & 0 & 0 & -2 \end{pmatrix}.$$The Jordan form of this (carefully chosen) matrix has two blocks of size 1 corresponding to the eigenvalues $\lambda_1 = +i$ and $\lambda_2 = -i$ and one block of size 2 corresponding to the repeated eigenvalue $\lambda_3 = -2$. We can compute the similarity transformation matrix $P$ and check that the decomposition is valid:
>>> J, P = A.jordan_form(transform=True) >>> P * J * P.inv() ca_mat of size 4 x 4 [ 20, 77, 59, 40] [ 0, -2, -3, -2] [-10, -35, -23, -15] [ 2, 7, 3, 1]
Calcium's ability to compute the Jordan decomposition and evaluate matrix functions rests on its ability to compute roots of the characteristic polynomial and to work with expressions involving these roots. It will always work as long as $A$ has rational entries, at least if the algebraic number computations do not become too expensive to run to completion. For example, here we obtain the eigenvalues of the five by five Hilbert matrix:
>>> J = ca_mat([[ca(1)/(i+j+1) for i in range(5)] for j in range(5)]).jordan_form()
>>> for e in J:
... if e != 0:
... print(e)
...
3.28793e-6 {a where a = 3.28793e-6 [266716800000*a^5-476703360000*a^4+92708406000*a^3-1022881200*a^2+307505*a-1=0]}
0.000305898 {a where a = 0.000305898 [266716800000*a^5-476703360000*a^4+92708406000*a^3-1022881200*a^2+307505*a-1=0]}
0.0114075 {a where a = 0.0114075 [266716800000*a^5-476703360000*a^4+92708406000*a^3-1022881200*a^2+307505*a-1=0]}
0.208534 {a where a = 0.208534 [266716800000*a^5-476703360000*a^4+92708406000*a^3-1022881200*a^2+307505*a-1=0]}
1.56705 {a where a = 1.56705 [266716800000*a^5-476703360000*a^4+92708406000*a^3-1022881200*a^2+307505*a-1=0]}
The following works (and looks deceptively easy, because the intermediate expressions are huge and simplification to rational entries implicitly involves computing in an algebraic number field of degree $5! = 120$):
>>> J, P = ca_mat([[ca(1)/(i+j+1) for i in range(5)] for j in range(5)]).jordan_form(transform=True)
>>> P * J * P.inv()
ca_mat of size 5 x 5
[ 1, 0.500000 {1/2}, 0.333333 {1/3}, 0.250000 {1/4}, 0.200000 {1/5}]
[0.500000 {1/2}, 0.333333 {1/3}, 0.250000 {1/4}, 0.200000 {1/5}, 0.166667 {1/6}]
[0.333333 {1/3}, 0.250000 {1/4}, 0.200000 {1/5}, 0.166667 {1/6}, 0.142857 {1/7}]
[0.250000 {1/4}, 0.200000 {1/5}, 0.166667 {1/6}, 0.142857 {1/7}, 0.125000 {1/8}]
[0.200000 {1/5}, 0.166667 {1/6}, 0.142857 {1/7}, 0.125000 {1/8}, 0.111111 {1/9}]
Calcium is not yet able to represent roots of arbitrary polynomials, and frequently fails to usefully simplify radical expressions where solutions in terms of radicals are possible, so something like the following currently succeeds only in very special cases:
>>> A = ca_mat([[log(2), log(3)], [log(4), log(5)]])
>>> J, P = A.jordan_form(transform=True)
>>> J
ca_mat of size 2 x 2
[2.46769 {(a+b+e)/2 where a = 2.63279 [Sqrt(6.93159 {b^2-2*b*e+8*d*e+e^2})], b = 1.60944 [Log(5)], c = 1.38629 [Log(4)], d = 1.09861 [Log(3)], e = 0.693147 [Log(2)]}, 0]
[0, -0.165103 {(-a+b+e)/2 where a = 2.63279 [Sqrt(6.93159 {b^2-2*b*e+8*d*e+e^2})], b = 1.60944 [Log(5)], c = 1.38629 [Log(4)], d = 1.09861 [Log(3)], e = 0.693147 [Log(2)]}]
>>> P * J * P.inv() == A
True
Matrix exponentials
Here are a handful of examples of evaluating the matrix exponential:
>>> ca_mat([[1,0,0],[0,2,0],[0,0,3]]).exp()
ca_mat of size 3 x 3
[2.71828 {a where a = 2.71828 [Exp(1)]}, 0, 0]
[ 0, 7.38906 {a where a = 7.38906 [Exp(2)]}, 0]
[ 0, 0, 20.0855 {a where a = 20.0855 [Exp(3)]}]
>>> ca_mat([[1,2],[0,3]]).exp()
ca_mat of size 2 x 2
[2.71828 {a where a = 2.71828 [Exp(1)]}, 17.3673 {b^3-b where a = 20.0855 [Exp(3)], b = 2.71828 [Exp(1)]}]
[ 0, 20.0855 {a where a = 20.0855 [Exp(3)]}]
>>> ca_mat([[1,2],[3,4]]).exp()[0,0]
51.9690 {(-a*c+11*a+b*c+11*b)/22 where a = 215.354 [Exp(5.37228 {(c+5)/2})], b = 0.689160 [Exp(-0.372281 {(-c+5)/2})], c = 5.74456 [c^2-33=0]}
>>> ca_mat([[1,2],[3,4]]).exp().det() == exp(5)
True
>>> ca_mat([[0,0,1],[1,0,0],[0,1,0]]).exp().det()
1
>>> ca_mat([[0,1,0,0,0],[0,0,2,0,0],[0,0,0,3,0],[0,0,0,0,4],[0,0,0,0,0]]).exp()
ca_mat of size 5 x 5
[1, 1, 1, 1, 1]
[0, 1, 2, 3, 4]
[0, 0, 1, 3, 6]
[0, 0, 0, 1, 4]
[0, 0, 0, 0, 1]
The exponential of a nilpotent matrix $A$ equals $P(A)$ where $P$ is a polynomial (the truncated Taylor series):
>>> A = ca_mat([[10,32,3,-13], [-8,-30,-4,11], [25,90,11,-34], [-6,-28,-5,9]])
>>> A.exp()
ca_mat of size 4 x 4
[ 9.50000 {19/2}, 29, 3, -11.5000 {-23/2}]
[-10.5000 {-21/2}, -40, -5.50000 {-11/2}, 15]
[ 28.5000 {57/2}, 109, 15, -40.5000 {-81/2}]
[-12.5000 {-25/2}, -53, -8, 19.5000 {39/2}]
>>> A**0 + A**1 + A**2/2 + A**3/6
ca_mat of size 4 x 4
[ 9.50000 {19/2}, 29, 3, -11.5000 {-23/2}]
[-10.5000 {-21/2}, -40, -5.50000 {-11/2}, 15]
[ 28.5000 {57/2}, 109, 15, -40.5000 {-81/2}]
[-12.5000 {-25/2}, -53, -8, 19.5000 {39/2}]
The exact forms of matrix functions can be extremely complicated (and unpredictably so) for moderately large $n$. To illustrate, let us generate two random $4 \times 4$ matrices with $\{-1,0,+1\}$ entries.
>>> A = ca_mat(4, 4, [randint(-1,1) for i in range(16)]) >>> B = ca_mat(4, 4, [randint(-1,1) for i in range(16)]) >>> A ca_mat of size 4 x 4 [ 1, 1, -1, 1] [ 0, -1, -1, -1] [ 1, -1, 0, 0] [-1, -1, 1, 1] >>> B ca_mat of size 4 x 4 [ 1, 0, 0, -1] [-1, 1, 0, 1] [ 0, -1, -1, -1] [-1, -1, -1, -1]
We print the top left entries of $e^A$ and $e^B$:
>>> A.exp()[0,0]
1.00868 - 0e-34*I {(a*c+a+b*c-b)/(2*c) where a = 0.404578 + 1.59831*I [Exp(0.500000 + 1.32288*I {(c+1)/2})], b = 0.404578 - 1.59831*I [Exp(0.500000 - 1.32288*I {(-c+1)/2})], c = 2.64575*I [c^2+7=0], d = 1.73205 [d^2-3=0]}
>>> B.exp()[0,0]
3.27965 + 0e-31*I {(-a*e^3*f*g^2+a*e^3*f*g-a*e^3*f*h^3+a*e^3*f*h^2+2*a*e^3*f*h-2*a*e^3*f+3*a*e^3*g^2*h^3+a*e^3*g^2*h^2-5*a*e^3*g^2*h+2*a*e^3*g^2+2*a*e^3*g*h^3+3*a*e^3*g*h^2-6*a*e^3*g*h-a*e^3*g-2*a*e^3*h^3-2*a*e^3*h^2+4*a*e^3*h-3*a*e^3-a*e^2*f*g^2+a*e^2*f*g-a*e^2*f*h^3+a*e^2*f*h^2+2*a*e^2*f*h-2*a*e^2*f+3*a*e^2*g^2*h^3+a*e^2*g^2*h^2-5*a*e^2*g^2*h+2*a*e^2*g^2+2*a*e^2*g*h^3+3*a*e^2*g*h^2-5*a*e^2*g*h-2*a*e^2*g-2*a*e^2*h^3+4*a*e^2*h-6*a*e^2+2*a*e*f*g^2-2*a*e*f*g+2*a*e*f*h^3-2*a*e*f*h^2-4*a*e*f*h+4*a*e*f-6*a*e*g^2*h^3-2*a*e*g^2*h^2+10*a*e*g^2*h-4*a*e*g^2-4*a*e*g*h^3-6*a*e*g*h^2+13*a*e*g*h+a*e*g+4*a*e*h^3+6*a*e*h^2-8*a*e*h+3*a*e-2*a*g*h+2*a*g-4*a*h^2+6*a-b*f^3*g^2*h^3+b*f^3*g^2*h-2*b*f^3*g^2+2*b*f^3*g*h^2+b*f^3*g*h-4*b*f^3*g+3*b*f^3*h^3+2*b*f^3*h^2-4*b*f^3*h+3*b*f^3-b*f^2*g^2*h^3+b*f^2*g^2*h-2*b*f^2*g^2+2*b*f^2*g*h^2+b*f^2*g*h-4*b*f^2*g+3*b*f^2*h^3+2*b*f^2*h^2-4*b*f^2*h+3*b*f^2+2*b*f*g^2*h^3-2*b*f*g^2*h+4*b*f*g^2-4*b*f*g*h^2-2*b*f*g*h+8*b*f*g-6*b*f*h^3-4*b*f*h^2+8*b*f*h-6*b*f-b*g^2*h^2+b*g^2*h+b*g^2-2*b*g*h^3+3*b*g*h-b*g+b*h^3-b*h^2-2*b*h+2*b-c*f*g^3*h^3-4*c*f*g^3*h^2+3*c*f*g^3-c*f*g^2*h^3-4*c*f*g^2*h^2+3*c*f*g^2+2*c*f*g*h^3+7*c*f*g*h^2+c*f*g*h-5*c*f*g+2*c*f*h^3+c*f*h^2-3*c*f*h+3*c*f-2*c*g^3*h^3-3*c*g^3*h^2+2*c*g^3*h+2*c*g^3-2*c*g^2*h^3-3*c*g^2*h^2+2*c*g^2*h+2*c*g^2+5*c*g*h^3+7*c*g*h^2-6*c*g*h-3*c*g+c*h^3+3*c*h^2-c*h-2*c+3*d*f*g^2*h^3+3*d*f*g^2*h^2-3*d*f*g^2*h+2*d*f*g*h^3+4*d*f*g*h^2-4*d*f*g*h-2*d*f*g-d*f*h^3-2*d*f*h^2-d*f+d*g^2*h^3+2*d*g^2*h^2-2*d*g^2*h-d*g^2+3*d*g*h^3+3*d*g*h^2-3*d*g*h-d*h^3+2*d*h-d)/(12*f*g^2*h^3-18*f*g^2*h+6*f*g^2+12*f*g*h^2-12*f*g*h-8*f*g-6*f*h^3-6*f*h^2+14*f*h-6*f+6*g^2*h^2-6*g^2*h-4*g^2+12*g*h^3-18*g*h+6*g-6*h^3+4*h^2+12*h-12) where a = 2.76183 + 1.45566*I [Exp(1.13846 + 0.485062*I {e})], b = 2.76183 - 1.45566*I [Exp(1.13846 - 0.485062*I {f})], c = 0.714243 [Exp(-0.336532 {g})], d = 0.143648 [Exp(-1.94039 {h})], e = 1.13846 + 0.485062*I [e^4-3*e^2+2*e+1=0], f = 1.13846 - 0.485062*I [f^4-3*f^2+2*f+1=0], g = -0.336532 [g^4-3*g^2+2*g+1=0], h = -1.94039 [h^4-3*h^2+2*h+1=0]}
We got lucky with $A$, while the exact form of $e^B$ is horrible! I suspect that this expression can be simplified, though this is beyond what Calcium can do at present. Superficially, the results could be tidied up a bit by rewriting the complex exponentials in terms of trigonometric functions. The 0e-34 and 0e-31 imaginary parts appear in the printed numerical values since Calcium does not know that this combination of complex exponentials is actually exactly real.
Matrix logarithms
Here are a handful of examples of evaluating the matrix logarithm:
>>> ca_mat([[1,0,0],[0,2,0],[0,0,3]]).log()
ca_mat of size 3 x 3
[0, 0, 0]
[0, 0.693147 {a where a = 0.693147 [Log(2)]}, 0]
[0, 0, 1.09861 {a where a = 1.09861 [Log(3)]}]
>>> ca_mat([[4,2],[2,4]]).log()
ca_mat of size 2 x 2
[ 1.24245 {(a+b)/2 where a = 1.79176 [Log(6)], b = 0.693147 [Log(2)]}, 0.549306 {(a-b)/2 where a = 1.79176 [Log(6)], b = 0.693147 [Log(2)]}]
[0.549306 {(a-b)/2 where a = 1.79176 [Log(6)], b = 0.693147 [Log(2)]}, 1.24245 {(a+b)/2 where a = 1.79176 [Log(6)], b = 0.693147 [Log(2)]}]
>>> ca_mat([[4,2],[2,4]]).log().det() == log(2)*(log(2)+log(3))
True
>>> ca_mat([[1,1],[0,1]]).log()
ca_mat of size 2 x 2
[0, 1]
[0, 0]
Singular matrices do not have a logarithm. Calcium detects this:
>>> ca_mat([[0,1],[0,0]]).log() Traceback (most recent call last): ... ZeroDivisionError: singular matrix
Checking that $e^{\log(A)} = A$ will work in very simple cases:
>>> ca_mat([[0,0,1],[0,1,0],[1,0,0]]).log() / (pi*I)
ca_mat of size 3 x 3
[ 0.500000 {1/2}, 0, -0.500000 {-1/2}]
[ 0, 0, 0]
[-0.500000 {-1/2}, 0, 0.500000 {1/2}]
>>> ca_mat([[0,0,1],[0,1,0],[1,0,0]]).log().exp()
ca_mat of size 3 x 3
[0, 0, 1]
[0, 1, 0]
[1, 0, 0]
Here is an example that does not work yet:
>>> ca_mat([[0,0,1],[0,2,0],[-1,0,0]]).log()
ca_mat of size 3 x 3
[ 0, 0, 1.57080 {(a)/2 where a = 3.14159 [Pi], b = I [b^2+1=0]}]
[ 0, 0.693147 {a where a = 0.693147 [Log(2)]}, 0]
[-1.57080 {(-a)/2 where a = 3.14159 [Pi], b = I [b^2+1=0]}, 0, 0]
>>> ca_mat([[0,0,1],[0,2,0],[-1,0,0]]).log().exp()
Traceback (most recent call last):
...
NotImplementedError: unable to compute matrix exponential
fredrikj.net |
Blog index |
RSS feed |
Follow me on Mastodon |
Become a sponsor