\documentclass{nature}

\usepackage[ruled,vlined]{algorithm2e}
\usepackage{amssymb}
\usepackage{xcolor}
\usepackage[utf8]{inputenc}
\usepackage{amsmath}
\usepackage{endfloat} %

\date{\today}

\newcommand{\citeasnoun}[1]{Ref.~\citenum{#1}}
\newcommand{\secref}[1]{Sec.~\ref{#1}}
\newcommand{\Secref}[1]{Section~\ref{#1}}
\renewcommand{\eqref}[1]{Eq.~(\ref{eq:#1})}
\newcommand{\Eqref}[1]{Equation~(\ref{eq:#1})}
\newcommand{\figref}[1]{Fig.~\ref{#1}}

\newcommand{\edit}[1]{{#1}}
\newcommand{\markup}[1]{{#1}}

\usepackage{graphicx}
\usepackage{url}
\title{Physics-enhanced deep surrogates for PDEs} 
\author{Rapha{\"e}l~Pestourie$^{1,\ast}$, Youssef~Mroueh$^{2,3}$, Chris~Rackauckas$^{1}$, Payel~Das$^{2,\ast}$ \& Steven~G.~Johnson$^1$}
\date{}
\date{\today}

\begin{document}

\maketitle

\noindent \normalsize{$^{1}$ MIT, 77 Massachusetts Ave, Cambridge, MA 02139, USA}\\
\normalsize{$^{2}$ IBM Research AI, IBM Thomas J Watson Research Center,  Yorktown Heights, NY 10598, USA}\\
\normalsize{$^{3}$ MIT-IBM Watson AI Lab, Cambridge, MA 02139, USA}\\
\normalsize{$^\ast$Correspondence to: rpestour@mit.edu; daspa@us.ibm.com.}

\begin{abstract}
  We present a ``physics-enhanced deep-surrogate'' (``PEDS'') approach towards developing fast surrogate models for complex physical systems, which is described by partial differential equations (PDEs) and similar models. Specifically, a unique combination of %
a low-fidelity, explainable physics simulator and %
a neural network generator is proposed, which is trained end-to-end to globally match the output of an expensive high-fidelity numerical solver.  We consider low-fidelity models derived from coarser discretizations and/or by simplifying the physical equations, which are several orders of magnitude faster than a high-fidelity ``brute-force'' PDE solver. The neural network generates an approximate input, which is adaptively mixed with a downsampled guess and fed into the low-fidelity simulator. In this way, by incorporating the limited physical knowledge from the differentiable low-fidelity model ``layer'', we ensure that the conservation laws and symmetries governing the system are respected by the design of our hybrid system. Experiments on three test problems---diffusion, reaction--diffusion, and electromagnetic scattering models---show that a PEDS surrogate 
can be \edit{up to} 3$\times$ more accurate than a ``black-box'' neural network with limited data  ($\approx 10^3$ training points), and reduces the data needed by at least a factor of 100 for \edit{a target error of 5\%, comparable to fabrication uncertainty}. PEDS even appears to learn with a steeper asymptotic power law than black-box surrogates. In summary, PEDS provides a general, data-driven strategy to bridge the gap between a vast array of simplified physical models with corresponding
brute-force numerical solvers, offering accuracy, speed, data efficiency, as well as physical insights into the process.
\end{abstract}
\section{Introduction}
In mechanics, optics, thermal transport, fluid dynamics, physical chemistry, climate models, crumpling theory, and many other fields, data-driven surrogate models---such as polynomial fits, radial basis functions, or neural networks---are widely used as an efficient solution to replace repetitive calls to slow numerical solvers~\cite{baker2019workshop, benner2015survey, willard2020integrating, hoffmann2019machine, pant2021deep, pestourie2018inverse}. 
However the reuse benefit of surrogate models comes at a significant training cost, in which a costly high-fidelity numerical solver must be evaluated many times to provide an adequate training set, and this cost rapidly increases with the number of model parameters (the ``curse of dimensionality'')~\cite{boyd2007chebyshev}.  %
In this paper, we explore one promising route to increasing training-data efficiency: incorporating \emph{some} knowledge of the underlying physics into the surrogate by training a generative neural network (NN) ``end-to-end'' with an \emph{approximate} physics model.  We call this hybrid system a ``physics-enhanced deep surrogate'' (PEDS). 
\markup{We demonstrate multiple-order-of-magnitude improvements in sample and time complexity on three different test problems involving the diffusion equation's flux, the reaction-diffusion equation's flux}, and Maxwell's-equations' complex transmission coefficient for optical metamaterials---composite materials whose properties are designed via microstructured geometries~\cite{pestourie2020active}. 
In inverse design (large-scale optimization) of nanostructured thermal materials, chemical reactors, or optical metamaterials, the same surrogate model capturing important geometric aspects of the system may be re-used thousands or millions of time~\cite{lu2022multifidelity,pestourie2018inverse, pestourie2020assume}, making surrogate models especially attractive to accelerate computational design~\cite{bayati2021inverse, li2021inverse}.
To obtain an accurate surrogate of a PDE, we apply a deep NN to \emph{generate a low-fidelity geometry, optimally mixed with the downsampled geometry}, which is then used as an input into an approximate low-fidelity solver and trained end-to-end to minimize the overall error, as depicted in Fig.~\ref{fig:PEDS_diagram} (Sec.~\ref{sec:results}).  The low-fidelity solver may simply be the same numerical method as the high-fidelity PDE solver except at a lower spatial resolution, or it may have additional simplifications in the physics (as in the reaction--diffusion example below, where the low-fidelity model discards the nonlinear term of the PDE).  By design, this low-fidelity solver
yields unacceptably large errors in the target output (perhaps $> 100\%$), but it is orders of magnitude faster than the high-fidelity model while qualitatively preserving at least some of the underlying physics. The NN is trained to nonlinearly correct for these errors in the low-fidelity model, but the low-fidelity model ``builds in'' some knowledge of the physics and geometry that improves the data efficiency of the training. For example, the low-fidelity diffusion model enforces conservation of mass, while the low-fidelity Maxwell model automatically respects conservation of energy and reciprocity~\cite{potton2004reciprocity},  and we can also enforce geometric symmetries; all of these augment the ``trustworthiness''~\cite{li2021trustworthy} of the model. \markup{Compared to a NN-only baseline model (SI, Implementation details of PEDS and baseline), 
\edit{we find that, with a very small dataset of $\approx 1000$ points, PEDS consistently increases the accuracy by up to 3$\times$ compared to the baseline, and reduces the need for training data by an order of magnitude. For the number of parameters of the surrogate models we tested, it amounts to a Cartesian product of less than two points in each input direction. To obtain a $\approx5$\% error, comparable to fabrication uncertainty, PEDS reduces the data need by a factor of at least 100 compared to competing approaches. }}In the more challenging case of our surrogate of the complex optical transmission, PEDS seems to improve the asymptotic \emph{rate} of learning ($\approx 5\times$ larger power law), so that the benefits increase as accuracy tolerance is lowered (Fig.~\ref{fig:resultfigure} and \secref{sec:results}). We show through an ablation study of the surrogate for Maxwell's equations that adding information from the downsampled structure increases the accuracy by 15\% in a low-data regime. %
Furthermore, when the low-fidelity solver layer is very inaccurate, we find that PEDS gains significant additional benefits by combining it with active-learning techniques from our earlier work~\cite{pestourie2020active}, and in fact the benefits of active learning (AL) seem to be even greater for PEDS than for competing approaches.  Although the resulting PEDS surrogate is more expensive to evaluate than a NN by itself due to the low-fidelity solver, it is still much faster than the high-fidelity solver with two to four orders of magnitude speedup. Furthermore, since the NN generates a downsampled version of the geometry, this output can be further examined to gain insight into the fundamental nonlinear physical processes captured by the low-fidelity solver.
\section{Results}
\markup{\subsection{PEDS Framework}
\label{sec:results}
In this work, we illustrate PEDS with three well-known PDEs, as shown in Table~\ref{tab:fourier}, which are implicated in wide varieties of important applications. First, we study the linear diffusion equation, which has applications in materials science, information theory, biophysics and probability, among others. In particular, we train a surrogate model for the thermal flux, which is a useful design property for thermoelectrics. Second, we build a surrogate model for the nonlinear reaction-diffusion equation. This PDE is used in chemistry and its surrogates can influence the design of chemical reactors. Third, we model the complex transmission of Maxwell's equations through a parameterized structure, which is typically used in the design of optical metamaterials~\cite{pestourie2020active, pestourie2018inverse, pestourie2020assume}.}

\begin{table}[h!]
\centering
\begin{tabular}{lll}
\hline
Equation name             & Equation formula & Model(\textit{input dimension})\\
\hline
Diffusion          & $\nabla\cdot D\nabla \textbf{u}= \textbf{s}_0$ & Fourier($d$)\\ %
Reaction-diffusion &  $\nabla\cdot D\nabla \textbf{u}= -k\textbf{u}(1-\textbf{u})+\textbf{s}_0$ & Fisher($d$)\\ 
2D Maxwell (Helmholtz)          & $\nabla^2\textbf{u}-\omega^2\varepsilon\textbf{u}=\textbf{s}_1$ & Maxwell($d$)\\
\hline
\end{tabular}
\caption{Governing equations of the surrogate models for our example problems. $d$ is the input dimension, i.e. the number of input variables in the surrogate model, which ranges from $10$ to $25$.}
\label{tab:fourier}
\end{table}
\begin{figure}[h!]
    \centering
    \includegraphics[width=\textwidth]{PEDS_diagram5.png}
    \caption{Diagram of PEDS: (Main) From the geometry parameterization, the surrogate generates a low-fidelity structure that is combined with a  downsampled geometry (e.g. downsampled by pixel averaging) to be fed into a low-fidelity solver (symbolized by a cartoon picture of James Clerk Maxwell). (Inset) The training data is generated by solving more costly simulations directly on a high-fidelity solver (symbolized by a photograph of James Clerk Maxwell).}
    \label{fig:PEDS_diagram}
\end{figure}
Before delving into implementation details and results, we present the core principles of PEDS which are common between all surrogates.

\subsubsection{Model and Methods}

The PEDS surrogate model $\tilde{f}(p)$ aims to predict $f^{hf}(\mathrm{hf}(p))$---an output property of interest as it would be computed from a computationally intensive high-fidelity (hf) solver $f^{hf}$. The hf solver computes the PDE solution for a high-fidelity geometry $\mathrm{hf}(p)$, with $p$ being some parameterization of the geometry (or other system parameters). PEDS is depicted schematically in~\figref{fig:PEDS_diagram}, and is implemented in the following stages:  %

\begin{enumerate}
    \item Given the parameters $p$ of the geometry,
a deep generative NN model yields a grid of pixels describing a %
low-fidelity geometry.  We call this function $\mathrm{generator}_\mathrm{NN}(p)$.

    \item We also compute a low-fidelity downsampling (e.g. via sub-pixel averaging~\cite{oskooi2009accurate}) of the geometry, denoted $\mathrm{downsample}(p)$; other prior knowledge could also be incorporated here as well.
    
    \item We define $G$ as a weighted combination $G(p) = w\cdot \mathrm{generator}_\mathrm{NN}(p) + (1-w)\cdot \mathrm{downsample}(p)$, with a weight $w\in[0,1]$ (independent of $p$) that is another learned parameter.
    
    \item If there are any additional constraints/symmetries that the physical problem imposes on the geometry, they can be applied as projections $P[G]$.  For example, mirror symmetry could be enforced by averaging $G$ with its mirror image.
    
    \item Finally, given the low-fidelity geometry $P[G(p)]$, we evaluate the low-fidelity solver $f^\mathrm{lf}$ to obtain the property of interest:  $\tilde{f}(p) = f^\mathrm{lf}(P[G(p)])$.
    
\end{enumerate}
In summary, the PEDS model $\tilde{f}(p)$ is
\begin{equation}
\tilde{f}(p) = f^\mathrm{lf}\left(P\left[w\cdot\mathrm{generator}_\mathrm{NN}(p) + (1-w)\cdot \mathrm{downsample}(p)\right]\right) \, .
    \label{eq:model}
\end{equation}

\paragraph{Dataset acquisition } PEDS is a supervised model that is trained on a labeled dataset.  We build the training set by querying the high-fidelity solver with parameterized geometries $S=\{ (p_i, t^{hf}_i) , i=1 ... N\}$, where $p_i$ are parameterized geometries in the training set and $t^{hf}_i=f^{hf}(p_i)$. The upfront cost of building the training dataset is the most time-consuming part of developing a supervised surrogate model $\tilde{f}(p)$. By building some approximate low-fidelity physics knowledge into the surrogate, we will show that PEDS greatly reduces the number $N$ of queries to expensive simulations.

\paragraph{Training loss }A basic PEDS training strategy could simply minimize the mean squared error $\sum_{(p,t^\mathrm{hf})\in S}|\tilde{f}(p) - t^\mathrm{hf}|^2$ (for a training set $S$) with respect to the parameters of the NN and the weight~$w$. When the data may have outliers, we use a Huber loss~\cite{huber1992robust}.

\begin{equation}\label{eq:huber}
    L_\delta (a) = \begin{cases}
 \frac{1}{2}{a^2}                   & \text{for } |a| \le \delta, \\
 \delta \cdot \left(|a| - \frac{1}{2}\delta\right), & \text{otherwise.}
\end{cases}
\end{equation}

We also employ a  more complicated loss function that allows us to easily incorporate active-learning strategies~\cite{pestourie2020active}.  We optimize the Gaussian negative log-likelihood of a Bayesian model~\cite{lakshminarayanan2016simple}
\begin{equation}\label{eq:loglikelihood}
    -\sum_{(p_i, t^{hf}_i)\in S} \log{\mathrm{P}_\Theta(t^{hf}_i|p_i)} \propto \sum_{(p_i, t^{hf}_i)\in S} \left[ \log{\sigma(p_i)} + \frac{(t^{hf}_i-\tilde{f}(p_i))^2}{2 \sigma(p_i)^2} \right]
\end{equation}
where $\mathrm{P}_\Theta$ is a Gaussian likelihood defined by $\Theta$ which includes the parameters of the generator model parameters and the combination weight $w$, and the heteroskedastic ``standard deviation'' $\sigma(p) > 0$ is the output of another NN (trained along with our surrogate model). 

\paragraph{Ensemble model} We also train surrogates that are an \emph{ensemble} of 5 independent surrogates. The prediction of the ensemble is the average of the predictions of each individual model.

\paragraph{Stochastic gradient descent }In practice, rather than examining the entire training set $S$ at each training step, we follow the standard ``batch'' approach~\cite{goodfellow2016deep} of sampling a random subset of $S$ and minimizing the expected loss with the Adam stochastic gradient-descent algorithm~\cite{kingma2014adam} (via the Flux.jl~\cite{innes:2018} software in the Julia language).

\paragraph{Adjoint method} The low-fidelity solver is a layer of the PEDS model, which is trained end-to-end, so we must backpropagate its gradient $\nabla_g f^\mathrm{lf}$ with respect to the low-fidelity geometry input $g$ through the other layers to obtain the overall sensitivities of the loss function. This is accomplished efficiently using the known ``adjoint'' methods~\cite{molesky2018inverse}. Such methods  yield a vector-Jacobian product that is then automatically composed with the other layers using automatic differentiation~(AD) (via the Zygote.jl~\cite{innes2018don} software). 

In particular, the low-fidelity solver layer is differentiable because each pixel of the low-fidelity geometry is assigned to a sub-pixel average of the infinite-resolution structure, which increases accuracy~\cite{oskooi2009accurate} and makes $\mathrm{downsample}(p)$ piecewise differentiable. In the same way, $\mathrm{hf}(p)$ is differentiable for the high-fidelity geometry.

\label{sec:model}

\markup{\paragraph{PEDS for diffusion equation} Our first two surrogate models are for the diffusion equation from Table~\ref{tab:fourier}. They are called Fourier($16$) and Fourier($25$), and they predict the thermal flux $\kappa(p)$ from the diffusion equation for 16 and 25 input variables, respectively. As showed in Fig.~\ref{fig:fffigure}~(left),  the 2D nanostructured material defines the coefficient matrix $D(p)$ where the parameter vector $p$ contains the 25 (resp. 16) independent side lengths of a five by five (resp. four by four) grid of air holes etched in the medium. The thermal conductivity coefficients in $D$ are set to 1 in the medium and 0.1 in the holes. The boundary conditions are periodic in $x$-direction and Dirichlet boundary conditions in the $y$ direction, fixing the temperature to $1$ at the bottom and to $0$ at the top, as illustrated by thick red and blue lines in Fig.~\ref{fig:fffigure}~(left). The Dirichlet boundary conditions are equivalent to the source term $\textbf{s}_0$ in Table~\ref{tab:fourier}.
Both the high-fidelity and the low-fidelity solvers employ a finite-difference solver that represents the geometry by a grid of discretized thermal conductivity. Sub-pixel averaging is employed at the boundary between the holes and the medium. For both Fourier($16$) and Fourier($25$), the high-fidelity solver has a resolution of 100. The low-fidelity solver has a resolution of 4 or 5, which corresponds to a single pixel per hole position. Each high-fidelity data point acquisition requires $\approx 35$~ms, and each low-fidelity data point acquisition requires $\approx 65~\mu$m and $\approx 75~\mu$m, respectively, which represents a speed-up of $\approx 500\times$ (Table~\ref{tab:lowfidresult}, Speedup). We compute the low-fidelity solver baseline error, by computing the solution with the low-fidelity solver and the geometry $\mathrm{downsample}(p)$, where $p$ is the geometry parameterization (i.e. without mixing with a neural generator output). Despite the much lower resolution, the low fidelity solvers have a fairly low error of 13.5\% and 8.5\%, respectively. This good performance of an averaged structure comes from the fact that the diffusion equation is a smoothing equation. Nonetheless, such errors would still be dominant compared to typical experimental uncertainties of $\approx$5\%.
Fourier($16$) and Fourier($25$) were trained to predict the flux through a plane as in Fig.~\ref{fig:fffigure}~(middle) by minimizing Huber loss in Eq.~\ref{eq:huber} with $\delta=10^{-3}$ to lower the sensitivity to outliers.}

\markup{\paragraph{PEDS for reaction--diffusion equation} Our next two surrogate models solve the reaction--diffusion equation from Table~\ref{tab:fourier}, and are called Fisher($16$) and Fisher($25$). They predict the flux $\kappa(p)$ through the same geometry as Fourier($16$) and Fourier($25$), respectively. As can be seen in Table~\ref{tab:fourier} the reaction--diffusion equation has an additional nonlinear term $k \textbf{u}(1-\textbf{u})$ compared to the diffusion equation. $k$ is a coefficient that controls the amount of nonlinearity in the PDE. In Fig.~\ref{fig:fffigure}~(middle and right), we see how much the nonlinearity impacts the PDE solution. The high-fidelity nonlinear solver is using finite different and Newton's method in conjunction with a continuation method that increases k from 0.1 to 10 in 5 multiplicative steps. The low-fidelity solvers of Fisher($16$) and Fisher($25$) are identical to that of Fourier($16$) and Fourier($25$), respectively. Importantly, the low-fidelity solver not only has a coarse resolution, but also uses an approximate physics that neglects the nonlinear term from the reaction--diffusion equation. Each high-fidelity data point requires $\approx700$~ms that is around $10^4\times$ slower than the low-fidelity solver (Table~\ref{tab:lowfidresult}, Speedup). The low-fidelity solvers have error of 38.1\% and 36.7\% respectively. Fisher($16$) and Fisher($25$)  were trained to predict the flux through a plane as in Huber loss in Eq.~\ref{eq:huber} with $\delta=10^{-3}$ to lower the sensitivity to outliers.}

\begin{figure}[h!]
    \centering
    \includegraphics[width=\textwidth]{figurefourierfisher2.png}
    \caption{(Left) Geometry with 5 by 5 air holes with varying widths. There are Dirichlet boundary conditions on top (blue line) forcing the temperature to 0 and at the bottom (red line) forcing to 1, and periodic boundary conditions on the sides. (Middle and Right) Temperature field for the diffusion equation and the reaction diffusion equation, respectively. The orange dotted line is where the flux is evaluated to compute $\kappa$.}
    \label{fig:fffigure}
\end{figure}

\paragraph{PEDS for Maxwell's equations} Similarly to~\citeasnoun{pestourie2020active}, our third surrogate model Maxwell($10$) predicts the complex transmission $t^{hf}(p)$ of a 2D ``meta-atom'' unit cell with a parameterized geometry $p$, which consists of ten layers of air holes with independent widths etched in a substrate (of dielectric constant $\varepsilon = 2.1$ corresponding to silica), with periodic boundary conditions in $x$ and outgoing radiation boundary conditions in the $y$ direction and an incoming normal-incident planewave from below, as shown in Fig.~\ref{fig:resultfigure}~(right). 
In terms of the vacuum wavelength $\lambda$ of the incident wave (for the largest $\lambda$ considered below), the period in $x$ is $0.95\lambda$ and the total thickness is $11\lambda$ (with hole heights of 0.75$\lambda$ and interstices of 0.35$\lambda$); the fact that the structure is several wavelengths in diameter causes the transmission $t^{hf}(p)$ to be a complicated oscillatory function that makes the surrogate training challenging~\cite{pestourie2020active}.  A ``metasurface'' consists of a collection of many of these meta-atoms, designed to perform some optical function such as focusing~\cite{li2021inverse}.  The full solution for a metasurface can be approximated in terms of the transmissions of individual periodic `unit cells via a local periodic approximation~\cite{pestourie2018inverse, pestourie2020assume}. A schematic unit cell with 3~holes is showed in Fig.~\ref{fig:PEDS_diagram}, and an example 10-hole structure from the training set is shown in Fig.~\ref{fig:resultfigure}~(right).

Both the high-fidelity and low-fidelity solvers
for Maxwell($10$) employ finite-difference frequency-domain (FDFD) discretizations of Maxwell's equations~\cite{champagne2001fdfd}, using perfectly matched layers (PMLs)~\cite{sacks1995perfectly} to implement outgoing boundary conditions.  Similarly to the solvers of the two previous equations,  FDFD represents the geometry by a grid of discretized $\varepsilon$ ``pixels,'' which is a function of the parameters (hole widths) $p$, $\mathrm{hf}(p)$, and $\mathrm{downsample}(p)$ for the high-fidelity solver and the baseline coarse solver, respectively. An FDFD resolution of 40 pixels per wavelength is used as our high-fidelity solver.
This resolution is typical for high-fidelity solvers in electromagnetism, because it is comparable to the manufacturing accuracy in nanophotonics and hence suffices for practical metalens design~\cite{li2021inverse, bayati2021inverse} within fabrication uncertainty. (Sharp/narrowband resonances can shift if one refines the resolution further, but the positions and the bandwidths of the resonances are accurate to within a few percent.) Each high-fidelity-solver data point required $\approx 1$~s (on a 3.5 GHz 6-Core Intel Xeon E5); an analogous simulation in 3D takes several hours.  Our PEDS surrogate uses an FDFD solver at a coarser resolution of 10 pixels per wavelength, which is about $100\times$ faster in 2D and $> 10^4\times$ faster in 3D, but has much worse accuracy. 
It differs from the high-fidelity solver's transmission by $124$\% on our test set, which is significantly more than the four other surrogates presented in this article. Maxwell($10$) model was trained to predict the complex transmission for 3 frequencies by minimizing the negative Gaussian likelihood loss function to enable comparison with and without using AL~\cite{pestourie2020active}. The input of the model $p$ is the concatenation of the 10 widths and the one-hot encoding of the frequency. 
\markup{\subsection{Overall benefits of PEDS}
Most importantly, in a low-data regime ($\approx10^3$ data points for $10$ to $25$ input parameters), we report that PEDS consistently increases the accuracy by \edit{up to $3\times$ and reduces the data needed by at least an order of magnitude. All PEDS surrogates reduce the need for training data by a factor of $>100$ to attain an error level of 5\% comparable to uncertainties in experiments (Table~\ref{tab:nnonlyresult}, Fig.~\ref{fig:resultfigure}), which is sufficient for design purposes}. In the case of Fourier($16$) and Fourier($25$), the mixing weight  $w$ of the neural generated structures is around $0.1$, whereas for Fisher($16$) and Fisher($25$), the mixing weight $w$ is around $0.45$. Since the low-fidelity solver is more inaccurate for the nonlinear reaction--diffusion equation where the linear relaxation results in errors $>0.35\%$, the neural generator has approximately a $5\times$ larger weight, indicating it has the stronger impact of including the nonlinear effects in PEDS. We report the exact optimal combining weights in (SI, Table 1) for Fourier($16$), Fourier($25$), Fisher($16$), and Fisher($25$).
Performance in a low-data regime are summarized in Table~\ref{tab:nnonlyresult} for accuracy improvement, computed as the fractional error (FE) on a test set (SI, fractional error). For Fourier($16$), Fourier($25$), Fisher($16$), Fisher($25$), and Maxwell($10$), the error of PEDS goes down to typical levels of experimental uncertainties of 3.7\%, 3.8\%, 4.5\%, 5.5\%, and 19\% respectively.

We compared Fourier($16$), Fourier($25$), Fisher($16$), Fisher($25$), and Maxwell($10$) against a NN-only baseline, which consists of an ensemble of neural networks with the same number of parameters as PEDS generators with an additional fully connected layer to replace PEDS low-fidelity solver layer (Table~\ref{tab:nnonlyresult}). \edit{With 1000 training points}, PEDS is an improvement compared to the neural network baseline of \edit{up to 3$\times$ (Table~\ref{tab:nnonlyresult}, PEDS ($\approx 10^3$) and NN-only ($\approx 10^3$)). Furthermore, the neural network baseline still cannot reach the reported PEDS accuracies when given an order of magnitude more data, which means that PEDS saves at least an order of magnitude in data (Table~\ref{tab:nnonlyresult}, NN-only ($\approx 10^4$)). Except Maxwell(10), the NN-only baselines cannot reach PEDS error with two orders of magnitude more data (Table~\ref{tab:nnonlyresult}, NN-only ($\approx 10^5$)). In particular for Fourier surrogates, going from $10^4$ to $10^5$ points reduces the error by less that $0.1\%$. Except Maxwell(10), which is further discussed in Section~\ref{sec:AL}, PEDS achieves error of $5\%$ in low-data regime (1000 training points), and reduces the data need by a factor of at least 100.}}

\begin{table}[h!]
\begin{tabular}{lllll}
\hline
Model(\textit{input dim})  & PEDS ($\approx 10^3$) & NN-only ($\approx 10^3$) & NN-only ($\approx 10^4$) & NN-only ($\approx 10^5$) \\ 
\hline
Fourier(16) & 3.7\% & 5.1\% & 4.8\% & 4.8\% \\
Fourier(25) & 3.8\% & 4.7\% & 4.4\% & 4.4\% \\
Fisher(16)  & 4.5\% & 10.1\% & 9.9\% & 9.5\% \\
Fisher(25)  & 5.5\% & 14.4\% & 14.0\% & 12.7\% \\
Maxwell(10)  & 19\% (AL) & 56\% & 19\% & 15\% \\
\hline
\end{tabular}
\caption{PEDS error versus NN-only baselines' errors (mean fractional error on the test set). We report the orders of magnitude of training points in parenthesis. With more than an order of magnitude extra data, NN-only baseline still has much higher error than PEDS.  Except Maxwell(10), all baselines still cannot achieve PEDS error with two orders of magnitude extra data. The improvement when going from $10^4$ to $10^5$ points with Fourier surrogates are smaller than $0.1\%$. In the Maxwell case, we show in section 3.3 that it is crucial to include active learning (AL) in addition to PEDS.}
\label{tab:nnonlyresult}
\end{table}

\markup{We further compared PEDS to a low-fidelity solver baseline, which uses the low-fidelity solver with $\mathrm{downsample}(p)$ as input, without mixing with the low-fidelity geometry generated by the neural network (Table~\ref{tab:lowfidresult}).  PEDS also boosts the accuracy of the low-fidelity solver by $3.6\times$, $2.2\times$, $8.5\times$, $6.7\times$, and $6.5\times$, respectively (Table~\ref{tab:lowfidresult}, Improvement). For the reaction--diffusion equation, the low-fidelity solver has a coarser resolution and a linear approximation of the physics (neglecting the nonlinear term of reaction--diffusion equation), but the neural network generator captures the necessary nonlinearity to get improvement $> 5\times$ (Table~\ref{tab:lowfidresult}, Improvement).  The speedups vary between two and four orders of magnitude (Table~\ref{tab:lowfidresult}, Speedup). For Maxwell($10$), using a coarser low-fidelity solver generally gains two orders of magnitude in 2D, which should translate into a four orders of magnitude speedup for three-dimensional problems. We see the biggest speedups when the low-fidelity solver is not only coarser than the high-fidelity solver, but also when it is a linear relaxation of the physics (reaction--diffusion equation). In that case, the speedup is four orders of magnitudes.}
\begin{table}[h!]
\begin{tabular}{llllll}
\hline
Model(\textit{input dim}) & PEDS error ($\approx 10^3$) & Low-fidelity error & Improvement & Speedup \\
\hline
Fourier(16) & 3.7\% & 13.5\% & $3.6\times$ & 500$\times$ 
\\%35ms/65µs=538 
Fourier(25) & 3.8\% & 8.5\% & $2.2\times$ & 500$\times$ 
\\%35ms/75µs=466
Fisher(16) & 4.5\% & 38.1\% & $8.5\times$ & $10^4\times$ 
\\% 700ms/65µs=10.8k     \\ 
Fisher(25) & 5.5\% & 36.7\% & $6.7\times$ & $10^4\times$ 
\\ %
Maxwell(10) & 19\% (AL) & 124\%  & $6.5\times$ & $10^2\times$ / $10^4\times$\\
\hline
\end{tabular}
\caption{With $\approx 10^3$ training points, PEDS consistently improves error (mean fractional error on the test set) by 2--8$\times$ compared to the low-fidelity solver. ``Improvement'' is the reduction in error by PEDS compared to the low-fidelity.  Speedups are shown for 2D simulations, and speedup for 3D simulations is also reported for Maxwell($10$)}
\label{tab:lowfidresult}
\end{table}
\subsection{Detailed analysis of Maxwell(10) case study}\label{sec:AL}
In previous section, we showed the general performance of PEDS in the low-data regime. For Maxwell($10$), where the low-accuracy solver has a very large error ($>100\%$), we study the training curve asymptotically and when combining with AL~\cite{pestourie2020active}. In contrast to the previous section, where we performed static training that takes a training set sampled at random, here we discuss results from AL experiments by dynamic Bayesian training, where the training set is iteratively expanded using an AL algorithm~\cite{pestourie2020active}. Essentially, AL attempts to sample training points where the model uncertainty is highest, thereby reducing the number of costly point acquisitions by querying the high-fidelity solver. 
Our previous work showed an order of magnitude improvement in terms of data efficiency by using AL, when compared to a black-box NN~\cite{pestourie2020active}. Consistently, in this study, we also report substantial improvements from active learning for PEDS.

The active-learning algorithm iteratively builds a training set by filtering randomly generated points with respect to a trained measure of uncertainty~\cite{pestourie2020active}. The hyperparameters of this algorithm are (i) $n_\mathrm{init}$, which is the number of points the surrogate models is initially trained with; (ii) $T$, the number of exploration iteration; (iii) $M$ and $K$, which are such that $M\times K$ points are randomly generated at each iteration and only $K$ points with highest uncertainty $\sigma(p)$ are explored (SI, Active learning implementation details). We run the expensive high-fidelity solver to get the PDE solutions of the explored points. %
We have trained surrogates as well as an \emph{ensemble} of 5 independent surrogates. We found that models optimizing the negative log-likelihood perform similarly to models optimizing the mean squared error in the case static training. This is not surprising, because the mean squared error is part of the negative log-likelihood objective. 
\label{sec:accuracy}

\begin{figure}[h!]
    \centering
    \includegraphics[width=\textwidth]{resultfigure12.png}
    \caption{(Left) Fractional error (FE) on the test set: PEDS outperforms the other baseline models significantly when combined with active learning (AL). 
    (Right) Geometry of the unit cell of the surrogate model. Each of the 10 air holes have independent widths, the simulation is performed with periodic boundary conditions on the long sides, the incident light comes from the bottom and the complex transmission is measured at the top of the geometry.}
    \label{fig:resultfigure}
\end{figure}

We compared PEDS to a NN-only baseline 
using the fractional error as an evaluation metric~(SI, Implementation details of PEDS and baselines).
In Fig.~\ref{fig:resultfigure}, we show that PEDS clearly outperforms all other models when combined with active learning. In low-data regime, it is $2.9\times$ more accurate than the baseline. Asymptotically, in high-data regime, it converges to the true value with a power law exponent $5\times$ better, with a slope of -0.5, in contrast to -0.1, for the baseline on the loglog plot. %

From a data-efficiency perspective, the PEDS+AL solver achieves 20\% error on the test set, while using only about $5\%$ of the training data needed to train the NN-only baseline, and $12.5\%$ of the training data needed to train the NN-only baseline with AL (Fig.~\ref{fig:resultfigure}).
Only PEDS+AL reaches a low $3.5$\% error with a training data size of $\approx500k$ (Fig.~\ref{fig:resultfigure}). However, if we extrapolate the other curves in Fig.~\ref{fig:resultfigure}, it is clear that they would require at \emph{least} two orders of magnitude more data to achieve similar low error. \edit{This completes the claim that PEDS saves at least two orders of magnitude in training data to achieve and error comparable to fabrication uncertainty.}

Evaluating the baseline (with an ensemble of neural networks) takes 500~$\mu s$, while PEDS evaluates in $5$~ms, which is about a ten times slower. However the high-fidelity solver is about a hundred times slower, evaluating at $\approx1$~s. In order to simulate the data set quickly, and without loss of generality, we showed results for PEDS in 2D (Fig.~\ref{fig:resultfigure}~(right). As PEDS is already faster than the high-fidelity model by two orders of magnitude, this difference will be even starker for 3D simulations. The simulation of the equivalent structure in 3D evaluates in about $100$~ms with the low-fidelity model, and in $2462$~s with the high-fidelity model. In this occurrence, PEDS would represent a speed-up by at least four orders of magnitude.  

\subsubsection{Ablation study} 

Next, we show results of ablation experiments in order to understand the effect of mixing the generated structure with a downsampled structure. Specifically, we performed an ablation study on an AL ensemble model in the low-data regime (1280 training points); results are shown in Table~\ref{tab:ablation}. The edge cases of using only the downsampled structure with the low-fidelity solver (Table~\ref{tab:ablation}, coarsified only) performs the worst (124\% error with respect to the high-fidelity solver), corresponding to $w=0.0$ in \eqref{model}. Conversely, using the NN generator only (Table~\ref{tab:ablation}, generator only), corresponding to $w=1.0$ in \eqref{model}, is still about 15\% worse (0.20 error) than using adaptive mixing $0 < w < 1$ (Table~\ref{tab:ablation}, PEDS). Imposing mirror symmetry, via $P[G] = (G + \mbox{mirror image})/2$ in \eqref{model} (Table~\ref{tab:ablation}, PEDS with symmetry), did not improve the accuracy of the model in this case (but is a useful option in general, since symmetry may have a larger effect on the physics in other applications).
\begin{table}[h!]
\centering
\begin{tabular}{lll}
\hline
Generative model for low-fidelity geometry & FE on test set & PEDS improvement \\\hline
$w = 0.0$ (coarsified only)                      & 1.24      & 86\%             \\
$w = 1.0$ (generator only)                   & 0.20      &           15\%  \\
PEDS with symmetry
&  0.18            &     5\%    \\
PEDS
&  0.17            &     ---    \\
\hline
\end{tabular}
\caption{Ablation study of PEDS with ensembling and active learning for 1280 training points, showing the impact of mixing generated and coarsified geometries, as well as imposing symmetry.}
\label{tab:ablation}
\end{table}

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{generatedstudy3.png}
    \caption{ (Left) First 9 principal components which explain most of the variation in the complex transmission. (Right) Coordinate of randomly generated structures on the two first principal components. Clusters can clearly discriminate the input geometries ($f=0.5$ in blue, $f=0.75$ in orange, $f=1.0$ in green). (Insets) Example generated geometries corresponding to the three frequencies of the surrogate model. The generated geometry is smoothest for the smallest frequency.}
    \label{fig:generatedstudy}
\end{figure}

\subsubsection{Analysis of generated geometries} Because the trained PEDS model includes a NN that generates ``equivalent'' coarse-grained geometries to the input structure, it is interesting to analyze these geometries and potentially extract physical insights.

\paragraph{Frequency dependence }
The neural network generates structures that are qualitatively different as a function of the input frequency (Fig.~\ref{fig:generatedstudy}, right insets). As might be expected on physical grounds (e.g. effective-medium theory~\cite{holloway2011characterizing}), the lowest frequency (longer wavelengths) corresponds to the smoothest generated structures, because the wavelength sets the minimum relevant lengthscale for wave scattering. To help quantify this, we performed a principal components analysis (PCA) of $\mathrm{generator}_\mathrm{NN}(p)$ for $10^5$ uniform random $p$ values (including random frequency). We show the first few principal components in Fig.~\ref{fig:generatedstudy}~(left). The first and second components explain 67\% and 13\% of the variation, respectively. We show in Fig.~\ref{fig:generatedstudy}~(right) that the coordinates of the first two components are sufficient to classify generated geometries according to the input frequency.

\paragraph{Scattering richness } To explore the effect of additional scattering physics produced by multiple layers of holes, we generated coarse geometries for different numbers of layers (equivalently, fixing the parameters of the ``deleted'' layers to zero).  We then decomposed the resulting $\mathrm{generator}_\mathrm{NN}(p)$ into the PCA components from above. As we increase the number of layers, the average coordinates of some principal components monotonically increase in magnitude. Since we know that more layers contain more scattering richness, the corresponding principal components geometries provides some geometrical insight into how scattering richness translates into the generated structure. From our analysis of generated structures for the smallest frequency, the first principal component geometry clearly contributes to scattering richness, with an average coordinate (across $10^3$ generated structures) increasing -11 to 26 as the number of layers goes from 1 to 9.
\section{Discussion}
\label{sec:discussion}

The significance of the PEDS approach is that it can easily be applied to a wide variety of physical systems.  It is common across many disciplines to have models at varying levels of fidelity, whether they simply differ in spatial resolution (as in Fourier($16$), Fourier($25$), and Maxwell($10$)) or in the types of physical processes they incorporate (as in Fisher($16$) and Fisher($25$)). %
For example, in fluid mechanics the low-fidelity model could be Stokes flow (neglecting inertia), while the high-fidelity model might be a full Navier--Stokes model (vastly more expensive to simulate)~\cite{ferziger2002computational}, with generator NN correcting for the deficiencies of the simpler model.  As another example, we are currently investigating a PEDS approach to construct a surrogate for complex Boltzmann-transport models~\cite{romano2021openbte} where the low-fidelity heat-transport equation can simply be a diffusion equation. Knowledge of priors can also be introduced in the low-fidelity geometry that is mixed with the neural generator output. PEDS provides a data-driven strategy to connect a vast array of simplified physical models with the accuracy of brute-force numerical solvers, offering both more insight and more data efficiency than physics-independent black-box surrogates.

When compared to related works, PEDS should not be confused with physics-informed neural networks~(PINNs), which solve the full PDE (imposed pointwise throughout the domain) for the entire PDE solution (\emph{not} a surrogate for a finite set of outputs like the complex transmission or the thermal flux)~\cite{karniadakis2021physics, lu2021physics}, and which do not employ any pre-existing solver. Current PINNs tend to be slower than conventional high-fidelity PDE solvers (e.g. based on finite elements)~\cite{shin2020convergence}, but offer potentially greater flexibility. Universal ordinary differential equations (UODEs)~\cite{rackauckas2020universal} also tackle a different problem from PEDS: they identify unknown dynamics in an ODE by replacing the unknown terms with neural networks trained on data. In contrast to DeepONet~\cite{lu2021learning, lu2022multifidelity} and Fourier neural operators~\cite{li2020fourier}, PEDS includes a numerical solver layer.  Our approach has some similarities with input space mapping (SM)~\cite{koziel2008space}, especially neural SM~\cite{bakr2000neural} and coarse mesh/fine SM~\cite{feng2019coarse}, where the input of a fine solver is mapped into the input of a coarse solver. However SM uses the same parameterization for the fine solver and the coarse solver, rather than mapping to ``downsampled'' resolution, and does not mix the generated input with a downsampled guess adaptively. We show that PEDS substantially outperforms SM in the SI (SM baseline). Finally, in contrast to error-correction techniques at the output level of the surrogate~\cite{lu2020extraction, koziel2006space}, PEDS includes the solver in an end-to-end fashion during the training process. In PEDS, the output of the low-fidelity solver layer is not further transformed, which preserves key properties of the low-fidelity solver such as conservation of energy or mass. Mappings between coarse and fine descriptions of a system is also leveraged in the renormalization group technique in physics~\cite{weinberg1995quantum}, but in the latter context this is accompanied by a change of scale---often to investigate self-similar phenomena---and not necessarily a change in the number of degrees of freedom.

In addition to applying the PEDS approach to additional physical systems, there are a number of other possible technical refinements.   For example, one could easily extend the PEDS NN to take an image of the high-fidelity-structure geometry rather than its parameterization, perhaps employing  convolutional neural networks to represent a translation-independent ``coarsification'' and/or a multiresolution architecture.    This type of surrogate could then be employed for topology optimization in which ``every pixel'' is a degree of freedom~\cite{molesky2018inverse}.  Another interesting direction might be to develop new low-fidelity
physics models that admit ultra-fast solvers but are too inaccurate to be used \emph{except} with PEDS;
for instance, mapping Maxwell's equations in 3D onto a simpler
(scalar-like) wave equation or mapping the materials into
objects that admit especially efficient solvers (such
as impedance surfaces~\cite{perez2018sideways} or compact objects for surface-integral equation methods~\cite{jin2015finite}).
\section*{Data Availability Statement}

The data that support the findings of this study is available from the corresponding
author upon reasonable request.

\section*{\edit{Code Availability statement}}
\edit{The code used for these findings is available upon reasonable request.}
\section*{Acknowledgements}

R.P. was supported by the U.S. Army Research Office through the Institute for Soldier Nanotechnologies (Award No. W911NF-18-2-0048) and the MIT-IBM Watson AI Laboratory (Challenge No. 2415). The authors thank Meredith Dost for her suggestions in proof reading.

\section*{Competing interests}
The authors declare no competing financial or non-financial interests.

\section*{Author contributions}
R.P., Y.M., C.R., P.D., and S.G.J. designed the study, contributed to the machine-learning approach, and analyzed results; R.P. led the code development, software implementation, and numerical experiments; R.P. and S.G.J. were responsible for the physical ideas and interpretation. All authors contributed to the algorithmic ideas and writing.

\section*{References}
\bibliographystyle{naturemag}
\bibliography{refs.bib}
\end{document}


Peer review:
Nice new method for developing deep surrogate PDE solvers in a data-efficient manner
The paper proposes a natural approach for improving sample (data) efficiency in deep surrogate PDE solvers. The high level idea is to use two spatial scales (first train a "low-fidelity" model and use its outputs as side information for a final, high fidelity model); this in combination with an active learning strategy leads to data reductions of upto 10x.

The paper is nicely written and the topic is timely.

Points of feedback: a) consider validating on more challenging complex systems. b) why only 2 scales? is there room for an extension to a hierarchy of coarse-to-fine mappings (analogous to multiresolution analysis).
Rating: 7: Good paper, accept
Confidence: 4: The reviewer is confident but not absolutely certain that the evaluation is correct