montecarlo_nystrom package¶
- montecarlo_nystrom.montecarlo_nystrom(*, random_samples: Callable[[int], Array], kernel: Callable[[Array, Array], Array], rhs: Callable[[Array], Array], n: int, n_mean: int = 1, solve: Callable[[Array, Array], Array] | None = None) Callable[[Array], Array][source]¶
Solve integral equations of the second kind of the following form.
\(\forall d \in \mathbb{N}.\) \(\forall \Omega \in \mathbb{R}^d [\Omega \text{ is bounded Lipschitz}].\) \(\forall p \in L^\infty(\Omega, {\mathbb{R}}_{\geq 0}) [\int_\Omega p(y) dy = 1].\) \(\forall f, z \in L^2(\Omega,\mathbb{C}).\) \(\forall k \in L^2(\Omega,L^2(\Omega,\mathbb{C}))\) \([z(y) + \int_\Omega k(y, y') z(y') p(y') dy' = f(y)].\)
Let \(N \in \mathbb{N}\), \((y_i)_{i=1}^N\) be i.i.d. samples drawn from \(p\).
Let \((z_{N,i})_{i=1}^N\) be the solution of the linear system
\[ z_{N,i} + \frac{1}{N} \sum_{j=1}^N k(y_i, y_j) z_{N,j} = f(y_i) \quad i \in \{1, \ldots, N\} \]and
\[ z_N(y) := f(y) - \frac{1}{N} \sum_{i=1}^N k(y, y_i) z_{N,i} \quad y \in \Omega \]Then \(z_N\) would approximate \(z\) as \(N \to \infty\).
- Parameters:
random_samples (Callable[[int], Array]) – A function that takes an integer n and returns n i.i.d. samples from the distribution p of shape (…, n, d).
kernel (Callable[[Array, Array], Array]) –
Kernel function k of (…, :, :, d), (…, :, :, d) -> (…, :, :).
The diagonal part will be automatically dropped and can be of any value.
rhs (Callable[[Array], Array]) – Right-hand side function f of (…, :, d) -> (…, :).
n (int) – The number of random samples to draw from the distribution p.
n_mean (int) – The number of independent runs to average over, by default 1.
solve (Callable[[Array, Array], Array] | None, optional) – A function that takes a matrix A of shape (…, n, n) and a right-hand side b of shape (…, n, 1) and returns the solution of shape (…, n, 1). If None, uses xp.linalg.solve, by default None.
- Returns:
Approximate solution function z_N that takes an array of shape (…(x), …, d) and returns an array of shape (…(x), …).
- Return type:
Callable[[Array], Array]
Examples
>>> import numpy as np >>> from montecarlo_nystrom import montecarlo_nystrom >>> rng = np.random.default_rng(0) >>> def random_samples(n): ... return rng.uniform(0, 1, size=(n, 1)) >>> def kernel(x, y): ... return np.linalg.vector_norm(x - y, axis=-1) ** -0.4 >>> def rhs(x): ... x0 = x[..., 0] ... return np.ones_like(x0) >>> z_N = montecarlo_nystrom( ... random_samples=random_samples, ... kernel=kernel, ... rhs=rhs, ... n=100, ... n_mean=10, ... ) >>> np.round(z_N(np.asarray((0.5,))), 6) # Evaluate at x=0.5 np.float64(0.272957)