Skip to content

Commit

Permalink
chore: fix formatting
Browse files Browse the repository at this point in the history
  • Loading branch information
ashutosh-b-b committed Jul 2, 2024
1 parent e533217 commit a3c2c36
Show file tree
Hide file tree
Showing 31 changed files with 510 additions and 380 deletions.
1 change: 0 additions & 1 deletion LICENSE.md
Original file line number Diff line number Diff line change
Expand Up @@ -19,4 +19,3 @@ The HighDimPDE.jl package is licensed under the MIT "Expat" License:
> LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
> OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
> SOFTWARE.
>
15 changes: 10 additions & 5 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@
[![codecov](https://codecov.io/gh/SciML/HighDimPDE.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/SciML/HighDimPDE.jl)
[![Build Status](https://github.com/SciML/HighDimPDE.jl/workflows/CI/badge.svg)](https://github.com/SciML/HighDimPDE.jl/actions?query=workflow%3ACI)

[![ColPrac: Contributor's Guide on Collaborative Practices for Community Packages](https://img.shields.io/badge/ColPrac-Contributor's%20Guide-blueviolet)](https://github.com/SciML/ColPrac)
[![ColPrac: Contributor's Guide on Collaborative Practices for Community Packages](https://img.shields.io/badge/ColPrac-Contributor%27s%20Guide-blueviolet)](https://github.com/SciML/ColPrac)
[![SciML Code Style](https://img.shields.io/static/v1?label=code%20style&message=SciML&color=9558b2&labelColor=389826)](https://github.com/SciML/SciMLStyle)

# HighDimPDE.jl
Expand All @@ -13,12 +13,13 @@

$$
\begin{aligned}
(\partial_t u)(t,x) &= \int_{\Omega} f\big(t,x,{\bf x}, u(t,x),u(t,{\bf x}), ( \nabla_x u )(t,x ),( \nabla_x u )(t,{\bf x} ) \big) d{\bf x} \\
& \quad + \big\langle \mu(t,x), ( \nabla_x u )( t,x ) \big\rangle + \tfrac{1}{2} \text{Trace} \big(\sigma(t,x) [ \sigma(t,x) ]^* ( \text{Hess}_x u)(t, x ) \big).
(\partial_t u)(t,x) &= \int_{\Omega} f\big(t,x,{\bf x}, u(t,x),u(t,{\bf x}), ( \nabla_x u )(t,x ),( \nabla_x u )(t,{\bf x} ) \big) d{\bf x} \\
& \quad + \big\langle \mu(t,x), ( \nabla_x u )( t,x ) \big\rangle + \tfrac{1}{2} \text{Trace} \big(\sigma(t,x) [ \sigma(t,x) ]^* ( \text{Hess}_x u)(t, x ) \big).
\end{aligned}
$$

where $u \colon [0,T] \times \Omega \to \mathbb{R}, \Omega \subseteq \mathbb{R}^{d}$ is subject to initial and boundary conditions, and where $d$ is large.

## Tutorials and Documentation

For information on using the package,
Expand All @@ -27,22 +28,26 @@ For information on using the package,
the documentation, which contains the unreleased features.

## Installation

Open Julia and type the following

```julia
using Pkg;
Pkg.add("HighDimPDE.jl")
```

This will download the latest version from the git repo and download all dependencies.

## Getting started

See documentation and `test` folders.

## Reference
- Boussange, V., Becker, S., Jentzen, A., Kuckuck, B., Pellissier, L., Deep learning approximations for non-local nonlinear PDEs with Neumann boundary conditions. [arXiv](https://arxiv.org/abs/2205.03672) (2022)

- Boussange, V., Becker, S., Jentzen, A., Kuckuck, B., Pellissier, L., Deep learning approximations for non-local nonlinear PDEs with Neumann boundary conditions. [arXiv](https://arxiv.org/abs/2205.03672) (2022)

<!-- - Becker, S., Braunwarth, R., Hutzenthaler, M., Jentzen, A., von Wurstemberger, P., Numerical simulations for full history recursive multilevel Picard approximations for systems of high-dimensional partial differential equations. [arXiv](https://arxiv.org/abs/2005.10206) (2020)
- Beck, C., Becker, S., Cheridito, P., Jentzen, A., Neufeld, A., Deep splitting method for parabolic PDEs. [arXiv](https://arxiv.org/abs/1907.03452) (2019)
- Han, J., Jentzen, A., E, W., Solving high-dimensional partial differential equations using deep learning. [arXiv](https://arxiv.org/abs/1707.02568) (2018) -->

<!-- ## Acknowledgements
`HighDimPDE.jl` is inspired from Sebastian Becker's scripts in Python, TensorFlow, and C++. Pr. Arnulf Jentzen largely contributed to the theoretical developments of the solver algorithms implemented. -->
11 changes: 5 additions & 6 deletions docs/pages.jl
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@ pages = [
"Solver Algorithms" => ["MLP.md",
"DeepSplitting.md",
"DeepBSDE.md",
"NNStopping.md",
"NNStopping.md",
"NNKolmogorov.md",
"NNParamKolmogorov.md"],
"Tutorials" => [
Expand All @@ -14,14 +14,13 @@ pages = [
"tutorials/mlp.md",
"tutorials/nnstopping.md",
"tutorials/nnkolmogorov.md",
"tutorials/nnparamkolmogorov.md",
"tutorials/nnparamkolmogorov.md"
],
"Extended Examples" => [
"examples/blackscholes.md",
"examples/blackscholes.md"
],
"Solver Algorithms" =>
["MLP.md",
"Solver Algorithms" => ["MLP.md",
"DeepSplitting.md",
"DeepBSDE.md"],
"Feynman Kac formula" => "Feynman_Kac.md",
"Feynman Kac formula" => "Feynman_Kac.md"
]
9 changes: 6 additions & 3 deletions docs/src/DeepBSDE.md
Original file line number Diff line number Diff line change
@@ -1,16 +1,19 @@
# [The `DeepBSDE` algorithm](@id deepbsde)

### Problems Supported:
1. [`ParabolicPDEProblem`](@ref)

1. [`ParabolicPDEProblem`](@ref)

```@autodocs
Modules = [HighDimPDE]
Pages = ["DeepBSDE.jl", "DeepBSDE_Han.jl"]
```

## The general idea 💡
The `DeepBSDE` algorithm is similar in essence to the `DeepSplitting` algorithm, with the difference that

The `DeepBSDE` algorithm is similar in essence to the `DeepSplitting` algorithm, with the difference that
it uses two neural networks to approximate both the the solution and its gradient.

## References
- Han, J., Jentzen, A., E, W., Solving high-dimensional partial differential equations using deep learning. [arXiv](https://arxiv.org/abs/1707.02568) (2018)

- Han, J., Jentzen, A., E, W., Solving high-dimensional partial differential equations using deep learning. [arXiv](https://arxiv.org/abs/1707.02568) (2018)
63 changes: 40 additions & 23 deletions docs/src/DeepSplitting.md
Original file line number Diff line number Diff line change
@@ -1,8 +1,9 @@
# [The `DeepSplitting` algorithm](@id deepsplitting)

### Problems Supported:
1. [`PIDEProblem`](@ref)
2. [`ParabolicPDEProblem`](@ref)

1. [`PIDEProblem`](@ref)
2. [`ParabolicPDEProblem`](@ref)

```@autodocs
Modules = [HighDimPDE]
Expand All @@ -13,76 +14,93 @@ The `DeepSplitting` algorithm reformulates the PDE as a stochastic learning prob

The algorithm relies on two main ideas:

- The approximation of the solution $u$ by a parametric function $\bf u^\theta$. This function is generally chosen as a (Feedforward) Neural Network, as it is a [universal approximator](https://en.wikipedia.org/wiki/Universal_approximation_theorem).
- The approximation of the solution $u$ by a parametric function $\bf u^\theta$. This function is generally chosen as a (Feedforward) Neural Network, as it is a [universal approximator](https://en.wikipedia.org/wiki/Universal_approximation_theorem).

- The training of $\bf u^\theta$ by simulated stochastic trajectories of particles, through the link between linear PDEs and the expected trajectory of associated Stochastic Differential Equations (SDEs), explicitly stated by the [Feynman Kac formula](https://en.wikipedia.org/wiki/Feynman–Kac_formula).
- The training of $\bf u^\theta$ by simulated stochastic trajectories of particles, through the link between linear PDEs and the expected trajectory of associated Stochastic Differential Equations (SDEs), explicitly stated by the [Feynman Kac formula](https://en.wikipedia.org/wiki/Feynman%E2%80%93Kac_formula).

## The general idea 💡

Consider the PDE

```math
\partial_t u(t,x) = \mu(t, x) \nabla_x u(t,x) + \frac{1}{2} \sigma^2(t, x) \Delta_x u(t,x) + f(x, u(t,x)) \tag{1}
```
with initial conditions $u(0, x) = g(x)$, where $u \colon \R^d \to \R$.

with initial conditions $u(0, x) = g(x)$, where $u \colon \R^d \to \R$.

### Local Feynman Kac formula

`DeepSplitting` solves the PDE iteratively over small time intervals by using an approximate [Feynman-Kac representation](@ref feynmankac) locally.

More specifically, considering a small time step $dt = t_{n+1} - t_n$ one has that

```math
u(t_{n+1}, X_{T - t_{n+1}}) \approx \mathbb{E} \left[ f(t, X_{T - t_{n}}, u(t_{n},X_{T - t_{n}}))(t_{n+1} - t_n) + u(t_{n}, X_{T - t_{n}}) | X_{T - t_{n+1}}\right] \tag{3}.
```

One can therefore use Monte Carlo integrations to approximate the expectations

```math
u(t_{n+1}, X_{T - t_{n+1}}) \approx \frac{1}{\text{batch\_size}}\sum_{j=1}^{\text{batch\_size}} \left[ u(t_{n}, X_{T - t_{n}}^{(j)}) + (t_{n+1} - t_n)\sum_{k=1}^{K} \big[ f(t_n, X_{T - t_{n}}^{(j)}, u(t_{n},X_{T - t_{n}}^{(j)})) \big] \right]
```


### Reformulation as a learning problem

The `DeepSplitting` algorithm approximates $u(t_{n+1}, x)$ by a parametric function ${\bf u}^\theta_n(x)$. It is advised to let this function be a neural network ${\bf u}_\theta \equiv NN_\theta$ as they are universal approximators.

For each time step $t_n$, the `DeepSplitting` algorithm
For each time step $t_n$, the `DeepSplitting` algorithm

1. Generates the particle trajectories $X^{x, (j)}$ satisfying [Eq. (2)](@ref feynmankac) over the whole interval $[0,T]$.
1. Generates the particle trajectories $X^{x, (j)}$ satisfying [Eq. (2)](@ref feynmankac) over the whole interval $[0,T]$.

2. Seeks ${\bf u}_{n+1}^{\theta}$ by minimizing the loss function
2. Seeks ${\bf u}_{n+1}^{\theta}$ by minimizing the loss function

```math
L(\theta) = ||{\bf u}^\theta_{n+1}(X_{T - t_n}) - \left[ f(t, X_{T - t_{n-1}}, {\bf u}_{n-1}(X_{T - t_{n-1}}))(t_{n} - t_{n-1}) + {\bf u}_{n-1}(X_{T - t_{n-1}}) \right] ||
```

This way, the PDE approximation problem is decomposed into a sequence of separate learning problems.
In `HighDimPDE.jl` the right parameter combination $\theta$ is found by iteratively minimizing $L$ using **stochastic gradient descent**.

!!! tip

To solve with `DeepSplitting`, one needs to provide to `solve`
- `dt`
- `batch_size`
- `maxiters`: the number of iterations for minimizing the loss function
- `abstol`: the absolute tolerance for the loss function
- `use_cuda`: if you have a Nvidia GPU, recommended.

- `dt`
- `batch_size`
- `maxiters`: the number of iterations for minimizing the loss function
- `abstol`: the absolute tolerance for the loss function
- `use_cuda`: if you have a Nvidia GPU, recommended.

## Solving point-wise or on a hypercube

### Pointwise

`DeepSplitting` allows obtaining $u(t,x)$ on a single point $x \in \Omega$ with the keyword $x$.

```julia
prob = PIDEProblem(μ, σ, x, tspan, g, f,)
prob = PIDEProblem(μ, σ, x, tspan, g, f)
```

### Hypercube

Yet more generally, one wants to solve Eq. (1) on a $d$-dimensional cube $[a,b]^d$. This is offered by `HighDimPDE.jl` with the keyword `x0_sample`.

```julia
prob = PIDEProblem(μ, σ, x, tspan, g, f; x0_sample = x0_sample)
```

Internally, this is handled by assigning a random variable as the initial point of the particles, i.e.

```math
X_t^\xi = \int_0^t \mu(X_s^x)ds + \int_0^t\sigma(X_s^x)dB_s + \xi,
```

where $\xi$ a random variable uniformly distributed over $[a,b]^d$. This way, the neural network is trained on the whole interval $[a,b]^d$ instead of a single point.

## Non-local PDEs

`DeepSplitting` can solve for non-local reaction diffusion equations of the type

```math
\partial_t u = \mu(x) \nabla_x u + \frac{1}{2} \sigma^2(x) \Delta u + \int_{\Omega}f(x,y, u(t,x), u(t,y))dy
```
Expand All @@ -93,13 +111,12 @@ The non-localness is handled by a Monte Carlo integration.
u(t_{n+1}, X_{T - t_{n+1}}) \approx \sum_{j=1}^{\text{batch\_size}} \left[ u(t_{n}, X_{T - t_{n}}^{(j)}) + \frac{(t_{n+1} - t_n)}{K}\sum_{k=1}^{K} \big[ f(t, X_{T - t_{n}}^{(j)}, Y_{X_{T - t_{n}}^{(j)}}^{(k)}, u(t_{n},X_{T - t_{n}}^{(j)}), u(t_{n},Y_{X_{T - t_{n}}^{(j)}}^{(k)})) \big] \right]
```

!!! tip
In practice, if you have a non-local model, you need to provide the sampling method and the number $K$ of MC integration through the keywords `mc_sample` and `K`.
```julia
alg = DeepSplitting(nn, opt = opt, mc_sample = mc_sample, K = 1)
```
`mc_sample` can be whether `UniformSampling(a, b)` or ` NormalSampling(σ_sampling, shifted)`.
!!! tip
In practice, if you have a non-local model, you need to provide the sampling method and the number $K$ of MC integration through the keywords `mc_sample` and `K`.
`julia alg = DeepSplitting(nn, opt = opt, mc_sample = mc_sample, K = 1) `
`mc_sample` can be whether `UniformSampling(a, b)` or ` NormalSampling(σ_sampling, shifted)`.

## References
- Boussange, V., Becker, S., Jentzen, A., Kuckuck, B., Pellissier, L., Deep learning approximations for non-local nonlinear PDEs with Neumann boundary conditions. [arXiv](https://arxiv.org/abs/2205.03672) (2022)
- Beck, C., Becker, S., Cheridito, P., Jentzen, A., Neufeld, A., Deep splitting method for parabolic PDEs. [arXiv](https://arxiv.org/abs/1907.03452) (2019)

- Boussange, V., Becker, S., Jentzen, A., Kuckuck, B., Pellissier, L., Deep learning approximations for non-local nonlinear PDEs with Neumann boundary conditions. [arXiv](https://arxiv.org/abs/2205.03672) (2022)
- Beck, C., Becker, S., Cheridito, P., Jentzen, A., Neufeld, A., Deep splitting method for parabolic PDEs. [arXiv](https://arxiv.org/abs/1907.03452) (2019)
43 changes: 31 additions & 12 deletions docs/src/Feynman_Kac.md
Original file line number Diff line number Diff line change
@@ -1,33 +1,39 @@
# [Feynman Kac formula](@id feynmankac)

The Feynman Kac formula is generally stated for terminal condition problems (see e.g. [Wikipedia](https://en.wikipedia.org/wiki/Feynman–Kac_formula)), where
The Feynman Kac formula is generally stated for terminal condition problems (see e.g. [Wikipedia](https://en.wikipedia.org/wiki/Feynman%E2%80%93Kac_formula)), where

```math
\partial_t u(t,x) + \mu(x) \nabla_x u(t,x) + \frac{1}{2} \sigma^2(x) \Delta_x u(t,x) + f(x, u(t,x)) = 0 \tag{1}
```
with terminal condition $u(T, x) = g(x)$, and $u \colon \R^d \to \R$.

with terminal condition $u(T, x) = g(x)$, and $u \colon \R^d \to \R$.

In this case, the FK formula states that for all $t \in (0,T)$ it holds that

```math
u(t, x) = \int_t^T \mathbb{E} \left[ f(X^x_{s-t}, u(s, X^x_{s-t}))ds \right] + \mathbb{E} \left[ u(0, X^x_{T-t}) \right] \tag{2}
```
where

where

```math
X_t^x = \int_0^t \mu(X_s^x)ds + \int_0^t\sigma(X_s^x)dB_s + x,
```

and $B_t$ is a [Brownian motion](https://en.wikipedia.org/wiki/Wiener_process).

![Brownian motion - Wikipedia](https://upload.wikimedia.org/wikipedia/commons/f/f8/Wiener_process_3d.png)

Intuitively, this formula is motivated by the fact that [the density of Brownian particles (motion) satisfies the diffusion equation](https://en.wikipedia.org/wiki/Brownian_motion#Einstein's_theory).

Intuitively, this formula is motivated by the fact that [the density of Brownian particles (motion) satisfies the diffusion equation](https://en.wikipedia.org/wiki/Brownian_motion#Einstein%27s_theory).

The equivalence between the average trajectory of particles and PDEs given by the Feynman-Kac formula allows overcoming the curse of dimensionality that standard numerical methods suffer from, because the expectations can be approximated [Monte Carlo integrations](https://en.wikipedia.org/wiki/Monte_Carlo_integration), which approximation error decreases as $1/\sqrt{N}$ and is therefore not dependent on the dimensions. On the other hand, the computational complexity of traditional deterministic techniques grows exponentially in the number of dimensions.
The equivalence between the average trajectory of particles and PDEs given by the Feynman-Kac formula allows overcoming the curse of dimensionality that standard numerical methods suffer from, because the expectations can be approximated [Monte Carlo integrations](https://en.wikipedia.org/wiki/Monte_Carlo_integration), which approximation error decreases as $1/\sqrt{N}$ and is therefore not dependent on the dimensions. On the other hand, the computational complexity of traditional deterministic techniques grows exponentially in the number of dimensions.

## Forward non-linear Feynman-Kac

> How to transform previous equation to an initial value problem?
Define $v(\tau, x) = u(T-\tau, x)$. Observe that $v(0,x) = u(T,x)$. Further, observe that by the chain rule

```math
\begin{aligned}
\partial_\tau v(\tau, x) &= \partial_\tau u(T-\tau,x)\\
Expand All @@ -36,20 +42,26 @@ Define $v(\tau, x) = u(T-\tau, x)$. Observe that $v(0,x) = u(T,x)$. Further, obs
\end{aligned}
```

From Eq. (1) we get that
From Eq. (1) we get that

```math
- \partial_t u(T - \tau,x) = \mu(x) \nabla_x u(T - \tau,x) + \frac{1}{2} \sigma^2(x) \Delta_x u(T - \tau,x) + f(x, u(T - \tau,x)).
```

Replacing $u(T-\tau, x)$ by $v(\tau, x)$ we get that $v$ satisfies

```math
\partial_\tau v(\tau, x) = \mu(x) \nabla_x v(\tau,x) + \frac{1}{2} \sigma^2(x) \Delta_x v(\tau,x) + f(x, v(\tau,x))
```

and from Eq. (2) we obtain

```math
v(\tau, x) = \int_{T-\tau}^T \mathbb{E} \left[ f(X^x_{s- T + \tau}, v(s, X^x_{s-T + \tau}))ds \right] + \mathbb{E} \left[ v(0, X^x_{\tau}) \right].
```
By using the substitution rule with $\tau \to \tau -T$ (shifting by T) and $\tau \to - \tau$ (inversing), and finally inversing the integral bound we get that

By using the substitution rule with $\tau \to \tau -T$ (shifting by T) and $\tau \to - \tau$ (inversing), and finally inversing the integral bound we get that

```math
\begin{aligned}
v(\tau, x) &= \int_{-\tau}^0 \mathbb{E} \left[ f(X^x_{s + \tau}, v(s + T, X^x_{s + \tau}))ds \right] + \mathbb{E} \left[ v(0, X^x_{\tau}) \right]\\
Expand All @@ -58,18 +70,25 @@ v(\tau, x) &= \int_{-\tau}^0 \mathbb{E} \left[ f(X^x_{s + \tau}, v(s + T, X^x_{s
\end{aligned}
```

This leads to the
This leads to the

!!! info "Non-linear Feynman Kac for initial value problems"

Consider the PDE

```math
\partial_t u(t,x) = \mu(t, x) \nabla_x u(t,x) + \frac{1}{2} \sigma^2(t, x) \Delta_x u(t,x) + f(x, u(t,x))
```
with initial conditions $u(0, x) = g(x)$, where $u \colon \R^d \to \R$.

with initial conditions $u(0, x) = g(x)$, where $u \colon \R^d \to \R$.
Then

```math
u(t, x) = \int_0^t \mathbb{E} \left[ f(X^x_{t - s}, u(T-s, X^x_{t - s}))ds \right] + \mathbb{E} \left[ u(0, X^x_t) \right] \tag{3}
```
with

with

```math
X_t^x = \int_0^t \mu(X_s^x)ds + \int_0^t\sigma(X_s^x)dB_s + x.
```
```
Loading

0 comments on commit a3c2c36

Please sign in to comment.