# Neural operator preconditioning for accelerating the solution of the parametric Helmholtz equations

Yanfei Xiang, Luc Giraud, Paul Mycek, Maksym Shpakovych, Carola Kruse

<figure><img src="image_1.png" /></figure>

To cite this version:

Yanfei Xiang, Luc Giraud, Paul Mycek, Maksym Shpakovych, Carola Kruse. Neural operator preconditioning for accelerating the solution of the parametric Helmholtz equations. 2024 SIAM Conference on Applied Linear Algebra (SIAM LA24), May 2024, Paris, France. pp.2. hal-04581208v4

HAL Id: hal-04581208

https://inria.hal.science/hal-04581208v4

Submitted on 10 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Neural operator preconditioning for accelerating the solution of the parametric Helmholtz equations

Yanfei Xiang, Luc Giraud, Paul Mycek, Maksym Shpakovych, Carola Kruse

► To cite this version:

Yanfei Xiang, Luc Giraud, Paul Mycek, Maksym Shpakovych, Carola Kruse. Neural operator preconditioning for accelerating the solution of the parametric Helmholtz equations. 2024 SIAM Conference on Applied Linear Algebra (SIAM LA24), Laura Grigori, Daniel Kressner, May 2024, Paris, France. pp.2. [hal-04581208v4](https://hal HAL Id:)

Id: hal-04581208

https://inria.hal.science/hal-04581208v4

Submitted on 10 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
<figure><img src="image_3.png" /></figure>

# Neural operator preconditioning for accelerating the solution of the parametric Helmholtz equations

SIAM LA24, Paris, France

May 15, 2024

Yanfei Xiang (yanfei.xiang@inria.fr)

Concace team, Centre Inria de l'Université de Bordeaux
<figure><img src="image_2.png" /></figure>

# Neural operator preconditioning for accelerating the solution of the parametric Helmholtz equations

SIAM LA24, Paris, France

May 15, 2024

Yanfei Xiang (yanfei.xiang@inria.fr)

Concace team, Centre Inria de l'Université de Bordeaux
# Parametric Helmholtz equations

## Helmholtz equation:

Parametric PDEs:

$$ 
\left\{
\begin{array}{l}
[\nabla^2 + (-\frac{\omega}{c(x)})^2] u(x) = \rho(x) \\
\text{(subject to the Sommerfeld radiation condition at infinity)} \\
+ \text{ Perfectly Matched Layer (PML) BCs on domain } \Omega
\end{array}
\right.
\qquad (1)
$$

* let $x \in \mathbb{R}^d$ be the grid points in the $d$-dimensional $\Omega$ ($d = 1, 2, 3$);

* $u(x) : \mathbb{R}^d \to \mathbb{C}$ is the complex acoustic wavefield to be computed;

* $\rho(x) : \mathbb{R}^d \to \mathbb{C}$ is the source function (could be fixed or not);

* $c(x) : \mathbb{R}^d \to \mathbb{R}_+$ is a speed of sound distribution function (could be fixed or not) and $\omega = 1 \in \mathbb{R}_+$ is an angular frequency of the source.

## Part of the Physical Applications:

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>
# Parametric Helmholtz equations

## Helmholtz equation:

Parametric PDEs:
$$ 
\left\{
\begin{array}{l}
[\nabla^2 + (-\frac{\omega}{c(x)})^2]u(x) = \rho(x) \\
\text{(subject to the Sommerfeld radiation condition at infinity)} \\
+ \text{ Perfectly Matched Layer (PML) BCs on domain } \Omega
\end{array}
\right.
\qquad (1)
$$

* let $x \in \mathbb{R}^d$ be the grid points in the $d$-dimensional $\Omega$ ($d = 1, 2, 3$);

* $u(x) : \mathbb{R}^d \to \mathbb{C}$ is the complex acoustic wavefield to be computed;

* $\rho(x) : \mathbb{R}^d \to \mathbb{C}$ is the source function (could be fixed or not);

* $c(x) : \mathbb{R}^d \to \mathbb{R}_+$ is a speed of sound distribution function (could be fixed or not) and $\omega = 1 \in \mathbb{R}_+$ is an angular frequency of the source.

## Part of the Physical Applications:

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>
Learning parametric Helmholtz operators

Transfer the differential expression to the discrete form:

$$
\text{Parametric PDEs:} \quad \left\{
\begin{array}{l}
\displaystyle \left[ \nabla^2 + \left(-\frac{\omega}{c(x)}\right)^2 \right] u(x) = \rho(x) \\[1.5ex]
\displaystyle \quad + \text{PML BCs on domain } \Omega \\[1.5ex]
\displaystyle \quad \qquad \text{Discrete derivatives } \Downarrow \text{ FEM, FDM, FFT, ...} \\[1.5ex]
\displaystyle \quad A(c)u = b, \quad \text{with LinearOp } A(c) \in \mathbb{C}^{n \times n}, u, b \in \mathbb{C}^n
\end{array}
\right.
\qquad (2)
$$

After discretization, numerical linear algebra methods, like subspace
methods, can be used to solve the parametric Eq. (2). However,

- simply using subspace methods without preconditioning is much less effective for solving the Helmholtz Eqs.ª;

- generate a properly algebraic preconditioner could be possible (if the n is not too big) but is as challenging as solve the system directly;

- generally, the algebraic preconditioning needs to be re-generated for each of the parametric Helmholtz Eqs. (2).

<figure><img src="image_8.png" /></figure>

<footnote>ªErnst and Gander. Why it is Difficult to Solve Helmholtz Problems with Classical Iterative Methods. 2012</footnote>
Learning parametric Helmholtz operators

Transfer the differential expression to the discrete form:

$$
\text{Parametric PDEs:} \quad \left\{
\begin{array}{l}
[\nabla^2 + \left(-\frac{\omega}{c(x)}\right)^2] u(x) = \rho(x) \\
\quad + \text{PML BCs on domain } \Omega \\
\quad \text{Discrete derivatives } \Downarrow \text{ FEM, FDM, FFT, ...} \\
A(c)u = b, \text{ with LinearOp } A(c) \in \mathbb{C}^{n \times n}, u, b \in \mathbb{C}^n
\end{array}
\right.
\qquad (2)
$$

After discretization, numerical linear algebra methods, like subspace
methods, can be used to solve the parametric Eq. (2). However,

- simply using subspace methods without preconditioning is much less effective for solving the Helmholtz Eqs.ª;

- generate a properly algebraic preconditioner could be possible (if the n is not too big) but is as challenging as solve the system directly;

- generally, the algebraic preconditioning needs to be re-generated for each of the parametric Helmholtz Eqs. (2).

<figure><img src="image_7.png" /></figure>

<footnote>ªErnst and Gander. Why it is Difficult to Solve Helmholtz Problems with Classical Iterative Methods. 2012</footnote>
Learning parametric Helmholtz operators

Transfer the differential expression to the discrete form:

$$ \text{Parametric PDEs:} \quad \left\{ \begin{array}{l} \left[ \nabla^2 + \left( -\frac{\omega}{c(x)} \right)^2 \right] u(x) = \rho(x) \\ \quad + \text{PML BCs on domain } \Omega \\ \text{Discrete derivatives } \Downarrow \text{ FEM, FDM, FFT, } \dots \\ A(c)u = b, \text{ with LinearOp } A(c) \in \mathbb{C}^{n \times n}, u, b \in \mathbb{C}^n \end{array} \right. \qquad (2) $$

In recent decade, the thriving neural networks (NNs) solvers, like the physics-informed neural networks (PINNs)ª, is used to solve the parametric Eq. (2) without discretization. However, these NNs solvers

- are usually costly in trainingᵇ, and may fail to solve challenging PDEs if without finely tuning of the hyper-parameters of the NNs;

- solve the PDEs without the theoretical convergence guarantee;

- generally reach limited accuracy and exhibit limited or NO network generalizability, thus re-training is required even it is costly.

<footnote>ªLu et al., Physics-Informed Neural Networks with Hard Constraints for Inverse Design. SISC. 2021</footnote>

<footnote>ᵇStrubell et al., Energy and policy considerations for deep learning in NLP. ACL meeting, Italy. 2019</footnote>
Learning parametric Helmholtz operators

Transfer the differential expression to the discrete form:

$$
\text{Parametric PDEs:} \quad \left\{
\begin{array}{l}
\displaystyle \left[ \nabla^2 + \left( -\frac{\omega}{c(x)} \right)^2 \right] u(x) = \rho(x) \\[1.5ex]
\displaystyle \quad + \text{PML BCs on domain } \Omega \\[1.5ex]
\displaystyle \quad \quad \text{Discrete derivatives } \Downarrow \text{ FEM, FDM, FFT, ...} \\[1.5ex]
\displaystyle \quad A(c)u = b, \text{ with LinearOp } A(c) \in \mathbb{C}^{n \times n}, u, b \in \mathbb{C}^n
\end{array}
\right.
\qquad (2)
$$

In recent decade, the thrived **neural networks (NNs)** solvers, like the
**physics-informed neural networks (PINNs)**ª, is used to solve the
parametric Eq. (2) without discretization. However, these NNs solvers

- are usually **costly in training**ᵇ, and may **fail to solve** challenging PDEs if without finely tuning of the hyper-parameters of the NNs;

- solve the PDEs **without the theoretical convergence guarantee**;

- generally reach **limited accuracy** and exhibit **limited or NO network generalizability**, thus **re-training** is required even it is costly.

<footnote>ªLu et al., Physics-Informed Neural Networks with Hard Constraints for Inverse Design. SISC. 2021</footnote>

<footnote>Inria-bStrubell et al., Energy and policy considerations for deep learning in NLP. ACL meeting, Italy. 2019</footnote>
Learning parametric Helmholtz operators

Transfer the differential expression to the discrete form:

$$
\text{Parametric PDEs:} \quad \left\{
\begin{array}{l}
\displaystyle \left[ \nabla^2 + \left(-\frac{\omega}{c(x)}\right)^2 \right] u(x) = \rho(x) \\[1.5ex]
\displaystyle \quad + \text{PML BCs on domain } \Omega \\[1.5ex]
\displaystyle \quad \text{Discrete derivatives } \Downarrow \text{ FEM, FDM, FFT, ...} \\[1.5ex]
\displaystyle \quad A(c)u = b, \text{ with LinearOp } A(c) \in \mathbb{C}^{n \times n}, u, b \in \mathbb{C}^n
\end{array}
\right.
\qquad (2)
$$

* Goal of this work $\Rightarrow$ To learn neural operator $\mathcal{F}_{\theta}$ that approximates

$$
A(c)^{-1} \text{ by } \mathcal{F}_{\theta}([b, c, (\text{BCs})]) \rightarrow u_{\theta} \sim A(c)^{-1}b
$$

The learned $\mathcal{F}_{\theta}$ can be used as a **flexible preconditioner** for the subspace methods (like FGMRES) to accelerate the solution of parametric Eq. (2) with varying **c**, varying **b** and varying domain $\Omega$ but **without re-training**.

* Loss-function:

$$
\min_{\theta} \frac{\|A(c)u_{\theta} - b\|_{2}^{2}}{\|b\|_{2}^{2}} \quad (3)
$$
Learning parametric Helmholtz operators

Transfer the differential expression to the discrete form:

$$
\text{Parametric PDEs:} \quad \left\{
\begin{array}{l}
\displaystyle \left[ \nabla^2 + \left(-\frac{\omega}{c(x)}\right)^2 \right] u(x) = \rho(x) \\[1.5ex]
\displaystyle \quad + \text{PML BCs on domain } \Omega \\[1.5ex]
\displaystyle \quad \text{Discrete derivatives } \Downarrow \text{ FEM, FDM, FFT, ...} \\[1.5ex]
\displaystyle \quad A(\mathbf{c})u = \mathbf{b}, \text{ with LinearOp } A(\mathbf{c}) \in \mathbb{C}^{n \times n}, \mathbf{u}, \mathbf{b} \in \mathbb{C}^n
\end{array}
\right.
\qquad (2)
$$

* Goal of this work $\Rightarrow$ To learn neural operator $\mathcal{F}_{\theta}$ that approximates

$$
A(c)^{-1} \text{ by } \mathcal{F}_{\theta}([b, c, (BCs)]) \rightarrow u_{\theta} \sim A(c)^{-1}b
$$

The learned $\mathcal{F}_{\theta}$ can be used as a **flexible preconditioner** for the subspace methods (like FGMRES) to accelerate the solution of parametric Eq. (2) with varying $c$, varying $b$ and varying domain $\Omega$ but **without re-training**.

* Loss-function:

$$
\min_{\theta} \frac{\|A(\mathbf{c})u_{\theta} - \mathbf{b}\|_{2}^{2}}{\|\mathbf{b}\|_{2}^{2}} \quad (3)
$$
The generalized Arnoldi relation

Originally, the FGMRES methodª have the generalized Arnoldi relation

$$AZ_m = V_{m+1} \bar{H}_m, \quad (4)$$

where $V_{m+1} = [v_1, ..., v_{m+1}]$ is the Krylov basis, and $Z_m = [z_1, ..., z_m]$ is the preconditioned Krylov basis.

In our FGMRES case, we have

$$z_i = F_{\theta}(v_i), \quad i = 1, \dots, m, \quad (5)$$

where $F_{\theta}$ is the trained non-linear neural operator satisfies $F_{\theta} \sim A^{-1}$. This part is computed with single machine precision (float 32), the one used in the training process, and other parts are in double precision.

⇒ Given there is no information about the data structure of the Krylov basis, the neural operator $F_{\theta}$ is trained with randomly generated datasets.

<footnote>ªY. Saad, A Flexible Inner-Outer Preconditioned GMRES Algorithm. SISC. 1993</footnote>
The generalized Arnoldi relation

Originally, the FGMRES methodª have the generalized Arnoldi relation

$$AZ_m = V_{m+1} \bar{H}_m, \quad (4)$$

where $V_{m+1} = [v_1, ..., v_{m+1}]$ is the Krylov basis, and $Z_m = [z_1, ..., z_m]$ is the preconditioned Krylov basis.

In our FGMRES case, we have

$$z_i = \mathcal{F}_{\theta}(v_i), \quad i = 1, \dots, m, \quad (5)$$

where $\mathcal{F}_{\theta}$ is the trained non-linear neural operator satisfies $\mathcal{F}_{\theta} \sim A^{-1}$. This part is computed with single machine precision (float 32), the one used in the training process, and other parts are in double precision.

⇒ Given there is no information about the data structure of the Krylov basis, the neural operator $\mathcal{F}_{\theta}$ is trained with randomly generated datasets.

<footnote>ªY. Saad, A Flexible Inner-Outer Preconditioned GMRES Algorithm. SISC. 1993</footnote>
U-Net architecture

U-Netª architecture with 4 depth:

<figure><img src="image_9.png" /></figure>

Each blue box corresponds to a multi-channel feature map. The number of channels is denoted on top of the box. The x-y-size is provided at the lower left edge of the box. White boxes represent copied feature maps. The arrows denote the different operations.

<figure><img src="image_10.png" /></figure>

<footnote>ª[Figure 1] Ronneberger et al., U-Net: Convolutional Networks for Biomedical Image Segmentation. 2015</footnote>
U-Netª architecture with 4 depth:

<figure><img src="image_8.png" /></figure>

Each blue box corresponds to a multi-channel feature map. The number of channels is denoted on top of the box. The x-y-size is provided at the lower left edge of the box. White boxes represent copied feature maps. The arrows denote the different operations.

<figure><img src="image_9.png" /></figure>

<footnote>ª[Figure 1] Ronneberger et al., U-Net: Convolutional Networks for Biomedical Image Segmentation. 2015</footnote>
# U-Net with meshes in 2D (domain 64 × 64)

<figure><img src="image_11.png" /></figure>

## Training (on 2 V100 GPUs, depth = 4, train. time: 49.83min):

* Fixed grid point $x \in \mathbb{R}^{n^2}$ in the 2-dimensional domain ($n = 64$ for 2D case)

* Random source term $b \in \mathbb{C}^{n^2} \sim \mathcal{N}(0, 1)$ satisfying normally distr.

* Random speed of sound $c \in \mathbb{R}^{n^2+}$ uniformly distributed on the interval $[1, 2]$ and fixed frequency of the source $\omega = 1 \in \mathbb{R}_+$

* Function $\sigma(x) : \mathbb{R}^{n^2} \to \mathbb{C}^{n^2}$ used in the definition of the PML boundary condition

* U-Net 2d: Training with single machine precision (float 32); Trainable params. - 832 K

## Testing:

* Trained U-Net preconditioner can accelerate the solution of Eq. (2) with varying **b**, **c**, domain size $\Omega$ (bcf. the discretisation invariance property from the convolution property)
# U-Net with meshes in 2D (domain 64 × 64)

<figure><img src="image_10.png" /></figure>

## Training (on 2 V100 GPUs, depth = 4, train. time: 49.83min):

* Fixed grid point $x \in \mathbb{R}^{n^2}$ in the 2-dimensional domain ($n = 64$ for 2D case)

* Random source term $b \in \mathbb{C}^{n^2} \sim N(0, 1)$ satisfying normally distr.

* Random speed of sound $c \in \mathbb{R}^{n^2+}$ uniformly distributed on the interval $[1, 2]$ and fixed frequency of the source $\omega = 1 \in \mathbb{R}_+$

* Function $\sigma(x) : \mathbb{R}^{n^2} \to \mathbb{C}^{n^2}$ used in the definition of the PML boundary condition

* U-Net 2d: Training with single machine precision (float 32); Trainable params. - 832 K

## Testing:

* Trained U-Net preconditioner can accelerate the solution of Eq. (2) with varying **b**, **c**, domain size $\Omega$ (bcf. the discretisation invariance property from the convolution property)
Compare trained U-Net to algebraic preconditioner

Visualization of the first 150 sequential elements of a $64^2 \times 64^2$ matrix:

<figure><img src="image_12.png" /><figcaption>(a). Coefficient matrix A with
$\frac{nnz}{n^2} \times \% = 33.38\%$</figcaption></figure>

A with elements $a_{ij}$ to be placed in matrix $M$.

$$ (1). M(i, j) = \begin{cases} a_{ii}, & \text{if } \|a_{ii}\| \ge \text{thresh for } i = j, \\ a_{ij}, & \text{if } \frac{\|a_{ij}\|}{\|a_{ii}\|\|a_{jj}\|} \ge \text{thresh for } i \ne j. \end{cases} $$

$$ (2). \frac{\|A - M\|_F}{\|A\|_F} = \begin{cases} 5.41 \times 10^{-2} & \text{if thresh} = 10^{-3}, \\ 2.04 \times 10^{-3} & \text{if thresh} = 10^{-4}, \\ 1.81 \times 10^{-16} & \text{if thresh} = 10^{-5}. \end{cases} $$

<figure><img src="image_13.png" /><figcaption>(b). $M$ with thresh= $10^{-3}$,
$\frac{nnz}{n^2} \times \% = 1.16\%$</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(c). $M$ with thresh= $10^{-4}$,
$\frac{nnz}{n^2} \times \% = 3.00\%$</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(d). $M$ with thresh= $10^{-5}$,
$\frac{nnz}{n^2} \times \% = 3.10\%$</figcaption></figure>
## Compare trained U-Net to algebraic preconditioner

Visualization of the first 150 sequential elements of a $64^2 \times 64^2$ matrix:

<figure><img src="image_11.png" /><figcaption>(a). Coefficient matrix A with
$\frac{nnz}{n^2} \times \% = 33.38\%$</figcaption></figure>

A with elements $a_{ij}$ to be placed in matrix $M$.

$$ (1). M(i, j) = \begin{cases} a_{ii}, & \text{if } \|a_{ii}\| \ge \text{thresh for } i = j, \\ a_{ij}, & \text{if } \frac{\|a_{ij}\|}{\|a_{ii}\|\|a_{jj}\|} \ge \text{thresh for } i \ne j. \end{cases} $$

$$ (2). \frac{\|A - M\|_F}{\|A\|_F} = \begin{cases} 5.41 \times 10^{-2} & \text{if thresh} = 10^{-3}, \\ 2.04 \times 10^{-3} & \text{if thresh} = 10^{-4}, \\ 1.81 \times 10^{-16} & \text{if thresh} = 10^{-5}. \end{cases} $$

<figure><img src="image_12.png" /><figcaption>(b). $M$ with thresh$= 10^{-3}$,
$\frac{nnz}{n^2} \times \% = 1.16\%$</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(c). $M$ with thresh$= 10^{-4}$,
$\frac{nnz}{n^2} \times \% = 3.00\%$</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(d). $M$ with thresh$= 10^{-5}$,
$\frac{nnz}{n^2} \times \% = 3.10\%$</figcaption></figure>
## Compare trained U-Net to algebraic preconditioner

Solve $A(c)u = b$ on 2D domain $64 \times 64$ (i.e., $A(c) \in \mathbb{C}^{64^2 \times 64^2}$):

* All involved algorithms are run on the same CPUs/GPUs device with Python prototype.

* $\eta_b = \frac{\|A(c)u-b\|}{\|b\|}$. We stop the iteration when $\eta_b \le \varepsilon$ with $\varepsilon = 10^{-13}$ by default.

<figure><img src="image_16.png" /><figcaption>PGMRES with varying thresh</figcaption></figure>

$$ \frac{\|A - M\|_F}{\|A\|_F} = \begin{cases} 5.41 \times 10^{-2} & \text{if thresh} = 10^{-3}, \\ 2.04 \times 10^{-3} & \text{if thresh} = 10^{-4}, \\ 1.81 \times 10^{-16} & \text{if thresh} = 10^{-5}. \end{cases} $$

* PGMRES: GMRES preconditioned by spilu($M$), which is realized by applying sparse ilu to the sparse matrix $M$ with selected elements from coefficient matrix $A$
## Compare trained U-Net to algebraic preconditioner

Solve $A(c)u = b$ on 2D domain $64 \times 64$ (i.e., $A(c) \in \mathbb{C}^{64^2 \times 64^2}$):

* All involved algorithms are run on the same CPUs/GPUs device with Python prototype.

* $\eta_b = \frac{\|A(c)u-b\|}{\|b\|}$. We stop the iteration when $\eta_b \le \varepsilon$ with $\varepsilon = 10^{-13}$ by default.

<figure><img src="image_15.png" /><figcaption>PGMRES with varying thresh</figcaption></figure>

$$ \frac{\|A - M\|_F}{\|A\|_F} = \begin{cases} 5.41 \times 10^{-2} & \text{if thresh} = 10^{-3}, \\ 2.04 \times 10^{-3} & \text{if thresh} = 10^{-4}, \\ 1.81 \times 10^{-16} & \text{if thresh} = 10^{-5}. \end{cases} $$

* PGMRES: GMRES preconditioned by spilu($M$), which is realized by applying sparse ilu to the sparse matrix $M$ with selected elements from coefficient matrix $A$
## Compare trained U-Net to algebraic preconditioner

Solve $A(c)u = b$ on 2D domain $64 \times 64$ (i.e., $A(c) \in \mathbb{C}^{64^2 \times 64^2}$):

* FGMRES is GMRES preconditioned by the trained U-Net with depth = 4.

* FGMRES is implemented by the mixed-precision calculation.

<figure><img src="image_17.png" /></figure>

* The effectiveness of trained U-Net preconditioner is close to (could be the same as) spilu($M$) with $\text{thresh} = 10^{-5}$ (the best algebra preconditioner for this example)
## Compare trained U-Net to algebraic preconditioner

Solve $A(c)u = b$ on 2D domain $64 \times 64$ (i.e., $A(c) \in \mathbb{C}^{64^2 \times 64^2}$):

* FGMRES is GMRES preconditioned by the trained U-Net with depth = 4.

* FGMRES is implemented by the mixed-precision calculation.

<figure><img src="image_16.png" /></figure>

* The effectiveness of trained U-Net preconditioner is close to (could be the same as) spilu($M$) with $\text{thresh} = 10^{-5}$ (the best algebra preconditioner for this example)
Network generalizability: varying the source function b

Solve $A(c)u = b$ on 2D domain $64 \times 64$ (i.e., $A(c) \in \mathbb{C}^{64^2 \times 64^2}$):

<figure><img src="image_18.png" /></figure>

* Speed of sound $c \in \mathbb{R}^{64^2}$ is fixed and satisfying uniform distribution on interval $[1, 2]$

* Performance of trained U-Net preconditioner is independent from varying the datatype of the $b \in \mathbb{C}^{64^2}$ or varying the mean and median value of the normal setting of b
Network generalizability: varying the source function b

Solve $A(c)u = b$ on 2D domain $64 \times 64$ (i.e., $A(c) \in \mathbb{C}^{64^2 \times 64^2}$):

<figure><img src="image_17.png" /></figure>

* Speed of sound $c \in \mathbb{R}^{64^2}$ is fixed and satisfying uniform distribution on interval $[1, 2]$

* Performance of trained U-Net preconditioner is independent from varying the datatype of the $b \in \mathbb{C}^{64^2}$ or varying the mean and median value of the normal setting of b
Network generalizability: varying the source function b

Visualize wavefiled $u \in C^{64^2}$ solved by 1 dirac for $b \in C^{64^2}$ (with fixed c):

<figure><img src="image_19.png" /></figure>

Visualize wavefiled $u \in C^{64^2}$ solved by 4 dirac for $b \in C^{64^2}$ (with fixed c):

<figure><img src="image_20.png" /></figure>
Network generalizability: varying the source function b

Visualize wavefiled $u \in C^{64^2}$ solved by 1 dirac for $b \in C^{64^2}$ (with fixed c):

<figure><img src="image_18.png" /></figure>

Visualize wavefiled $u \in C^{64^2}$ solved by 4 dirac for $b \in C^{64^2}$ (with fixed c):

<figure><img src="image_19.png" /></figure>
Network generalizability: varying the speed of sound c

Solve $A(c)u = b$ on 2D domain $64 \times 64$ (i.e., $A(c) \in \mathbb{C}^{64^2 \times 64^2}$):

<figure><img src="image_21.png" /><figcaption>(a). shrink the interval</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>(b). enlarge the interval & special structure</figcaption></figure>

■ Source function $b \in \mathbb{C}^{64^2}$ is fixed with 1 dirac setting, 1 in grid position (32, 32), the center of domain $64 \times 64$, and 0 elsewhere in domain

■ Square shape of $c \in \mathbb{R}^{64^2}$: domain $[20:44, 20:44] = [2:1]$ and $[44:20, 44:20] = [2:44, 20:44]$ elsewhere

■ Varying the range of the c satisfying uniform distribution and varying the its shape can effect the performance of the trained U-Net preconditioner
Network generalizability: varying the speed of sound c

Solve $A(c)u = b$ on 2D domain $64 \times 64$ (i.e., $A(c) \in \mathbb{C}^{64^2 \times 64^2}$):

<figure><img src="image_20.png" /><figcaption>(a). shrink the interval</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(b). enlarge the interval & special structure</figcaption></figure>

■ Source function $b \in \mathbb{C}^{64^2}$ is fixed with 1 dirac setting, 1 in grid position (32, 32), the center of domain $64 \times 64$, and 0 elsewhere in domain

■ Square shape of $c \in \mathbb{R}^{64^2}$: domain $[20:44, 20:44] = [2:1]$ and $[44:20, 44:20] = [2:44, 20:44]$ elsewhere

■ Varying the range of the c satisfying uniform distribution and varying the its shape can effect the performance of the trained U-Net preconditioner
## Network generalizability: varying the speed of sound c

Visualize $c \in \mathbb{R}^{64^2}$ with square shape:

<figure><img src="image_23.png" /></figure>

Visualize complex wavefiled $u \in \mathbb{C}^{64^2}$ solved by $c \in \mathbb{R}^{64^2}$ in square shape:

<figure><img src="image_24.png" /></figure>
Network generalizability: varying the speed of sound c

Visualize $c \in \mathbb{R}^{64^2}$ with square shape:

<figure><img src="image_22.png" /></figure>

Visualize complex wavefiled $u \in \mathbb{C}^{64^2}$ solved by $c \in \mathbb{R}^{64^2}$ in square shape:

<figure><img src="image_23.png" /></figure>
Network generalizability: varying the domain size

Solve $A(c)u = b$ from 2D domain $64 \times 64$ until $512 \times 512$ (i.e., from linear operator $A(c) \in \mathbb{C}^{64^2 \times 64^2}$ to $A(c) \in \mathbb{C}^{512^2 \times 512^2}$):

<figure><img src="image_25.png" /><figcaption>Numerical results on different domain size</figcaption></figure>

* U-Net trained on 2D domain $64 \times 64$ (with normal distr. $N(0, 1)$ $b \in \mathbb{C}^{64^2}$ and uniform distr. $c \in \mathbb{R}^{64^2}$ in interval $[1,2]$) exhibits excellent network generalizability from the 4x larger domain $128 \times 128$ until the **64x larger** domain $512 \times 512$
Network generalizability: varying the domain size

Solve $A(c)u = b$ from 2D domain $64 \times 64$ until $512 \times 512$ (i.e., from linear operator $A(c) \in \mathbb{C}^{64^2 \times 64^2}$ to $A(c) \in \mathbb{C}^{512^2 \times 512^2}$):

<figure><img src="image_24.png" /><figcaption>Numerical results on different domain size</figcaption></figure>

- U-Net trained on 2D domain $64 \times 64$ (with normal distr. $N(0, 1)$ $b \in \mathbb{C}^{64^2}$ and uniform distr. $c \in \mathbb{R}^{64^2}$ in interval $[1,2]$) exhibits excellent network generalizability from the 4x larger domain $128 \times 128$ until the **64x larger** domain $512 \times 512$
## Network generalizability: Adult human skull example

Solve $A(c)u = b$ with the **structured dataset** for both $c \in \mathbb{R}^{512^2}$ and $b \in \mathbb{R}^{512^2}$ on a 64x larger 2D domain $512 \times 512$ ($A(c) \in \mathbb{C}^{512^2 \times 512^2}$):

<figure><img src="image_26.png" /></figure>

<figure><img src="image_27.png" /></figure>

This example is downloaded from the qure.ai dataset<sup>a</sup>

Recall part of the training settings:

■ Random training datasets on domain 64 × 64: $b \in \mathbb{C}^{64^2}$ satisfies standard Normal distr. $\mathcal{N}(0, 1)$; $c \in \mathbb{R}^{64^2}$ is Uniform distr. on interval $[1, 2]$

<footnote><sup>a</sup>http://headctstudy.qure.ai/dataset</footnote>
Network generalizability: Adult human skull example

Solve $A(c)u = b$ with the **structured dataset** for both $c \in \mathbb{R}^{512^2}$ and $b \in \mathbb{R}^{512^2}$ on a 64x larger 2D domain $512 \times 512$ ($A(c) \in \mathbb{C}^{512^2 \times 512^2}$):

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /></figure>

This example is downloaded from the **qure.ai** dataset<sup>a</sup>

Recall part of the training settings:

■ Random training datasets on domain 64 × 64: $b \in \mathbb{C}^{64^2}$ satisfies standard Normal distr. $N(0, 1)$; $c \in \mathbb{R}^{64^2}$ is Uniform distr. on interval $[1, 2]$

<footnote><sup>a</sup>http://headctstudy.qure.ai/dataset</footnote>
## Network generalizability: Adult human skull example

Solve $A(c)u = b$ with the **structured dataset** for both c and b on 2D domain $512 \times 512$ (i.e., $A(c) \in \mathbb{C}^{512^2 \times 512^2}$):

### Adult human skull example

<figure><img src="image_28.png" /></figure>

- Comparison of FGMRES (NNs trained on domain $64 \times 64$) to GMRES on solving the skull example (on $512 \times 512$ from CQ500 with 491 scans) to reaching accuracy at $10^{-3}$

- Consumed its & CPU time (Training time of U-Net on 2 V100 GPUs is 49.83min):
  FGMRES = 4097 & 4.24h; GMRES = 14849 & 13.16h
## Network generalizability: Adult human skull example

Solve $A(c)u = b$ with the **structured dataset** for both c and b on 2D domain $512 \times 512$ (i.e., $A(c) \in \mathbb{C}^{512^2 \times 512^2}$):

<figure><img src="image_27.png" /></figure>

- Comparison of FGMRES (NNs trained on domain $64 \times 64$) to GMRES on solving the skull example (on $512 \times 512$ from CQ500 with 491 scans) to reaching accuracy at $10^{-3}$

- Consumed its & CPU time (Training time of U-Net on 2 V100 GPUs is 49.83min):
  FGMRES = 4097 & 4.24h; GMRES = 14849 & 13.16h
# Network generalizability: Adult human skull example

Visualization of the structured dataset of the skull example on 2D domain $512 \times 512$:

<figure><img src="image_29.png" /></figure>

Visualize the complex wavefiled $c \in \mathbb{R}^{512^2}$ of human skull example solved by FGMRES:

<figure><img src="image_30.png" /></figure>
# Network generalizability: Adult human skull example

Visualization of the structured dataset of the skull example on 2D domain $512 \times 512$:

<figure><img src="image_28.png" /></figure>

Visualize the complex wavefiled $c \in \mathbb{R}^{512^2}$ of human skull example solved by FGMRES:

<figure><img src="image_29.png" /></figure>
PINNs & DeepONet & FNO

Other neural networks (NNs) architectures for the parametric Helmholtz equations:

■ NNs solver:

■ PINNs: Lu et al., Physics-Informed Neural Networks with Hard Constraints for Inverse Design. SISC. 2021

■ NNs operators:

■ DeepONet: Lu et al., Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators. 2021;

■ FNO: Kovachki et al., Neural Operator: Learning Maps Between Function Spaces With Applications to PDEs. 2023

⇒ Compare to above NNs architectures, **U-Net** is the optimal one that meets all our goal and exhibits the robust preconditioning behavior

Three random setting for the source function b:

■ Normal distr. $N(0,1)$: satisfying the standard normal distribution $N(0,1)$;

■ Uniform distr. $[-1, 1]:$ satisfying the uniformly distribution on interval $[-1, 1]$;

■ Dirac setting: with one non-zero value in the random position in the 2D domain;
PINNs & DeepONet & FNO

Other neural networks (NNs) architectures for the parametric Helmholtz equations:

■ NNs solver:

■ PINNs: Lu et al., Physics-Informed Neural Networks with Hard Constraints for Inverse Design. SISC. 2021

■ NNs operators:

■ DeepONet: Lu et al., Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators. 2021;

■ FNO: Kovachki et al., Neural Operator: Learning Maps Between Function Spaces With Applications to PDEs. 2023

⇒ Compare to above NNs architectures, **U-Net** is the optimal one that meets all our goal and exhibits the robust preconditioning behavior

Three random setting for the source function b:

■ Normal distr. $N(0,1)$: satisfying the standard normal distribution $N(0,1)$;

■ Uniform distr. $[-1, 1]:$ satisfying the uniformly distribution on interval $[-1, 1]$;

■ Dirac setting: with one non-zero value in the random position in the 2D domain;
Take home message of this work

**Goal:** To learn the optimal neural operator $F_\theta$ for accelerating the solution of the parametric Helmholtz Eqs.

Training different neural networks (NNs) ⇒ choose the optimal NNs architecture → **U-Net is the optimal one for our case**

Training U-Net on 2D with different settings for b ⇒ choose the optimal training dataset → **Standard Normal distr. $N(0,1)$ for b**

⇒ The performance of trained neural operator depends on:

- the choosing of NNs architecture;

- the setting of training datasets;

- the tuning of hyper-parameters of NNs;

- the utilities of training tools like pytorch-lightning.

© Research report version of this work will soon be accessible online at HAL (https://inria.hal.science/) !

→ Yanfei Xiang (yanfei.xiang@inria.fr)
Take home message of this work

**Goal:** To learn the optimal neural operator $F_\theta$ for accelerating the solution of the parametric Helmholtz Eqs.

Training different neural networks (NNs) ⇒ choose the optimal NNs architecture → **U-Net is the optimal one for our case**

Training U-Net on 2D with different settings for b ⇒ choose the optimal training dataset → **Standard Normal distr. N(0,1) for b**

⇒ The performance of trained neural operator depends on:

■ the choosing of NNs architecture;

■ the setting of training datasets;

■ the tuning of hyper-parameters of NNs;

■ the utilities of training tools like pytorch-lightning.

© Research report version of this work will soon be accessible online at HAL (https://inria.hal.science/) !

→ Yanfei Xiang (yanfei.xiang@inria.fr)
# My co-workers from Inria and Cerfacs

- Inria: Luc Giraud and Maksym Shpakovych

- Cerfacs: Paul Mycek and Carola Kruse

<figure><img src="image_31.png" /><figcaption>Luc Giraud</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>Paul Mycek</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>Maksym Shpakovych</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>Carola Kruse</figcaption></figure>

We are **Concace** joint Inria team with Airbus CR&T and Cerfacs.
# My co-workers from Inria and Cerfacs

- Inria: Luc Giraud and Maksym Shpakovych

- Cerfacs: Paul Mycek and Carola Kruse

<figure><img src="image_30.png" /><figcaption>Luc Giraud</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>Paul Mycek</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>Maksym Shpakovych</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>Carola Kruse</figcaption></figure>

We are **Concace** joint Inria team with Airbus CR&T and Cerfacs.
# Concace

## Numerical & Parallel Composability for High Performance Computing

### Joint Inria-Industry Project Team with Airbus Central R&T and Cerfacs

To deal with large size and large dimension problems coming from model- and data-driven applications, Concace will take advantage of modern development tools and languages to design high-level expressions of complex parallel algorithms.

While the traditional approach to HPC is to fully exploit hardware, Concace's complementary approach will enable a richer composability of numerical methods, allowing to fully exploit existing and new numerical algorithms.

<figure><img src="image_35.png" /></figure>

**AIRBUS**

CERFACS

Inria

**Thank you for your attention!**

**Questions ?**
# Concace

## Numerical & Parallel Composability for High Performance Computing

### Joint Inria-Industry Project Team with Airbus Central R&T and Cerfacs

To deal with large size and large dimension problems coming from model- and data-driven applications, Concace will take advantage of modern development tools and languages to design high-level expressions of complex parallel algorithms.

While the traditional approach to HPC is to fully exploit hardware, Concace's complementary approach will enable a richer composability of numerical methods, allowing to fully exploit existing and new numerical algorithms.

<figure><img src="image_34.png" /></figure>

**AIRBUS**

CERFACS

Inria

**Thank you for your attention!**

**Questions ?**
Consuming time

Solve $A(c)u = b$ on 2D domain $64 \times 64$ ($A(c) \in \mathbb{C}^{64^2 \times 64^2}$):

<figure><img src="image_36.png" /><figcaption>(a). PGMRES with varying thresh</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>(b). FGMRES and PGMRES</figcaption></figure>

* PGMRES: GMRES preconditioned by spilu($M$), which is realized by applying sparse ilu to the sparse matrix $M$ with selected elements from coefficient matrix $A$

* The effectiveness of trained U-Net preconditioner is close to (could be the same as) spilu($M$) with thresh= $10^{-5}$ (the best algebra preconditioner for this example)

* Iterations & GPUs time: GMRES: 1978 & 141.8737s; FGMRES: 29 & 4.3403s; PGMRES (thresh= $10^{-3}$): 522 & 43.0996s; PGMRES (thresh= $10^{-4}$): 111 & 3.2276s; PGMRES (thresh= $10^{-5}$): 23 & 0.6582s
Consuming time

Solve $A(c)u = b$ on 2D domain $64 \times 64$ ($A(c) \in \mathbb{C}^{64^2 \times 64^2}$):

<figure><img src="image_35.png" /><figcaption>(a). PGMRES with varying thresh</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>(b). FGMRES and PGMRES</figcaption></figure>

* PGMRES: GMRES preconditioned by spilu($M$), which is realized by applying sparse ilu to the sparse matrix $M$ with selected elements from coefficient matrix $A$

* The effectiveness of trained U-Net preconditioner is close to (could be the same as) spilu($M$) with thresh$= 10^{-5}$ (the best algebra preconditioner for this example)

* Iterations & GPUs time: GMRES: 1978 & 141.8737s; FGMRES: 29 & 4.3403s; PGMRES (thresh$= 10^{-3}$): 522 & 43.0996s; PGMRES (thresh$= 10^{-4}$): 111 & 3.2276s; PGMRES (thresh$= 10^{-5}$): 23 & 0.6582s