Coverage for /usr/lib/python3/dist-packages/mpmath/matrices/calculus.py: 6%
159 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1from ..libmp.backend import xrange
3# TODO: should use diagonalization-based algorithms
5class MatrixCalculusMethods(object):
7 def _exp_pade(ctx, a):
8 """
9 Exponential of a matrix using Pade approximants.
11 See G. H. Golub, C. F. van Loan 'Matrix Computations',
12 third Ed., page 572
14 TODO:
15 - find a good estimate for q
16 - reduce the number of matrix multiplications to improve
17 performance
18 """
19 def eps_pade(p):
20 return ctx.mpf(2)**(3-2*p) * \
21 ctx.factorial(p)**2/(ctx.factorial(2*p)**2 * (2*p + 1))
22 q = 4
23 extraq = 8
24 while 1:
25 if eps_pade(q) < ctx.eps:
26 break
27 q += 1
28 q += extraq
29 j = int(max(1, ctx.mag(ctx.mnorm(a,'inf'))))
30 extra = q
31 prec = ctx.prec
32 ctx.dps += extra + 3
33 try:
34 a = a/2**j
35 na = a.rows
36 den = ctx.eye(na)
37 num = ctx.eye(na)
38 x = ctx.eye(na)
39 c = ctx.mpf(1)
40 for k in range(1, q+1):
41 c *= ctx.mpf(q - k + 1)/((2*q - k + 1) * k)
42 x = a*x
43 cx = c*x
44 num += cx
45 den += (-1)**k * cx
46 f = ctx.lu_solve_mat(den, num)
47 for k in range(j):
48 f = f*f
49 finally:
50 ctx.prec = prec
51 return f*1
53 def expm(ctx, A, method='taylor'):
54 r"""
55 Computes the matrix exponential of a square matrix `A`, which is defined
56 by the power series
58 .. math ::
60 \exp(A) = I + A + \frac{A^2}{2!} + \frac{A^3}{3!} + \ldots
62 With method='taylor', the matrix exponential is computed
63 using the Taylor series. With method='pade', Pade approximants
64 are used instead.
66 **Examples**
68 Basic examples::
70 >>> from mpmath import *
71 >>> mp.dps = 15; mp.pretty = True
72 >>> expm(zeros(3))
73 [1.0 0.0 0.0]
74 [0.0 1.0 0.0]
75 [0.0 0.0 1.0]
76 >>> expm(eye(3))
77 [2.71828182845905 0.0 0.0]
78 [ 0.0 2.71828182845905 0.0]
79 [ 0.0 0.0 2.71828182845905]
80 >>> expm([[1,1,0],[1,0,1],[0,1,0]])
81 [ 3.86814500615414 2.26812870852145 0.841130841230196]
82 [ 2.26812870852145 2.44114713886289 1.42699786729125]
83 [0.841130841230196 1.42699786729125 1.6000162976327]
84 >>> expm([[1,1,0],[1,0,1],[0,1,0]], method='pade')
85 [ 3.86814500615414 2.26812870852145 0.841130841230196]
86 [ 2.26812870852145 2.44114713886289 1.42699786729125]
87 [0.841130841230196 1.42699786729125 1.6000162976327]
88 >>> expm([[1+j, 0], [1+j,1]])
89 [(1.46869393991589 + 2.28735528717884j) 0.0]
90 [ (1.03776739863568 + 3.536943175722j) (2.71828182845905 + 0.0j)]
92 Matrices with large entries are allowed::
94 >>> expm(matrix([[1,2],[2,3]])**25)
95 [5.65024064048415e+2050488462815550 9.14228140091932e+2050488462815550]
96 [9.14228140091932e+2050488462815550 1.47925220414035e+2050488462815551]
98 The identity `\exp(A+B) = \exp(A) \exp(B)` does not hold for
99 noncommuting matrices::
101 >>> A = hilbert(3)
102 >>> B = A + eye(3)
103 >>> chop(mnorm(A*B - B*A))
104 0.0
105 >>> chop(mnorm(expm(A+B) - expm(A)*expm(B)))
106 0.0
107 >>> B = A + ones(3)
108 >>> mnorm(A*B - B*A)
109 1.8
110 >>> mnorm(expm(A+B) - expm(A)*expm(B))
111 42.0927851137247
113 """
114 if method == 'pade':
115 prec = ctx.prec
116 try:
117 A = ctx.matrix(A)
118 ctx.prec += 2*A.rows
119 res = ctx._exp_pade(A)
120 finally:
121 ctx.prec = prec
122 return res
123 A = ctx.matrix(A)
124 prec = ctx.prec
125 j = int(max(1, ctx.mag(ctx.mnorm(A,'inf'))))
126 j += int(0.5*prec**0.5)
127 try:
128 ctx.prec += 10 + 2*j
129 tol = +ctx.eps
130 A = A/2**j
131 T = A
132 Y = A**0 + A
133 k = 2
134 while 1:
135 T *= A * (1/ctx.mpf(k))
136 if ctx.mnorm(T, 'inf') < tol:
137 break
138 Y += T
139 k += 1
140 for k in xrange(j):
141 Y = Y*Y
142 finally:
143 ctx.prec = prec
144 Y *= 1
145 return Y
147 def cosm(ctx, A):
148 r"""
149 Gives the cosine of a square matrix `A`, defined in analogy
150 with the matrix exponential.
152 Examples::
154 >>> from mpmath import *
155 >>> mp.dps = 15; mp.pretty = True
156 >>> X = eye(3)
157 >>> cosm(X)
158 [0.54030230586814 0.0 0.0]
159 [ 0.0 0.54030230586814 0.0]
160 [ 0.0 0.0 0.54030230586814]
161 >>> X = hilbert(3)
162 >>> cosm(X)
163 [ 0.424403834569555 -0.316643413047167 -0.221474945949293]
164 [-0.316643413047167 0.820646708837824 -0.127183694770039]
165 [-0.221474945949293 -0.127183694770039 0.909236687217541]
166 >>> X = matrix([[1+j,-2],[0,-j]])
167 >>> cosm(X)
168 [(0.833730025131149 - 0.988897705762865j) (1.07485840848393 - 0.17192140544213j)]
169 [ 0.0 (1.54308063481524 + 0.0j)]
170 """
171 B = 0.5 * (ctx.expm(A*ctx.j) + ctx.expm(A*(-ctx.j)))
172 if not sum(A.apply(ctx.im).apply(abs)):
173 B = B.apply(ctx.re)
174 return B
176 def sinm(ctx, A):
177 r"""
178 Gives the sine of a square matrix `A`, defined in analogy
179 with the matrix exponential.
181 Examples::
183 >>> from mpmath import *
184 >>> mp.dps = 15; mp.pretty = True
185 >>> X = eye(3)
186 >>> sinm(X)
187 [0.841470984807897 0.0 0.0]
188 [ 0.0 0.841470984807897 0.0]
189 [ 0.0 0.0 0.841470984807897]
190 >>> X = hilbert(3)
191 >>> sinm(X)
192 [0.711608512150994 0.339783913247439 0.220742837314741]
193 [0.339783913247439 0.244113865695532 0.187231271174372]
194 [0.220742837314741 0.187231271174372 0.155816730769635]
195 >>> X = matrix([[1+j,-2],[0,-j]])
196 >>> sinm(X)
197 [(1.29845758141598 + 0.634963914784736j) (-1.96751511930922 + 0.314700021761367j)]
198 [ 0.0 (0.0 - 1.1752011936438j)]
199 """
200 B = (-0.5j) * (ctx.expm(A*ctx.j) - ctx.expm(A*(-ctx.j)))
201 if not sum(A.apply(ctx.im).apply(abs)):
202 B = B.apply(ctx.re)
203 return B
205 def _sqrtm_rot(ctx, A, _may_rotate):
206 # If the iteration fails to converge, cheat by performing
207 # a rotation by a complex number
208 u = ctx.j**0.3
209 return ctx.sqrtm(u*A, _may_rotate) / ctx.sqrt(u)
211 def sqrtm(ctx, A, _may_rotate=2):
212 r"""
213 Computes a square root of the square matrix `A`, i.e. returns
214 a matrix `B = A^{1/2}` such that `B^2 = A`. The square root
215 of a matrix, if it exists, is not unique.
217 **Examples**
219 Square roots of some simple matrices::
221 >>> from mpmath import *
222 >>> mp.dps = 15; mp.pretty = True
223 >>> sqrtm([[1,0], [0,1]])
224 [1.0 0.0]
225 [0.0 1.0]
226 >>> sqrtm([[0,0], [0,0]])
227 [0.0 0.0]
228 [0.0 0.0]
229 >>> sqrtm([[2,0],[0,1]])
230 [1.4142135623731 0.0]
231 [ 0.0 1.0]
232 >>> sqrtm([[1,1],[1,0]])
233 [ (0.920442065259926 - 0.21728689675164j) (0.568864481005783 + 0.351577584254143j)]
234 [(0.568864481005783 + 0.351577584254143j) (0.351577584254143 - 0.568864481005783j)]
235 >>> sqrtm([[1,0],[0,1]])
236 [1.0 0.0]
237 [0.0 1.0]
238 >>> sqrtm([[-1,0],[0,1]])
239 [(0.0 - 1.0j) 0.0]
240 [ 0.0 (1.0 + 0.0j)]
241 >>> sqrtm([[j,0],[0,j]])
242 [(0.707106781186547 + 0.707106781186547j) 0.0]
243 [ 0.0 (0.707106781186547 + 0.707106781186547j)]
245 A square root of a rotation matrix, giving the corresponding
246 half-angle rotation matrix::
248 >>> t1 = 0.75
249 >>> t2 = t1 * 0.5
250 >>> A1 = matrix([[cos(t1), -sin(t1)], [sin(t1), cos(t1)]])
251 >>> A2 = matrix([[cos(t2), -sin(t2)], [sin(t2), cos(t2)]])
252 >>> sqrtm(A1)
253 [0.930507621912314 -0.366272529086048]
254 [0.366272529086048 0.930507621912314]
255 >>> A2
256 [0.930507621912314 -0.366272529086048]
257 [0.366272529086048 0.930507621912314]
259 The identity `(A^2)^{1/2} = A` does not necessarily hold::
261 >>> A = matrix([[4,1,4],[7,8,9],[10,2,11]])
262 >>> sqrtm(A**2)
263 [ 4.0 1.0 4.0]
264 [ 7.0 8.0 9.0]
265 [10.0 2.0 11.0]
266 >>> sqrtm(A)**2
267 [ 4.0 1.0 4.0]
268 [ 7.0 8.0 9.0]
269 [10.0 2.0 11.0]
270 >>> A = matrix([[-4,1,4],[7,-8,9],[10,2,11]])
271 >>> sqrtm(A**2)
272 [ 7.43715112194995 -0.324127569985474 1.8481718827526]
273 [-0.251549715716942 9.32699765900402 2.48221180985147]
274 [ 4.11609388833616 0.775751877098258 13.017955697342]
275 >>> chop(sqrtm(A)**2)
276 [-4.0 1.0 4.0]
277 [ 7.0 -8.0 9.0]
278 [10.0 2.0 11.0]
280 For some matrices, a square root does not exist::
282 >>> sqrtm([[0,1], [0,0]])
283 Traceback (most recent call last):
284 ...
285 ZeroDivisionError: matrix is numerically singular
287 Two examples from the documentation for Matlab's ``sqrtm``::
289 >>> mp.dps = 15; mp.pretty = True
290 >>> sqrtm([[7,10],[15,22]])
291 [1.56669890360128 1.74077655955698]
292 [2.61116483933547 4.17786374293675]
293 >>>
294 >>> X = matrix(\
295 ... [[5,-4,1,0,0],
296 ... [-4,6,-4,1,0],
297 ... [1,-4,6,-4,1],
298 ... [0,1,-4,6,-4],
299 ... [0,0,1,-4,5]])
300 >>> Y = matrix(\
301 ... [[2,-1,-0,-0,-0],
302 ... [-1,2,-1,0,-0],
303 ... [0,-1,2,-1,0],
304 ... [-0,0,-1,2,-1],
305 ... [-0,-0,-0,-1,2]])
306 >>> mnorm(sqrtm(X) - Y)
307 4.53155328326114e-19
309 """
310 A = ctx.matrix(A)
311 # Trivial
312 if A*0 == A:
313 return A
314 prec = ctx.prec
315 if _may_rotate:
316 d = ctx.det(A)
317 if abs(ctx.im(d)) < 16*ctx.eps and ctx.re(d) < 0:
318 return ctx._sqrtm_rot(A, _may_rotate-1)
319 try:
320 ctx.prec += 10
321 tol = ctx.eps * 128
322 Y = A
323 Z = I = A**0
324 k = 0
325 # Denman-Beavers iteration
326 while 1:
327 Yprev = Y
328 try:
329 Y, Z = 0.5*(Y+ctx.inverse(Z)), 0.5*(Z+ctx.inverse(Y))
330 except ZeroDivisionError:
331 if _may_rotate:
332 Y = ctx._sqrtm_rot(A, _may_rotate-1)
333 break
334 else:
335 raise
336 mag1 = ctx.mnorm(Y-Yprev, 'inf')
337 mag2 = ctx.mnorm(Y, 'inf')
338 if mag1 <= mag2*tol:
339 break
340 if _may_rotate and k > 6 and not mag1 < mag2 * 0.001:
341 return ctx._sqrtm_rot(A, _may_rotate-1)
342 k += 1
343 if k > ctx.prec:
344 raise ctx.NoConvergence
345 finally:
346 ctx.prec = prec
347 Y *= 1
348 return Y
350 def logm(ctx, A):
351 r"""
352 Computes a logarithm of the square matrix `A`, i.e. returns
353 a matrix `B = \log(A)` such that `\exp(B) = A`. The logarithm
354 of a matrix, if it exists, is not unique.
356 **Examples**
358 Logarithms of some simple matrices::
360 >>> from mpmath import *
361 >>> mp.dps = 15; mp.pretty = True
362 >>> X = eye(3)
363 >>> logm(X)
364 [0.0 0.0 0.0]
365 [0.0 0.0 0.0]
366 [0.0 0.0 0.0]
367 >>> logm(2*X)
368 [0.693147180559945 0.0 0.0]
369 [ 0.0 0.693147180559945 0.0]
370 [ 0.0 0.0 0.693147180559945]
371 >>> logm(expm(X))
372 [1.0 0.0 0.0]
373 [0.0 1.0 0.0]
374 [0.0 0.0 1.0]
376 A logarithm of a complex matrix::
378 >>> X = matrix([[2+j, 1, 3], [1-j, 1-2*j, 1], [-4, -5, j]])
379 >>> B = logm(X)
380 >>> nprint(B)
381 [ (0.808757 + 0.107759j) (2.20752 + 0.202762j) (1.07376 - 0.773874j)]
382 [ (0.905709 - 0.107795j) (0.0287395 - 0.824993j) (0.111619 + 0.514272j)]
383 [(-0.930151 + 0.399512j) (-2.06266 - 0.674397j) (0.791552 + 0.519839j)]
384 >>> chop(expm(B))
385 [(2.0 + 1.0j) 1.0 3.0]
386 [(1.0 - 1.0j) (1.0 - 2.0j) 1.0]
387 [ -4.0 -5.0 (0.0 + 1.0j)]
389 A matrix `X` close to the identity matrix, for which
390 `\log(\exp(X)) = \exp(\log(X)) = X` holds::
392 >>> X = eye(3) + hilbert(3)/4
393 >>> X
394 [ 1.25 0.125 0.0833333333333333]
395 [ 0.125 1.08333333333333 0.0625]
396 [0.0833333333333333 0.0625 1.05]
397 >>> logm(expm(X))
398 [ 1.25 0.125 0.0833333333333333]
399 [ 0.125 1.08333333333333 0.0625]
400 [0.0833333333333333 0.0625 1.05]
401 >>> expm(logm(X))
402 [ 1.25 0.125 0.0833333333333333]
403 [ 0.125 1.08333333333333 0.0625]
404 [0.0833333333333333 0.0625 1.05]
406 A logarithm of a rotation matrix, giving back the angle of
407 the rotation::
409 >>> t = 3.7
410 >>> A = matrix([[cos(t),sin(t)],[-sin(t),cos(t)]])
411 >>> chop(logm(A))
412 [ 0.0 -2.58318530717959]
413 [2.58318530717959 0.0]
414 >>> (2*pi-t)
415 2.58318530717959
417 For some matrices, a logarithm does not exist::
419 >>> logm([[1,0], [0,0]])
420 Traceback (most recent call last):
421 ...
422 ZeroDivisionError: matrix is numerically singular
424 Logarithm of a matrix with large entries::
426 >>> logm(hilbert(3) * 10**20).apply(re)
427 [ 45.5597513593433 1.27721006042799 0.317662687717978]
428 [ 1.27721006042799 42.5222778973542 2.24003708791604]
429 [0.317662687717978 2.24003708791604 42.395212822267]
431 """
432 A = ctx.matrix(A)
433 prec = ctx.prec
434 try:
435 ctx.prec += 10
436 tol = ctx.eps * 128
437 I = A**0
438 B = A
439 n = 0
440 while 1:
441 B = ctx.sqrtm(B)
442 n += 1
443 if ctx.mnorm(B-I, 'inf') < 0.125:
444 break
445 T = X = B-I
446 L = X*0
447 k = 1
448 while 1:
449 if k & 1:
450 L += T / k
451 else:
452 L -= T / k
453 T *= X
454 if ctx.mnorm(T, 'inf') < tol:
455 break
456 k += 1
457 if k > ctx.prec:
458 raise ctx.NoConvergence
459 finally:
460 ctx.prec = prec
461 L *= 2**n
462 return L
464 def powm(ctx, A, r):
465 r"""
466 Computes `A^r = \exp(A \log r)` for a matrix `A` and complex
467 number `r`.
469 **Examples**
471 Powers and inverse powers of a matrix::
473 >>> from mpmath import *
474 >>> mp.dps = 15; mp.pretty = True
475 >>> A = matrix([[4,1,4],[7,8,9],[10,2,11]])
476 >>> powm(A, 2)
477 [ 63.0 20.0 69.0]
478 [174.0 89.0 199.0]
479 [164.0 48.0 179.0]
480 >>> chop(powm(powm(A, 4), 1/4.))
481 [ 4.0 1.0 4.0]
482 [ 7.0 8.0 9.0]
483 [10.0 2.0 11.0]
484 >>> powm(extraprec(20)(powm)(A, -4), -1/4.)
485 [ 4.0 1.0 4.0]
486 [ 7.0 8.0 9.0]
487 [10.0 2.0 11.0]
488 >>> chop(powm(powm(A, 1+0.5j), 1/(1+0.5j)))
489 [ 4.0 1.0 4.0]
490 [ 7.0 8.0 9.0]
491 [10.0 2.0 11.0]
492 >>> powm(extraprec(5)(powm)(A, -1.5), -1/(1.5))
493 [ 4.0 1.0 4.0]
494 [ 7.0 8.0 9.0]
495 [10.0 2.0 11.0]
497 A Fibonacci-generating matrix::
499 >>> powm([[1,1],[1,0]], 10)
500 [89.0 55.0]
501 [55.0 34.0]
502 >>> fib(10)
503 55.0
504 >>> powm([[1,1],[1,0]], 6.5)
505 [(16.5166626964253 - 0.0121089837381789j) (10.2078589271083 + 0.0195927472575932j)]
506 [(10.2078589271083 + 0.0195927472575932j) (6.30880376931698 - 0.0317017309957721j)]
507 >>> (phi**6.5 - (1-phi)**6.5)/sqrt(5)
508 (10.2078589271083 - 0.0195927472575932j)
509 >>> powm([[1,1],[1,0]], 6.2)
510 [ (14.3076953002666 - 0.008222855781077j) (8.81733464837593 + 0.0133048601383712j)]
511 [(8.81733464837593 + 0.0133048601383712j) (5.49036065189071 - 0.0215277159194482j)]
512 >>> (phi**6.2 - (1-phi)**6.2)/sqrt(5)
513 (8.81733464837593 - 0.0133048601383712j)
515 """
516 A = ctx.matrix(A)
517 r = ctx.convert(r)
518 prec = ctx.prec
519 try:
520 ctx.prec += 10
521 if ctx.isint(r):
522 v = A ** int(r)
523 elif ctx.isint(r*2):
524 y = int(r*2)
525 v = ctx.sqrtm(A) ** y
526 else:
527 v = ctx.expm(r*ctx.logm(A))
528 finally:
529 ctx.prec = prec
530 v *= 1
531 return v