Coverage for /home/martinb/.local/share/virtualenvs/camcops/lib/python3.6/site-packages/scipy/linalg/_expm_frechet.py : 9%

Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
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.
39 For ``compute_expm = False``, only `expm_frechet_AE` is returned.
41 See also
42 --------
43 expm : Compute the exponential of a matrix.
45 Notes
46 -----
47 This section describes the available implementations that can be selected
48 by the `method` parameter. The default method is *SPS*.
50 Method *blockEnlarge* is a naive algorithm.
52 Method *SPS* is Scaling-Pade-Squaring [1]_.
53 It is a sophisticated implementation which should take
54 only about 3/8 as much time as the naive implementation.
55 The asymptotics are the same.
57 .. versionadded:: 0.13.0
59 References
60 ----------
61 .. [1] Awad H. Al-Mohy and Nicholas J. Higham (2009)
62 Computing the Frechet Derivative of the Matrix Exponential,
63 with an application to Condition Number Estimation.
64 SIAM Journal On Matrix Analysis and Applications.,
65 30 (4). pp. 1639-1657. ISSN 1095-7162
67 Examples
68 --------
69 >>> import scipy.linalg
70 >>> A = np.random.randn(3, 3)
71 >>> E = np.random.randn(3, 3)
72 >>> expm_A, expm_frechet_AE = scipy.linalg.expm_frechet(A, E)
73 >>> expm_A.shape, expm_frechet_AE.shape
74 ((3, 3), (3, 3))
76 >>> import scipy.linalg
77 >>> A = np.random.randn(3, 3)
78 >>> E = np.random.randn(3, 3)
79 >>> expm_A, expm_frechet_AE = scipy.linalg.expm_frechet(A, E)
80 >>> M = np.zeros((6, 6))
81 >>> M[:3, :3] = A; M[:3, 3:] = E; M[3:, 3:] = A
82 >>> expm_M = scipy.linalg.expm(M)
83 >>> np.allclose(expm_A, expm_M[:3, :3])
84 True
85 >>> np.allclose(expm_frechet_AE, expm_M[:3, 3:])
86 True
88 """
89 if check_finite:
90 A = np.asarray_chkfinite(A)
91 E = np.asarray_chkfinite(E)
92 else:
93 A = np.asarray(A)
94 E = np.asarray(E)
95 if A.ndim != 2 or A.shape[0] != A.shape[1]:
96 raise ValueError('expected A to be a square matrix')
97 if E.ndim != 2 or E.shape[0] != E.shape[1]:
98 raise ValueError('expected E to be a square matrix')
99 if A.shape != E.shape:
100 raise ValueError('expected A and E to be the same shape')
101 if method is None:
102 method = 'SPS'
103 if method == 'SPS':
104 expm_A, expm_frechet_AE = expm_frechet_algo_64(A, E)
105 elif method == 'blockEnlarge':
106 expm_A, expm_frechet_AE = expm_frechet_block_enlarge(A, E)
107 else:
108 raise ValueError('Unknown implementation %s' % method)
109 if compute_expm:
110 return expm_A, expm_frechet_AE
111 else:
112 return expm_frechet_AE
115def expm_frechet_block_enlarge(A, E):
116 """
117 This is a helper function, mostly for testing and profiling.
118 Return expm(A), frechet(A, E)
119 """
120 n = A.shape[0]
121 M = np.vstack([
122 np.hstack([A, E]),
123 np.hstack([np.zeros_like(A), A])])
124 expm_M = scipy.linalg.expm(M)
125 return expm_M[:n, :n], expm_M[:n, n:]
128"""
129Maximal values ell_m of ||2**-s A|| such that the backward error bound
130does not exceed 2**-53.
131"""
132ell_table_61 = (
133 None,
134 # 1
135 2.11e-8,
136 3.56e-4,
137 1.08e-2,
138 6.49e-2,
139 2.00e-1,
140 4.37e-1,
141 7.83e-1,
142 1.23e0,
143 1.78e0,
144 2.42e0,
145 # 11
146 3.13e0,
147 3.90e0,
148 4.74e0,
149 5.63e0,
150 6.56e0,
151 7.52e0,
152 8.53e0,
153 9.56e0,
154 1.06e1,
155 1.17e1,
156 )
159# The b vectors and U and V are copypasted
160# from scipy.sparse.linalg.matfuncs.py.
161# M, Lu, Lv follow (6.11), (6.12), (6.13), (3.3)
163def _diff_pade3(A, E, ident):
164 b = (120., 60., 12., 1.)
165 A2 = A.dot(A)
166 M2 = np.dot(A, E) + np.dot(E, A)
167 U = A.dot(b[3]*A2 + b[1]*ident)
168 V = b[2]*A2 + b[0]*ident
169 Lu = A.dot(b[3]*M2) + E.dot(b[3]*A2 + b[1]*ident)
170 Lv = b[2]*M2
171 return U, V, Lu, Lv
174def _diff_pade5(A, E, ident):
175 b = (30240., 15120., 3360., 420., 30., 1.)
176 A2 = A.dot(A)
177 M2 = np.dot(A, E) + np.dot(E, A)
178 A4 = np.dot(A2, A2)
179 M4 = np.dot(A2, M2) + np.dot(M2, A2)
180 U = A.dot(b[5]*A4 + b[3]*A2 + b[1]*ident)
181 V = b[4]*A4 + b[2]*A2 + b[0]*ident
182 Lu = (A.dot(b[5]*M4 + b[3]*M2) +
183 E.dot(b[5]*A4 + b[3]*A2 + b[1]*ident))
184 Lv = b[4]*M4 + b[2]*M2
185 return U, V, Lu, Lv
188def _diff_pade7(A, E, ident):
189 b = (17297280., 8648640., 1995840., 277200., 25200., 1512., 56., 1.)
190 A2 = A.dot(A)
191 M2 = np.dot(A, E) + np.dot(E, A)
192 A4 = np.dot(A2, A2)
193 M4 = np.dot(A2, M2) + np.dot(M2, A2)
194 A6 = np.dot(A2, A4)
195 M6 = np.dot(A4, M2) + np.dot(M4, A2)
196 U = A.dot(b[7]*A6 + b[5]*A4 + b[3]*A2 + b[1]*ident)
197 V = b[6]*A6 + b[4]*A4 + b[2]*A2 + b[0]*ident
198 Lu = (A.dot(b[7]*M6 + b[5]*M4 + b[3]*M2) +
199 E.dot(b[7]*A6 + b[5]*A4 + b[3]*A2 + b[1]*ident))
200 Lv = b[6]*M6 + b[4]*M4 + b[2]*M2
201 return U, V, Lu, Lv
204def _diff_pade9(A, E, ident):
205 b = (17643225600., 8821612800., 2075673600., 302702400., 30270240.,
206 2162160., 110880., 3960., 90., 1.)
207 A2 = A.dot(A)
208 M2 = np.dot(A, E) + np.dot(E, A)
209 A4 = np.dot(A2, A2)
210 M4 = np.dot(A2, M2) + np.dot(M2, A2)
211 A6 = np.dot(A2, A4)
212 M6 = np.dot(A4, M2) + np.dot(M4, A2)
213 A8 = np.dot(A4, A4)
214 M8 = np.dot(A4, M4) + np.dot(M4, A4)
215 U = A.dot(b[9]*A8 + b[7]*A6 + b[5]*A4 + b[3]*A2 + b[1]*ident)
216 V = b[8]*A8 + b[6]*A6 + b[4]*A4 + b[2]*A2 + b[0]*ident
217 Lu = (A.dot(b[9]*M8 + b[7]*M6 + b[5]*M4 + b[3]*M2) +
218 E.dot(b[9]*A8 + b[7]*A6 + b[5]*A4 + b[3]*A2 + b[1]*ident))
219 Lv = b[8]*M8 + b[6]*M6 + b[4]*M4 + b[2]*M2
220 return U, V, Lu, Lv
223def expm_frechet_algo_64(A, E):
224 n = A.shape[0]
225 s = None
226 ident = np.identity(n)
227 A_norm_1 = scipy.linalg.norm(A, 1)
228 m_pade_pairs = (
229 (3, _diff_pade3),
230 (5, _diff_pade5),
231 (7, _diff_pade7),
232 (9, _diff_pade9))
233 for m, pade in m_pade_pairs:
234 if A_norm_1 <= ell_table_61[m]:
235 U, V, Lu, Lv = pade(A, E, ident)
236 s = 0
237 break
238 if s is None:
239 # scaling
240 s = max(0, int(np.ceil(np.log2(A_norm_1 / ell_table_61[13]))))
241 A = A * 2.0**-s
242 E = E * 2.0**-s
243 # pade order 13
244 A2 = np.dot(A, A)
245 M2 = np.dot(A, E) + np.dot(E, A)
246 A4 = np.dot(A2, A2)
247 M4 = np.dot(A2, M2) + np.dot(M2, A2)
248 A6 = np.dot(A2, A4)
249 M6 = np.dot(A4, M2) + np.dot(M4, A2)
250 b = (64764752532480000., 32382376266240000., 7771770303897600.,
251 1187353796428800., 129060195264000., 10559470521600.,
252 670442572800., 33522128640., 1323241920., 40840800., 960960.,
253 16380., 182., 1.)
254 W1 = b[13]*A6 + b[11]*A4 + b[9]*A2
255 W2 = b[7]*A6 + b[5]*A4 + b[3]*A2 + b[1]*ident
256 Z1 = b[12]*A6 + b[10]*A4 + b[8]*A2
257 Z2 = b[6]*A6 + b[4]*A4 + b[2]*A2 + b[0]*ident
258 W = np.dot(A6, W1) + W2
259 U = np.dot(A, W)
260 V = np.dot(A6, Z1) + Z2
261 Lw1 = b[13]*M6 + b[11]*M4 + b[9]*M2
262 Lw2 = b[7]*M6 + b[5]*M4 + b[3]*M2
263 Lz1 = b[12]*M6 + b[10]*M4 + b[8]*M2
264 Lz2 = b[6]*M6 + b[4]*M4 + b[2]*M2
265 Lw = np.dot(A6, Lw1) + np.dot(M6, W1) + Lw2
266 Lu = np.dot(A, Lw) + np.dot(E, W)
267 Lv = np.dot(A6, Lz1) + np.dot(M6, Z1) + Lz2
268 # factor once and solve twice
269 lu_piv = scipy.linalg.lu_factor(-U + V)
270 R = scipy.linalg.lu_solve(lu_piv, U + V)
271 L = scipy.linalg.lu_solve(lu_piv, Lu + Lv + np.dot((Lu - Lv), R))
272 # squaring
273 for k in range(s):
274 L = np.dot(R, L) + np.dot(L, R)
275 R = np.dot(R, R)
276 return R, L
279def vec(M):
280 """
281 Stack columns of M to construct a single vector.
283 This is somewhat standard notation in linear algebra.
285 Parameters
286 ----------
287 M : 2-D array_like
288 Input matrix
290 Returns
291 -------
292 v : 1-D ndarray
293 Output vector
295 """
296 return M.T.ravel()
299def expm_frechet_kronform(A, method=None, check_finite=True):
300 """
301 Construct the Kronecker form of the Frechet derivative of expm.
303 Parameters
304 ----------
305 A : array_like with shape (N, N)
306 Matrix to be expm'd.
307 method : str, optional
308 Extra keyword to be passed to expm_frechet.
309 check_finite : bool, optional
310 Whether to check that the input matrix contains only finite numbers.
311 Disabling may give a performance gain, but may result in problems
312 (crashes, non-termination) if the inputs do contain infinities or NaNs.
314 Returns
315 -------
316 K : 2-D ndarray with shape (N*N, N*N)
317 Kronecker form of the Frechet derivative of the matrix exponential.
319 Notes
320 -----
321 This function is used to help compute the condition number
322 of the matrix exponential.
324 See also
325 --------
326 expm : Compute a matrix exponential.
327 expm_frechet : Compute the Frechet derivative of the matrix exponential.
328 expm_cond : Compute the relative condition number of the matrix exponential
329 in the Frobenius norm.
331 """
332 if check_finite:
333 A = np.asarray_chkfinite(A)
334 else:
335 A = np.asarray(A)
336 if len(A.shape) != 2 or A.shape[0] != A.shape[1]:
337 raise ValueError('expected a square matrix')
339 n = A.shape[0]
340 ident = np.identity(n)
341 cols = []
342 for i in range(n):
343 for j in range(n):
344 E = np.outer(ident[i], ident[j])
345 F = expm_frechet(A, E,
346 method=method, compute_expm=False, check_finite=False)
347 cols.append(vec(F))
348 return np.vstack(cols).T
351def expm_cond(A, check_finite=True):
352 """
353 Relative condition number of the matrix exponential in the Frobenius norm.
355 Parameters
356 ----------
357 A : 2-D array_like
358 Square input matrix with shape (N, N).
359 check_finite : bool, optional
360 Whether to check that the input matrix contains only finite numbers.
361 Disabling may give a performance gain, but may result in problems
362 (crashes, non-termination) if the inputs do contain infinities or NaNs.
364 Returns
365 -------
366 kappa : float
367 The relative condition number of the matrix exponential
368 in the Frobenius norm
370 Notes
371 -----
372 A faster estimate for the condition number in the 1-norm
373 has been published but is not yet implemented in SciPy.
375 .. versionadded:: 0.14.0
377 See also
378 --------
379 expm : Compute the exponential of a matrix.
380 expm_frechet : Compute the Frechet derivative of the matrix exponential.
382 Examples
383 --------
384 >>> from scipy.linalg import expm_cond
385 >>> A = np.array([[-0.3, 0.2, 0.6], [0.6, 0.3, -0.1], [-0.7, 1.2, 0.9]])
386 >>> k = expm_cond(A)
387 >>> k
388 1.7787805864469866
390 """
391 if check_finite:
392 A = np.asarray_chkfinite(A)
393 else:
394 A = np.asarray(A)
395 if len(A.shape) != 2 or A.shape[0] != A.shape[1]:
396 raise ValueError('expected a square matrix')
398 X = scipy.linalg.expm(A)
399 K = expm_frechet_kronform(A, check_finite=False)
401 # The following norm choices are deliberate.
402 # The norms of A and X are Frobenius norms,
403 # and the norm of K is the induced 2-norm.
404 A_norm = scipy.linalg.norm(A, 'fro')
405 X_norm = scipy.linalg.norm(X, 'fro')
406 K_norm = scipy.linalg.norm(K, 2)
408 kappa = (K_norm * A_norm) / X_norm
409 return kappa