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

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"""Compute the action of the matrix exponential.
2"""
4import numpy as np
6import scipy.linalg
7import scipy.sparse.linalg
8from scipy.sparse.linalg import aslinearoperator
9from scipy.sparse.sputils import is_pydata_spmatrix
11__all__ = ['expm_multiply']
14def _exact_inf_norm(A):
15 # A compatibility function which should eventually disappear.
16 if scipy.sparse.isspmatrix(A):
17 return max(abs(A).sum(axis=1).flat)
18 elif is_pydata_spmatrix(A):
19 return max(abs(A).sum(axis=1))
20 else:
21 return np.linalg.norm(A, np.inf)
24def _exact_1_norm(A):
25 # A compatibility function which should eventually disappear.
26 if scipy.sparse.isspmatrix(A):
27 return max(abs(A).sum(axis=0).flat)
28 elif is_pydata_spmatrix(A):
29 return max(abs(A).sum(axis=0))
30 else:
31 return np.linalg.norm(A, 1)
34def _trace(A):
35 # A compatibility function which should eventually disappear.
36 if scipy.sparse.isspmatrix(A):
37 return A.diagonal().sum()
38 elif is_pydata_spmatrix(A):
39 return A.to_scipy_sparse().diagonal().sum()
40 else:
41 return np.trace(A)
44def _ident_like(A):
45 # A compatibility function which should eventually disappear.
46 if scipy.sparse.isspmatrix(A):
47 return scipy.sparse.construct.eye(A.shape[0], A.shape[1],
48 dtype=A.dtype, format=A.format)
49 elif is_pydata_spmatrix(A):
50 import sparse
51 return sparse.eye(A.shape[0], A.shape[1], dtype=A.dtype)
52 else:
53 return np.eye(A.shape[0], A.shape[1], dtype=A.dtype)
56def expm_multiply(A, B, start=None, stop=None, num=None, endpoint=None):
57 """
58 Compute the action of the matrix exponential of A on B.
60 Parameters
61 ----------
62 A : transposable linear operator
63 The operator whose exponential is of interest.
64 B : ndarray
65 The matrix or vector to be multiplied by the matrix exponential of A.
66 start : scalar, optional
67 The starting time point of the sequence.
68 stop : scalar, optional
69 The end time point of the sequence, unless `endpoint` is set to False.
70 In that case, the sequence consists of all but the last of ``num + 1``
71 evenly spaced time points, so that `stop` is excluded.
72 Note that the step size changes when `endpoint` is False.
73 num : int, optional
74 Number of time points to use.
75 endpoint : bool, optional
76 If True, `stop` is the last time point. Otherwise, it is not included.
78 Returns
79 -------
80 expm_A_B : ndarray
81 The result of the action :math:`e^{t_k A} B`.
83 Notes
84 -----
85 The optional arguments defining the sequence of evenly spaced time points
86 are compatible with the arguments of `numpy.linspace`.
88 The output ndarray shape is somewhat complicated so I explain it here.
89 The ndim of the output could be either 1, 2, or 3.
90 It would be 1 if you are computing the expm action on a single vector
91 at a single time point.
92 It would be 2 if you are computing the expm action on a vector
93 at multiple time points, or if you are computing the expm action
94 on a matrix at a single time point.
95 It would be 3 if you want the action on a matrix with multiple
96 columns at multiple time points.
97 If multiple time points are requested, expm_A_B[0] will always
98 be the action of the expm at the first time point,
99 regardless of whether the action is on a vector or a matrix.
101 References
102 ----------
103 .. [1] Awad H. Al-Mohy and Nicholas J. Higham (2011)
104 "Computing the Action of the Matrix Exponential,
105 with an Application to Exponential Integrators."
106 SIAM Journal on Scientific Computing,
107 33 (2). pp. 488-511. ISSN 1064-8275
108 http://eprints.ma.man.ac.uk/1591/
110 .. [2] Nicholas J. Higham and Awad H. Al-Mohy (2010)
111 "Computing Matrix Functions."
112 Acta Numerica,
113 19. 159-208. ISSN 0962-4929
114 http://eprints.ma.man.ac.uk/1451/
116 Examples
117 --------
118 >>> from scipy.sparse import csc_matrix
119 >>> from scipy.sparse.linalg import expm, expm_multiply
120 >>> A = csc_matrix([[1, 0], [0, 1]])
121 >>> A.todense()
122 matrix([[1, 0],
123 [0, 1]], dtype=int64)
124 >>> B = np.array([np.exp(-1.), np.exp(-2.)])
125 >>> B
126 array([ 0.36787944, 0.13533528])
127 >>> expm_multiply(A, B, start=1, stop=2, num=3, endpoint=True)
128 array([[ 1. , 0.36787944],
129 [ 1.64872127, 0.60653066],
130 [ 2.71828183, 1. ]])
131 >>> expm(A).dot(B) # Verify 1st timestep
132 array([ 1. , 0.36787944])
133 >>> expm(1.5*A).dot(B) # Verify 2nd timestep
134 array([ 1.64872127, 0.60653066])
135 >>> expm(2*A).dot(B) # Verify 3rd timestep
136 array([ 2.71828183, 1. ])
137 """
138 if all(arg is None for arg in (start, stop, num, endpoint)):
139 X = _expm_multiply_simple(A, B)
140 else:
141 X, status = _expm_multiply_interval(A, B, start, stop, num, endpoint)
142 return X
145def _expm_multiply_simple(A, B, t=1.0, balance=False):
146 """
147 Compute the action of the matrix exponential at a single time point.
149 Parameters
150 ----------
151 A : transposable linear operator
152 The operator whose exponential is of interest.
153 B : ndarray
154 The matrix to be multiplied by the matrix exponential of A.
155 t : float
156 A time point.
157 balance : bool
158 Indicates whether or not to apply balancing.
160 Returns
161 -------
162 F : ndarray
163 :math:`e^{t A} B`
165 Notes
166 -----
167 This is algorithm (3.2) in Al-Mohy and Higham (2011).
169 """
170 if balance:
171 raise NotImplementedError
172 if len(A.shape) != 2 or A.shape[0] != A.shape[1]:
173 raise ValueError('expected A to be like a square matrix')
174 if A.shape[1] != B.shape[0]:
175 raise ValueError('the matrices A and B have incompatible shapes')
176 ident = _ident_like(A)
177 n = A.shape[0]
178 if len(B.shape) == 1:
179 n0 = 1
180 elif len(B.shape) == 2:
181 n0 = B.shape[1]
182 else:
183 raise ValueError('expected B to be like a matrix or a vector')
184 u_d = 2**-53
185 tol = u_d
186 mu = _trace(A) / float(n)
187 A = A - mu * ident
188 A_1_norm = _exact_1_norm(A)
189 if t*A_1_norm == 0:
190 m_star, s = 0, 1
191 else:
192 ell = 2
193 norm_info = LazyOperatorNormInfo(t*A, A_1_norm=t*A_1_norm, ell=ell)
194 m_star, s = _fragment_3_1(norm_info, n0, tol, ell=ell)
195 return _expm_multiply_simple_core(A, B, t, mu, m_star, s, tol, balance)
198def _expm_multiply_simple_core(A, B, t, mu, m_star, s, tol=None, balance=False):
199 """
200 A helper function.
201 """
202 if balance:
203 raise NotImplementedError
204 if tol is None:
205 u_d = 2 ** -53
206 tol = u_d
207 F = B
208 eta = np.exp(t*mu / float(s))
209 for i in range(s):
210 c1 = _exact_inf_norm(B)
211 for j in range(m_star):
212 coeff = t / float(s*(j+1))
213 B = coeff * A.dot(B)
214 c2 = _exact_inf_norm(B)
215 F = F + B
216 if c1 + c2 <= tol * _exact_inf_norm(F):
217 break
218 c1 = c2
219 F = eta * F
220 B = F
221 return F
224# This table helps to compute bounds.
225# They seem to have been difficult to calculate, involving symbolic
226# manipulation of equations, followed by numerical root finding.
227_theta = {
228 # The first 30 values are from table A.3 of Computing Matrix Functions.
229 1: 2.29e-16,
230 2: 2.58e-8,
231 3: 1.39e-5,
232 4: 3.40e-4,
233 5: 2.40e-3,
234 6: 9.07e-3,
235 7: 2.38e-2,
236 8: 5.00e-2,
237 9: 8.96e-2,
238 10: 1.44e-1,
239 # 11
240 11: 2.14e-1,
241 12: 3.00e-1,
242 13: 4.00e-1,
243 14: 5.14e-1,
244 15: 6.41e-1,
245 16: 7.81e-1,
246 17: 9.31e-1,
247 18: 1.09,
248 19: 1.26,
249 20: 1.44,
250 # 21
251 21: 1.62,
252 22: 1.82,
253 23: 2.01,
254 24: 2.22,
255 25: 2.43,
256 26: 2.64,
257 27: 2.86,
258 28: 3.08,
259 29: 3.31,
260 30: 3.54,
261 # The rest are from table 3.1 of
262 # Computing the Action of the Matrix Exponential.
263 35: 4.7,
264 40: 6.0,
265 45: 7.2,
266 50: 8.5,
267 55: 9.9,
268 }
271def _onenormest_matrix_power(A, p,
272 t=2, itmax=5, compute_v=False, compute_w=False):
273 """
274 Efficiently estimate the 1-norm of A^p.
276 Parameters
277 ----------
278 A : ndarray
279 Matrix whose 1-norm of a power is to be computed.
280 p : int
281 Non-negative integer power.
282 t : int, optional
283 A positive parameter controlling the tradeoff between
284 accuracy versus time and memory usage.
285 Larger values take longer and use more memory
286 but give more accurate output.
287 itmax : int, optional
288 Use at most this many iterations.
289 compute_v : bool, optional
290 Request a norm-maximizing linear operator input vector if True.
291 compute_w : bool, optional
292 Request a norm-maximizing linear operator output vector if True.
294 Returns
295 -------
296 est : float
297 An underestimate of the 1-norm of the sparse matrix.
298 v : ndarray, optional
299 The vector such that ||Av||_1 == est*||v||_1.
300 It can be thought of as an input to the linear operator
301 that gives an output with particularly large norm.
302 w : ndarray, optional
303 The vector Av which has relatively large 1-norm.
304 It can be thought of as an output of the linear operator
305 that is relatively large in norm compared to the input.
307 """
308 #XXX Eventually turn this into an API function in the _onenormest module,
309 #XXX and remove its underscore,
310 #XXX but wait until expm_multiply goes into scipy.
311 return scipy.sparse.linalg.onenormest(aslinearoperator(A) ** p)
313class LazyOperatorNormInfo:
314 """
315 Information about an operator is lazily computed.
317 The information includes the exact 1-norm of the operator,
318 in addition to estimates of 1-norms of powers of the operator.
319 This uses the notation of Computing the Action (2011).
320 This class is specialized enough to probably not be of general interest
321 outside of this module.
323 """
324 def __init__(self, A, A_1_norm=None, ell=2, scale=1):
325 """
326 Provide the operator and some norm-related information.
328 Parameters
329 ----------
330 A : linear operator
331 The operator of interest.
332 A_1_norm : float, optional
333 The exact 1-norm of A.
334 ell : int, optional
335 A technical parameter controlling norm estimation quality.
336 scale : int, optional
337 If specified, return the norms of scale*A instead of A.
339 """
340 self._A = A
341 self._A_1_norm = A_1_norm
342 self._ell = ell
343 self._d = {}
344 self._scale = scale
346 def set_scale(self,scale):
347 """
348 Set the scale parameter.
349 """
350 self._scale = scale
352 def onenorm(self):
353 """
354 Compute the exact 1-norm.
355 """
356 if self._A_1_norm is None:
357 self._A_1_norm = _exact_1_norm(self._A)
358 return self._scale*self._A_1_norm
360 def d(self, p):
361 """
362 Lazily estimate d_p(A) ~= || A^p ||^(1/p) where ||.|| is the 1-norm.
363 """
364 if p not in self._d:
365 est = _onenormest_matrix_power(self._A, p, self._ell)
366 self._d[p] = est ** (1.0 / p)
367 return self._scale*self._d[p]
369 def alpha(self, p):
370 """
371 Lazily compute max(d(p), d(p+1)).
372 """
373 return max(self.d(p), self.d(p+1))
375def _compute_cost_div_m(m, p, norm_info):
376 """
377 A helper function for computing bounds.
379 This is equation (3.10).
380 It measures cost in terms of the number of required matrix products.
382 Parameters
383 ----------
384 m : int
385 A valid key of _theta.
386 p : int
387 A matrix power.
388 norm_info : LazyOperatorNormInfo
389 Information about 1-norms of related operators.
391 Returns
392 -------
393 cost_div_m : int
394 Required number of matrix products divided by m.
396 """
397 return int(np.ceil(norm_info.alpha(p) / _theta[m]))
400def _compute_p_max(m_max):
401 """
402 Compute the largest positive integer p such that p*(p-1) <= m_max + 1.
404 Do this in a slightly dumb way, but safe and not too slow.
406 Parameters
407 ----------
408 m_max : int
409 A count related to bounds.
411 """
412 sqrt_m_max = np.sqrt(m_max)
413 p_low = int(np.floor(sqrt_m_max))
414 p_high = int(np.ceil(sqrt_m_max + 1))
415 return max(p for p in range(p_low, p_high+1) if p*(p-1) <= m_max + 1)
418def _fragment_3_1(norm_info, n0, tol, m_max=55, ell=2):
419 """
420 A helper function for the _expm_multiply_* functions.
422 Parameters
423 ----------
424 norm_info : LazyOperatorNormInfo
425 Information about norms of certain linear operators of interest.
426 n0 : int
427 Number of columns in the _expm_multiply_* B matrix.
428 tol : float
429 Expected to be
430 :math:`2^{-24}` for single precision or
431 :math:`2^{-53}` for double precision.
432 m_max : int
433 A value related to a bound.
434 ell : int
435 The number of columns used in the 1-norm approximation.
436 This is usually taken to be small, maybe between 1 and 5.
438 Returns
439 -------
440 best_m : int
441 Related to bounds for error control.
442 best_s : int
443 Amount of scaling.
445 Notes
446 -----
447 This is code fragment (3.1) in Al-Mohy and Higham (2011).
448 The discussion of default values for m_max and ell
449 is given between the definitions of equation (3.11)
450 and the definition of equation (3.12).
452 """
453 if ell < 1:
454 raise ValueError('expected ell to be a positive integer')
455 best_m = None
456 best_s = None
457 if _condition_3_13(norm_info.onenorm(), n0, m_max, ell):
458 for m, theta in _theta.items():
459 s = int(np.ceil(norm_info.onenorm() / theta))
460 if best_m is None or m * s < best_m * best_s:
461 best_m = m
462 best_s = s
463 else:
464 # Equation (3.11).
465 for p in range(2, _compute_p_max(m_max) + 1):
466 for m in range(p*(p-1)-1, m_max+1):
467 if m in _theta:
468 s = _compute_cost_div_m(m, p, norm_info)
469 if best_m is None or m * s < best_m * best_s:
470 best_m = m
471 best_s = s
472 best_s = max(best_s, 1)
473 return best_m, best_s
476def _condition_3_13(A_1_norm, n0, m_max, ell):
477 """
478 A helper function for the _expm_multiply_* functions.
480 Parameters
481 ----------
482 A_1_norm : float
483 The precomputed 1-norm of A.
484 n0 : int
485 Number of columns in the _expm_multiply_* B matrix.
486 m_max : int
487 A value related to a bound.
488 ell : int
489 The number of columns used in the 1-norm approximation.
490 This is usually taken to be small, maybe between 1 and 5.
492 Returns
493 -------
494 value : bool
495 Indicates whether or not the condition has been met.
497 Notes
498 -----
499 This is condition (3.13) in Al-Mohy and Higham (2011).
501 """
503 # This is the rhs of equation (3.12).
504 p_max = _compute_p_max(m_max)
505 a = 2 * ell * p_max * (p_max + 3)
507 # Evaluate the condition (3.13).
508 b = _theta[m_max] / float(n0 * m_max)
509 return A_1_norm <= a * b
512def _expm_multiply_interval(A, B, start=None, stop=None,
513 num=None, endpoint=None, balance=False, status_only=False):
514 """
515 Compute the action of the matrix exponential at multiple time points.
517 Parameters
518 ----------
519 A : transposable linear operator
520 The operator whose exponential is of interest.
521 B : ndarray
522 The matrix to be multiplied by the matrix exponential of A.
523 start : scalar, optional
524 The starting time point of the sequence.
525 stop : scalar, optional
526 The end time point of the sequence, unless `endpoint` is set to False.
527 In that case, the sequence consists of all but the last of ``num + 1``
528 evenly spaced time points, so that `stop` is excluded.
529 Note that the step size changes when `endpoint` is False.
530 num : int, optional
531 Number of time points to use.
532 endpoint : bool, optional
533 If True, `stop` is the last time point. Otherwise, it is not included.
534 balance : bool
535 Indicates whether or not to apply balancing.
536 status_only : bool
537 A flag that is set to True for some debugging and testing operations.
539 Returns
540 -------
541 F : ndarray
542 :math:`e^{t_k A} B`
543 status : int
544 An integer status for testing and debugging.
546 Notes
547 -----
548 This is algorithm (5.2) in Al-Mohy and Higham (2011).
550 There seems to be a typo, where line 15 of the algorithm should be
551 moved to line 6.5 (between lines 6 and 7).
553 """
554 if balance:
555 raise NotImplementedError
556 if len(A.shape) != 2 or A.shape[0] != A.shape[1]:
557 raise ValueError('expected A to be like a square matrix')
558 if A.shape[1] != B.shape[0]:
559 raise ValueError('the matrices A and B have incompatible shapes')
560 ident = _ident_like(A)
561 n = A.shape[0]
562 if len(B.shape) == 1:
563 n0 = 1
564 elif len(B.shape) == 2:
565 n0 = B.shape[1]
566 else:
567 raise ValueError('expected B to be like a matrix or a vector')
568 u_d = 2**-53
569 tol = u_d
570 mu = _trace(A) / float(n)
572 # Get the linspace samples, attempting to preserve the linspace defaults.
573 linspace_kwargs = {'retstep': True}
574 if num is not None:
575 linspace_kwargs['num'] = num
576 if endpoint is not None:
577 linspace_kwargs['endpoint'] = endpoint
578 samples, step = np.linspace(start, stop, **linspace_kwargs)
580 # Convert the linspace output to the notation used by the publication.
581 nsamples = len(samples)
582 if nsamples < 2:
583 raise ValueError('at least two time points are required')
584 q = nsamples - 1
585 h = step
586 t_0 = samples[0]
587 t_q = samples[q]
589 # Define the output ndarray.
590 # Use an ndim=3 shape, such that the last two indices
591 # are the ones that may be involved in level 3 BLAS operations.
592 X_shape = (nsamples,) + B.shape
593 X = np.empty(X_shape, dtype=np.result_type(A.dtype, B.dtype, float))
594 t = t_q - t_0
595 A = A - mu * ident
596 A_1_norm = _exact_1_norm(A)
597 ell = 2
598 norm_info = LazyOperatorNormInfo(t*A, A_1_norm=t*A_1_norm, ell=ell)
599 if t*A_1_norm == 0:
600 m_star, s = 0, 1
601 else:
602 m_star, s = _fragment_3_1(norm_info, n0, tol, ell=ell)
604 # Compute the expm action up to the initial time point.
605 X[0] = _expm_multiply_simple_core(A, B, t_0, mu, m_star, s)
607 # Compute the expm action at the rest of the time points.
608 if q <= s:
609 if status_only:
610 return 0
611 else:
612 return _expm_multiply_interval_core_0(A, X,
613 h, mu, q, norm_info, tol, ell,n0)
614 elif not (q % s):
615 if status_only:
616 return 1
617 else:
618 return _expm_multiply_interval_core_1(A, X,
619 h, mu, m_star, s, q, tol)
620 elif (q % s):
621 if status_only:
622 return 2
623 else:
624 return _expm_multiply_interval_core_2(A, X,
625 h, mu, m_star, s, q, tol)
626 else:
627 raise Exception('internal error')
630def _expm_multiply_interval_core_0(A, X, h, mu, q, norm_info, tol, ell, n0):
631 """
632 A helper function, for the case q <= s.
633 """
635 # Compute the new values of m_star and s which should be applied
636 # over intervals of size t/q
637 if norm_info.onenorm() == 0:
638 m_star, s = 0, 1
639 else:
640 norm_info.set_scale(1./q)
641 m_star, s = _fragment_3_1(norm_info, n0, tol, ell=ell)
642 norm_info.set_scale(1)
644 for k in range(q):
645 X[k+1] = _expm_multiply_simple_core(A, X[k], h, mu, m_star, s)
646 return X, 0
649def _expm_multiply_interval_core_1(A, X, h, mu, m_star, s, q, tol):
650 """
651 A helper function, for the case q > s and q % s == 0.
652 """
653 d = q // s
654 input_shape = X.shape[1:]
655 K_shape = (m_star + 1, ) + input_shape
656 K = np.empty(K_shape, dtype=X.dtype)
657 for i in range(s):
658 Z = X[i*d]
659 K[0] = Z
660 high_p = 0
661 for k in range(1, d+1):
662 F = K[0]
663 c1 = _exact_inf_norm(F)
664 for p in range(1, m_star+1):
665 if p > high_p:
666 K[p] = h * A.dot(K[p-1]) / float(p)
667 coeff = float(pow(k, p))
668 F = F + coeff * K[p]
669 inf_norm_K_p_1 = _exact_inf_norm(K[p])
670 c2 = coeff * inf_norm_K_p_1
671 if c1 + c2 <= tol * _exact_inf_norm(F):
672 break
673 c1 = c2
674 X[k + i*d] = np.exp(k*h*mu) * F
675 return X, 1
678def _expm_multiply_interval_core_2(A, X, h, mu, m_star, s, q, tol):
679 """
680 A helper function, for the case q > s and q % s > 0.
681 """
682 d = q // s
683 j = q // d
684 r = q - d * j
685 input_shape = X.shape[1:]
686 K_shape = (m_star + 1, ) + input_shape
687 K = np.empty(K_shape, dtype=X.dtype)
688 for i in range(j + 1):
689 Z = X[i*d]
690 K[0] = Z
691 high_p = 0
692 if i < j:
693 effective_d = d
694 else:
695 effective_d = r
696 for k in range(1, effective_d+1):
697 F = K[0]
698 c1 = _exact_inf_norm(F)
699 for p in range(1, m_star+1):
700 if p == high_p + 1:
701 K[p] = h * A.dot(K[p-1]) / float(p)
702 high_p = p
703 coeff = float(pow(k, p))
704 F = F + coeff * K[p]
705 inf_norm_K_p_1 = _exact_inf_norm(K[p])
706 c2 = coeff * inf_norm_K_p_1
707 if c1 + c2 <= tol * _exact_inf_norm(F):
708 break
709 c1 = c2
710 X[k + i*d] = np.exp(k*h*mu) * F
711 return X, 2