Modeling Linear PDE Systems with Probabilistic Machine Learning
Markus Lange-Hegermann
Institut für industrielle Informationstechnik - inIT
Department of Electrical Engineering and Computer Science, TH OWL
Bogdan Raiță
Department of Mathematics and Statistics
Georgetown University
With many, many contributions by Andreas Besginow, Marc Härkönen, Jianlei Huang, Xin Li, and Daniel Robertz
Differential Equations and Data: Why?
- Combine expert knowledge and data.
- Want intepretable and reliable models.
- Focus on (irregularly spaced) data rather than IC/BC.
- Numerical simulations might be computationally too costly.
Example
Damped spring-mass system
\[\frac{\mathrm{d}^2y}{\mathrm{d}t^2} + 2 \frac{\mathrm{d}y}{\mathrm{d}t} + 10y = 0\] Sample 5 noisy points
Example
Neural Network
- ⚠️ ~300,000 parameters
- ✅ Fast convergence
- ❌ Bad inter/extrapolation
Example
Physics Informed Neural Network
- ✅ Good interpolation
- ✅ Fast inference
- ⚠️ ~300,000 parameters
- ⚠️ Requires extra points
- ⚠️ Not an exact solution
- ❌ Slow convergence
- ❌ Hard to optimize
A more constrained approach
How to solve? \[\frac{\mathrm{d}^2y}{\mathrm{d}t^2} + 2 \frac{\mathrm{d}y}{\mathrm{d}t} + 10y = 0\]
Algebraic preprocessing: characteristic frequencies \[z^2 + 2z + 10 = 0 \Rightarrow z = -1 \pm 3\sqrt{-1}\]
Solution Space \[y(t) = (\textcolor{#b51963}{c_1} \cdot e^{-t} \cos 3t + \textcolor{#b51963}{c_2} \cdot e^{-t} \sin 3t)\]
A more constrained approach
How to solve? \[\frac{\mathrm{d}^2y}{\mathrm{d}t^2} + 2 \frac{\mathrm{d}y}{\mathrm{d}t} + 10y = 0\]
Algebraic preprocessing: characteristic frequencies \[z^2 + 2z + 10 = 0 \Rightarrow z = -1 \pm 3\sqrt{-1}\]
Solution Space \[y(t) = (\textcolor{#b51963}{c_1} \cdot e^{-t} \cos 3t + \textcolor{#b51963}{c_2} \cdot e^{-t} \sin 3t)\]
Gaussian process prior \[Y(t) \sim (\textcolor{#b51963}{C_1} \cdot e^{-t} \cos 3t + \textcolor{#b51963}{C_2} \cdot e^{-t} \sin 3t) + \textcolor{#b51963}{\epsilon}\] where \[ \textcolor{#b51963}{C_1 \sim \mathcal{N}(0, \sigma_1)}, \textcolor{#b51963}{C_2 \sim \mathcal{N}(0, \sigma_2)}, \textcolor{#b51963}{\epsilon \sim \mathcal{N}(0, \sigma_0)} \]
A more constrained approach
How to solve? \[\frac{\mathrm{d}^2y}{\mathrm{d}t^2} + 2 \frac{\mathrm{d}y}{\mathrm{d}t} + 10y = 0\]
Algebraic preprocessing: characteristic frequencies \[z^2 + 2z + 10 = 0 \Rightarrow z = -1 \pm 3\sqrt{-1}\]
Solution Space \[y(t) = (\textcolor{#b51963}{c_1} \cdot e^{-t} \cos 3t + \textcolor{#b51963}{c_2} \cdot e^{-t} \sin 3t)\]
Gaussian process prior \[Y(t) \sim (\textcolor{#b51963}{C_1} \cdot e^{-t} \cos 3t + \textcolor{#b51963}{C_2} \cdot e^{-t} \sin 3t) + \textcolor{#b51963}{\epsilon}\] where \[ \textcolor{#b51963}{C_1 \sim \mathcal{N}(0, \sigma_1)}, \textcolor{#b51963}{C_2 \sim \mathcal{N}(0, \sigma_2)}, \textcolor{#b51963}{\epsilon \sim \mathcal{N}(0, \sigma_0)} \]
Algebra: determine suitable frequencies
Stochastic: weigh frequencies
Gaussian Process
(or ridge regression)
- ✅ 3 parameters
- ✅ Fast convergence
- ✅ Good inter/extrapolation
- ✅ Exact solution
- ✅ No extra points
Overview
| EPGP (ours) | S-EPGP(ours) | Vanilla PINN | Vanilla Num. Solver | |
|---|---|---|---|---|
| Differential equations | linear, c.c. | linear, c.c. | any | well-posed |
| Additional Information | data | data | data | IC/BC |
| Evaluation time | $\approx$ data | $\approx$ sparsity | neural net | mediocre |
| Solutions? | yes | yes | approx. near data | yes |
| Prerequisites | integral | - | - | - |
| Sampling | yes | yes | no | not applicable |
Another ODE example
How to solve? \[\frac{\mathrm{d}^4y}{\mathrm{d}t^4} + 2 \frac{\mathrm{d}^2y}{\mathrm{d}t^2} + y = 0\]
Another ODE example
How to solve? \[\frac{\mathrm{d}^4y}{\mathrm{d}t^4} + 2 \frac{\mathrm{d}^2y}{\mathrm{d}t^2} + y = 0\]
Algebraic preprocessing: characteristic frequencies \[z^4 + 2z^2 + 1 = 0 \Rightarrow z = \pm i\text{ with multiplicity 2}\]
Another ODE example
How to solve? \[\frac{\mathrm{d}^4y}{\mathrm{d}t^4} + 2 \frac{\mathrm{d}^2y}{\mathrm{d}t^2} + y = 0\]
Algebraic preprocessing: characteristic frequencies \[z^4 + 2z^2 + 1 = 0 \Rightarrow z = \pm i\text{ with multiplicity 2}\]
Solution Space \[y(t) = \textcolor{#b51963}{c_1}\textcolor{#0073e6}{t}e^{\textcolor{#0073e6}{-i}t} + \textcolor{#b51963}{c_2}e^{\textcolor{#0073e6}{-i}t}+\textcolor{#b51963}{c_3}\textcolor{#0073e6}{t}e^{\textcolor{#0073e6}{i}t}+\textcolor{#b51963}{c_4}e^{\textcolor{#0073e6}{i}t}\]
Another ODE example
How to solve? \[\frac{\mathrm{d}^4y}{\mathrm{d}t^4} + 2 \frac{\mathrm{d}^2y}{\mathrm{d}t^2} + y = 0\]
Algebraic preprocessing: characteristic frequencies \[z^4 + 2z^2 + 1 = 0 \Rightarrow z = \pm i\text{ with multiplicity 2}\]
Solution Space \[y(t) = \textcolor{#b51963}{c_1}\textcolor{#0073e6}{t}e^{\textcolor{#0073e6}{-i}t} + \textcolor{#b51963}{c_2}e^{\textcolor{#0073e6}{-i}t}+\textcolor{#b51963}{c_3}\textcolor{#0073e6}{t}e^{\textcolor{#0073e6}{i}t}+\textcolor{#b51963}{c_4}e^{\textcolor{#0073e6}{i}t}\]
Gaussian process prior \[Y(t) \sim \textcolor{#b51963}{C_1}\textcolor{#0073e6}{t}e^{\textcolor{#0073e6}{-i}t} + \textcolor{#b51963}{C_2}e^{\textcolor{#0073e6}{-i}t}+\textcolor{#b51963}{C_3}\textcolor{#0073e6}{t}e^{\textcolor{#0073e6}{i}t}+\textcolor{#b51963}{C_4}e^{\textcolor{#0073e6}{i}t} + \textcolor{#b51963}{\epsilon}\] where \[ \textcolor{#b51963}{C_j \sim \mathcal{N}(0, \sigma_j)}, \textcolor{#b51963}{\epsilon \sim \mathcal{N}(0, \sigma_0)} \]
Principle for ODEs
How to solve? \[\sum_{j=0}^k \alpha_j\frac{\mathrm d^jy}{\mathrm d t^j} = 0\]
Principle for ODEs
How to solve? \[\sum_{j=0}^k \alpha_j\frac{\mathrm d^jy}{\mathrm d t^j} = 0\]
Algebraic preprocessing: characteristic frequencies \[\sum_{j=0}^k \alpha_jz^j=0 \Rightarrow z\in V = \text{ complex roots with multiplicities}\]
Principle for ODEs
How to solve? \[\sum_{j=0}^k \alpha_j\frac{\mathrm d^jy}{\mathrm d t^j} = 0\]
Algebraic preprocessing: characteristic frequencies \[\sum_{j=0}^k \alpha_jz^k=0 \Rightarrow z\in V = \text{ complex roots with multiplicities}\]
Solution Space \[y(t) = \sum_{\textcolor{#0073e6}{z}\in V}\sum_{\text{multipliers for z}}\textcolor{#b51963}{c_{h}}\textcolor{#0073e6}{B_h(t)}e^{\textcolor{#0073e6}{z}t} \]
Principle for ODEs
How to solve? \[\sum_{j=0}^k \alpha_j\frac{\mathrm d^jy}{\mathrm d t^j} = 0\]
Algebraic preprocessing: characteristic frequencies \[\sum_{j=0}^k \alpha_jz^k=0 \Rightarrow z\in V = \text{ complex roots with multiplicities}\]
Solution Space \[y(t) = \sum_{\textcolor{#0073e6}{z}\in V}\sum_{\text{multipliers for z}}\textcolor{#b51963}{c_{h}}\textcolor{#0073e6}{B_h(t)}e^{\textcolor{#0073e6}{z}t} \]
Principle for ODEs:
All solutions are linear combination of exponential-polynomial solutions.
Fourier & Heat Equation
\[\frac{\partial T}{\partial t}-\frac{\partial^2 T}{\partial x^2}=0\]
Fourier & Heat Equation
\[\frac{\partial T}{\partial t}-\frac{\partial^2 T}{\partial x^2}=0\]
Fourier (1807):
- Noticed that $T_z(x,t)=e^{\textcolor{#0073e6}{-z^2}t+\textcolor{#0073e6}{iz}x}=e^{-z^2t}(\cos zx+i\sin zx)$ are solutions.
- Noticed linearity: $T^c(x,z)=\sum_{z=-N}^N \textcolor{#b51963}{c_z}e^{\textcolor{#0073e6}{-z^2}t+\textcolor{#0073e6}{iz}x}$ are solutions.
- Postulated density: All solutions $T$ can be approximated with solutions $T^c$.
- Solved from initial condition: If $T(x,0)=g(x)$ and $T\approx T^c$ then \[g(x)= \sum_{z\in\mathbb Z} \textcolor{#b51963}{c_z}e^{\textcolor{#0073e6}{iz}x}.\]
- Invented Fourier series.
Our goals
- Use additional knowledge of linear PDE systems with constant coefficients.
- Construct a computationally nice probabily distribution on such PDE systems.
- Restricting to this important special case allows more suitable methods.
- Introduce two such methods: EPGP and S-EPGP.
Principle for PDEs
All solutions are linear combination of exponential-polynomial solutions up to approximation.
Let $A \in R^{\ell' \times {\ell}}$ for $R=\R[\partial_{x_1},\ldots,\partial_{x_n}]$
Want to solve $Af=0$ for $f=(f_1,f_2,\ldots,f_\ell)$.
Principle for PDEs
All solutions are linear combination of exponential-polynomial solutions up to approximation.
Let $A \in R^{\ell' \times {\ell}}$ for $R=\R[\partial_{x_1},\ldots,\partial_{x_n}]$
Want to solve $Af=0$ for $f=(f_1,f_2,\ldots,f_\ell)$.
Algebraic preprocessing: characteristic frequencies \[\ker A(z)=\{0\} \Rightarrow z\in V = \text{ characteristic variety}\]
Principle for PDEs
All solutions are linear combination of exponential-polynomial solutions up to approximation.
Let $A \in R^{\ell' \times {\ell}}$ for $R=\R[\partial_{x_1},\ldots,\partial_{x_n}]$
Want to solve $Af=0$ for $f=(f_1,f_2,\ldots,f_\ell)$.
Algebraic preprocessing: characteristic frequencies \[\ker A(z)=\{0\} \Rightarrow z\in V = \text{ characteristic variety}\]
Solution Space is the closure of \[f(x) = \sum_{j=1}^M\sum_{h=1}^m\textcolor{#b51963}{c_{j}}\textcolor{#0073e6}{B_h(x,z_j)}e^{\textcolor{#0073e6}{z_j}\cdot x} \] where $\textcolor{#0073e6}{z_j}\in V$.
Principle for PDEs
All solutions are linear combination of exponential-polynomial solutions up to approximation.
Let $A \in R^{\ell' \times {\ell}}$ for $R=\R[\partial_{x_1},\ldots,\partial_{x_n}]$
Want to solve $Af=0$ for $f=(f_1,f_2,\ldots,f_\ell)$.
Algebraic preprocessing: characteristic frequencies \[\ker A(z)=\{0\} \Rightarrow z\in V = \text{ characteristic variety}\]
Solution Space is the closure of \[f(x) = \sum_{j=1}^M\sum_{h=1}^m\textcolor{#b51963}{c_{j}}\textcolor{#0073e6}{B_h(x,z_j)}e^{\textcolor{#0073e6}{z_j}\cdot x} \] where $\textcolor{#0073e6}{z_j}\in V$.
This is the FUNDAMENTAL PRINCIPLE of Ehrenpreis-Palamodov ('70).
Fundamental Principle for PDEs
Explain Ehrenpreis-Palamodov in more detail:
Let $A \in R^{\ell' \times {\ell}}$ for $R=\R[\partial_{x_1},\ldots,\partial_{x_n}]$ and let $\Omega \subseteq \R^n$ be a convex, compact set. There exist
- characteristic varieties $\textcolor{#0073e6}{\{V_1,\dotsc,V_s\}}$ and
- $\ell\times1$ polynomial operators $\textcolor{#0073e6}{\{B_{i,1}(x, z), \dotsc, B_{i,m_i}(x,z)\}_{i=1,\dotsc,s}}$ in $2n$ variables,
such that any smooth solution $f \colon \Omega \to \R^{\ell}$ to the equation $Af = 0$ can be approximated by \[ \begin{align*} f(x) &= \sum_{i=1}^s \sum_{j=1}^{m_i}\sum_h \textcolor{#b51963}{c_h^{ij}} \textcolor{#0073e6}{B_{i,j}}(x, \textcolor{#0073e6}{z^{ij}_h}) e^{\langle x, \textcolor{#0073e6}{z^{ij}_h} \rangle} \end{align*} \] where $\textcolor{#0073e6}{z^{ij}_h}\in V_i$. This is a NONLINEAR FOURIER SERIES with frequencies NOT from a lattice.
S-EPGP
Sparse version of EPGP.
- Define a GP prior with realizations of the form \[ \begin{align*} f(x) = \sum_{j=1}^m \sum_{i=1}^r \textcolor{#b51963}{C_{i,j}} \textcolor{#0073e6}{B_j}(x, \textcolor{#0073e6}{z_{i,j}}) e^{\langle x, \textcolor{#0073e6}{z_{i,j}} \rangle} =: \textcolor{#b51963}{C^T}\cdot \textcolor{#0073e6}{\phi(x)}, \end{align*} \]
- Turn $f(x)$ into a GP: $\textcolor{#b51963}{C_{i,j} \sim \mathcal{N}\left(0, \Sigma\right)}$ for $\textcolor{#b51963}{\Sigma=\operatorname{diag}(\sigma_{i}^2)}$
S-EPGP
Sparse version of EPGP. Suitable, when integration is not possible.
- Define a GP prior with realizations of the form \[ \begin{align*} f(x) = \sum_{j=1}^m \sum_{i=1}^r \textcolor{#b51963}{C_{i,j}} \textcolor{#0073e6}{B_j}(x, \textcolor{#0073e6}{z_{i,j}}) e^{\langle x, \textcolor{#0073e6}{z_{i,j}} \rangle} =: \textcolor{#b51963}{C^T}\cdot \textcolor{#0073e6}{\phi(x)}, \end{align*} \]
- Turn $f(x)$ into a GP: $\textcolor{#b51963}{C_{i,j} \sim \mathcal{N}\left(0, \Sigma\right)}$ for $\textcolor{#b51963}{\Sigma=\operatorname{diag}(\sigma_{i}^2)}$
Wave equation 2D: Initial Condition
\[\frac{\partial^2 u}{\partial t^2} = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\]
Wave equation 2D: Initial Speed
\[\frac{\partial^2 u}{\partial t^2} = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\]
Wave equation 2D
\[\frac{\partial^2 u}{\partial t^2} = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\]
Numerical solution:
Mean of our S-EPGP:
PINN:
First three frames of the numerical solution are the training data, with a $21\times21$ spatial grid.
Wave equation 2D
\[\frac{\partial^2 u}{\partial t^2} = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\]
Heat Equation
Heat Equation
Heat Equation
Heat Equation
Heat Equation
Heat Equation
Heat Equation
Cuomo, Di Cola, Giampaolo, Rozza, Raissi, Piccialli; Scientific machine learning through physics-informed neural networks: Where we are and what’s next. Journal of Scientific Computing; arXiv:2201.05624
Heat Equation
Maxwell's equations
\[\mathbf{E}(x,y,z,t)\]
3D wave equation
\[ \mathbf{B}(x,y,z,t) \]
Wave equation 2D: Boundary Conditions
\[\frac{\partial^2 u}{\partial t^2} = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\]
Wave equation 2D: BC, old method
\[\frac{\partial^2 u}{\partial t^2} = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\]
Wave equation 2D: BC, new method
\[\frac{\partial^2 u}{\partial t^2} = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\]
Wave equation 2D
\[\frac{\partial^2 u}{\partial t^2} = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\]
Wave equation 2D
\[\frac{\partial^2 u}{\partial t^2} = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\]
Heat equation 2D: IBVP
\[\frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\]
Conclusion
-
EPGP: A class of GP priors for solving PDE
- S-EPGP: sparse version
- Fully algorithmic
- No approximations anywhere $\Rightarrow$ exact solutions
- Improved accuracy and convergence speed over PINN
- Compatible with standard GP techniques
- Interpretable hyperparameters (optimizable)
- Only assumptions: linear and constant coefficients
- Very little data necessary
- We can say more for parametrizable/controllable systems
\[ \frac{\partial^2 u}{\partial t^2} = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \]
Conclusion
\[ \frac{\partial^2 u}{\partial t^2} = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \]
Open Problems
- Many
- Large scale problems and comparison to numerical solvers
- Extensions to neural networks
- Bringing these methods to industry
- Boundary conditions for S-EPGP
- Inverse problems
- Convergence and density questions
- Topologically interesting domains and monodromy phenomena
Conclusion
Thank you for your attention!
Literature
- arXiv:2411.16663
- arXiv:2212.14319 (ICML 2023)
- arXiv:2208.12515 (NeurIPS 2022)
- arXiv:2205.03185 (CASC 2022)
- arXiv:2002.00818 (AISTATS 2021)
- arXiv:1801.09197 (NeurIPS 2018)
\[ \frac{\partial^2 u}{\partial t^2} = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \]