Skip to content

ultrasphere-dev/montecarlo-nystrom

Folders and files

NameName
Last commit message
Last commit date

Latest commit

 

History

100 Commits
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

Monte-Carlo Nystrom

CI Status Documentation Status Test coverage percentage

uv Ruff pre-commit

PyPI Version Supported Python versions License


Documentation: https://montecarlo-nystrom.readthedocs.io

Source Code: https://github.com/ultrasphere-dev/montecarlo-nystrom


Monte-Carlo Nystrom method in NumPy / PyTorch

Installation

Install this via pip (or your favourite package manager):

pip install montecarlo-nystrom

Usage

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$.

The below example solves the case where $d = 1$, $\Omega = [0, 1]$, $p$ (random_samples) is the uniform distribution on $[0, 1]$, $k(x, y) = |x - y|^{-0.4}$ (kernel), and $f(x) = 1$ (rhs), and evaluates the solution at $x = (0.5,)$.

>>> 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)

References

Contributors ✨

Thanks goes to these wonderful people (emoji key):

ultrasphere-dev
ultrasphere-dev

💻 🤔 📖

This project follows the all-contributors specification. Contributions of any kind welcome!

Credits

Copier

This package was created with Copier and the browniebroke/pypackage-template project template.

About

Monte-Carlo Nystrom method in NumPy / PyTorch

Topics

Resources

License

Code of conduct

Contributing

Stars

Watchers

Forks

Sponsor this project

Packages

 
 
 

Contributors