papyri
{ } Raw JSON

bundles / scipy 1.17.1 / scipy / sparse / linalg / _funm_multiply_krylov / funm_multiply_krylov

function

scipy.sparse.linalg._funm_multiply_krylov:funm_multiply_krylov

source: /scipy/sparse/linalg/_funm_multiply_krylov.py :149

Signature

def   funm_multiply_krylov ( f ,  A ,  b ,  * ,  assume_a = general ,  t = 1.0 ,  atol = 0.0 ,  rtol = 1e-06 ,  restart_every_m = None ,  max_restarts = 20 )

Summary

A restarted Krylov method for evaluating y = f(tA) b from [1] [2].

Parameters

f : callable

Callable object that computes the matrix function F = f(X).

A : {sparse array, ndarray, LinearOperator}

A real or complex N-by-N matrix. Alternatively, A can be a linear operator which can produce Ax using, e.g., scipy.sparse.linalg.LinearOperator.

b : ndarray

A vector to multiply the f(tA) with.

assume_a : string, optional

Indicate the structure of A. The algorithm will use this information to select the appropriated code path. The available options are 'hermitian'/'her' and 'general'/'gen'. If ommited, then it is assumed that A has a 'general' structure.

t : float, optional

The value to scale the matrix A with. The default is t = 1.0

atol, rtol : float, optional

Parameters for the convergence test. For convergence, norm(||y_k - y_k-1||) <= max(rtol*norm(b), atol) should be satisfied. The default is atol=0. and rtol=1e-6.

restart_every_m : integer

If the iteration number reaches this value a restart is triggered. Larger values increase iteration cost but may be necessary for convergence. If omitted, min(20, n) is used.

max_restarts : int, optional

Maximum number of restart cycles. The algorithm will stop after max_restarts cycles even if the specified tolerance has not been achieved. The default is max_restarts=20

Returns

y : ndarray

The result of f(tA) b.

Notes

The convergence of the Krylov method heavily depends on the spectrum of A and the function f. With restarting, there are only formal proofs for functions of order 1 (e.g., exp, sin, cos) and Stieltjes functions [2] [3], while the general case remains an open problem.

Examples

import numpy as np
from scipy.sparse import csr_array
from scipy.sparse.linalg import funm_multiply_krylov
from scipy.linalg import expm, solve
A = csr_array([[3, 2, 0], [1, -1, 0], [0, 5, 1]], dtype=float)
b = np.array([2, 4, -1], dtype=float)
t = 0.1
✓
Compute ``y = exp(tA) b``. 
y = funm_multiply_krylov(expm, A, b, t = t)
✓
y
✗
ref = expm(t * A.todense()) @ b
err = y - ref
✓
err
✗
Compute :math:`y = (A^3 - A) b`. 
poly = lambda X : X @ X @ X - X
y = funm_multiply_krylov(poly, A, b)
✓
y
✗
ref = poly(A.todense()) @ b
err = y - ref
✓
err
✗
Compute :math:`y = f(tA) b`, where :math:`f(X) = X^{-1}(e^{X} - I)`. This is known as the "phi function" from the exponential integrator literature. 
phim_1 = lambda X : solve(X, expm(X) - np.eye(X.shape[0]))
y = funm_multiply_krylov(phim_1, A, b, t = t)
✓
y
✗
ref = phim_1(t * A.todense()) @ b
err = y - ref
✓
err
✗

Aliases

  • scipy.sparse.linalg.funm_multiply_krylov

Referenced by

This package