Eigen: Functions of matrix argument#
Matrix Exponential#
- matA.Expm()#
Computes the matrix exponential of a square matrix \(A\), which is defined by the power series \(\displaystyle \exp(A) = I + A + \frac{A^2}{2!} + \frac{A^3}{3!} + \ldots\)
See also: Wikipedia [1546], Wikipedia [1556].
See also: Eigen [216], Eigen [215].
Basic examples:
>>> from mpmath import * >>> mp.dps = 15; mp.pretty = True >>> expm(zeros(3)) [1.0 0.0 0.0] [0.0 1.0 0.0] [0.0 0.0 1.0] >>> expm(eye(3)) [2.71828182845905 0.0 0.0] [ 0.0 2.71828182845905 0.0] [ 0.0 0.0 2.71828182845905] >>> expm([[1,1,0],[1,0,1],[0,1,0]]) [ 3.86814500615414 2.26812870852145 0.841130841230196] [ 2.26812870852145 2.44114713886289 1.42699786729125] [0.841130841230196 1.42699786729125 1.6000162976327] >>> expm([[1,1,0],[1,0,1],[0,1,0]], method='pade') [ 3.86814500615414 2.26812870852145 0.841130841230196] [ 2.26812870852145 2.44114713886289 1.42699786729125] [0.841130841230196 1.42699786729125 1.6000162976327] >>> expm([[1+j, 0], [1+j,1]]) [(1.46869393991589 + 2.28735528717884j) 0.0] [ (1.03776739863568 + 3.536943175722j) (2.71828182845905 + 0.0j)]
Matrices with large entries are allowed:
>>> expm(matrix([[1,2],[2,3]])**25) [5.65024064048415e+2050488462815550 9.14228140091932e+2050488462815550] [9.14228140091932e+2050488462815550 1.47925220414035e+2050488462815551]
The identity \(\exp(A+B) = \exp(A) \exp(B)\) does not hold for noncommuting matrices:
>>> A = hilbert(3) >>> B = A + eye(3) >>> chop(mnorm(A*B - B*A)) 0.0 >>> chop(mnorm(expm(A+B) - expm(A)*expm(B))) 0.0 >>> B = A + ones(3) >>> mnorm(A*B - B*A) 1.8 >>> mnorm(expm(A+B) - expm(A)*expm(B)) 42.0927851137247
Matrix Sine#
- matA.Sinm(A)#
Calculates the sine function of the matrix.
See also: Wikipedia [1546], Wikipedia [1564].
See also: Eigen [216], Eigen [221].
The cosine of a square matrix \(A\) is defined in analogy with the matrix exponential.
\[\cos(A) = \frac{e^{iA} + e^{-iA}}{2}\]\[\sin(A) = \frac{e^{iA} - e^{-iA}}{2i}\]\[\cos^2(A) + \sin^2(A) = I\]For real \(A\), we can write \(\cos(A) = \Re ( e^{iA} )\) and \(\sin(A) = \Im ( e^{iA} )\).
\[\cosh(A) = \frac{e^{A} + e^{-A}}{2}\]\[\sinh(A) = \frac{e^{A} - e^{-A}}{2}\]Examples:
>>> from mpmath import * >>> mp.dps = 15; mp.pretty = True >>> X = eye(3) >>> cosm(X) [0.54030230586814 0.0 0.0] [ 0.0 0.54030230586814 0.0] [ 0.0 0.0 0.54030230586814] >>> X = hilbert(3) >>> cosm(X) [ 0.424403834569555 -0.316643413047167 -0.221474945949293] [-0.316643413047167 0.820646708837824 -0.127183694770039] [-0.221474945949293 -0.127183694770039 0.909236687217541] >>> X = matrix([[1+j,-2],[0,-j]]) >>> cosm(X) [(0.833730025131149 - 0.988897705762865j) (1.07485840848393 - 0.17192140544213j)] [ 0.0 (1.54308063481524 + 0.0j)]
Matrix Cosine#
- matA.Cosm()#
Calculates the cosine function of the matrix.
See also: Wikipedia [1546], Wikipedia [1564].
See also: Eigen [216], Eigen [214].
The cosine of a square matrix \(A\) is defined in analogy with the matrix exponential.
\[\cos(A) = \frac{e^{iA} + e^{-iA}}{2}\]\[\sin(A) = \frac{e^{iA} - e^{-iA}}{2i}\]\[\cos^2(A) + \sin^2(A) = I\]For real \(A\), we can write \(\cos(A) = \Re ( e^{iA} )\) and \(\sin(A) = \Im ( e^{iA} )\).
\[\cosh(A) = \frac{e^{A} + e^{-A}}{2}\]\[\sinh(A) = \frac{e^{A} - e^{-A}}{2}\]Examples:
>>> from mpmath import * >>> mp.dps = 15; mp.pretty = True >>> X = eye(3) >>> cosm(X) [0.54030230586814 0.0 0.0] [ 0.0 0.54030230586814 0.0] [ 0.0 0.0 0.54030230586814] >>> X = hilbert(3) >>> cosm(X) [ 0.424403834569555 -0.316643413047167 -0.221474945949293] [-0.316643413047167 0.820646708837824 -0.127183694770039] [-0.221474945949293 -0.127183694770039 0.909236687217541] >>> X = matrix([[1+j,-2],[0,-j]]) >>> cosm(X) [(0.833730025131149 - 0.988897705762865j) (1.07485840848393 - 0.17192140544213j)] [ 0.0 (1.54308063481524 + 0.0j)]
Matrix Hyperbolic Sine#
- matA.Sinhm()#
Calculates the hyperbolic sine function of the matrix.
See also: Wikipedia [1546], Wikipedia [1564].
See also: Eigen [216], Eigen [218].
The hyperbolic sine of a square matrix \(A\) is defined in analogy with the matrix exponential.
\[\sinh(A) = \frac{e^{A} - e^{-A}}{2}\]\[\cosh(A) = \frac{e^{A} + e^{-A}}{2}\]\[\cosh^2(A) + \sinh^2(A) = I\]Sub DemoCplxMatrixFunctions() Console.WriteLine("Hello CplxMatrixFunctions!") Dim digits = 15 Dim n As Int32 = 4 Dim A, B, C, D, E, F As New cplx_mat_t A.RandomSymmetric(n) A.Print("A: ") B = A.SinhMat() B.Print("B = Sinh(A): ") C = A.CoshMat() C.Print("C = Cosh(A): ") D = C * C - B * B D.Print("C * C - B * B: ") End Sub
Hello CplxMatrixFunctions! A: 0.254341-1.994995j, 1.110691+0.032289j, -0.311350-0.647176j, 0.547685-1.281594j; 1.110691+0.032289j, 1.435774-1.303568j, -0.641133-0.284371j, -0.827662-0.328684j; -0.311350-0.647176j, -0.641133-0.284371j, -0.217231+1.954100j, -0.776330-0.558306j; 0.547685-1.281594j, -0.827662-0.328684j, -0.776330-0.558306j, 0.652181-1.335063j; B = Sinh(A): -0.344029-0.627502j, 0.956056-0.480516j, -0.134878-0.034162j, -0.590098+0.338244j; 0.956056-0.480516j, 1.739419-3.149695j, -0.793048-0.141610j, -1.597165+0.360539j; -0.134878-0.034162j, -0.793048-0.141610j, -0.152880+1.211517j, -0.351965+0.026714j; -0.590098+0.338244j, -1.597165+0.360539j, -0.351965+0.026714j, 0.171917-0.584919j; C = Cosh(A): -0.039723+0.009230j, 0.532073-1.170311j, -0.476465+0.181642j, -1.605643-0.168838j; 0.532073-1.170311j, 2.022858-2.326341j, -0.294411-0.149198j, -1.095355+1.019715j; -0.476465+0.181642j, -0.294411-0.149198j, -0.216192+0.186437j, 0.187127+0.027892j; -1.605643-0.168838j, -1.095355+1.019715j, 0.187127+0.027892j, 0.851060-0.052593j; C * C - B * B: 1.000000+0.000000j, 0.000000+0.000000j, 0.000000+0.000000j, 0.000000+0.000000j; 0.000000+0.000000j, 1.000000+0.000000j, 0.000000+0.000000j, 0.000000+0.000000j; 0.000000+0.000000j, 0.000000+0.000000j, 1.000000+0.000000j, 0.000000+0.000000j; 0.000000+0.000000j, 0.000000+0.000000j, 0.000000+0.000000j, 1.000000+0.000000j;
Matrix Hyperbolic Cosine#
- matA.Coshm()#
See also: Wikipedia [1546], Wikipedia [1564].
See also: Eigen [216], Eigen [217].
Calculates the hyperbolic cosine function of the matrix.
The hyperbolic cosine of a square matrix \(A\) is defined in analogy with the matrix exponential.
\[\cosh(A) = \frac{e^{A} + e^{-A}}{2}\]\[\sinh(A) = \frac{e^{A} - e^{-A}}{2}\]\[\cosh^2(A) + \sinh^2(A) = I\]Sub DemoCplxMatrixFunctions() Console.WriteLine("Hello CplxMatrixFunctions!") Dim digits = 15 Dim n As Int32 = 4 Dim A, B, C, D, E, F As New cplx_mat_t A.RandomSymmetric(n) A.Print("A: ") B = A.ExpMat() B.Print("B = Exp(A): ") C = B.LogMat() C.Print("C = Log(B): ") D = B.SqrtMat() D.Print("D = Sqrt(B): ") E = D * D E.Print("E = D * D: ") B = A.SinMat() B.Print("B = Sin(A): ") C = A.CosMat() C.Print("C = Cos(A): ") D = B * B + C * C D.Print("B * B + C * C: ") B = A.SinhMat() B.Print("B = Sinh(A): ") C = A.CoshMat() C.Print("C = Cosh(A): ") D = C * C - B * B D.Print("C * C - B * B: ") End Sub
Hello CplxMatrixFunctions! A: 0.254341-1.994995j, 1.110691+0.032289j, -0.311350-0.647176j, 0.547685-1.281594j; 1.110691+0.032289j, 1.435774-1.303568j, -0.641133-0.284371j, -0.827662-0.328684j; -0.311350-0.647176j, -0.641133-0.284371j, -0.217231+1.954100j, -0.776330-0.558306j; 0.547685-1.281594j, -0.827662-0.328684j, -0.776330-0.558306j, 0.652181-1.335063j; B = Exp(A): -0.383752-0.618272j, 1.488129-1.650828j, -0.611343+0.147480j, -2.195741+0.169406j; 1.488129-1.650828j, 3.762277-5.476036j, -1.087459-0.290808j, -2.692520+1.380254j; -0.611343+0.147480j, -1.087459-0.290808j, -0.369072+1.397953j, -0.164837+0.054606j; -2.195741+0.169406j, -2.692520+1.380254j, -0.164837+0.054606j, 1.022977-0.637512j; C = Log(B): 0.254341-1.994995j, 1.110691+0.032289j, -0.311350-0.647176j, 0.547685-1.281594j; 1.110691+0.032289j, 1.435774-1.303568j, -0.641133-0.284371j, -0.827662-0.328684j; -0.311350-0.647176j, -0.641133-0.284371j, -0.217231+1.954100j, -0.776330-0.558306j; 0.547685-1.281594j, -0.827662-0.328684j, -0.776330-0.558306j, 0.652181-1.335063j; D = Sqrt(B): 0.508620-0.931943j, 0.661969-0.319593j, -0.337003-0.172068j, -0.485698-0.655636j; 0.661969-0.319593j, 2.077099-1.250104j, -0.401242-0.208424j, -0.746596+0.018077j; -0.337003-0.172068j, -0.401242-0.208424j, 0.502520+0.963369j, -0.377613-0.193978j; -0.485698-0.655636j, -0.746596+0.018077j, -0.377613-0.193978j, 1.013700-0.687531j; E = D * D: -0.383752-0.618272j, 1.488129-1.650828j, -0.611343+0.147480j, -2.195741+0.169406j; 1.488129-1.650828j, 3.762277-5.476036j, -1.087459-0.290808j, -2.692520+1.380254j; -0.611343+0.147480j, -1.087459-0.290808j, -0.369072+1.397953j, -0.164837+0.054606j; -2.195741+0.169406j, -2.692520+1.380254j, -0.164837+0.054606j, 1.022977-0.637512j; B = Sin(A): 3.185996-4.423204j, 1.572018+0.576777j, -0.236161-2.007119j, 3.744383-3.392271j; 1.572018+0.576777j, 0.829215-0.073722j, -0.452703-0.362702j, 0.213357-0.343298j; -0.236161-2.007119j, -0.452703-0.362702j, -0.274184+2.749772j, -0.899405-1.828027j; 3.744383-3.392271j, 0.213357-0.343298j, -0.899405-1.828027j, 3.157659-2.589512j; C = Cos(A): 4.524926+3.229556j, -0.735268+1.275592j, 1.749846+0.069384j, 3.245352+3.579443j; -0.735268+1.275592j, -0.328697+0.628146j, 0.172580+0.683175j, 0.290239+0.309901j; 1.749846+0.069384j, 0.172580+0.683175j, 3.469410-0.664084j, 0.996030+0.538410j; 3.245352+3.579443j, 0.290239+0.309901j, 0.996030+0.538410j, 2.603737+2.713445j; B * B + C * C: 1.000000+0.000000j, 0.000000+0.000000j, 0.000000+0.000000j, 0.000000+0.000000j; 0.000000+0.000000j, 1.000000+0.000000j, 0.000000+0.000000j, 0.000000+0.000000j; 0.000000+0.000000j, 0.000000+0.000000j, 1.000000+0.000000j, 0.000000+0.000000j; 0.000000+0.000000j, 0.000000+0.000000j, 0.000000+0.000000j, 1.000000+0.000000j; B = Sinh(A): -0.344029-0.627502j, 0.956056-0.480516j, -0.134878-0.034162j, -0.590098+0.338244j; 0.956056-0.480516j, 1.739419-3.149695j, -0.793048-0.141610j, -1.597165+0.360539j; -0.134878-0.034162j, -0.793048-0.141610j, -0.152880+1.211517j, -0.351965+0.026714j; -0.590098+0.338244j, -1.597165+0.360539j, -0.351965+0.026714j, 0.171917-0.584919j; C = Cosh(A): -0.039723+0.009230j, 0.532073-1.170311j, -0.476465+0.181642j, -1.605643-0.168838j; 0.532073-1.170311j, 2.022858-2.326341j, -0.294411-0.149198j, -1.095355+1.019715j; -0.476465+0.181642j, -0.294411-0.149198j, -0.216192+0.186437j, 0.187127+0.027892j; -1.605643-0.168838j, -1.095355+1.019715j, 0.187127+0.027892j, 0.851060-0.052593j; C * C - B * B: 1.000000+0.000000j, 0.000000+0.000000j, 0.000000+0.000000j, 0.000000+0.000000j; 0.000000+0.000000j, 1.000000+0.000000j, 0.000000+0.000000j, 0.000000+0.000000j; 0.000000+0.000000j, 0.000000+0.000000j, 1.000000+0.000000j, 0.000000+0.000000j; 0.000000+0.000000j, 0.000000+0.000000j, 0.000000+0.000000j, 1.000000+0.000000j;
Matrix Square Root#
- matA.Sqrtm()#
Computes a square root of the square matrix \(A\), i.e. returns a matrix \(B = A^{1/2}\) such that \(B^2 = A\). The square root of a matrix, if it exists, is not unique.
See also: Wikipedia [1546], Wikipedia [1563].
See also: Eigen [216], Eigen [222].
Examples
Square roots of some simple matrices:
>>> from mpmath import * >>> mp.dps = 15; mp.pretty = True >>> sqrtm([[1,0], [0,1]]) [1.0 0.0] [0.0 1.0] >>> sqrtm([[0,0], [0,0]]) [0.0 0.0] [0.0 0.0] >>> sqrtm([[2,0],[0,1]]) [1.4142135623731 0.0] [ 0.0 1.0] >>> sqrtm([[1,1],[1,0]]) [ (0.920442065259926 - 0.21728689675164j) (0.568864481005783 + 0.351577584254143j)] [(0.568864481005783 + 0.351577584254143j) (0.351577584254143 - 0.568864481005783j)] >>> sqrtm([[1,0],[0,1]]) [1.0 0.0] [0.0 1.0] >>> sqrtm([[-1,0],[0,1]]) [(0.0 - 1.0j) 0.0] [ 0.0 (1.0 + 0.0j)] >>> sqrtm([[j,0],[0,j]]) [(0.707106781186547 + 0.707106781186547j) 0.0] [ 0.0 (0.707106781186547 + 0.707106781186547j)]
A square root of a rotation matrix, giving the corresponding half-angle rotation matrix:
>>> t1 = 0.75 >>> t2 = t1 * 0.5 >>> A1 = matrix([[cos(t1), -sin(t1)], [sin(t1), cos(t1)]]) >>> A2 = matrix([[cos(t2), -sin(t2)], [sin(t2), cos(t2)]]) >>> sqrtm(A1) [0.930507621912314 -0.366272529086048] [0.366272529086048 0.930507621912314] >>> A2 [0.930507621912314 -0.366272529086048] [0.366272529086048 0.930507621912314]
The identity \((A^2)^{1/2} = A\) does not necessarily hold:
>>> A = matrix([[4,1,4],[7,8,9],[10,2,11]]) >>> sqrtm(A**2) [ 4.0 1.0 4.0] [ 7.0 8.0 9.0] [10.0 2.0 11.0] >>> sqrtm(A)**2 [ 4.0 1.0 4.0] [ 7.0 8.0 9.0] [10.0 2.0 11.0] >>> A = matrix([[-4,1,4],[7,-8,9],[10,2,11]]) >>> sqrtm(A**2) [ 7.43715112194995 -0.324127569985474 1.8481718827526] [-0.251549715716942 9.32699765900402 2.48221180985147] [ 4.11609388833616 0.775751877098258 13.017955697342] >>> chop(sqrtm(A)**2) [-4.0 1.0 4.0] [ 7.0 -8.0 9.0] [10.0 2.0 11.0]
For some matrices, a square root does not exist:
>>> sqrtm([[0,1], [0,0]]) Traceback (most recent call last): ... ZeroDivisionError: matrix is numerically singular
Two examples from the documentation for Matlab’s
sqrtm:>>> mp.dps = 15; mp.pretty = True >>> sqrtm([[7,10],[15,22]]) [1.56669890360128 1.74077655955698] [2.61116483933547 4.17786374293675] >>> >>> X = matrix(\ ... [[5,-4,1,0,0], ... [-4,6,-4,1,0], ... [1,-4,6,-4,1], ... [0,1,-4,6,-4], ... [0,0,1,-4,5]]) >>> Y = matrix(\ ... [[2,-1,-0,-0,-0], ... [-1,2,-1,0,-0], ... [0,-1,2,-1,0], ... [-0,0,-1,2,-1], ... [-0,-0,-0,-1,2]]) >>> mnorm(sqrtm(X) - Y) 4.53155328326114e-19
Matrix Logarithm#
- matA.Logm()#
Calculates the logarithm of the matrix.
See also: Wikipedia [1546], Wikipedia [1554].
See also: Eigen [216], Eigen [219].
Computes a logarithm of the square matrix \(A\), i.e. returns a matrix \(B = \log(A)\) such that \(\exp(B) = A\). The logarithm of a matrix, if it exists, is not unique.
Examples
Logarithms of some simple matrices:
>>> from mpmath import * >>> mp.dps = 15; mp.pretty = True >>> X = eye(3) >>> logm(X) [0.0 0.0 0.0] [0.0 0.0 0.0] [0.0 0.0 0.0] >>> logm(2*X) [0.693147180559945 0.0 0.0] [ 0.0 0.693147180559945 0.0] [ 0.0 0.0 0.693147180559945] >>> logm(expm(X)) [1.0 0.0 0.0] [0.0 1.0 0.0] [0.0 0.0 1.0]
A logarithm of a complex matrix:
>>> X = matrix([[2+j, 1, 3], [1-j, 1-2*j, 1], [-4, -5, j]]) >>> B = logm(X) >>> nprint(B) [ (0.808757 + 0.107759j) (2.20752 + 0.202762j) (1.07376 - 0.773874j)] [ (0.905709 - 0.107795j) (0.0287395 - 0.824993j) (0.111619 + 0.514272j)] [(-0.930151 + 0.399512j) (-2.06266 - 0.674397j) (0.791552 + 0.519839j)] >>> chop(expm(B)) [(2.0 + 1.0j) 1.0 3.0] [(1.0 - 1.0j) (1.0 - 2.0j) 1.0] [ -4.0 -5.0 (0.0 + 1.0j)]
A matrix \(X\) close to the identity matrix, for which \(\log(\exp(X)) = \exp(\log(X)) = X\) holds:
>>> X = eye(3) + hilbert(3)/4 >>> X [ 1.25 0.125 0.0833333333333333] [ 0.125 1.08333333333333 0.0625] [0.0833333333333333 0.0625 1.05] >>> logm(expm(X)) [ 1.25 0.125 0.0833333333333333] [ 0.125 1.08333333333333 0.0625] [0.0833333333333333 0.0625 1.05] >>> expm(logm(X)) [ 1.25 0.125 0.0833333333333333] [ 0.125 1.08333333333333 0.0625] [0.0833333333333333 0.0625 1.05]
A logarithm of a rotation matrix, giving back the angle of the rotation:
>>> t = 3.7 >>> A = matrix([[cos(t),sin(t)],[-sin(t),cos(t)]]) >>> chop(logm(A)) [ 0.0 -2.58318530717959] [2.58318530717959 0.0] >>> (2*pi-t) 2.58318530717959
For some matrices, a logarithm does not exist:
>>> logm([[1,0], [0,0]]) Traceback (most recent call last): ... ZeroDivisionError: matrix is numerically singular
Logarithm of a matrix with large entries:
>>> logm(hilbert(3) * 10**20).apply(re) [ 45.5597513593433 1.27721006042799 0.317662687717978] [ 1.27721006042799 42.5222778973542 2.24003708791604] [0.317662687717978 2.24003708791604 42.395212822267]
Matrix power#
- mat.Powm(r)#
Computes \(A^r = \exp(A \log r)\) for a matrix \(A\) and complex number \(r\).
See also: Wikipedia [1546], Wikipedia [1556], Wikipedia [1554].
See also: Eigen [216], Eigen [220].
Examples
Powers and inverse powers of a matrix:
>>> from mpmath import * >>> mp.dps = 15; mp.pretty = True >>> A = matrix([[4,1,4],[7,8,9],[10,2,11]]) >>> powm(A, 2) [ 63.0 20.0 69.0] [174.0 89.0 199.0] [164.0 48.0 179.0] >>> chop(powm(powm(A, 4), 1/4.)) [ 4.0 1.0 4.0] [ 7.0 8.0 9.0] [10.0 2.0 11.0] >>> powm(extraprec(20)(powm)(A, -4), -1/4.) [ 4.0 1.0 4.0] [ 7.0 8.0 9.0] [10.0 2.0 11.0] >>> chop(powm(powm(A, 1+0.5j), 1/(1+0.5j))) [ 4.0 1.0 4.0] [ 7.0 8.0 9.0] [10.0 2.0 11.0] >>> powm(extraprec(5)(powm)(A, -1.5), -1/(1.5)) [ 4.0 1.0 4.0] [ 7.0 8.0 9.0] [10.0 2.0 11.0]
A Fibonacci-generating matrix:
>>> powm([[1,1],[1,0]], 10) [89.0 55.0] [55.0 34.0] >>> fib(10) 55.0 >>> powm([[1,1],[1,0]], 6.5) [(16.5166626964253 - 0.0121089837381789j) (10.2078589271083 + 0.0195927472575932j)] [(10.2078589271083 + 0.0195927472575932j) (6.30880376931698 - 0.0317017309957721j)] >>> (phi**6.5 - (1-phi)**6.5)/sqrt(5) (10.2078589271083 - 0.0195927472575932j) >>> powm([[1,1],[1,0]], 6.2) [ (14.3076953002666 - 0.008222855781077j) (8.81733464837593 + 0.0133048601383712j)] [(8.81733464837593 + 0.0133048601383712j) (5.49036065189071 - 0.0215277159194482j)] >>> (phi**6.2 - (1-phi)**6.2)/sqrt(5) (8.81733464837593 - 0.0133048601383712j)