Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/linalg/_expm_frechet.py: 9%
152 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1"""Frechet derivative of the matrix exponential."""
2import numpy as np
3import scipy.linalg
5__all__ = ['expm_frechet', 'expm_cond']
8def expm_frechet(A, E, method=None, compute_expm=True, check_finite=True):
9 """
10 Frechet derivative of the matrix exponential of A in the direction E.
12 Parameters
13 ----------
14 A : (N, N) array_like
15 Matrix of which to take the matrix exponential.
16 E : (N, N) array_like
17 Matrix direction in which to take the Frechet derivative.
18 method : str, optional
19 Choice of algorithm. Should be one of
21 - `SPS` (default)
22 - `blockEnlarge`
24 compute_expm : bool, optional
25 Whether to compute also `expm_A` in addition to `expm_frechet_AE`.
26 Default is True.
27 check_finite : bool, optional
28 Whether to check that the input matrix contains only finite numbers.
29 Disabling may give a performance gain, but may result in problems
30 (crashes, non-termination) if the inputs do contain infinities or NaNs.
32 Returns
33 -------
34 expm_A : ndarray
35 Matrix exponential of A.
36 expm_frechet_AE : ndarray
37 Frechet derivative of the matrix exponential of A in the direction E.
38 For ``compute_expm = False``, only `expm_frechet_AE` is returned.
40 See Also
41 --------
42 expm : Compute the exponential of a matrix.
44 Notes
45 -----
46 This section describes the available implementations that can be selected
47 by the `method` parameter. The default method is *SPS*.
49 Method *blockEnlarge* is a naive algorithm.
51 Method *SPS* is Scaling-Pade-Squaring [1]_.
52 It is a sophisticated implementation which should take
53 only about 3/8 as much time as the naive implementation.
54 The asymptotics are the same.
56 .. versionadded:: 0.13.0
58 References
59 ----------
60 .. [1] Awad H. Al-Mohy and Nicholas J. Higham (2009)
61 Computing the Frechet Derivative of the Matrix Exponential,
62 with an application to Condition Number Estimation.
63 SIAM Journal On Matrix Analysis and Applications.,
64 30 (4). pp. 1639-1657. ISSN 1095-7162
66 Examples
67 --------
68 >>> import numpy as np
69 >>> from scipy import linalg
70 >>> rng = np.random.default_rng()
72 >>> A = rng.standard_normal((3, 3))
73 >>> E = rng.standard_normal((3, 3))
74 >>> expm_A, expm_frechet_AE = linalg.expm_frechet(A, E)
75 >>> expm_A.shape, expm_frechet_AE.shape
76 ((3, 3), (3, 3))
78 Create a 6x6 matrix containg [[A, E], [0, A]]:
80 >>> M = np.zeros((6, 6))
81 >>> M[:3, :3] = A
82 >>> M[:3, 3:] = E
83 >>> M[3:, 3:] = A
85 >>> expm_M = linalg.expm(M)
86 >>> np.allclose(expm_A, expm_M[:3, :3])
87 True
88 >>> np.allclose(expm_frechet_AE, expm_M[:3, 3:])
89 True
91 """
92 if check_finite:
93 A = np.asarray_chkfinite(A)
94 E = np.asarray_chkfinite(E)
95 else:
96 A = np.asarray(A)
97 E = np.asarray(E)
98 if A.ndim != 2 or A.shape[0] != A.shape[1]:
99 raise ValueError('expected A to be a square matrix')
100 if E.ndim != 2 or E.shape[0] != E.shape[1]:
101 raise ValueError('expected E to be a square matrix')
102 if A.shape != E.shape:
103 raise ValueError('expected A and E to be the same shape')
104 if method is None:
105 method = 'SPS'
106 if method == 'SPS':
107 expm_A, expm_frechet_AE = expm_frechet_algo_64(A, E)
108 elif method == 'blockEnlarge':
109 expm_A, expm_frechet_AE = expm_frechet_block_enlarge(A, E)
110 else:
111 raise ValueError('Unknown implementation %s' % method)
112 if compute_expm:
113 return expm_A, expm_frechet_AE
114 else:
115 return expm_frechet_AE
118def expm_frechet_block_enlarge(A, E):
119 """
120 This is a helper function, mostly for testing and profiling.
121 Return expm(A), frechet(A, E)
122 """
123 n = A.shape[0]
124 M = np.vstack([
125 np.hstack([A, E]),
126 np.hstack([np.zeros_like(A), A])])
127 expm_M = scipy.linalg.expm(M)
128 return expm_M[:n, :n], expm_M[:n, n:]
131"""
132Maximal values ell_m of ||2**-s A|| such that the backward error bound
133does not exceed 2**-53.
134"""
135ell_table_61 = (
136 None,
137 # 1
138 2.11e-8,
139 3.56e-4,
140 1.08e-2,
141 6.49e-2,
142 2.00e-1,
143 4.37e-1,
144 7.83e-1,
145 1.23e0,
146 1.78e0,
147 2.42e0,
148 # 11
149 3.13e0,
150 3.90e0,
151 4.74e0,
152 5.63e0,
153 6.56e0,
154 7.52e0,
155 8.53e0,
156 9.56e0,
157 1.06e1,
158 1.17e1,
159 )
162# The b vectors and U and V are copypasted
163# from scipy.sparse.linalg.matfuncs.py.
164# M, Lu, Lv follow (6.11), (6.12), (6.13), (3.3)
166def _diff_pade3(A, E, ident):
167 b = (120., 60., 12., 1.)
168 A2 = A.dot(A)
169 M2 = np.dot(A, E) + np.dot(E, A)
170 U = A.dot(b[3]*A2 + b[1]*ident)
171 V = b[2]*A2 + b[0]*ident
172 Lu = A.dot(b[3]*M2) + E.dot(b[3]*A2 + b[1]*ident)
173 Lv = b[2]*M2
174 return U, V, Lu, Lv
177def _diff_pade5(A, E, ident):
178 b = (30240., 15120., 3360., 420., 30., 1.)
179 A2 = A.dot(A)
180 M2 = np.dot(A, E) + np.dot(E, A)
181 A4 = np.dot(A2, A2)
182 M4 = np.dot(A2, M2) + np.dot(M2, A2)
183 U = A.dot(b[5]*A4 + b[3]*A2 + b[1]*ident)
184 V = b[4]*A4 + b[2]*A2 + b[0]*ident
185 Lu = (A.dot(b[5]*M4 + b[3]*M2) +
186 E.dot(b[5]*A4 + b[3]*A2 + b[1]*ident))
187 Lv = b[4]*M4 + b[2]*M2
188 return U, V, Lu, Lv
191def _diff_pade7(A, E, ident):
192 b = (17297280., 8648640., 1995840., 277200., 25200., 1512., 56., 1.)
193 A2 = A.dot(A)
194 M2 = np.dot(A, E) + np.dot(E, A)
195 A4 = np.dot(A2, A2)
196 M4 = np.dot(A2, M2) + np.dot(M2, A2)
197 A6 = np.dot(A2, A4)
198 M6 = np.dot(A4, M2) + np.dot(M4, A2)
199 U = A.dot(b[7]*A6 + b[5]*A4 + b[3]*A2 + b[1]*ident)
200 V = b[6]*A6 + b[4]*A4 + b[2]*A2 + b[0]*ident
201 Lu = (A.dot(b[7]*M6 + b[5]*M4 + b[3]*M2) +
202 E.dot(b[7]*A6 + b[5]*A4 + b[3]*A2 + b[1]*ident))
203 Lv = b[6]*M6 + b[4]*M4 + b[2]*M2
204 return U, V, Lu, Lv
207def _diff_pade9(A, E, ident):
208 b = (17643225600., 8821612800., 2075673600., 302702400., 30270240.,
209 2162160., 110880., 3960., 90., 1.)
210 A2 = A.dot(A)
211 M2 = np.dot(A, E) + np.dot(E, A)
212 A4 = np.dot(A2, A2)
213 M4 = np.dot(A2, M2) + np.dot(M2, A2)
214 A6 = np.dot(A2, A4)
215 M6 = np.dot(A4, M2) + np.dot(M4, A2)
216 A8 = np.dot(A4, A4)
217 M8 = np.dot(A4, M4) + np.dot(M4, A4)
218 U = A.dot(b[9]*A8 + b[7]*A6 + b[5]*A4 + b[3]*A2 + b[1]*ident)
219 V = b[8]*A8 + b[6]*A6 + b[4]*A4 + b[2]*A2 + b[0]*ident
220 Lu = (A.dot(b[9]*M8 + b[7]*M6 + b[5]*M4 + b[3]*M2) +
221 E.dot(b[9]*A8 + b[7]*A6 + b[5]*A4 + b[3]*A2 + b[1]*ident))
222 Lv = b[8]*M8 + b[6]*M6 + b[4]*M4 + b[2]*M2
223 return U, V, Lu, Lv
226def expm_frechet_algo_64(A, E):
227 n = A.shape[0]
228 s = None
229 ident = np.identity(n)
230 A_norm_1 = scipy.linalg.norm(A, 1)
231 m_pade_pairs = (
232 (3, _diff_pade3),
233 (5, _diff_pade5),
234 (7, _diff_pade7),
235 (9, _diff_pade9))
236 for m, pade in m_pade_pairs:
237 if A_norm_1 <= ell_table_61[m]:
238 U, V, Lu, Lv = pade(A, E, ident)
239 s = 0
240 break
241 if s is None:
242 # scaling
243 s = max(0, int(np.ceil(np.log2(A_norm_1 / ell_table_61[13]))))
244 A = A * 2.0**-s
245 E = E * 2.0**-s
246 # pade order 13
247 A2 = np.dot(A, A)
248 M2 = np.dot(A, E) + np.dot(E, A)
249 A4 = np.dot(A2, A2)
250 M4 = np.dot(A2, M2) + np.dot(M2, A2)
251 A6 = np.dot(A2, A4)
252 M6 = np.dot(A4, M2) + np.dot(M4, A2)
253 b = (64764752532480000., 32382376266240000., 7771770303897600.,
254 1187353796428800., 129060195264000., 10559470521600.,
255 670442572800., 33522128640., 1323241920., 40840800., 960960.,
256 16380., 182., 1.)
257 W1 = b[13]*A6 + b[11]*A4 + b[9]*A2
258 W2 = b[7]*A6 + b[5]*A4 + b[3]*A2 + b[1]*ident
259 Z1 = b[12]*A6 + b[10]*A4 + b[8]*A2
260 Z2 = b[6]*A6 + b[4]*A4 + b[2]*A2 + b[0]*ident
261 W = np.dot(A6, W1) + W2
262 U = np.dot(A, W)
263 V = np.dot(A6, Z1) + Z2
264 Lw1 = b[13]*M6 + b[11]*M4 + b[9]*M2
265 Lw2 = b[7]*M6 + b[5]*M4 + b[3]*M2
266 Lz1 = b[12]*M6 + b[10]*M4 + b[8]*M2
267 Lz2 = b[6]*M6 + b[4]*M4 + b[2]*M2
268 Lw = np.dot(A6, Lw1) + np.dot(M6, W1) + Lw2
269 Lu = np.dot(A, Lw) + np.dot(E, W)
270 Lv = np.dot(A6, Lz1) + np.dot(M6, Z1) + Lz2
271 # factor once and solve twice
272 lu_piv = scipy.linalg.lu_factor(-U + V)
273 R = scipy.linalg.lu_solve(lu_piv, U + V)
274 L = scipy.linalg.lu_solve(lu_piv, Lu + Lv + np.dot((Lu - Lv), R))
275 # squaring
276 for k in range(s):
277 L = np.dot(R, L) + np.dot(L, R)
278 R = np.dot(R, R)
279 return R, L
282def vec(M):
283 """
284 Stack columns of M to construct a single vector.
286 This is somewhat standard notation in linear algebra.
288 Parameters
289 ----------
290 M : 2-D array_like
291 Input matrix
293 Returns
294 -------
295 v : 1-D ndarray
296 Output vector
298 """
299 return M.T.ravel()
302def expm_frechet_kronform(A, method=None, check_finite=True):
303 """
304 Construct the Kronecker form of the Frechet derivative of expm.
306 Parameters
307 ----------
308 A : array_like with shape (N, N)
309 Matrix to be expm'd.
310 method : str, optional
311 Extra keyword to be passed to expm_frechet.
312 check_finite : bool, optional
313 Whether to check that the input matrix contains only finite numbers.
314 Disabling may give a performance gain, but may result in problems
315 (crashes, non-termination) if the inputs do contain infinities or NaNs.
317 Returns
318 -------
319 K : 2-D ndarray with shape (N*N, N*N)
320 Kronecker form of the Frechet derivative of the matrix exponential.
322 Notes
323 -----
324 This function is used to help compute the condition number
325 of the matrix exponential.
327 See Also
328 --------
329 expm : Compute a matrix exponential.
330 expm_frechet : Compute the Frechet derivative of the matrix exponential.
331 expm_cond : Compute the relative condition number of the matrix exponential
332 in the Frobenius norm.
334 """
335 if check_finite:
336 A = np.asarray_chkfinite(A)
337 else:
338 A = np.asarray(A)
339 if len(A.shape) != 2 or A.shape[0] != A.shape[1]:
340 raise ValueError('expected a square matrix')
342 n = A.shape[0]
343 ident = np.identity(n)
344 cols = []
345 for i in range(n):
346 for j in range(n):
347 E = np.outer(ident[i], ident[j])
348 F = expm_frechet(A, E,
349 method=method, compute_expm=False, check_finite=False)
350 cols.append(vec(F))
351 return np.vstack(cols).T
354def expm_cond(A, check_finite=True):
355 """
356 Relative condition number of the matrix exponential in the Frobenius norm.
358 Parameters
359 ----------
360 A : 2-D array_like
361 Square input matrix with shape (N, N).
362 check_finite : bool, optional
363 Whether to check that the input matrix contains only finite numbers.
364 Disabling may give a performance gain, but may result in problems
365 (crashes, non-termination) if the inputs do contain infinities or NaNs.
367 Returns
368 -------
369 kappa : float
370 The relative condition number of the matrix exponential
371 in the Frobenius norm
373 See Also
374 --------
375 expm : Compute the exponential of a matrix.
376 expm_frechet : Compute the Frechet derivative of the matrix exponential.
378 Notes
379 -----
380 A faster estimate for the condition number in the 1-norm
381 has been published but is not yet implemented in SciPy.
383 .. versionadded:: 0.14.0
385 Examples
386 --------
387 >>> import numpy as np
388 >>> from scipy.linalg import expm_cond
389 >>> A = np.array([[-0.3, 0.2, 0.6], [0.6, 0.3, -0.1], [-0.7, 1.2, 0.9]])
390 >>> k = expm_cond(A)
391 >>> k
392 1.7787805864469866
394 """
395 if check_finite:
396 A = np.asarray_chkfinite(A)
397 else:
398 A = np.asarray(A)
399 if len(A.shape) != 2 or A.shape[0] != A.shape[1]:
400 raise ValueError('expected a square matrix')
402 X = scipy.linalg.expm(A)
403 K = expm_frechet_kronform(A, check_finite=False)
405 # The following norm choices are deliberate.
406 # The norms of A and X are Frobenius norms,
407 # and the norm of K is the induced 2-norm.
408 A_norm = scipy.linalg.norm(A, 'fro')
409 X_norm = scipy.linalg.norm(X, 'fro')
410 K_norm = scipy.linalg.norm(K, 2)
412 kappa = (K_norm * A_norm) / X_norm
413 return kappa