mathx”17

# Physics-informed deep learning and compressive collocation for high-
dimensional diffusion-reaction equations:
practical existence theory and numerics

Simone Brugiapaglia<EMAIL_ADDRESS>
Department of Mathematics and Statistics
Concordia University
Montréal, QC, Canada Nick Dexter<EMAIL_ADDRESS>
Department of Scientific Computing
Florida State University
Tallahassee, FL, USA Samir Karam<EMAIL_ADDRESS>
Department of Mathematics and Statistics
Concordia University
Montréal, QC, Canada Weiqi Wang<EMAIL_ADDRESS>
Department of Physics & Astronomy
University of Victoria
Victoria, BC, Canada

###### Abstract

On the forefront of scientific computing, Deep Learning (DL), i.e., machine
learning with Deep Neural Networks (DNNs), has emerged a powerful new tool for
solving Partial Differential Equations (PDEs). It has been observed that DNNs
are particularly well suited to weakening the effect of the _curse of
dimensionality_ , a term coined by Richard E. Bellman in the late ‘50s to
describe challenges such as the exponential dependence of the sample
complexity, i.e., the number of samples required to solve an approximation
problem, on the dimension of the ambient space. However, although DNNs have
been used to solve PDEs since the ‘90s, the literature underpinning their
mathematical efficiency in terms of numerical analysis (i.e., stability,
accuracy, and sample complexity), is only recently beginning to emerge. In
this paper, we leverage recent advancements in function approximation using
sparsity-based techniques and random sampling to develop and analyze an
efficient high-dimensional PDE solver based on DL. We show, both theoretically
and numerically, that it can compete with a novel stable and accurate
compressive spectral collocation method. In particular, we demonstrate a new
_practical existence theorem_ , which establishes the existence of a class of
trainable DNNs with suitable bounds on the network architecture and a
sufficient condition on the sample complexity, with logarithmic or, at worst,
linear scaling in dimension, such that the resulting networks stably and
accurately approximate a diffusion-reaction PDE with high probability.

Keywords: physics-informed neural networks, compressive Fourier collocation,
numerical methods for high-dimensional PDEs, practical existence theorem

##  1 Introduction

PDEs over high-dimensional domains are a powerful mathematical modelling tool
adopted in a variety of applications including molecular dynamics,
computational finance, optimal control, and statistical mechanics. Important
high-dimensional Partial Differential Equation (PDE) models in these areas are
the many-electron Schrödinger equation, the Hamilton–Jacobi–Bellman equation,
the Fokker-Planck equation and the Black-Scholes model. Analytic solutions to
these equations are in general not available and, hence, it is necessary to
design efficient numerical PDE solvers to approximate their solutions. A
crucial challenge that immediately arises in this context is the _curse of
dimensionality_ , see Bellman (1957, 1961). This refers to the tendency of
numerical methods for solving high-dimensional problems to exhibit a
computational cost or require an amount of data that scales exponentially with
the problem’s dimension.

Recent work has shown that _compressive sensing_ and _Deep Learning (DL)_ are
promising techniques to develop efficient high-dimensional PDE solvers and
lessen the curse. This success is part of a larger research trend in the area
of _scientific machine learning_ (Baker et al., 2019), where state-of-the-art
techniques from machine learning are applied to solve challenging scientific
computing problems, including the numerical solution of PDEs. One of the most
popular recent examples in this area are _Physics-Informed Neural Networks
(PINNs)_ , see Raissi et al. (2019) and earlier studies on the topic by
Lagaris et al. (1998), which recently gained an impressive amount of attention
in the scientific computing community. In particular, DL based methods have
shown great promise for high-dimensional PDEs, see Han et al. (2018) and the
recent review paper by E et al. (2021), and PDEs on domains with complex
geometries, see, e.g., Chen et al. (2022).

Concurrently, recent advancements involving compressive sensing in scientific
computing include the adoption of sparsity-based techniques for function
approximation from random samples, see Rauhut and Ward (2012), whose initial
success was due to their application in the field of _Uncertainty
Quantification (UQ)_ of parametric PDEs by Doostan and Owhadi (2011), see also
Adcock et al. (2022b) for a comprehensive review of the topic. Here we focus
on _Compressive Fourier Collocation_ (_CFC_), a method proposed in Wang and
Brugiapaglia (2024) as an improvement of the compressive spectral collocation
method from Brugiapaglia (2020) and able to lessen the curse of dimensionality
in the number of collocation points. A detailed literature review on
compressive sensing and DL methods for PDEs can be found in §1.2.

Motivated by these recent advances, in this paper we study and compare
numerical solvers for high-dimensional PDEs based on compressive sensing and
DL from both the theoretical and the numerical viewpoint. Our methodological
approach is inspired by the recent paper by Adcock and Dexter (2021), where a
similar practical and theoretical study was made in the context of high-
dimensional function approximation from pointwise samples.

###  1.1 Main contributions

Our main contributions, of both theoretical and computational nature, are
summarized below.

##### i) New convergence theorem for high-dimensional periodic PINNs.

The main theoretical contribution of the paper is a new convergence theorem
for periodic PINNs applied to (possibly high-dimensional) diffusion-reaction
problems. Our theoretical guarantee is a convergence result in the form of
_practical existence theorem_ , a theoretical approach recently developed in
the context of scalar-valued function approximation via DL in Adcock and
Dexter (2021) and Hilbert-valued function approximation relevant to parametric
PDEs and UQ in Adcock et al. (2022a). Our result, stated in Theorem 1 and
proved in Appendix 5, establishes the existence of a class of trainable
periodic PINNs with explicit architecture bounds and shows that networks in
this class can achieve near-optimal approximation rates for PDE solutions that
are sparse with respect to the Fourier basis through training using a number
of samples that is only mildly affected by the curse of dimensionality (i.e.,
that scales logarithmically or, at worst, linearly with the domain’s dimension
$d$).

##### ii) Numerical study of CFC and periodic PINNs in high dimensions.

Our second contribution is of computational nature and it is the
implementation and numerical study of CFC and periodic PINNs for high-
dimensional diffusion-reaction problems on the torus $\mathbb{T}^{d}$, up to
dimension $d=30$. We also compare these two methods by studying their accuracy
as a function of the training set size (i.e., the number of collocation
points). The code to reproduce our experiments can be found in the GitHub
repository https://github.com/WeiqiWangMath/PINN_high_dimensional_PDE.

##### iii) Improved implementation of CFC.

We propose a new variant of CFC, based on _adaptive lower Orthogonal Matching
Pursuit (OMP)_ recovery (see Algorithm 1). This variant improves the accuracy
of the CFC method from Wang and Brugiapaglia (2024) under suitable structural
assumptions on the PDE solution and allows the method to scale in higher
dimensions. This is due to the way it constructs the approximation support set
by adaptively exploring the set of possible candidates to be added in an
iterative procedure, thereby avoiding dealing with a large _a priori_
truncation set.

##### iv) CFC convergence theorem for diffusion-reaction problems.

We show the convergence of CFC for diffusion-reaction problems in Theorem 8.
Although this is mainly an auxiliary step needed to prove our main theoretical
result (Theorem 1) it is an extension of the analysis in Wang and Brugiapaglia
(2024) of independent interest.

###  1.2 Contextualization of our contributions

Numerical methods for PDEs based on compressive sensing and sparse recovery
have been considered in tandem with different discretization techniques. This
includes Galerkin (Jokar et al., 2010), Petrov-Galerkin (Brugiapaglia et al.,
2021c, 2015), Fourier-Galerkin (Gross and Iwen, 2023), isogeometric analysis
(Brugiapaglia et al., 2020), spectral collocation approaches (Brugiapaglia,
2020; Wang and Brugiapaglia, 2024) and methods based on the sparse Fourier
transform (Daubechies et al., 2007). In this paper, we are interested in
methods that can be applied to high-dimensional domains. In particular, we
focus on _Compressive Fourier Collocation_ (_CFC_) (see Wang and Brugiapaglia,
2024), which aims to compute a sparse approximation to the PDE solution with
respect to the Fourier basis from random collocation points. This method will
be reviewed in detail in §3. The Wavelet-Fourier CORSING method (Brugiapaglia
et al., 2021c) can in principle be implemented in domains in arbitrary
dimension, although making practical numerical implementations scale in high
dimensions is nontrivial. The sparse spectral method proposed in Gross and
Iwen (2023), based on Fourier-Galerkin discretization and sublinear time
algorithms, can scale to extremely high-dimensional problems.

The design of numerical PDE solvers based on neural networks dates back to the
1990s (Lagaris et al., 1998; Lee and Kang, 1990). More recently, this field
became extremely popular thanks to the introduction of _Physics-Informed
Neural Networks (PINNs)_ (see Raissi et al., 2019; Karniadakis et al., 2021).
Solvers based on DL have shown great promise specifically in the case of high-
dimensional PDEs. In this direction, approaches proposed in the literature
include the Deep Galerkin Method (Sirignano and Spiliopoulos, 2018), methods
based on the reformulation of high-dimensional PDEs as Stochastic Differential
Equations (SDEs) (Han and Jentzen, 2017; Han et al., 2018), PINNs for high-
dimensional problems (Hu et al., 2023; Zeng et al., 2022), and deep genetic
algorithms (Putri et al., 2024). In this paper, we will consider PINNs
combined with a periodic layer (Dong and Ni, 2021) to solve high-dimensional
PDEs on the $d$-dimensional torus. This approach will be presented in detail
in §2.1.

Other approaches for the numerical solution of PDEs able to scale to
moderately high dimensions include sparse grid methods (see, e.g., Shen and Yu
(2010, 2012)), methods based on tensor-based approximation and low rank
structures (Bachmayr and Dahmen, 2015; Bachmayr et al., 2016; Dahmen et al.,
2016), and sparse grid spectral methods (Kupka, 1997).

Our main theoretical result (Theorem 1) is a convergence theorem for PINNs
over high-dimensional periodic domains. Currently, the convergence analysis of
PINNs is an active research area and several studies have appeared in the
literature. However, the theory is arguably far from being fully developed.
The analysis in Shin et al. (2020) shows the convergence of PINNs for linear
second-order elliptic problems, but it is based on Hölder-type regularization
that is in general not implementable. The study in Shin et al. (2023) provides
_asymptotic_ convergence results (i.e., with training set size $m\to\infty$)
relying on assumptions in terms of Bernstein-type inequalities or Rademacher’s
complexity that could be challenging to verify in practice. The analysis in
Doumèche et al. (2023) shows convergence results for PINNs of asymptotic type
and nonasymptotic convergence rates for the expected squared $L^{2}$-error of
the form $\log^{2}(m)/m^{1/2}$ for Sobolev-type regularized loss functions,
where $m$ is the number of training data points used to collocate the PDE in
the PINN’s loss function. The work (De Ryck and Mishra, 2022) proves
generalization bounds for PINNs for Kolmogorov equations showing that PINNs
can lessen the curse of dimensionality, under the assumption that the
network’s weights are bounded. Compared to these results, our practical
existence theorem (Theorem 1) has the advantage of being of nonasymptotic type
and the corresponding error bound leads to fast convergence rates with respect
to the size of the training set when the PDE solution is sparse or
compressible with respect to the Fourier basis.

Finally, the framework of practical existence theorems developed in Adcock and
Dexter (2021); Adcock et al. (2022a) for scalar- and Hilbert-valued high-
dimensional function approximation was recently extended to Banach-valued
function approximation in Adcock et al. (2023) and reduced-order modelling of
parametric PDEs based on convolutional autoencoders in Franco and Brugiapaglia
(2024). This framework relies on recently proposed convergence results for
compressive sensing-based approximation in high dimensions by Adcock et al.
(2022b). For a review of practical existence theory, see Adcock et al. (2024).
All the practical existence theorems developed so far rely on the emulation of
orthogonal polynomials with neural networks (see Daws and Webster, 2019; De
Ryck et al., 2021; Opschoor et al., 2022) and are only applied to
approximating Hilbert or Banach-valued functions relevant to parametric PDEs
or scalar-valued functions. In this paper, we extend the scope of practical
existence theorems by emulating Fourier basis functions and considering the
case of (non-parametric) PDE solvers.

###  1.3 Outline of the paper

We briefly outline the organization of the paper. We start by illustrating the
model problem (a high-dimensional periodic diffusion-reaction equation),
periodic PINNs and our main theoretical result in §2. Then, we present CFC and
adaptive lower OMP in §3. §4 contains an extensive numerical study of CFC and
periodic PINNs. The proof of our convergence result is presented in §5.
Finally, we draw some conclusions and describe possible directions of future
work in §6.

##  2 Problem setting

##### Notation.

We start by recalling some standard mathematical notation that will be
employed throughout the paper. We denote the $d$-dimensional torus by
$\mathbb{T}^{d}$, with $d\in\mathbb{N}$, where $\mathbb{T}:=[0,1]/\sim$ and
$\sim$ is the equivalence relation on $[0,1]$ defined by $x\sim y$ if and only
if $x-y\in\mathbb{Z}$. $L^{2}(\mathbb{T}^{d})$ denotes the space of square-
integrable functions, equipped with inner product $\langle
v,w\rangle:=\int_{\mathbb{T}^{d}}v(\bm{x})\overline{w(\bm{x})}\,\mathrm{d}\bm{x}$
and norm $\|v\|_{L^{2}}=\langle v,v\rangle^{1/2}$. Moreover,
$L^{\infty}(\mathbb{T}^{d})$ denotes the space of functions such that
$\|v\|_{L^{\infty}}:=\text{ess}\sup_{\bm{x}\in\mathbb{T}^{d}}|v(\bm{x})|<\infty$.
We also consider the Sobolev spaces $H^{k}(\mathbb{T}^{d})$, $k=1,2$ equipped
with norms $\|v\|_{H^{1}}:=(\|v\|_{L^{2}}^{2}+\|\nabla v\|_{L^{2}}^{2})^{1/2}$
and $\|v\|_{H^{2}}:=(\|v\|_{H^{1}}^{2}+\|\nabla^{2}v\|_{L^{2}}^{2})^{1/2}$,
respectively. Here $\nabla$ and $\nabla^{2}$ denote the gradient and the
Hessian operators, respectively. Moreover, $\|\nabla
v\|_{L^{2}}^{2}=\int_{\mathbb{T}^{d}}\|\nabla
v(\bm{x})\|_{2}^{2}\,\mathrm{d}\bm{x}$ and
$\|\nabla^{2}v\|_{L^{2}}^{2}=\int_{\mathbb{T}^{d}}\|\nabla^{2}v(\bm{x})\|_{F}^{2}\,\mathrm{d}\bm{x}$,
where $\|\cdot\|_{2}$ is the discrete $2$-norm and $\|\cdot\|_{F}$ is the
Frobenius norm. We will also consider the Sobolev spaces
$W^{k,\infty}(\mathbb{T}^{d})$, $k=1,2$, equipped with norms
$\|v\|_{W^{1,\infty}}:=\max\\{\|v\|_{L^{\infty}},\|\nabla v\|_{L^{\infty}}\\}$
and
$\|v\|_{W^{2,\infty}}:=\max\\{\|v\|_{W^{1,\infty}},\|\nabla^{2}v\|_{L^{\infty}}\\}$,
respectively. Here $\|\nabla
v\|_{L^{\infty}}=\mathrm{ess}\sup_{\bm{x}\in\mathbb{T}^{d}}\|\nabla
v(\bm{x})\|_{\infty}$ and
$\|\nabla^{2}v\|_{L^{\infty}}=\mathrm{ess}\sup_{\bm{x}\in\mathbb{T}^{d}}\max_{i,j}|(\nabla^{2}v(\bm{x}))_{ij}|$.
$C^{k}(\mathbb{T}^{d})$ denotes the space of $k$-times continuously
differentiable functions over $\mathbb{T}^{d}$. We denote the set of first $n$
positive integers by $[n]:=\\{1,\ldots,n\\}$. The cardinality of a set $X$ is
denoted as $|X|$. For a vector $\bm{z}\in\mathbb{C}^{N}$, we define its
support as $\mathrm{supp}(\bm{z})=\\{j\in[N]:z_{j}\neq 0\\}$.

##### Model problem.

Our model problem is a periodic diffusion-reaction equation over the
$d$-dimensional torus $\mathbb{T}^{d}$. In this paper, we are interested in
the scenario where $d\gg 1$. Moreover, we consider the following periodic
diffusion-reaction equation as our model problem:

$-\nabla\cdot(a(\bm{x})\nabla u(\bm{x}))+\rho
u(\bm{x})=f(\bm{x}),\quad\forall\bm{x}\in\mathbb{T}^{d},$ (1)

where $u:\mathbb{T}^{d}\to\mathbb{R}$ is the PDE solution, and the diffusion
coefficient $a:\mathbb{T}^{d}\to\mathbb{R}$, the reaction term
$\rho\in\mathbb{R}$ and the forcing term $f:\mathbb{T}^{d}\to\mathbb{R}$ are
assumed to satisfy

$a\in C^{1}(\mathbb{T}^{d}),\quad\min_{\bm{x}\in\mathbb{T}^{d}}a(\bm{x})\geq
a_{\min}>0,\quad\rho>0,\quad\text{and}\quad f\in L^{2}(\mathbb{T}^{d}).$ (2)

These conditions are sufficient for the problem (1) to be well-posed (see,
e.g., (Brugiapaglia et al., 2021c, Proposition 2.1)) and guarantee that its
weak solutions belong to the Sobolev space $H^{2}(\mathbb{T}^{d})$ (see, e.g.,
(Evans, 2010, §6.3)). Throughout the paper, we also assume $f$ to be regular
enough for pointwise evaluations to be well-defined. Moreover, we define the
PDE operator

$\mathscr{L}[u]:=-\nabla\cdot(a(\bm{x})\nabla u(\bm{x}))+\rho
u(\bm{x}),\quad\forall u\in H^{2}(\mathbb{T}^{d}).$ (3)

Despite its simplicity, (1) is an interesting model problem since it shares
the same second-order diffusion term with more complex PDEs such as the Black-
Scholes, Schrödinger and Fokker-Planck models.

###  2.1 Physics-Informed Neural Networks (PINNs) with periodic layer

We start by illustrating the PINN setting adopted to solve the periodic high-
dimensional diffusion-reaction problem (1). First, we describe the Deep Neural
Network (DNN) architecture employed in the method, then illustrate the
training strategy. The framework presented here will encompass both our main
theoretical result (presented in §2.2) and the numerical experiments in §4.

To enforce periodic boundary conditions, we consider an approach proposed in
Dong and Ni (2021). This is achieved by adding a $C^{\infty}$ _periodic layer_
to the DNN as the first layer. The periodic layer contains, in turn, two
layers of width $dl$, for some $l\in\mathbb{N}$, denoted as
$\bm{q}^{(1)}=(q_{ij}^{(1)})_{i\in[d],j\in[l]}$ and
$\bm{v}^{(2)}=(v_{ij}^{(2)})_{i\in[d],j\in[l]}$, respectively. The neuron
$q_{ij}^{(1)}$ operates on the $i$th component of the input vector
$\bm{x}\in\mathbb{R}^{d}$ and applies a cosine transformation to enforce
periodicity while adding a phase shift parameter $\phi_{ij}$:

$q_{ij}^{(1)}(\bm{x})=\cos(2\pi
x_{i}+\phi_{ij}),\quad\forall\bm{x}\in\mathbb{R}^{d},\quad\forall
i\in[d],\forall j\in[l],$ (4)

where $\bm{\phi}=(\phi_{ij})_{i\in[d],j\in[l]}\in\mathbb{R}^{dl}$ are
trainable parameters and $l$ controls the number of neurons per dimension.
Then, the neurons in $\bm{v}^{(2)}$ collect the outputs of the corresponding
nodes in $\bm{q}^{(1)}$, and apply an affine transformation and (possibly
nonlinear) activation function
$\bm{\sigma}^{(2)}:\mathbb{R}^{dl}\to\mathbb{R}^{dl}$, i.e.,

$\bm{v}^{(2)}(\bm{x})=\bm{\sigma}^{(2)}\left(\text{diag}(\bm{w}^{(2)})\bm{x}+\bm{b}^{(2)}\right),\qquad\forall\bm{x}\in\mathbb{R}^{dl},$
(5)

where $\bm{w}^{(2)}$ and $\bm{b}^{(2)}$ are $dl$-dimensional vectors of
trainable parameters and $\text{diag}(\bm{w}^{(2)})$ is a $dl\times dl$ matrix
with $\bm{w}^{(2)}$ on the main diagonal and zeros elsewhere. After the
periodic layer, the DNN has $h\in\mathbb{N}$ traditional hidden layers having
width $w\in\mathbb{N}$. Each of these takes the output of the previous layer
as input and applies a trainable affine transformation and a componentwise
activation function $\sigma$ to it. These hidden layers are denoted by

$\bm{v}^{(k)}(\bm{x})=\bm{\sigma}^{(k)}(W^{(k)}\bm{x}+\bm{b}^{(k)}),\quad\forall\bm{x}\in\begin{cases}\mathbb{R}^{dl}&\text{if
}k=3\\\ \mathbb{R}^{w}&\text{otherwise}\end{cases},\quad\forall
k=3,\ldots,h+2,$ (6)

where the weight matrix $W^{(k)}$ is $w\times dl$ if $k=3$ and $w\times w$
otherwise, $\bm{b}^{(k)}\in\mathbb{R}^{w}$ and
$\bm{\sigma}^{(k)}:\mathbb{R}^{w}\to\mathbb{R}^{w}$. §4.2.3 discusses the
selection of the hyper-parameters $l$ (number of neurons per dimension in the
periodic layer), $h$ and $w$ defining the DNN architecture. The last hidden
layer of the network $\bm{v}^{(h+2)}$ activates linearly into one output
neuron so that the output of the network is scalar-valued. That is,

$v^{(h+3)}(\bm{x})=W^{(h+3)}\bm{x},\quad\forall\bm{x}\in\mathbb{R}^{w},$

where $W^{(h+3)}$ is $1\times w$. In this paper, we consider componentwise
activations $\bm{\sigma}^{(k)}$ of the form
$\bm{\sigma}^{(k)}(\bm{x})=(\sigma^{(k)}_{j}(x_{j}))$, where
$\sigma^{(k)}_{j}$ could be either a linear activation (i.e.,
$\sigma^{(k)}_{j}(x)=x$) or a nonlinear activation of the following two types:
the hyperbolic tangent (i.e., $\sigma^{(k)}_{j}(x)=\tanh(x)$) or the
_Rectified Power Unit (RePU)_ , defined by

$\textnormal{RePU}_{\ell}(x):=\max\\{0,x^{\ell}\\},\quad\ell\in\mathbb{N}.$
(7)

Note that when $\ell=1$ the function $\textnormal{RePU}_{\ell}$ is the
Rectified Linear Unit (ReLU). In our numerical experiments (see §4) we will
consider tanh activations. However, RePU and linear activations will be used
to derive our theoretical result (see §2.2).

In summary, we consider DNNs $\psi:\mathbb{R}^{d}\to\mathbb{R}$ of the form

$\psi=v^{(h+3)}\circ\bm{v}^{(h+2)}\cdots\circ\bm{v}^{(2)}\circ\bm{q}^{(1)}.$

Fig.1 depicts the architecture of a DNN with periodic layer in dimension
$d=2$.

Figure 1: Architecture of the neural network with the periodic layer ($d=2$).
The number in superscript represents the layer number. In this case, $h=3$ and
depth$(\psi)=7$.

We also define the depth and width of the network as

$\text{width}(\psi)=\max\\{dl,w\\}\quad\text{and}\quad\text{depth}(\psi)=h+4,$
(8)

where we included the input and the output layer in the depth count.

The DNN $\psi(\bm{x})$ is then trained to approximate the solution $u(\bm{x})$
of the high-dimensional PDE. Given collocation points drawn independently and
uniformly at random from $\mathbb{T}^{d}$, i.e.,

$\bm{x}_{1},\ldots,\bm{x}_{m}\stackrel{{\scriptstyle\text{i.i.d.}}}{{\sim}}\text{Uniform}(\mathbb{T}^{d}),$
(9)

the parameters $\bm{\phi}$, $\bm{w}^{(2)}$, $W^{(k)}$, $\bm{b}^{(k)}$ of the
DNN are learned by minimizing a regularized _Root Mean Squared Error (RMSE)_
loss, i.e.,

$\min_{\psi}\sqrt{\frac{1}{m}\sum_{i=1}^{m}|\mathscr{L}[\psi](\bm{x}_{i})-f(\bm{x}_{i})|^{2}}+\lambda\mathcal{R}(\psi),$
(10)

where $\mathscr{L}$ is the PDE operator defined in (3), $\lambda\geq 0$ is a
tuning parameter, and $\mathcal{R}$ is a regularization term that usually
involves the $\ell^{2}$\- or $\ell^{1}$-norm of the networks’ weights
(corresponding to the _weight decay_ or _sparse regularization_ strategies,
respectively). Minimizing the loss function (10) corresponds to finding the
network $\psi$ that minimizes the PDE residual at the collocation points in a
(regularized) least-squares sense. In the general PINN setting for solving
stationary PDEs, the loss function usually consists of the sum of two
components: the first one aims at minimizing the PDE residual (like in (10))
and the second component enforces boundary conditions. However, in our setting
the periodic layer forces $\psi$ to be periodic, hence automatically enforcing
boundary conditions. For this reason, the training loss does not contain a
boundary condition term. The loss is then minimized by stochastic gradient
descent methods. More technical details on the training procedure can be found
in §4.2. We conclude by noting that other losses different from the RMSE can
be considered for the PDE residual minimization. These include losses based
on, e.g., $L^{p}$-norms (Wang et al., 2022) or Sobolev norms (Son et al.,
2021). Here we limit our attention to the regularized RMSE loss for the sake
of simplicity and also because our main theoretical result holds for this
loss. The presence of the regularization term $\mathcal{R}$ will be crucial
for our convergence theorem, illustrated in the next subsection. In our
numerical results we will train by simply minimizing the unregularized RMSE
(or, equivalently, unregularized MSE) loss.

###  2.2 A practical existence theorem for periodic PINNs

Before presenting our main result, namely a convergence theorem for periodic
PINNs based on the framework of practical existence theory (Theorem 1), we
need to introduce some definitions and further technical ingredients. In a
nutshell, our main result shows that trained periodic PINNs are able to
achieve the same accuracy as a sparse Fourier approximation of the PDE
solution using a training set of collocation points whose size scales
logarithmically or, at worst, linearly with the dimension $d$. This is a
natural choice given the presence of periodic boundary conditions in (1). In
addition, we will require some technical conditions on the PDE coefficients
$a$ and $\rho$.

##### Target accuracy: sparse Fourier approximation.

The elements of the $L^{2}$-orthonormal Fourier basis are defined as

$F_{\bm{\nu}}(\bm{x})=\exp(2\pi\mathrm{i}\,\bm{\nu}\cdot\bm{x}),\quad\forall\bm{\nu}\in\mathbb{Z}^{d},\;\forall\bm{x}\in\mathbb{T}^{d}.$
(11)

In particular, we will focus on _sparse_ Fourier approximations supported on
_hyperbolic crosses_ (see, e.g., Dũng et al. (2018); Temlyakov (2018) and
references therein). The _hyperbolic cross_ of order $n$ is a multi-index set
of $\mathbb{Z}^{d}$ defined as

$\Lambda^{\mathrm{HC}}_{d,n}=\left\\{\bm{\nu}\in\mathbb{Z}^{d}:\prod_{k=1}^{d}(|\nu_{k}|+1)\leq
n\right\\}.$ (12)

The hyperbolic cross is a convenient choice in high-dimensional approximation
since its cardinality grows moderately with respect to $d$ and $n$ when
compared to other standard multi-index set choices such as the _tensor
product_ and the _total degree_ sets, see, e.g., (Adcock et al., 2022b, §2.3).

In order to leverage the CFC convergence theory needed for our practical
existence theorem, we consider a rescaled version of the system
$\\{F_{\bm{\nu}}\\}_{\bm{\nu}\in\mathbb{Z}^{d}}$, defined by

$\Psi_{\bm{\nu}}=\frac{1}{4\pi^{2}\|\bm{\nu}\|_{2}^{2}+\rho/a_{\bm{0}}}F_{\bm{\nu}},\quad\forall\bm{\nu}\in\mathbb{Z}^{d},$
(13)

where $a_{\bm{0}}=\langle
a,F_{\bm{0}}\rangle=\int_{\mathbb{T}^{d}}a(\bm{x})\,\mathrm{d}x$. As explained
in detail in §5.1, this rescaling ensures that
$\\{\mathscr{L}[\Psi_{\bm{\nu}}]\\}_{\bm{\nu}\in\mathbb{Z}^{d}}$ is a _bounded
Riesz system_ (Brugiapaglia et al., 2021a) of $L^{2}(\mathbb{T}^{d})$ under
sufficient conditions on $a$ and $\rho$, see (15) and (16) below, which is a
crucial property needed to show the convergence of the CFC method. To gain
some intuition about this rescaling, consider the simple case of constant
diffusion $a\equiv 1$. Applying the PDE operator to the rescaled Fourier
system yields
$\mathscr{L}[\Psi_{\bm{\nu}}]=(-\Delta+\rho)\Psi_{\bm{\nu}}=F_{\bm{\nu}}$ for
all $\bm{\nu}\in\mathbb{Z}^{d}$. In this simple scenario,
$\\{\mathscr{L}[\Psi_{\bm{\nu}}]\\}_{\bm{\nu}\in\mathbb{Z}^{d}}$ is a bounded
orthonormal system and, as such, is ideally suited for compressive sensing,
(see Foucart and Rauhut, 2013, §12).

In this setting, we consider a finite-dimensional truncation of the solution
$u$ to a finite multi-index set $\Lambda\subseteq\mathbb{Z}^{d}$ (expressed
with respect to the rescaled Fourier basis), i.e.,

$u_{\Lambda}=\sum_{\bm{\nu}\in\Lambda}c_{\bm{\nu}}\Psi_{\bm{\nu}},\quad\text{with
}c_{\bm{\nu}}=(4\pi^{2}\|\bm{\nu}\|_{2}^{2}+\rho/a_{\bm{0}})\cdot\langle
u,F_{\bm{\nu}}\rangle,$ (14)

and let
$\bm{c}_{\Lambda}=(c_{\bm{\nu}})_{\bm{\nu}\in\Lambda}\in\mathbb{C}^{N}$, where
$N=|\Lambda|$. An _$s$ -sparse_ approximation to $u$ is then obtained by only
keeping $s$ terms in the expansion (14). In general, we recall that a vector
is said to be $s$-sparse if it has at most $s$ nonzero entries. The best
possible accuracy of such an approximation is measured by the _best $s$-term
approximation error_, see (Cohen et al., 2009) and references therein, defined
as

$\sigma_{s}(\bm{c}_{\Lambda})_{p}=\min_{\bm{z}\in\mathbb{C}^{N}}\left\\{\|\bm{c}_{\Lambda}-\bm{z}\|_{p}:\bm{z}\text{
is $s$-sparse}\right\\}.$

##### Sufficient conditions on the PDE coefficients.

We also introduce a technical condition on the PDE coefficients $a$ and $\rho$
necessary for the convergence result. We consider diffusion coefficients
$a:\mathbb{T}^{d}\to\mathbb{R}$ having a sparse Fourier expansion, namely,

$\displaystyle a=a_{\bm{0}}+\sum_{\bm{\nu}\in
T}a_{\bm{\nu}}F_{\bm{\nu}},\quad\text{with }a_{\bm{\nu}}=\langle
a,F_{\bm{\nu}}\rangle,\quad\forall\bm{\nu}\in T\cup\\{\bm{0}\\},$ (15)

for some $T\subseteq\mathbb{Z}^{d}\setminus\\{\bm{0}\\}$. In our main result,
we will assume that $a$ and $\rho$ satisfy the following assumption:

$\displaystyle\sqrt{|T|}\cdot\|a-a_{\bm{0}}\|_{H^{1}}<\sqrt{\left(a_{\bm{0}}+\frac{\rho}{4\pi^{2}}\right)^{2}+a_{\bm{0}}^{2}}-\left(a_{\bm{0}}+\frac{\rho}{4\pi^{2}}\right).$
(16)

Note that in the case of constant diffusion (i.e., $a\equiv a_{\bm{0}}$) the
condition above is always satisfied under (2). In general, condition (16)
controls the oscillatory behaviour of $a$. We emphasize that the sparsity of
$a$ and condition (16) are sufficient for our convergence theorem to hold, but
we do not claim (nor believe) they are necessary.

We are now in a position to state our main result, which provides the
existence of a class of trained neural networks with architecture of the form
described in §2.1 able to approximate the solution to (1) with accuracy
comparable to that of an $s$-sparse Fourier approximation with high
probability.

######  Theorem 1 (Practical existence theorem for high-dimensional periodic
PINNs)

Given a dimension $d\in\mathbb{N}$, target sparsity $s\in\mathbb{N}$,
hyperbolic cross order $n\in\mathbb{N}$, RePU power $\ell\in\mathbb{N}$, with
$\ell\geq 2$, and probability of failure $\varepsilon\in(0,1)$, there exist:

  1. (i)

a class of neural networks $\mathcal{N}$ of the form described in §2.1 with
$d$-dimensional input and $1$-dimensional output layers,
$\textnormal{RePU}_{\ell}$ or linear activations, complex-valued weights and
biases, and such that, for all $\psi\in\mathcal{N}$,

$\displaystyle\textnormal{width}(\psi)$ $\displaystyle\leq
c_{\ell}^{(1)}\cdot\min\left\\{4n^{5}16^{d},e^{2}n^{2+\log_{2}d}\right\\}\cdot
d\cdot\min\\{2^{d},n\\},$ (17) $\displaystyle\textnormal{depth}(\psi)$
$\displaystyle\leq
c^{(2)}\cdot\left(\log_{2}(n)+\min\\{\log_{2}d,n\\}\right),$ (18)

  2. (ii)

a regularization function $\mathcal{R}:\mathcal{N}\to[0,\infty)$,

  3. (iii)

a choice of tuning parameter $\lambda$ depending only on $a$, $\rho$ and $s$,

such that the following holds with probability $1-\varepsilon$. For all
$a:\mathbb{T}^{d}\to\mathbb{R}$ and $\rho\in\mathbb{R}$ satisfying (2), (15)
and (16), let

$m\geq c^{(3)}_{a,\rho}\cdot s\cdot\log^{2}\left(c^{(4)}_{a,\rho}\cdot
s\right)\cdot\left(\min\\{\log(n)+d,\log(2n)\log(2d)\\}+\log(\varepsilon^{-1})\right),$
(19)

$\Lambda=\Lambda^{\mathrm{HC}}_{d,n}$ as in (12), and consider collocation
points $\bm{x}_{1},\ldots\bm{x}_{m}$ randomly and independently sampled from
the uniform measure on $\mathbb{T}^{d}$. Then, every minimizer $\hat{\psi}$ of
the training program

$\min_{\psi\in\mathcal{N}}\sqrt{\frac{1}{m}\sum_{i=1}^{m}\left|\mathscr{L}[\psi](\bm{x}_{i})-f(\bm{x}_{i})\right|^{2}}+\lambda\mathcal{R}(\psi)$
(20)

satisfies

$\displaystyle\|u-\hat{\psi}\|_{L^{2}}+\|(\Delta-\rho)(u-\hat{\psi})\|_{L^{2}}$
$\displaystyle\leq
C^{(1)}_{a,\rho}\cdot\frac{\sigma_{s}(\bm{c}_{\Lambda})_{1}}{\sqrt{s}}+C^{(2)}_{a,d,\rho}\cdot\left(\frac{\|u-u_{\Lambda}\|_{W^{2,\infty}}}{\sqrt{s}}+\|u-u_{\Lambda}\|_{H^{2}}\right).$
(21)

Moreover, if $\rho<1$, we also have

$\displaystyle\|u-\hat{\psi}\|_{H^{2}}\leq
C^{(3)}_{a,\rho}\cdot\frac{\sigma_{s}(\bm{c}_{\Lambda})_{1}}{\sqrt{s}}+C^{(4)}_{a,d,\rho}\cdot\left(\frac{\|u-u_{\Lambda}\|_{W^{2,\infty}}}{\sqrt{s}}+\|u-u_{\Lambda}\|_{H^{2}}\right).$
(22)

Here each constant depends only on the subscripted parameters. Moreover, the
dependence of each constant on $d$ is at most linear (when present).

Proof sketch. A complete proof of Theorem 1 is given in §5.2. It leverages the
convergence theory of CFC (presented and proved subsequently in §5.1) thanks
to an _ad hoc_ construction of the network class $\mathcal{N}$. The idea is to
construct networks $\psi\in\mathcal{N}$ so as to explicitly replicate linear
combinations of Fourier functions $\\{F_{\bm{\nu}}\\}_{\bm{\nu}\in\Lambda}$
supported on the hyperbolic cross $\Lambda$, recall (14). In this
construction, only the last layer, corresponding to the coefficients of the
linear combination is trained. The rest of the network is explicitly
constructed to replicate (a suitably rescaled version of) the Fourier basis
functions. In this setting, minimizing the regularized loss (20) where
$\mathcal{R}(\psi)$ is the $1$-norm of the weights in the last layer of the
network $\psi$ is equivalent to solving a sparse regularization problem,
specifically, a square-root LASSO problem (Adcock et al., 2019; Belloni et
al., 2011), which allows to rigorously connect periodic PINNs’ training with
the CFC convergence theory (see Theorem 8).

We conclude this section by highlighting some important features of Theorem 1.

  1. (i)

Theorem 1 is called a _practical_ existence theorem since, as opposed to
standard neural network existence results such as universal approximation
theorems, see, e.g., (Elbrächter et al., 2021) and references therein, it not
only guarantees the existence of neural networks (in this case, periodic
PINNs) with favorable approximation properties, but also establishes that such
networks can be computed by training a regularized RMSE loss and provides a
sufficient condition on the minimum number of samples (i.e., collocation
points) for the training process to be successful.

  2. (ii)

One of the key benefits of Theorem 1 is that the minimum number of training
samples $m$ needed to successfully train a periodic PINN from the class
$\mathcal{N}$ to solve a $d$-dimensional reaction-diffusion problem scales
logarithmically (when $d\gg n$) or, at worst, linearly (when $n\gg d$) in $d$.
This indicates that periodic PINNs are provably able to alleviate the curse of
dimensionality.

  3. (iii)

As established by the error bound (21), the periodic PINN approximation
accuracy guaranteed by Theorem 1 is controlled by the best $s$-term
approximation error $\sigma_{s}(\bm{c}_{\Lambda})_{1}$ and by the truncation
error $u-u_{\Lambda}$ (measured with respect to the $W^{2,\infty}$\- and the
$H^{2}$-norm). This accuracy is inherited by the CFC convergence result that
Theorem 1’s proof relies on (see Theorem 8). In this paper we do not assume
that $u$ belongs to a specific function class, but we note that the best
$s$-term approximation and the hyperbolic cross truncation error could be
bounded for functions satisfying suitable mixed regularity conditions, see
(Dũng et al., 2018; Temlyakov, 2018) and the discussion in (Wang et al., 2022,
§2.3). In general, $\sigma_{s}(\bm{c}_{\Lambda})_{1}$ can be estimated via
Stechkin’s inequality, see, e.g., (Adcock et al., 2022b, Lemma 3.5).

  4. (iv)

Theorem 1 also provides explicit bounds on the architecture of the periodic
PINNs from class $\mathcal{N}$. The networks’ depth scales logarithmically in
the dimension $d$ and the networks’ width scales polynomially in $d$ (note
that $n^{\log_{2}d}=d^{\log_{2}n}$), for sufficiently large $d$.

The main limitations of Theorem 1 and related open problems are discussed in
§6.

##  3 Compressive Fourier Collocation (CFC)

In this section, we illustrate the CFC method and its efficient numerical
implementation via adaptive lower OMP. For more details, we refer to Wang and
Brugiapaglia (2024). Similarly to §2.2, we consider a finite multi-index set
$\Lambda\subset\mathbb{Z}^{d}$, the rescaled Fourier basis
$\\{\Psi_{\bm{\nu}}\\}_{\bm{\nu}\in\mathbb{Z}^{d}}$ defined in (13) and a
finite-dimensional expansion $u_{\Lambda}$ of the form (14). Then, similarly
to the PINN approach, we collocate the diffusion-reaction equation (1) by
means of Monte Carlo sampling. Hence, we randomly generate $m$ i.i.d. uniform
points $\bm{x}_{1},\ldots,\bm{x}_{m}\in\mathbb{T}^{d}$, as in (9). Letting
$N=|\Lambda|$, we assume to have an ordering for the multi-indices in
$\Lambda=\\{\bm{\nu}_{1},\ldots,\bm{\nu}_{N}\\}$ (e.g., the lexicographic
ordering). The PDE collocation process leads to the linear system

$A\bm{z}=\bm{b},$ (23)

where $A\in\mathbb{C}^{m\times N}$ and $\bm{b}\in\mathbb{C}^{m}$ are defined
by

$A_{ij}=\frac{1}{\sqrt{m}}[-\nabla\cdot(a\nabla\Psi_{\bm{\nu}_{j}})+\rho\Psi_{\bm{\nu}_{j}}](\bm{x}_{i})\quad\text{and}\quad
b_{i}=\frac{1}{\sqrt{m}}f(\bm{x}_{i}),\quad\forall i\in[m],j\in[N].$ (24)

We refer to $A$ as the CFC matrix. This collocation method is _compressive_
since we choose $m\ll N$. This makes the linear system (23) underdetermined.

##### Computing the CFC solution.

Following the approach in Wang and Brugiapaglia (2024), a CFC approximation
$\hat{u}$ to $u$ can be computed by approximately solving the underdetermined
linear system (23) via sparse recovery techniques. Therefore, we need to
determine (i) a suitable truncation multi-index set $\Lambda$ and (ii) a
sparse recovery method to approximately solve the linear system (23). In this
paper, the truncation set $\Lambda$ is chosen as a _hyperbolic cross_ of order
$n$, i.e., $\Lambda=\Lambda^{\mathrm{HC}}_{d,n}$, recall (12). The hyperbolic
cross offers a twofold advantage. First, as discussed in §2.2, the cardinality
of $\Lambda^{\mathrm{HC}}_{d,n}$ grows moderately in $n$ and $d$ when compared
to other standard multi-index families (e.g., tensor product and total degree
set). Second, a hyperbolic cross can be characterized as the union of _lower_
(or, equivalently, _downward closed_ or _monotone_) sets of a given
cardinality (_cp_. Definition 2).

After fixing $\Lambda$, we compute an approximate solution
$\hat{\bm{c}}=(\hat{c}_{\bm{\nu}})_{\bm{\nu}\in\Lambda}\in\mathbb{C}^{N}$ to
the underdetermined linear system (23). To do so, we employ sparse recovery
techniques, such as _Orthogonal Matching Pursuit_ (_OMP_) or $\ell^{1}$
minimization (Foucart and Rauhut, 2013). Finally, we define the CFC
approximation as

$\hat{u}=\sum_{\bm{\nu}\in\Lambda}\hat{c}_{\bm{\nu}}\Psi_{\bm{\nu}}.$ (25)

In Wang and Brugiapaglia (2024) it was shown that $\hat{u}$ is an accurate and
stable approximation to $u$ for high-dimensional diffusion equations, under
sufficient conditions on the diffusion term and for a number of collocation
points that scales only logarithmically with the ambient dimension $d$. In
Theorem 8, we will extend the CFC convergence analysis from Wang and
Brugiapaglia (2024) to diffusion-reaction problems. From the computational
viewpoint, if one wants to compute a very sparse CFC approximation to the
solution $u$, then OMP typically offers a faster reconstruction than solving
an $\ell^{1}$ minimization problem via a convex optimization solver.

###  3.1 Adaptive lower Orthogonal Matching Pursuit (OMP)

For very high-dimensional domains (say, $d>20$), the cardinality of the
hyperbolic cross, despite being moderate with respect to other multi-index
choices, becomes considerably large. This can make storing $A$ and computing
$\hat{\bm{c}}$ genuinely challenging. To deal with higher dimensions, we need
a more efficient recovery procedure that does not rely on choosing a large _a
priori_ truncation set $\Lambda$, but constructs it iteratively. This can be
achieved by considering a reconstruction strategy called _adaptive lower
Orthogonal Matching Pursuit (OMP)_ , that we now illustrate. The adaptive
lower OMP algorithm presented here is based on analogous techniques employed
in adaptive high-dimensional approximation, in particular in sparse grids
(Gerstner and Griebel, 2003) and least squares methods (Migliorati, 2014,
2019).

First, we define the notion of lower set. Lower sets are an important class of
multi-index sets in approximation theory and we refer to, e.g., Cohen and
Migliorati (2018) or (Adcock et al., 2022b, §1.5) and references therein for
further reading. Note that for polynomial approximations, lower sets are
typically defined in $\mathbb{N}_{0}^{d}$. However, since we employ the
complex Fourier basis here we consider lower sets in $\mathbb{Z}^{d}$. Before
defining lower sets, we introduce a convenient notation to compare multi-
indices of $\mathbb{Z}^{d}$. For $\bm{\mu},\,\bm{\nu}\in\mathbb{Z}^{d}$, we
use $\bm{\mu}\preceq\bm{\nu}$ to indicate that $|\mu_{i}|\leq|\nu_{i}|$ for
every $i\in[d]$. Moreover, if $\bm{\mu}\preceq\bm{\nu}$ and there exists an
$i\in[d]$ such that $|\mu_{i}|<|\nu_{i}|$, then we say that
$\bm{\mu}\prec\bm{\nu}$.

######  Definition 2 (Lower set of $\mathbb{Z}^{d}$)

A multi-index set $\Lambda\subseteq\mathbb{Z}^{d}$ is said to be _lower_ if
the following holds for all $\bm{\nu}\in\Lambda$ and
$\bm{\mu}\in\mathbb{Z}^{d}$: if $\bm{\mu}\preceq\bm{\nu}$, then
$\bm{\mu}\in\Lambda$.

An example of a lower set is given in Fig. 2 (blue dots). An equivalent
condition for $\Lambda\subseteq\mathbb{Z}^{d}$ to be a lower set is the
following: if $\bm{\nu}\in\Lambda$, then the box
$\prod_{k=1}^{d}[-|\nu_{k}|,|\nu_{k}|]\subseteq\Lambda$. Note that lower sets
of $\mathbb{Z}^{d}$ are symmetric with respect to all coordinate hyperplanes
$\\{\bm{\nu}\in\mathbb{Z}^{d}:\nu_{k}=0\\}$, with $k\in[d]$. This symmetry can
be justified as follows. For real-valued solutions
$u:\mathbb{T}^{d}\to\mathbb{R}$ one could employ a real Fourier expansion with
respect to basis functions of the form
$\prod_{k=1}^{d}\xi_{k}(2\pi\nu_{k}x_{k})$, where $\xi_{k}\in\\{\sin,\cos\\}$
and $\bm{\nu}\in\mathbb{N}_{0}^{d}$. It is a simple exercise to verify that
the complex Fourier expansion of each of these real-valued trigonometric
functions is supported on a set symmetric with respect to every coordinate
hyperplane.

In order to create a mechanism able to iteratively enlarge a multi-index set
while preserving the lower set structure, we introduce the notion of reduced
margin.

######  Definition 3 (Reduced margin)

The _reduced margin_ of a (nonempty) lower set
$\Lambda\subseteq\mathbb{Z}^{d}$ is the multi-index set
$\mathcal{R}(\Lambda):=\\{\bm{\nu}\in\mathbb{Z}^{d}:\,\bm{\nu}\notin\Lambda\text{
and }\forall\bm{\mu}\prec\bm{\nu},\,\bm{\mu}\in\Lambda\\}.$ Moreover, we let
$\mathcal{R}(\emptyset):=\\{\bm{0}\\}$.

We provide an illustration of the reduced margin in Fig. 2 (red dots). A key
property of the reduced margin is that if $\Lambda\subseteq\mathbb{Z}^{d}$ is
a lower set and $\bm{\nu}\in\mathcal{R}(\Lambda)$ then
$\Lambda\cup\\{\bm{\mu}\in\mathbb{Z}^{d}:|\mu_{i}|=|\nu_{i}|,\;\forall
i\in[d]\\}$ is also a lower set.

The adaptive lower OMP algorithm is a structured variant of OMP, see, e.g.,
(Foucart and Rauhut, 2013, §3.2) and references therein, with the only
difference that the greedy search is restricted to iteratively enlarged lower
sets. The algorithm generates a nested sequence of lower sets
$\emptyset=\Lambda^{(0)}\subseteq\Lambda^{(1)}\subseteq\ldots\subseteq\Lambda^{(K)}$
and then computes a sequence of least-squares solutions $\bm{z}^{(n)}$ to the
linear system (23) such that
$\mathrm{supp}(\bm{z}^{(n)})\subseteq\Lambda^{(n)}$. At each iteration,
adaptive lower OMP extends the existing multi-index set $\Lambda_{n}$ by
picking an element $\bm{\nu}^{(n)}\in\mathcal{R}(\Lambda^{(n)})$ corresponding
to the largest absolute residual
$\left|(A^{*}(\bm{b}-A\bm{z}^{(n)}))_{\bm{\nu}}\right|$. Then, it also adds
all reflections of $\bm{\nu}^{(n)}$ with respect to the coordinate hyperplanes
$\\{\bm{\nu}\in\mathbb{Z}^{d}:\nu_{k}=0\\}$ to preserve the lower structure.
The greedy selection criterion of adaptive lower OMP (Line 5 in Algorithm 1)
is illustrated in Fig. 2.

Figure 2:  Illustration of the greedy selection criterion of adaptive lower
OMP (Line 5 in Algorithm 1). Considering a given lower set
$\Lambda^{(n)}\subset\mathbb{Z}^{2}$ (blue), its reduced margin
$\mathcal{R}(\Lambda^{(n)})$ is drawn in red. If the quantity
$|A^{*}(\bm{b}-A\bm{z}^{(n)})_{\bm{\nu}}|$ is minimized at $\bm{\nu}=(3,1)$,
then $\Lambda^{(n+1)}$ is constructed by adding the red circled dots to
$\Lambda^{(n)}$.

The method is summarized in Algorithm 1.

Algorithm 1 Adaptive lower Orthogonal Matching Pursuit (OMP)

0: $A\in\mathbb{C}^{m\times N}$, $\bm{b}\in\mathbb{C}^{m}$, and number of
iterations $K\in\mathbb{N}$

0: A $K$-sparse vector $\hat{\bm{c}}\in\mathbb{C}^{N}$, approximately solving
(23)

1: $\Lambda^{(0)}=\emptyset,\ \bm{z}^{(0)}=\bm{0},$

2: for $n=0,\dots,K-1$ do

3: Compute the columns of $A$ corresponding to $\mathcal{R}(\Lambda^{(n)})$
(see Definition 3)

4: $d_{\bm{\nu}}\leftarrow 1/\sqrt{\sum^{m}_{i=1}|A_{i\bm{\nu}}|^{2}}$,
$A_{i\bm{\nu}}\leftarrow A_{i\bm{\nu}}d_{\bm{\nu}}$, for all
$\bm{\nu}\in\mathcal{R}(\Lambda^{(n)})$ ($\ell^{2}$-normalize the new columns
of $A$)

5:
$\bm{\nu}^{(n)}\leftarrow\displaystyle\operatorname*{arg\,max}_{\bm{\nu}\in\mathcal{R}(\Lambda^{(n)})}\left|(A^{*}(\bm{b}-A\bm{z}^{(n)}))_{\bm{\nu}}\right|$
(greedy multi-index selection)

6: $\Lambda^{(n+1)}\leftarrow\
\Lambda^{(n)}\cup\\{\bm{\mu}\in\mathbb{Z}^{d}:|\mu_{i}|=|\nu_{i}^{(n)}|,\,\forall
i\in[d]\\}$

7:
$\bm{z}^{(n+1)}\leftarrow\displaystyle\arg\min_{\bm{z}\in\mathbb{C}^{N}}\left\\{\left\|\bm{b}-A\bm{z}\right\|_{2}^{2},\,\mbox{supp}(\bm{z})\subseteq\Lambda^{(n+1)}\right\\}$

8: end for

9: $\hat{\bm{c}}\leftarrow D\bm{z}^{(K)}$, with $D=\text{diag}(\bm{d})$

##  4 A numerical study of CFC and periodic PINNs

In this section, we present numerical experiments on the high-dimensional
periodic diffusion-reaction equation (1). After illustrating the numerical
setup in §4.1, we conduct tests on PINNs with periodic layer in §4.2 and the
(adaptive) CFC method in §4.3. Then, we compare the two approaches in §4.4.

###  4.1 Numerical setup

We start by introducing the setup of our numerical experiments.

##### Measurement of errors.

In all numerical experiments, we use the relative $L^{2}$-error to measure the
approximation error. It is defined as

$\text{relative}\
L^{2}\text{-error}=\frac{\|u-\hat{u}\|_{L^{2}}}{\|u\|_{L^{2}}},$

where $\hat{u}$ is the computed approximation to an exact solution $u$.
Computing the $L^{2}$-norms is a challenge in high dimension, therefore the
norms $\|\cdot\|_{L^{2}}$ are approximated using Monte Carlo integration,
i.e.,

$\|u\|_{L^{2}}\approx\sqrt{\frac{1}{M}\sum^{M}_{i=1}|u(\bm{x}_{i})|^{2}},$

where the $\bm{x}_{1},\ldots,\bm{x}_{M}$ are $M\gg 1$ random independent
points uniformly distributed in ${\mathbb{T}^{d}}$. Throughout this section,
we choose $M=10000$.

##### Diffusion coefficient.

As shown in Wang and Brugiapaglia (2024), the choice of the type of diffusion
coefficient (e.g., constant, sparse, or non-sparse) does not impact the
numerical results for the compressive Fourier collocation method solving high-
dimensional PDEs. Thus, in all experiments, we use the following diffusion
coefficient:

$a(\bm{x})=1+0.25\sin{(2\pi x_{1})}\sin{(2\pi x_{2})}.$

##### Exact solutions.

In our experiments, we will consider the problems having the following three
high-dimensional functions as their solution:

$\displaystyle u_{1}(x)$ $\displaystyle=\sin(4\pi x_{1})\sin(2\pi x_{2}),$
(Example 1) (26) $\displaystyle u_{2}(x)$ $\displaystyle=\exp\left(\sin(2\pi
x_{1})+\sin(2\pi x_{2})\right),$ (Example 2) (27) $\displaystyle u_{3}(x)$
$\displaystyle=\exp\left(\sum^{d}_{k=1}\frac{1}{k^{2}}\sin(2\pi
x_{k})\right).$ (Example 3) (28)

In the following, we will use Examples 1, 2, and 3 to refer to diffusion-
reaction equations with exact solutions given by (26), (27), and (28),
respectively. Note that for any given coefficients $a$ and $\rho$, we will
enforce the PDE to have a prescribed exact solution by suitable choice of the
forcing term $f$. Here, $u_{1}$ is defined in terms of only one real-valued
trigonometric function. This corresponds to a 4-sparse solution with respect
to the complex Fourier system. Hence, in this case we have
$\sigma_{s}(\bm{c}_{\Lambda})=0$ and $u-u_{\Lambda}=0$ as soon as the
truncation set $\Lambda$ contains the four multi-indices $(\pm 2,\pm
1,0,0,\dots)$. Note, though, that these four multi-indices do not form a lower
set. The solution $u_{2}$ is a smooth periodic function active only in the
variables $x_{1}$ and $x_{2}$. Hence, it exhibits a highly anisotropic
behaviour. Finally, $u_{3}$ is another smooth solution defined using all
variables from $x_{1},\ldots,x_{d}$. The behaviour of this solution is also
anisotropic and the solution is uniformly bounded by
$\exp(\sum_{k=1}^{\infty}1/{k^{2}})=\exp(\pi^{2}/6)$.

##### Randomization of experiments and visualization.

Due to the random nature of collocation points in the CFC and the periodic
PINN method, we consider 25 random runs for each test. However, for the
experiments in §4.2.3, where we examine the effect of the architecture on
performance, we use 10 random runs for each width-depth combination. In our
plots, the main curves represent the sample geometric mean of the error, and
the size of the lightly shaded areas corresponds to its corrected geometric
standard deviation. Note that calculating the geometric mean and standard
deviation of the errors is equivalent to compute the classical mean and
standard deviation of the log transformed errors. For more details on the
visualization strategy, we refer to (Adcock et al., 2022b, §A.1.3).

###  4.2 Physics-Informed Neural Networks

We run tests using PINNs to solve the PDE defined in (1) with the various
exact solutions from (26)-(28). For the experiments in §4.2.1 and §4.2.2, we
use DNNs with $l=11$ nodes in the periodic layer, depth $h=3$ hidden layers,
and width to depth ratio $r=10$ so that the networks have $30$ nodes per
hidden layer. Moreover, we consider a tanh activation in every layer. In
§4.2.3, we examine the impact of the architecture on performance by comparing
results obtained by varying the number of nodes on the periodic layer and on
the hidden layers of the fully connected part of the network, see Fig. 1. In
this section, we train periodic PINNs without regularization, i.e., we let
$\lambda=0$ in (10). Moreover, we simply minimize the MSE loss (which is
equivalent to minimizing the RMSE in the non-regularized setting).

####  4.2.1 Performance and the impact of the number of samples

We test the PINNs with the exact solutions defined in (26)-(28) and refer to
them as Example 1, 2, 3, respectively. Fig. 3 presents the performance of the
PINNs over 30000 epochs for each exact solution, with parameters $d=6$,
$\rho=0.5$.

Figure 3: (Impact of the number of samples) Relative $L^{2}$-error versus the
number of epochs for approximating the exact solutions defined in (26)-(28)
with $d=6$, where $m$ is the number of sample points.

The PINNs have 11 nodes on the periodic layer and three non-periodic layers
with 30 nodes each. We choose the commonly used $\tanh$ activation function
and the Adam optimizer (Kingma and Ba, 2017) with early stopping. We use this
neural network to conduct all numerical tests except those in §4.2.3. As shown
in Fig. 1, the PINN model approximates all three exact solutions with
$L^{2}$-error on the order of $10^{-2}$ to $10^{-3}$ within the $30000$ epoch
training budget. Notably, the error decreases rapidly in the first $10000$
epochs and can be seen to saturate for the remaining epochs. Also, in Examples
1 and 3, we see some small improvement in the error as the size of the
training set is increased, while for Example 2, we note that moving from
$m=2000$ to $m=3000$ points does not noticeably improve the error. These
differences in results obtained for our examples illustrate that the number of
samples needed to reach a minimum error varies depending on the complexity of
the exact solution.

####  4.2.2 Impact of the dimension

We now consider the performance of PINNs in higher-dimensional settings. We
present the results for exact solutions with $d=6,10,20$ in Fig. 4. With 3000
sample points, the PINNs approximate Example 1’s exact solution with relative
$L^{2}$-error on the order of $10^{-2}$ when $d=6$.

Figure 4: (Impact of the dimension) Relative $L^{2}$-error versus the number
of epochs for approximating the exact solutions defined in (26)-(28) with (top
left, bottom left, bottom right) $m=3000$ and (top right) $m=10000$ samples,
in $d=6,10,20$ dimensions.

However, we observe that the error increases with increasing dimension with
the error for $d=10$ above that for $d=6$, while for $d=20$ the PINNs do not
converge at all within the budget of epochs. Increasing the number of samples
to 10000 (result in the top right of Fig. 4), the PINNs can be seen to be
converging to the exact solution, saturating at an identical error around
$10^{-2}$ in dimensions $d=6,10$, and $20$ after 30000 epochs. This highlights
the dependence on the problem dimension in the number of samples required to
achieve a given accuracy. This can also be observed for Examples 2 and 3,
where we see that with fixed $m=3000$ samples the error increases as the
dimension increases.

As discussed in the introduction, we seek methods to solve high-dimensional
PDEs which overcome or substantially mitigate the curse of dimensionality. The
experiments in Fig. 5 aim to illuminate the relationship between the dimension
of the problem and the number of samples required to achieve a given error.

Figure 5: (Dimensionality and number of samples) Relative $L^{2}$-error after
30000 epochs versus number of sample points $m$ for the exact solutions
defined in Eq. (26)-(28), where $d=6,10,20$ are the dimension of the problem.
Left: Example 1. Middle: Example 2. Right: Example 3.

There we plot the relative $L^{2}$-error vs. the number of sample points in
$d=6,10$, and $20$ dimensions for each of our examples. We observe in all
three cases that the number of samples required to achieve a given error
increases with the dimension. However, focusing on the point at which the
PINNs begin to saturate in their relative $L^{2}$-error, this scaling does not
appear to be exponential in the dimension which would certainly be the case if
the method was not capable of mitigating the curse. Rather, the scaling in the
required number of samples appears to be approximately linear. Comparing the
results for all three examples when $d=6,10$ and $20$, we observe that the
number of samples required for the PINNs to reach saturation at relative
$L^{2}$-error around $10^{-2}$ roughly doubles moving from $6$ to $10$
dimensions and again moving from $10$ to $20$ dimensions. This empirical
observation is in accordance with Theorem 1, which establishes that $m$ scales
at worst linearly in $d$ (see (19)).

####  4.2.3 Impact of the architecture

We also examine the impact of the neural network’s architecture on
performance. Specifically, we consider varying the number of nodes in the
periodic layer and the width-depth ratio $r=w/h$ of the hidden layers. In Fig.
3 and 4, we see that the neural network approach produces accurate convergence
results in all three examples. We then naturally raise a question: can a
better choice of the hyper-parameters $l$, $h$, and $r=w/h$ further improve
the accuracy of the neural network approach? To answer this question, we
design experiments by changing these hyper-parameters and testing the
performance in Example 3. We change the number of nodes in the periodic layer
of the neural network in the left plot of Fig. 6.

Figure 6: (Impact of the architecture) Left: Relative $L^{2}$-error versus
number of nodes in the periodic layer. Right: Relative $L^{2}$-error versus
the number of hidden layers. Different curves represent the different width-
depth ratios of the hidden layers. For both figures, we test the neural
networks on Example 3, in $d=6$ dimensions using 5000 sample points, and the
relative $L^{2}$-error are measured after 30000 epochs.

The results show that increasing the number of nodes beyond 10 in the first
layer does not improve the $L^{2}$-error. The results for tuning the number of
layers and the number of nodes per layer for non-periodic layers are
illustrated in the right plot of Fig. 6. We observe that adding more hidden
layers to the fully-connected part of the network does not improve the
accuracy of the periodic-NN solution. However, we do observe larger standard
deviation for the choice of $r=3$ for the width-depth ratio in comparison to
the choices of $r=5,10,20$. We also observe diminishing returns in further
increasing the ratio to $r=20$ as the cost of training increases due to the
larger width of the network. Since we observe that the relative $L^{2}$-error
is not sensitive to these hyper-parameters, we choose hyper-parameters as
stated at the beginning of §4.2.

###  4.3 Adaptive lower OMP results

In this subsection, we represent numerical results for the adaptive lower OMP
method (Algorithm 1), including the convergence results and a comparison with
the non-adaptive Fourier compressive collocation method described in Wang and
Brugiapaglia (2024).

####  4.3.1 Number of iterations

The lower OMP method is an adaptive method. On each iteration, the method
adaptively changes the size of the lower index set and outputs an
approximation to the PDE solution. As in the previous sections, we use the
relative $L^{2}$-error to measure accuracy. Here the size of the index set
directly relates to the computational cost of the lower OMP method. Hence, we
study the relative $L^{2}$-error and the cardinality of the lower index set
over the number of iterations to measure the performance of the method. From
the numerical experiments in Fig. 7, we observe that for the different
examples considered, the convergence rate and the increments in the size of
the index set vary substantially.

Figure 7: (Number of iterations) Relative $L^{2}$-error (blue curve, left
y-label) and cardinality (red curve, right y-label) of the lower set versus
the number of iterations for the exact solutions defined in Eq. (26)-(28).
Other parameters include dimension $d=6$ and sample points $m=3000$.

In the case of Example 1, once the algorithm discovers the index for the
complex Fourier basis function $\sin{4\pi x_{1}}\sin{2\pi x_{2}}$, then the
$L^{2}$-error decays to less than $10^{-10}$, otherwise, the approach does not
converge. In this example, the visualization using sample geometric mean can
not fully describe the convergence behavior, so we also plot each run’s
relative $L^{2}$-error using black data points. As the number of iterations
increases, we observe the algorithm is more likely to discover the correct
Fourier basis function (i.e., more runs reach relative $L^{2}$-error below
$10^{-10}$). However, even after 150 iterations, not all trials have converged
to the exact solution as we see that some data points are still near the top
of the figure. This suboptimal behaviour is not surprising since the multi-
index set corresponding to the basis function $\sin{4\pi x_{1}}\sin{2\pi
x_{2}}$, i.e., $\\{\pm 2\bm{e}_{1}\pm\bm{e}_{2}\\}$, is not lower. Hence,
Example 1 is not ideally suited for the adaptive lower OMP setting.

The situation is different for the other two examples. Example 2 has only two
activated variables, and the coefficients satisfy a lower structure. The
adaptive lower OMP method successfully finds the appropriate two-dimensional
lower set and converges to the exact solution after 100 iterations for all
trials. The cardinality of the lower set increases at a slower rate in this
example because the lower set only expands in the direction of $x_{1}$ and
$x_{2}$. For Example 3, the $L^{2}$-error converges to approximately $10^{-3}$
within 200 iterations. In this example we also observe that the index set’s
cardinality significantly increases because the index set is being adaptively
and anisotropically extended in all active dimensions.

####  4.3.2 Comparison with traditional OMP method

We now compare the performance for both the lower OMP and traditional OMP
methods in moderate dimension $d=6$. We note that the traditional OMP method
is not applicable in higher dimensional problems due to the requirement of
using an enormous index set. On the other hand, the lower OMP method can
produce results in problems as high as $d=30$ dimensions because it adaptively
searches for the best index set (see §4.4 for more lower OMP results with
$d=30$). Fig. 8 compares the $L^{2}$-error of the traditional and the lower
OMP methods with the same solution sparsity. Traditional OMP has a better
performance when the Fourier expansion of the exact solution is fully captured
by the traditional OMP ambient space (see, e.g., the left plot for Example 1).
Otherwise, the lower OMP method adaptively extends the index set using the
existing residual data and achieves better accuracy after more iterations
(see, e.g., the middle and right plots for Examples 2 and 3).

Figure 8: (Comparison OMP) Relative $L^{2}$-error versus the sparsity of the
OMP solution for the exact solutions defined in Eq. (27)-(28), Here dimension
$d=6$, and we choose the hyperbolic cross multi-index set
$\Lambda^{\mathrm{HC}}_{d,n}$ with $n=18$ (cardinality $=3418$) as ambient
space for the traditional OMP.

###  4.4 Periodic PINNs vs. CFC

In this section we compare the performance of periodic PINNs and CFC with
adaptive lower OMP recovery on high-dimensional test problems. To make a fair
comparison, we use identical sets of sample points as the training data for
both methods. For the neural network, we use 30000 epochs to examine the
performance of the PINNs, which we empirically observe allows for the networks
to saturate on the training data. In the adaptive lower OMP algorithm, the
stopping criterion is set to be the size of the index set (which we recall is
also an indicator of the computational cost). We choose this size to be
greater than half of the number of sample points. As illustrated in Fig. 9,
the adaptive lower OMP method has advantages in discovering the underlying
anisotropy of the high-dimensional functions, see, e.g., the results in
Examples 2 and 3 where the function is inactive in most variables or when the
coefficients satisfy a lower structure. The adaptive lower OMP method captures
the most important terms in Example 2 as the number of samples reaches 1000.
In all three cases, the neural network accurately approximates the solution as
the number of sample points increases and reaches error approximately
$10^{-2}$ (Example 1) and $10^{-3}$ (Examples 2 and 3). In Example 3, the
error of the periodic PINN reaches a plateau after approximately 2500 samples
(see also §4.2). However, similar to Example 1 in Fig. 8, the lower OMP method
struggles to converge for Example 1, while the periodic neural network
provides slightly better average relative $L^{2}$-error at 8000 sample points
with tighter spread in the standard deviation over the trials. In summary, CFC
can achieve much higher accuracy than periodic PINNs (gaining from 1 to 8
orders of magnitude), but its performance depends on the sparsity properties
of $u$. On the other hand, periodic PINNs are able to achieve a consistent
accuracy level ($10^{-2}$ relative $L^{2}$-error) on all examples.

Figure 9: (Comparison) Relative $L^{2}$-error versus the number of sample
points $m$ for the exact solutions defined in Eq. (27)-(28), Here dimension
$d=30$, number of sample points $m=3000$. The neural networks have the same
structure as in §4.2 and the relative error is measured after 30000 epochs.
The cardinality of the index sets in the Lower adaptive OMP is constrained to
be less than or equal to $m/2$

##  5 Proof of the practical existence theorem

In order to prove Theorem 1, we first extend the CFC Theorem of Wang and
Brugiapaglia (2024) to handle periodic diffusion equations with a reaction
term in §5.1. This requires first outlining the theory of bounded Riesz
systems and forming a connection with the CFC matrix and its corresponding
Gram matrix. After illustrating this CFC convergence result, we will show a
proof of Theorem 1 in §5.2.

###  5.1 Compressive Fourier collocation for diffusion-reaction problems

In this subsection we extend the CFC convergence analysis of Wang and
Brugiapaglia (2024) from diffusion to diffusion-reaction problems of the form
(1). The general proof strategy is to present sufficient conditions for the
PDE coefficients $a$ and $\rho$ such that the CFC matrix (24) is a random
sampling matrix from a bounded Riesz system (Brugiapaglia et al., 2021b). This
step is the most technical aspect of the proof. The rest uses the existing
framework of sparse recovery in bounded Riesz systems whereby a sufficient
condition on the sampling complexity is chosen such that the CFC satisfies the
robust null space property with high-probability. From compressive sensing
theory, we then obtain recovery guarantees for the CFC approximation.

####  5.1.1 Some preliminary facts

We start by recalling the definition of the bounded Riesz system. We restrict
our attention to Riesz systems in $L^{2}(\mathbb{T}^{d})$, although the
definition can be extended to general Hilbert spaces. Note in the following
definition that $\ell^{2}(\Lambda;\mathbb{C})$ is the set of
$\ell^{2}$-integrable complex sequences, indexed by $\Lambda$, i.e.,
$\ell^{2}(\Lambda;\mathbb{C})=\left\\{\bm{z}=(z_{\bm{\nu}})_{\bm{\nu}\in\Lambda}:\sum_{\bm{\nu}\in\Lambda}|z_{\bm{\nu}}|^{2}<\infty\right\\}.$

######  Definition 4 (Bounded Riesz System)

Let $\Lambda\subseteq\mathbb{Z}^{d}$, $0<b_{\Phi}\leq B_{\Phi}<\infty$, and
let $\ell^{2}(\Lambda;\mathbb{C})$ denote the space of sequences
$\bm{z}=(z_{\bm{\nu}})_{\bm{\nu}\in\Lambda}$ with $\|\bm{z}\|_{2}<\infty$. A
set of functions $\\{\Phi_{\bm{\nu}}\\}_{\bm{\nu}\in\Lambda}\subset
L^{2}(\mathbb{T}^{d})$ is a _Riesz system_ with constants $b_{\Phi}$ and
$B_{\Phi}$ if

$b_{\Phi}\|\bm{z}\|^{2}_{2}\leq\left\|\sum_{\bm{\nu}\in\Lambda}z_{\bm{\nu}}\Phi_{\bm{\nu}}\right\|^{2}_{L^{2}}\leq
B_{\Phi}\|\bm{z}\|^{2}_{2},\quad\forall\bm{z}=(z_{\bm{\nu}})_{\bm{\nu}\in\Lambda}\in\ell^{2}(\Lambda;\mathbb{C}).$

The constants $b_{\Phi}$ and $B_{\Phi}$ are called _lower and upper Riesz
constants_ , respectively. Moreover, the system
$\\{\Phi_{\bm{\nu}}\\}_{\bm{\nu}\in\Lambda}$ is _bounded_ if there exists a
constant $0<K_{\Phi}<\infty$ such that

$\|\Phi_{\bm{\nu}}\|_{L^{\infty}}\leq
K_{\Phi},\quad\forall\bm{\nu}\in\Lambda.$

Note that any $L^{2}$-orthonormal system is a Riesz system with
$b_{\Phi}=B_{\Phi}=1$. In particular, the Fourier system
$\\{F_{\bm{\nu}}\\}_{\bm{\nu}\in\mathbb{Z}^{d}}$ is a bounded Riesz system
with $b_{\Phi}=B_{\Phi}=K_{\Phi}=1$

Recalling the definition (3) of $\mathscr{L}$, we define

$\Phi_{\bm{\nu}}=\mathscr{L}[\Psi_{\bm{\nu}}],\quad\forall\bm{\nu}\in\Lambda,$
(29)

where $\\{\Psi_{\bm{\nu}}\\}_{\bm{\nu}\in\mathbb{Z}^{d}}$ is the renormalized
Fourier system given in (13). To show that the system
$\\{\Phi_{\bm{\nu}}\\}_{\bm{\nu}\in\Lambda}$ defined in (29) is a Riesz
system, it is convenient to consider its _Gram matrix_
$G\in\mathbb{C}^{N\times N}$, where $N=|\Lambda|$, defined by

$G_{\bm{\nu}\bm{\mu}}=\left<\Phi_{\bm{\nu}},\Phi_{\bm{\mu}}\right>,\quad\forall\bm{\nu},\bm{\mu}\in\Lambda.$
(30)

We note in passing that, thanks to the normalization factor $1/\sqrt{m}$ in
(24), we have $\mathbb{E}[A^{*}A]=G$ (this follows from a direct computation
and the fact that the random collocation points $\bm{y}_{1},\ldots,\bm{y}_{m}$
are independently and uniformly distributed over $\mathbb{T}^{d}$). The
significance of the Gram matrix $G$ relies on the fact that it yields the
following norm equivalence:

$\left\|\sum_{\bm{\nu}\in\Lambda}c_{\bm{\nu}}\Phi_{\bm{\nu}}\right\|_{L^{2}}^{2}=\left<\sum_{\bm{\nu}\in\Lambda}c_{\bm{\nu}}\Phi_{\bm{\nu}},\sum_{\bm{\nu}\in\Lambda}c_{\bm{\nu}}\Phi_{\bm{\nu}}\right>=\bm{c}^{T}G\bm{c},\quad\forall\bm{c}\in\mathbb{C}^{N}.$
(31)

Note that $G$ is a Hermitian positive semidefinite matrix. Hence, it has only
real nonnegative eigenvalues. The Courant–Fischer–Weyl min-max principle
implies that, if $0<b_{\Phi}\leq B_{\Phi}<\infty$ are such that

$b_{\Phi}\leq\lambda_{\min}(G)\leq\lambda_{\max}(G)\leq B_{\Phi},$ (32)

then $\\{\Phi_{\bm{\nu}}\\}_{\bm{\nu}\in\Lambda}$ is a Riesz system with
constants $b_{\Phi}$ and $B_{\Phi}$. Hence, estimating the lower and upper
Riesz constants of $\\{\Phi_{\bm{\nu}}\\}_{\bm{\nu}\in\Lambda}$ corresponds to
finding two-sided spectral bounds for the Gram matrix $G$. To obtain this type
of spectral bounds, we will employ Gershgorin’s circle theorem, see, e.g.,
(Horn and Johnson, 2012, Theorem 6.1.1), of which Lemma 6 is a direct
consequence. Here we state a useful corollary of the Gershgorin circle
theorem.

######  Lemma 5 (Gershgorin’s circle theorem for Hermitian matrices)

Let $A\in\mathbb{C}^{N\times N}$ be a Hermitian matrix. Then, all eigenvalues
of $A$ lie in the real interval

$\left[\min_{i\in[N]}\left\\{A_{ii}-\sum_{j\neq
i}|A_{ij}|\right\\},\max_{i\in[N]}\left\\{A_{ii}+\sum_{j\neq
i}|A_{ij}|\right\\}\right].$

Using the Fourier expansion (15) of $a$, it is possible to compute an explicit
formula for the entries of the Gram matrix $G$.

######  Lemma 6 (Explicit formula for the Gram matrix)

Let $\Lambda\subseteq\mathbb{Z}^{d}$, $r_{\bm{\nu}}\in\mathbb{C}$, with
$\bm{\nu}\in\Lambda$, be generic rescaling constants for the Fourier system
such that $\Psi_{\bm{\nu}}=r_{\bm{\nu}}F_{\bm{\nu}}$, and consider a diffusion
coefficient $a\in C^{1}(\mathbb{T}^{d})$ having Fourier expansion (15) with
$T=\mathbb{Z}^{d}\setminus\\{\bm{0}\\}$. Then, the elements of $G$ admit the
following explicit formula in terms of the Fourier coefficients
$(a_{\bm{\tau}})_{\bm{\tau}\in\mathbb{Z}^{d}}$ of a and the reaction term
$\rho\in\mathbb{R}$:

$\displaystyle G_{\bm{\nu}\bm{\mu}}$
$\displaystyle=r_{\bm{\nu}}\bar{r}_{\bm{\mu}}\bigg{(}16\pi^{4}\sum_{\bm{\tau}\in\mathbb{Z}^{d}}\sum_{\bm{\tau^{\prime}}\in\mathbb{Z}^{d}}(\bm{\tau}\cdot\bm{\nu}+\|\bm{\nu}\|_{2}^{2})(\bm{\tau^{\prime}}\cdot\bm{\mu}+\|\bm{\mu}\|_{2}^{2})a_{\bm{\tau}}\bar{a}_{\bm{\tau^{\prime}}}\delta_{\bm{\tau}+\bm{\nu},\bm{\tau^{\prime}}+\bm{\mu}}+\rho^{2}\delta_{\bm{\nu},\bm{\mu}}$
$\displaystyle\quad+4\pi^{2}\rho\sum_{\bm{\tau}\in\mathbb{Z}^{d}}(\bm{\tau}\cdot\bm{\nu}+\|\bm{\nu}\|_{2}^{2})a_{\bm{\tau}}\delta_{\bm{\tau}+\bm{\nu},\bm{\mu}}+4\pi^{2}\rho\sum_{\bm{\tau^{\prime}}\in\mathbb{Z}^{d}}(\bm{\tau^{\prime}}\cdot\bm{\mu}+\|\bm{\mu}\|_{2}^{2})\bar{a}_{\bm{\tau^{\prime}}}\delta_{\bm{\nu},\bm{\tau^{\prime}}+\bm{\mu}}\bigg{)},$

for each $\bm{\nu},\bm{\mu}$ in $\Lambda\setminus\\{\bm{0}\\}$ and where
$\delta_{\bm{\nu},\bm{\mu}}$ denotes the Kronecker delta. Moreover,

$G_{\bm{0}\bm{0}}=r_{\bm{0}}^{2}\rho^{2}\quad\text{and}\quad
G_{\bm{\bm{\nu}}\bm{0}}=G_{\bm{\bm{0}\bm{\nu}}}=0,\quad\forall\bm{\nu}\in\Lambda\setminus\\{\bm{0}\\}.$
(33)

Proof  Before proving the identity, we note that gradients and Laplacians of
the Fourier basis functions defined in (11) can be easily computed as

$\nabla F_{\bm{\nu}}=(2\pi i\bm{\nu})F_{\bm{\nu}}\quad\text{and}\quad\Delta
F_{\bm{\nu}}=-4\pi^{2}\|\bm{\nu}\|_{2}^{2}F_{\bm{\nu}},\quad\forall\bm{\nu}\in\mathbb{Z}^{d}.$
(34)

Moreover,

$F_{\bm{\nu}}F_{\bm{\mu}}=F_{\bm{\nu}+\bm{\mu}},\quad\forall\bm{\nu},\bm{\mu}\in\mathbb{Z}^{d}.$
(35)

To prove the desired formula for $G_{\bm{\nu},\bm{\mu}}$, we expand the inner
product in (30). Using the above properties, the $L^{2}$-orthonormality of the
Fourier basis $\\{F_{\bm{\nu}}\\}_{\bm{\nu}\in\mathbb{Z}^{d}}$, and recalling
the expansion (15) of $a$, if $\rho$ is a constant, we see that

$\displaystyle G_{\bm{\nu}\bm{\mu}}$ $\displaystyle=\langle\nabla
a\cdot\nabla\Psi_{\bm{\nu}}+a\Delta\Psi_{\bm{\nu}}-\rho\Psi_{\bm{\nu}},\nabla
a\cdot\nabla\Psi_{\bm{\mu}}+a\Delta\Psi_{\bm{\mu}}-\rho\Psi_{\bm{\mu}}\rangle$
$\displaystyle=\langle\nabla
a\cdot\nabla\Psi_{\bm{\nu}}+a\Delta\Psi_{\bm{\nu}},\nabla
a\cdot\nabla\Psi_{\bm{\mu}}+a\Delta\Psi_{\bm{\mu}}\rangle+\rho^{2}\langle\Psi_{\bm{\nu}},\Psi_{\bm{\mu}}\rangle$
$\displaystyle\quad-\rho\langle\nabla
a\cdot\nabla\Psi_{\bm{\nu}}+a\Delta\Psi_{\bm{\nu}},\Psi_{\bm{\mu}}\rangle-\rho\langle\Psi_{\bm{\nu}},\nabla
a\cdot\nabla\Psi_{\bm{\mu}}+a\Delta\Psi_{\bm{\mu}}\rangle$
$\displaystyle=\sum_{\bm{\tau}}\sum_{\bm{\tau^{\prime}}}\langle((2i\pi\bm{\tau})\cdot(2i\pi\bm{\nu})-4\pi^{2}\|\bm{\nu}\|_{2}^{2})a_{\bm{\tau}}F_{\bm{\tau}}\Psi_{\bm{\nu}},((2i\pi\bm{\bm{\tau}}^{\prime})\cdot(2i\pi\bm{\mu})-4\pi^{2}\|\bm{\mu}\|_{2}^{2})a_{\bm{\tau^{\prime}}}F_{\bm{\tau^{\prime}}}\Psi_{\bm{\mu}}\rangle$
$\displaystyle\quad+\rho^{2}\langle\Psi_{\bm{\nu}},\Psi_{\bm{\mu}}\rangle-\rho\sum_{\bm{\tau}}\langle((2i\pi\bm{\tau})\cdot(2i\pi\bm{\nu})-4\pi^{2}\|\bm{\nu}\|_{2}^{2})a_{\bm{\tau}}F_{\bm{\tau}}\Psi_{\bm{\nu}},\Psi_{\bm{\mu}}\rangle$
$\displaystyle\quad-\rho\sum_{\bm{\tau^{\prime}}}\langle\Psi_{\bm{\nu}},((2i\pi\bm{\bm{\tau}}^{\prime})\cdot(2i\pi\bm{\mu})-4\pi^{2}\|\bm{\mu}\|_{2}^{2})a_{\bm{\tau^{\prime}}}F_{\bm{\tau^{\prime}}}\Psi_{\bm{\mu}}\rangle$
$\displaystyle=16\pi^{4}r_{\bm{\nu}}\bar{r}_{\bm{\mu}}\sum_{\bm{\tau}}\sum_{\bm{\tau^{\prime}}}(\bm{\tau}\cdot\bm{\nu}+\|\bm{\nu}\|_{2}^{2})(\bm{\tau^{\prime}}\cdot\bm{\mu}+\|\bm{\mu}\|_{2}^{2})a_{\bm{\tau}}\bar{a}_{\bm{\tau^{\prime}}}\langle
F_{\bm{\tau}}F_{\bm{\nu}},F_{\bm{\tau^{\prime}}}F_{\bm{\mu}}\rangle+\rho^{2}r_{\bm{\nu}}\bar{r}_{\bm{\mu}}\langle
F_{\bm{\nu}},F_{\bm{\mu}}\rangle$ $\displaystyle\quad+4\pi^{2}\rho
r_{\bm{\nu}}\bar{r}_{\bm{\mu}}\sum_{\bm{\tau}}(\bm{\tau}\cdot\bm{\nu}+\|\bm{\nu}\|_{2}^{2})a_{\bm{\tau}}\langle
F_{\bm{\tau}}F_{\bm{\nu}},F_{\bm{\mu}}\rangle+4\pi^{2}\rho
r_{\bm{\nu}}\bar{r}_{\bm{\mu}}\sum_{\bm{\tau^{\prime}}}(\bm{\tau^{\prime}}\cdot\bm{\mu}+\|\bm{\mu}\|_{2}^{2})\bar{a}_{\bm{\tau^{\prime}}}\langle
F_{\bm{\nu}},F_{\bm{\tau^{\prime}}}F_{\bm{\mu}}\rangle$
$\displaystyle=r_{\bm{\nu}}\bar{r}_{\bm{\mu}}\bigg{(}16\pi^{4}\sum_{\bm{\tau}}\sum_{\bm{\tau^{\prime}}}(\bm{\tau}\cdot\bm{\nu}+\|\bm{\nu}\|_{2}^{2})(\bm{\tau^{\prime}}\cdot\bm{\mu}+\|\bm{\mu}\|_{2}^{2})a_{\bm{\tau}}\bar{a}_{\bm{\tau^{\prime}}}\delta_{\bm{\tau}+\bm{\nu},\bm{\tau^{\prime}}+\bm{\mu}}+\rho^{2}\delta_{\bm{\nu},\bm{\mu}}$
$\displaystyle\quad+4\pi^{2}\rho\sum_{\bm{\tau}}(\bm{\tau}\cdot\bm{\nu}+\|\bm{\nu}\|_{2}^{2})a_{\bm{\tau}}\delta_{\bm{\tau}+\bm{\nu},\bm{\mu}}+4\pi^{2}\rho\sum_{\bm{\tau^{\prime}}}(\bm{\tau^{\prime}}\cdot\bm{\mu}+\|\bm{\mu}\|_{2}^{2})\bar{a}_{\bm{\tau^{\prime}}}\delta_{\bm{\nu},\bm{\tau^{\prime}}+\bm{\mu}}\bigg{)},$

where all the summations are over
$\bm{\tau},\bm{\tau}^{\prime}\in\mathbb{Z}^{d}$, and both $\bm{\nu},\bm{\mu}$
not equal to $\bm{0}$. Moreover, for any $\bm{\nu}\neq\bm{0}$, we have

$\displaystyle G_{\bm{\nu}\bm{0}}$ $\displaystyle=\langle\nabla
a\cdot\nabla\Psi_{\bm{\nu}}+a\Delta\Psi_{\bm{\nu}}-\rho\Psi_{\bm{\nu}},-\rho\Psi_{\bm{0}}\rangle$
$\displaystyle=\sum_{\bm{\tau}}\langle((2i\pi\bm{\tau})\cdot(2i\pi\bm{\nu})-4\pi^{2}\|\bm{\nu}\|_{2}^{2})a_{\bm{\tau}}F_{\bm{\tau}}\Psi_{\bm{\nu}},-\rho\Psi_{\bm{0}}\rangle$
$\displaystyle=4\pi^{2}\rho\sum_{\bm{\tau}}(\bm{\tau}\cdot\bm{\nu}+\|\bm{\nu}\|_{2}^{2})a_{\bm{\tau}}\rho\langle
F_{\bm{\tau}}F_{\bm{\nu}},F_{\bm{0}}\rangle$
$\displaystyle=4\pi^{2}\rho\sum_{\bm{\tau}}(\bm{\tau}\cdot\bm{\nu}+\|\bm{\nu}\|_{2}^{2})a_{\bm{\tau}}\rho\delta_{\bm{\tau}+\bm{\nu},\bm{0}}=0.$

Similarly, one can show that $G_{\bm{0}\bm{\mu}}=0$, and
$G_{\bm{0}\bm{0}}=\rho^{2}\langle
r_{\bm{0}}F_{\bm{0}},r_{\bm{0}}F_{\bm{0}}\rangle=r_{\bm{0}}^{2}\rho^{2}$.

We also need an auxiliary result about norm equivalencies in
$H^{2}(\mathbb{T}^{d})$ that will be necessary to derive the final error bound
of Theorem 8. Letting
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|u\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}:=\|u\|_{L^{2}}^{2}+\|\Delta
u\|_{L^{2}}^{2}$, the following lemma shows that this norm is equivalent to
the canonical $H^{2}$-norm.

###### Lemma 7

The norms $\vvvert\cdot\vvvert$ and $\|\cdot\|_{H^{2}}$ are equivalent.
Specifically, $\sqrt{2/3}\|u\|_{H^{2}}\leq\vvvert u\vvvert\leq\|u\|_{H^{2}}$
for every $u\in H^{2}(\mathbb{T}^{d})$.

Proof  Let $u(\bm{x})=\sum_{\bm{\nu}\in\mathbb{Z}^{d}}{c}_{\bm{\nu}}\exp(2\pi
i\bm{\nu}\cdot\bm{x})$ and recall that
$\|u\|_{H^{2}}^{2}=\|u\|_{L^{2}}^{2}+\|\nabla
u\|_{L^{2}}^{2}+\|\nabla^{2}u\|_{L^{2}}^{2}.$ Using Parseval’s identity, we
have

$\displaystyle\|\nabla
u\|_{L^{2}}^{2}=(2\pi)^{2}\int_{\mathbb{T}^{d}}\sum_{k=1}^{d}\left|\sum_{\bm{\nu}\in\mathbb{Z}^{d}}c_{\bm{\nu}}\nu_{k}\exp(2\pi
i\bm{\nu}\cdot\bm{x})\right|^{2}d\bm{x}$
$\displaystyle=(2\pi)^{2}\sum_{k=1}^{d}\sum_{\bm{\nu}\in\mathbb{Z}^{d}}|c_{\bm{\nu}}|^{2}|\nu_{k}|^{2}$
$\displaystyle=(2\pi)^{2}\sum_{\bm{\nu}\in\mathbb{Z}^{d}}|c_{\bm{\nu}}|^{2}\|\bm{\nu}\|_{2}^{2}.$

Moreover, recalling that
$\|\nabla^{2}u\|_{L^{2}}^{2}=\int_{\mathbb{T}^{d}}\|\nabla^{2}u(\bm{x})\|_{F}^{2}d\bm{x}$
and using Parseval’s identity again, we obtain

$\displaystyle\|\nabla^{2}u\|_{L^{2}}^{2}$
$\displaystyle=(2\pi)^{4}\int_{\mathbb{T}^{d}}\sum_{k=1}^{d}\sum_{l=1}^{d}\left|\sum_{\bm{\nu}\in\mathbb{Z}^{d}}c_{\bm{\nu}}\nu_{k}\nu_{l}\exp(2\pi
i\bm{\nu}\cdot\bm{x})\right|^{2}d\bm{x}$
$\displaystyle=(2\pi)^{4}\sum_{k=1}^{d}\sum_{l=1}^{d}\sum_{\bm{\nu}\in\mathbb{Z}^{d}}|c_{\bm{\nu}}|^{2}(\nu_{k}\nu_{l})^{2}=(2\pi)^{4}\sum_{\bm{\nu}\in\mathbb{Z}^{d}}|c_{\bm{\nu}}|^{2}\|\bm{\nu}\|_{2}^{4}.$

Now, observe that $\|u\|_{L^{2}}=\|\bm{c}\|_{2}$. Thus,

$\|u\|_{H^{2}}^{2}=\sum_{\bm{\nu}\in\mathbb{Z}^{d}}|c_{\bm{\nu}}|^{2}(1+(2\pi)^{2}\|\bm{\nu}\|_{2}^{2}+(2\pi)^{4}\|\bm{\nu}\|_{2}^{4})=\bm{c}^{*}D_{1}\bm{c},$

with
$D_{1}=\text{diag}((1+(2\pi)^{2}\|\bm{\nu}\|_{2}^{2}+(2\pi)^{4}\|\bm{\nu}\|_{2}^{4})_{\bm{\nu}\in\mathbb{Z}^{d}})$.
Moreover,

$\displaystyle\|\Delta u\|_{L^{2}}^{2}$
$\displaystyle=(2\pi)^{4}\int_{\mathbb{T}^{d}}\left|\sum_{\bm{\nu}\in\mathbb{Z}^{d}}c_{\bm{\nu}}\|\bm{\nu}\|_{2}^{2}\exp(2\pi
i\bm{\nu}\cdot\bm{x})\right|^{2}d\bm{x}=(2\pi)^{4}\sum_{\bm{\nu}\in\mathbb{Z}^{d}}|c_{\bm{\nu}}|^{2}\|\bm{\nu}\|_{2}^{4}.$

Therefore,

$\vvvert
u\vvvert^{2}=\sum_{\bm{\nu}\in\mathbb{Z}^{d}}|c_{\bm{\nu}}|^{2}(1+(2\pi)^{4}\|\bm{\nu}\|_{2}^{4})=\bm{c}^{*}D_{2}\bm{c},$

where
$D_{2}=\text{diag}((1+(2\pi)^{4}\|\bm{\nu}\|_{2}^{4})_{\bm{\nu}\in\mathbb{Z}^{d}})$.
Hence,

$\frac{\|u\|_{H^{2}}^{2}}{\vvvert
u\vvvert^{2}}=\frac{\bm{c}^{*}D_{1}\bm{c}}{\bm{c}^{*}D_{2}\bm{c}}=\frac{\bm{d}^{*}D_{3}\bm{d}}{\bm{d}^{*}\bm{d}}=\frac{\sum_{\bm{\nu}\in\mathbb{Z}^{d}}(D_{3})_{\bm{\nu},\bm{\nu}}|d_{\bm{\nu}}|^{2}}{\sum_{\bm{\nu}\in\mathbb{Z}^{d}}|d_{\bm{\nu}}|^{2}}\leq\sup_{\bm{\nu}\in\mathbb{Z}^{d}}((D_{3})_{\bm{\nu},\bm{\nu}}),$

where we made the change of variable $\bm{d}=D_{2}^{\frac{1}{2}}\bm{c}$ and
where

$D_{3}=D_{2}^{-\frac{1}{2}}D_{1}D_{2}^{-\frac{1}{2}}=\text{diag}\left(\left(\frac{1+(2\pi)^{2}\|\bm{\nu}\|_{2}^{2}+(2\pi)^{4}\|\bm{\nu}\|_{2}^{4}}{1+(2\pi)^{4}\|\bm{\nu}\|_{2}^{4}}\right)_{\bm{\nu}\in\mathbb{Z}^{d}}\right).$

In particular,

$\sup_{\bm{\nu}\in\mathbb{Z}^{d}}((D_{3})_{\bm{\nu},\bm{\nu}})=\sup_{\bm{\nu}\in\mathbb{Z}^{d}}\frac{1+(2\pi)^{2}\|\bm{\nu}\|_{2}^{2}+(2\pi)^{4}\|\bm{\nu}\|_{2}^{4}}{1+(2\pi)^{4}\|\bm{\nu}\|_{2}^{4}}\leq\frac{3}{2}$

The last bound is obtained by maximizing $x\mapsto(1+x^{2}+x^{4})/(1+x^{4})$
over $\mathbb{R}$. Similarly, we obtain the second inequality in the statement
of equivalence from

$\frac{\vvvert
u\vvvert^{2}}{\|u\|_{H^{2}}^{2}}\leq\sup_{\bm{\nu}\in\mathbb{Z}^{d}}\frac{1+(2\pi)^{4}\|\bm{\nu}\|_{2}^{4}}{1+(2\pi)^{2}\|\bm{\nu}\|_{2}^{2}+(2\pi)^{4}\|\bm{\nu}\|_{2}^{4}}\leq
1.$

This concludes the proof.

####  5.1.2 CFC convergence

We are now ready to illustrate our CFC convergence result. Keeping the CFC
setup of §3 in mind, in order to recover a compressible solution
$\hat{\bm{c}}$ to the linear system (23), we utilize the _Square-Root LASSO_
(in short, SR-LASSO), see (Adcock et al., 2019; Belloni et al., 2011) and
references therein, a modified version of the original LASSO (Least Absolute
Shrinkage and Selection Operator), which lacks a power of 2 on the data
fidelity term. The SR-LASSO is defined as

$\min_{\bm{z}\in\mathbb{C}^{N}}\left\\{\|A\bm{z}-\bm{b}\|_{2}+\lambda\|\bm{z}\|_{1}\right\\}.$
(36)

The main benefit of this formulation is that the optimal choice of the tuning
parameter $\lambda>0$ is independent of the noise level, see (Adcock et al.,
2022b, §6.6.2).

######  Theorem 8 (Convergence of CFC for diffusion-reaction problems)

Given a dimension $d\in\mathbb{N}$, target sparsity $s\in\mathbb{N}$,
hyperbolic cross order $n\in\mathbb{N}$ and probability of failure
$\varepsilon\in(0,1)$, let
$\Lambda=\Lambda^{\textnormal{HC}}_{d,n}\subset\mathbb{Z}^{d}$ and suppose
$a\in C^{1}(\mathbb{T}^{d})$ and $\rho\in\mathbb{R}$ satisfy (2) and
(15)–(16). Then, the system $\\{\Phi_{\bm{\nu}}\\}_{\bm{\nu}\in\Lambda}$
defined in (29) is a bounded Riesz system in the sense of Definition 4 with
constants

$\displaystyle b_{\Phi}$
$\displaystyle=a_{\bm{0}}^{2}-\left(2a_{\bm{0}}+\frac{\rho}{2\pi^{2}}\right)\beta-\beta^{2}>0,$
(37) $\displaystyle B_{\Phi}$
$\displaystyle=\left\|a\right\|_{H^{1}}^{2}+\frac{\rho^{2}}{16\pi^{4}}+\frac{a_{\bm{0}}\rho}{2\pi^{2}}+\left(2a_{\bm{0}}+\frac{\rho}{2\pi^{2}}\right)\beta+\beta^{2},$
(38)

where $\beta=\sqrt{|T|}\|a-a_{\bm{0}}\|_{H^{1}}$, and
$K_{\Phi}=a_{\bm{0}}+\beta+\frac{\rho}{4\pi^{2}}.$ Moreover, the CFC solution
$\hat{u}$ in (25) approximating the high-dimensional periodic diffusion-
reaction equation (1), whose coefficients $\hat{\bm{c}}$ are computed by
solving the SR-LASSO problem (36) with $A$ and $\bm{b}$ defined as in (24) and
with tuning parameter $d^{(1)}_{a,\rho}\sqrt{B_{\Phi}/s}<\lambda\leq
d^{(2)}_{a,\rho}\sqrt{B_{\Phi}/s}$, satisfies the following with probability
at least $1-\varepsilon$: if

$m\geq c^{(3)}_{a,\rho}\cdot s\cdot\log^{2}\left(c^{(4)}_{a,\rho}\cdot
s\right)\cdot\left(\min\\{\log(n)+d,\log(2n)\log(2d)\\}+\log(\varepsilon^{-1})\right),$
(39)

then

$\displaystyle\|u-\hat{u}\|_{L^{2}}+\|(\Delta-\rho)(u-\hat{u})\|_{L^{2}}$
$\displaystyle\leq
C^{(1)}_{a,\rho}\cdot\frac{\sigma_{s}(\bm{c}_{\Lambda})_{1}}{\sqrt{s}}+C^{(2)}_{a,d,\rho}\cdot\left(\frac{\|u-u_{\Lambda}\|_{W^{2,\infty}}}{\sqrt{s}}+\|u-u_{\Lambda}\|_{H^{2}}\right).$
(40)

Moreover, if $\rho<1$, we also have

$\displaystyle\|u-\hat{u}\|_{H^{2}}\leq
C^{(3)}_{a,\rho}\cdot\frac{\sigma_{s}(\bm{c}_{\Lambda})_{1}}{\sqrt{s}}+C^{(4)}_{a,d,\rho}\cdot\left(\frac{\|u-u_{\Lambda}\|_{W^{2,\infty}}}{\sqrt{s}}+\|u-u_{\Lambda}\|_{H^{2}}\right).$
(41)

Here each constant depends only on the subscripted parameters. Moreover, the
dependence of each constant on $d$ is at most linear (when present).

###### Remark 9

The above theorem holds for OMP recovery as well, under a sufficient condition
on the ratio $b_{\Phi}/B_{\Phi}$, see (Wang and Brugiapaglia, 2024, Theorem
3).

###### Remark 10

There is a gap between Theorem 8 and the corresponding result in the diffusion
setting (Wang and Brugiapaglia, 2024, Theorem 3.5). Namely, in the latter the
diffusion coefficient $a$ can be nonsparse. Extending Theorem 8 to allow for
more general diffusion coefficients is an open problem.

Proof  Given sparse diffusion and reaction terms $a$ and $\rho$ as in (2) and
(15)–(16), we determine that
$\left\\{\Phi_{\bm{\nu}}\right\\}_{\bm{\nu}\in\Lambda}$ is a bounded Riesz
system using the explicit form of the Gram matrix in Lemma 6.

##### Step 1: Riesz property.

We find lower and upper Riesz constants $b_{\Phi}$ and $B_{\Phi}$ by
establishing a two-sided spectral bound for the Gram matrix $G$, recall
equation (30). First, observing that the first row (and column) of $G$ only
has one nonzero entry (namely, $G_{\bm{00}}$), one eigenvalue of $G$ is
$G_{\bm{00}}=a_{\bm{0}}^{2}$. To estimate the remaining eigenvalues, we use a
special version of Gershgorin’s circle theorem for Hermitian matrices,
presented in Lemma 5, see, e.g., (Horn and Johnson, 2012, Theorem 6.1.1) for
the general statement. Using Lemma 6 with the generic normalization, the
diagonal entries $G_{\bm{\nu}\bm{\nu}}$ with
$\bm{\nu}\in\Lambda\setminus\\{\bm{0}\\}$ are given by

$\displaystyle G_{\bm{\nu}\bm{\nu}}$
$\displaystyle=16\pi^{4}|r_{\bm{\nu}}|^{2}\sum_{\bm{\tau}\in
T\cup\\{\bm{0}\\}}(\bm{\tau}\cdot\bm{\nu}+\|\bm{\nu}\|_{2}^{2})^{2}|a_{\bm{\tau}}|^{2}+\rho^{2}|r_{\bm{\nu}}|^{2}+8\pi^{2}|r_{\bm{\nu}}|^{2}\|\bm{\nu}\|_{2}^{2}\rho
a_{\bm{0}}$
$\displaystyle=16\pi^{4}\|\bm{\nu}\|_{2}^{4}|r_{\bm{\nu}}|^{2}\left[\left(a_{\bm{0}}+\frac{\rho}{4\pi^{2}\|\bm{\nu}\|_{2}^{2}}\right)^{2}+\sum_{\bm{\tau}\in
T}\left(\frac{\bm{\tau}\cdot\bm{\nu}+\|\bm{\nu}\|_{2}^{2}}{\|\bm{\nu}\|_{2}^{2}}\right)^{2}\left|a_{\bm{\tau}}\right|^{2}\right].$

Using the Cauchy-Schwarz inequality and inputting the actual normalization
used in (13), we see that, for all $\bm{\nu}\in\Lambda\setminus\\{\bm{0}\\}$,

$|r_{\bm{\nu}}|^{2}\left(a_{\bm{0}}+\frac{\rho}{4\pi^{2}\|\bm{\nu}\|_{2}^{2}}\right)^{2}\leq\frac{\lvert
G_{\bm{\nu}\bm{\nu}}\rvert}{16\pi^{4}\|\bm{\nu}\|_{2}^{4}}\leq|r_{\bm{\nu}}|^{2}\left[a_{\bm{0}}^{2}+\frac{\rho^{2}}{16\pi^{4}}+\frac{|a_{\bm{0}}|\rho}{2\pi^{2}}+\sum_{\bm{\tau}\in
T}\left(\|\bm{\tau}\|_{2}+1\right)^{2}\lvert a_{\bm{\tau}}\rvert^{2}\right]$

which, in turn, implies

$a_{\bm{0}}^{2}\leq\lvert
G_{\bm{\nu}\bm{\nu}}\rvert\leq\left\|a\right\|^{2}_{H^{1}}+\frac{\rho^{2}}{16\pi^{4}}+\frac{|a_{\bm{0}}|\rho}{2\pi^{2}}.$
(42)

The inequality defining the upper bound can be proved as follows. Using the
definition of $H^{1}$-norm, the differentiation properties (34), and the fact
that the Fourier system $\\{F_{\bm{\nu}}\\}_{\bm{\nu}\in\mathbb{Z}^{d}}$ is
$L^{2}$-orthonormal, we obtain

$\displaystyle\left\|a\right\|^{2}_{H^{1}}$
$\displaystyle=\left\|a\right\|^{2}_{L^{2}}+\sum^{d}_{l=1}\left\|\frac{\partial
a}{\partial x_{l}}\right\|^{2}_{L^{2}}=\sum_{\bm{\tau}\in
T\cup\\{\bm{0}\\}}|a_{\bm{\tau}}|^{2}+\sum^{d}_{l=1}\sum_{\bm{\tau}\in
T\cup\\{\bm{0}\\}}|2\pi a_{\bm{\tau}}\bm{\tau}_{l}|^{2}$
$\displaystyle=|a_{\bm{0}}|^{2}+\sum_{\bm{\tau}\in
T}\left(1+(2\pi)^{2}\left\|\bm{\tau}\right\|_{2}^{2}\right)|a_{\bm{\tau}}|^{2}\geq|a_{\bm{0}}|^{2}+\sum_{\bm{\tau}\in
T}\left(\|\bm{\tau}\|_{2}+1\right)^{2}\lvert a_{\bm{\tau}}\rvert^{2},$ (43)

which proves (42). To apply Gershgorin’s circle theorem, we now bound the sum
of all off-diagonal entries in the $\bm{\nu}$-th row of $G$. Using Lemma 6
again, the definition of the Kronecker delta, (33) and the fact that
$|r_{\bm{\nu}}|\leq 1/(4\pi^{2}\|\bm{\nu}\|_{2}^{2})\leq 1/4\pi^{2}$ for every
$\bm{\nu}\neq\bm{0}$, we obtain

$\displaystyle\sum_{\bm{\mu}\in\Lambda\setminus\\{\bm{\nu}\\}}\lvert
G_{\bm{\nu}\bm{\mu}}\rvert$
$\displaystyle=\sum_{\bm{\mu}\in\Lambda\setminus\\{\bm{\nu}\\}}\bigg{|}r_{\bm{\nu}}\bar{r}_{\bm{\mu}}\bigg{(}16\pi^{4}\sum_{\bm{\tau}\in\mathbb{Z}^{d}}(\bm{\tau}\cdot\bm{\nu}+\|\bm{\nu}\|_{2}^{2})\left((\bm{\tau}+\bm{\nu}-\bm{\mu})\cdot\bm{\mu}+\|\bm{\mu}\|_{2}^{2}\right)a_{\bm{\tau}}\bar{a}_{\bm{\tau}+\bm{\nu}-\bm{\mu}}$
$\displaystyle\quad+4\pi^{2}\rho\left((\bm{\mu}-\bm{\nu})\cdot\bm{\nu}+\|\bm{\nu}\|_{2}^{2}\right)a_{\bm{\mu}-\bm{\nu}}+4\pi^{2}\rho\left((\bm{\nu}-\bm{\mu})\cdot\bm{\mu}+\|\bm{\mu}\|_{2}^{2}\right)\bar{a}_{\bm{\nu}-\bm{\mu}}\bigg{)}\bigg{|}$
$\displaystyle\leq\sum_{\bm{\tau}\in\mathbb{Z}^{d}}\left(\frac{|\bm{\tau}\cdot\bm{\nu}|}{\|\bm{\nu}\|_{2}^{2}}+1\right)\lvert
a_{\bm{\tau}}\rvert\sum_{\bm{\mu}\in\Lambda\setminus\\{\bm{\nu}\\}}\left(\frac{|\left(\bm{\tau}+\bm{\nu}-\bm{\mu}\right)\cdot\bm{\mu}|}{\|\bm{\mu}\|_{2}^{2}}+1\right)\lvert
a_{\bm{\bm{\tau}+\bm{\nu}-\bm{\mu}}}\rvert$
$\displaystyle\quad+\frac{\rho}{4\pi^{2}}\sum_{\bm{\mu}\in\Lambda\setminus\\{\bm{\nu}\\}}\left(\frac{|(\bm{\mu}-\bm{\nu})\cdot\bm{\nu}|}{\|\bm{\nu}\|_{2}^{2}}+1\right)|a_{\bm{\mu}-\bm{\nu}}|$
$\displaystyle\quad+\frac{\rho}{4\pi^{2}}\sum_{\bm{\mu}\in\Lambda\setminus\\{\bm{\nu}\\}}\left(\frac{|(\bm{\nu}-\bm{\mu})\cdot\bm{\mu}|}{\|\bm{\mu}\|_{2}^{2}}+1\right)|\bar{a}_{\bm{\nu}-\bm{\mu}}|.$

Using the Cauchy-Schwarz inequality
$|\bm{\tau}\cdot\bm{\nu}|\leq\|\bm{\tau}\|_{2}\cdot\|\bm{\nu}\|_{2}$ and the
fact that $\|\bm{\mu}\|_{2}\geq 1$ and $\|\bm{\nu}\|_{2}\geq 1$ for all
$\bm{\mu},\bm{\nu}\neq\bm{0}$, we see that, for all
$\bm{\nu}\in\Lambda\setminus\\{\bm{0}\\}$,

$\displaystyle\sum_{\bm{\mu}\in\Lambda\setminus\\{\bm{\nu}\\}}\lvert
G_{\bm{\nu}\bm{\mu}}\rvert$ $\displaystyle\leq\sum_{\bm{\tau}\in
T\cup\\{\bm{0}\\}}\left(\|\bm{\tau}\|_{2}+1\right)\lvert
a_{\bm{\tau}}\rvert\sum_{\bm{\mu}\in\Lambda\setminus\\{\bm{\nu}\\}}\left(\|\bm{\tau}+\bm{\nu}-\bm{\mu}\|_{2}+1\right)\lvert
a_{\bm{\bm{\tau}+\bm{\nu}-\bm{\mu}}}\rvert$
$\displaystyle\quad+\frac{\rho}{4\pi^{2}}\sum_{\bm{\tau}\in
T}\left(\|\bm{\tau}\|_{2}+1\right)\lvert
a_{\bm{\tau}}\rvert+\frac{\rho}{4\pi^{2}}\sum_{\bm{\tau}\in
T}\left(\|\bm{\tau}\|_{2}+1\right)\lvert\bar{a}_{\bm{\tau}}\rvert.$

Substituting $\bm{\tau}^{\prime}=\bm{\tau}+\bm{\nu}-\bm{\mu}$ (which implies
$\bm{\tau}^{\prime}\neq\bm{\tau}$), recalling that
$(a_{\bm{\nu}})_{\bm{\nu}\in\mathbb{Z}^{d}}$ is supported on
$T\cup\\{\bm{0}\\}$, and separating the $a_{\bm{0}}$ term, we obtain

$\displaystyle\sum_{\bm{\mu}\in\Lambda\setminus\\{\bm{\nu}\\}}\lvert
G_{\bm{\nu}\bm{\mu}}\rvert$ $\displaystyle\leq\sum_{\bm{\tau}\in
T\cup\\{\bm{0}\\}}\left(\|\bm{\tau}\|_{2}+1\right)\lvert
a_{\bm{\tau}}\rvert\sum_{\bm{\tau^{\prime}}\in
T\cup\\{\bm{0}\\}\setminus\\{\bm{\tau}\\}}\left(\|\bm{\tau^{\prime}}\|_{2}+1\right)\lvert
a_{\bm{\tau^{\prime}}}\rvert+\frac{\rho}{2\pi^{2}}\sum_{\bm{\tau}\in
T}\left(\|\bm{\tau}\|_{2}+1\right)\lvert a_{\bm{\tau}}\rvert$
$\displaystyle=|a_{\bm{0}}|\sum_{\bm{\tau^{\prime}}\in
T}\left(\|\bm{\tau^{\prime}}\|_{2}+1\right)\lvert
a_{\bm{\tau^{\prime}}}\rvert+\sum_{\bm{\tau}\in
T}\left(\|\bm{\tau}\|_{2}+1\right)\lvert
a_{\bm{\tau}}\rvert\sum_{\bm{\tau^{\prime}}\in
T\cup\\{\bm{0}\\}\setminus\\{\bm{\tau}\\}}\left(\|\bm{\tau^{\prime}}\|_{2}+1\right)\lvert
a_{\bm{\tau^{\prime}}}\rvert$
$\displaystyle\quad+\frac{\rho}{2\pi^{2}}\sum_{\bm{\tau}\in
T}\left(\|\bm{\tau}\|_{2}+1\right)\lvert a_{\bm{\tau}}\rvert$
$\displaystyle\leq 2|a_{\bm{0}}|\sum_{\bm{\tau^{\prime}}\in
T}\left(\|\bm{\tau^{\prime}}\|_{2}+1\right)\lvert
a_{\bm{\tau^{\prime}}}\rvert+\sum_{\bm{\tau}\in
T}\left(\|\bm{\tau}\|_{2}+1\right)\lvert
a_{\bm{\tau}}\rvert\sum_{\bm{\tau^{\prime}}\in
T}\left(\|\bm{\tau^{\prime}}\|_{2}+1\right)\lvert
a_{\bm{\tau^{\prime}}}\rvert$
$\displaystyle\quad+\frac{\rho}{2\pi^{2}}\sum_{\bm{\tau}\in
T}\left(\|\bm{\tau}\|_{2}+1\right)\lvert a_{\bm{\tau}}\rvert$
$\displaystyle=\left(2a_{\bm{0}}+\frac{\rho}{2\pi^{2}}\right)\sum_{\bm{\tau^{\prime}}\in
T}\left(\|\bm{\tau^{\prime}}\|_{2}+1\right)\lvert
a_{\bm{\tau^{\prime}}}\rvert+\left(\sum_{\bm{\tau}\in
T}\left(\|\bm{\tau}\|_{2}+1\right)\lvert a_{\bm{\tau}}\rvert\right)^{2}.$

Applying the Cauchy–Schwarz inequality, denoting $t=|T|$, and separating the
$a_{\bm{0}}$ term, we obtain the bound

$\sum_{\bm{\tau}\in T}\left(\|\bm{\tau}\|_{2}+1\right)\lvert
a_{\bm{\tau}}\rvert\leq\sqrt{t}\left(\sum_{\bm{\tau}\in
T}\left(\|\bm{\tau}\|_{2}+1\right)^{2}\lvert
a_{\bm{\tau}}\rvert^{2}\right)^{1/2}\leq\sqrt{t}\left(\left\|a\right\|^{2}_{H^{1}}-a_{\bm{0}}^{2}\right)^{1/2}:=\beta.$
(44)

Combining the above inequalities yields

$\sum_{\bm{\mu}\in\Lambda\setminus\\{\bm{\nu}\\}}\lvert
G_{\bm{\nu}\bm{\mu}}\rvert\leq\left(2a_{\bm{0}}+\frac{\rho}{2\pi^{2}}\right)\beta+\beta^{2}.$
(45)

Finally, applying the Gershgorin circle theorem on $G$ combining (42) and (45)
and recalling that $\rho^{2}$ is an eigenvalue of $G$, we obtain the Riesz
constants

$\displaystyle b_{\Phi}$
$\displaystyle=a_{\bm{0}}^{2}-\left(2a_{\bm{0}}+\frac{\rho}{2\pi^{2}}\right)\beta-\beta^{2}>0,$
(46) $\displaystyle B_{\Phi}$
$\displaystyle=\left\|a\right\|_{H^{1}}^{2}+\frac{\rho^{2}}{16\pi^{4}}+\frac{a_{\bm{0}}\rho}{2\pi^{2}}+\left(2a_{\bm{0}}+\frac{\rho}{2\pi^{2}}\right)\beta+\beta^{2}.$
(47)

We note that the positivity of $b_{\Phi}$ is required for equation (39) to be
well-defined and this is assured by the sufficient condition (16).

##### Step 2: Boundedness.

To bound the essential supremum of the Riesz system, we use properties (34)
and (35)

$\displaystyle\|\Phi_{\bm{\nu}}\|_{L^{\infty}}$
$\displaystyle=\|\mathscr{L}\left[\Psi_{\bm{\nu}}\right]\|_{L^{\infty}}=\|-\nabla\cdot\left(a\nabla\Psi_{\bm{\nu}}\right)+\rho\Psi_{\bm{\nu}}\|_{L^{\infty}}$
$\displaystyle\leq\|\nabla
a\cdot\nabla\Psi_{\bm{\nu}}\|_{L^{\infty}}+\|a\Delta\Psi_{\bm{\nu}}\|_{L^{\infty}}+\|\rho\Psi_{\bm{\nu}}\|_{L^{\infty}}$
$\displaystyle\leq\left\|\sum_{\bm{\tau}\in T\cup\\{\bm{0}\\}}2\pi
i\bm{\tau}a_{\bm{\tau}}F_{\bm{\tau}}\right\|_{L^{\infty}}\cdot\left\|\frac{2\pi
i\bm{\nu}}{4\pi^{2}\|\bm{\nu}\|^{2}_{2}}F_{\bm{\nu}}\right\|_{L^{\infty}}+\left\|\sum_{\bm{\tau}\in
T\cup\\{\bm{0}\\}}a_{\bm{\tau}}F_{\bm{\tau}}\right\|_{L^{\infty}}+\frac{\rho}{4\pi^{2}\|\bm{\nu}\|^{2}_{2}}$
$\displaystyle\leq\sum_{\bm{\tau}\in
T\cup\\{\bm{0}\\}}\frac{\lvert\bm{\tau}\cdot\bm{\nu}\rvert}{\|\bm{\nu}\|^{2}_{2}}|a_{\bm{\tau}}|+\sum_{\bm{\tau}\in
T\cup\\{\bm{0}\\}}|a_{\bm{\tau}}|+\frac{\rho}{4\pi^{2}\|\bm{\nu}\|^{2}_{2}}$
$\displaystyle=\sum_{\bm{\tau}\in
T\cup\\{\bm{0}\\}}\left(\frac{\lvert\bm{\tau}\cdot\bm{\nu}\rvert}{\|\bm{\nu}\|^{2}_{2}}+1\right)|a_{\bm{\tau}}|+\frac{\rho}{4\pi^{2}\|\bm{\nu}\|^{2}_{2}}$

Finally, applying the Cauchy inequality, equation (5.1.2) and the fact that
$\|\bm{\nu}\|_{2}\geq 1$, we obtain

$\displaystyle\|\Phi_{\bm{\nu}}\|_{L^{\infty}}$
$\displaystyle\leq\sum_{\bm{\tau}\in
T\cup\\{\bm{0}\\}}\left(\frac{|\bm{\tau}\cdot\bm{\nu}|}{\|\bm{\nu}\|^{2}}+1\right)|a_{\bm{\tau}}|+\frac{\rho}{4\pi^{2}\|\bm{\nu}\|_{2}^{2}}$
$\displaystyle\leq a_{\bm{0}}+\sum_{\bm{\tau}\in
T}\left(\|\bm{\tau}\|+1\right)|a_{\bm{\tau}}|+\frac{\rho}{4\pi^{2}}\leq
a_{\bm{0}}+\beta+\frac{\rho}{4\pi^{2}}.$ (48)

This proves that $\\{\Phi_{\bm{\nu}}\\}_{\bm{\nu}\in\Lambda}$ is a bounded
Riesz system (Definition 4). Now, it remains to show accurate and stable
recovery guarantees for the problem (36). The machinery for this is enabled by
recent advances in sparse recovery for bounded Riesz systems (Brugiapaglia et
al., 2021a).

##### Step 3: Bounded Riesz property $\Longrightarrow$ error bound.

The strategy is to pick a lower bound on the sample complexity, such that $A$
satisfies the robust null-space property with high probability, and then
appeal to the recovery bounds already supplied by the literature surrounding
SR-LASSO. We begin with a simple re-normalization of the compressive Fourier
collocation matrix, by letting $\tilde{A}=A/\sqrt{B_{\phi}}$. We thus consider
a rescaled version of the SR-LASSO problem

$\hat{\bm{c}}\in\operatorname*{arg\,min}_{\bm{z}\in\mathbb{C}^{N}}\|A\bm{z}-\bm{b}\|_{2}+\lambda\|\bm{z}\|_{1},$
(49)

defined as

$\hat{\tilde{\bm{c}}}\in\operatorname*{arg\,min}_{\tilde{\bm{z}}\in\mathbb{C}^{N}}\|\tilde{A}\tilde{\bm{z}}-\bm{b}\|_{2}+\tilde{\lambda}\|\tilde{\bm{z}}\|_{1},$
(50)

where $\tilde{\bm{z}}=\sqrt{B_{\Phi}}\bm{z}$ and
$\tilde{\lambda}=\lambda/\sqrt{B_{\Phi}}$. The two minimizers are such that
$\sqrt{B_{\Phi}}\hat{\bm{c}}=\hat{\tilde{\bm{c}}}$. An inspection of the proof
of (Brugiapaglia et al., 2021a, Theorem 2.6) reveals that condition

$m\geq
c_{0}\left(\frac{\max\\{1,B_{\Phi}\\}}{b_{\Phi}}\right)^{2}K_{\Phi}^{2}s\log^{2}\left(sK_{\Phi}^{2}\frac{\max\\{1,B_{\Phi}\\}}{b_{\Phi}}\right)\log(eN),$
(51)

where $c_{0}>0$ is a universal constant, is sufficient for $\tilde{A}$ to
satisfy the rNSP with constants $\rho=1/2$ and $\gamma=2B_{\Phi}/b_{\Phi}$ and
with probability at least $1-\varepsilon/2$.111Observe that the probability of
failure (Brugiapaglia et al., 2021a, Theorem 2.6) is
$4\exp{-c_{0}^{\prime}(b_{\Phi}/\max\\{1,B_{\Phi}\\})^{2}m/(sK_{\Phi}^{2})}$,
for some constant $c_{0}^{\prime}>0$, and which is bounded by $\varepsilon$
thanks to (39).

Note that condition (51) is implied by (39) with

$\displaystyle
c^{(3)}_{a,\rho}=c_{0}\left(\frac{\max\\{1,B_{\Phi}\\}}{b_{\Phi}}\right)^{2}K_{\Phi}^{2}\quad\text{and}\quad
c^{(4)}_{a,\rho}=K_{\Phi}^{2}\frac{\max\\{1,B_{\Phi}\\}}{b_{\Phi}},$

and thanks to the following upper bound on the cardinality of the hyperbolic
cross obtained from Kühn et al. (2015); Chernov and Dũng (2016), see also Eq.
(2.10) from (Wang and Brugiapaglia, 2024):

$N=|\Lambda|\leq\min\\{4n^{5}16^{d},e^{2}n^{2+\log_{2}(d)}\\}.$ (52)

Hence, we can apply (Adcock et al., 2022b, Theorem 6.29) with $\bm{w}=1$ and
$\bm{b}=\tilde{A}\tilde{\bm{c}}_{\Lambda}+\bm{e}=A\bm{c}_{\Lambda}+\bm{e}$,
where

$\bm{e}=\frac{1}{\sqrt{m}}\left(\mathscr{L}[u-u_{\Lambda}](\bm{x}_{i})\right)_{i\in[m]}\in\mathbb{C}^{m},$

yielding the following coefficient recovery guarantee:

$\|\hat{\tilde{\bm{c}}}-\tilde{\bm{c}}_{\Lambda}\|_{2}\leq
c_{1}\frac{\sigma_{s}(\tilde{\bm{c}}_{\Lambda})_{1}}{\sqrt{s}}+\frac{1}{2}\left(\frac{c_{1}}{\sqrt{s}\tilde{\lambda}}+c_{2}\right)\|\bm{e}\|_{2},$
(53)

for constants $c_{1}=2(1+\rho)^{2}/(1-\rho)=9$ and
$c_{2}=2(3+\rho)\gamma/(1-\rho)=28B_{\Phi}/b_{\Phi}$, see also Eq. (6.15) from
(Adcock et al., 2022b). Note also that Eq. (6.47) (Adcock et al., 2022b)
allows us to pick a specific range for $\tilde{\lambda}$ and, consequently,
for $\lambda$:

$\sqrt{B_{\Phi}}\frac{d^{(1)}_{a,\rho}}{\sqrt{s}}\leq\lambda\leq\sqrt{B_{\Phi}}\frac{d^{(2)}_{a,\rho}}{\sqrt{s}},$
(54)

where $d^{(1)}_{a,\rho}$ and $d^{(2)}_{a,\rho}$ are such that
$0<d^{(1)}_{a,\rho}\leq
d^{(2)}_{a,\rho}\leq\frac{1+\rho}{(3+\rho)\gamma}=\frac{3b_{\Phi}}{14B_{\Phi}}$.
Finally, changing back to the original variables, and inserting values for
constants, we obtain

$\displaystyle\|\hat{\bm{c}}-\bm{c}_{\Lambda}\|_{2}$
$\displaystyle\leq\frac{c_{1}\sigma_{s}(\bm{c}_{\Lambda})_{1}}{\sqrt{s}}+\frac{1}{2\sqrt{B_{\Phi}}}\left(\frac{c_{1}\sqrt{B_{\Phi}}}{\sqrt{s}\lambda}+c_{2}\right)\|\bm{e}\|_{2}$
$\displaystyle=\frac{9\sigma_{s}(\bm{c}_{\Lambda})_{1}}{\sqrt{s}}+\frac{1}{2}\left(\frac{9}{\sqrt{s}\lambda}+\frac{28\sqrt{B_{\Phi}}}{b_{\Phi}}\right)\|\bm{e}\|_{2}$
$\displaystyle\leq\frac{9\sigma_{s}(\bm{c}_{\Lambda})_{1}}{\sqrt{s}}+\frac{1}{2}\left(\frac{9}{\sqrt{B_{\Phi}}d_{1}}+\frac{28\sqrt{B_{\Phi}}}{b_{\Phi}}\right)\|\bm{e}\|_{2}.$
(55)

Thanks to the rescaling used in (13), we have
$(-\Delta+\rho)\Psi_{\bm{\nu}}=F_{\bm{\nu}}$ for all
$\bm{\nu}\in\mathbb{Z}^{d}$. Hence, for any $v\in H^{2}(\mathbb{T}^{d})$ such
that $v=\sum_{\bm{\nu}\in\mathbb{Z}^{d}}d_{\bm{\nu}}\Psi_{\bm{\nu}}$, using
Parseval’s identity, we obtain

$\|(\Delta-\rho)v\|_{L^{2}}=\|\bm{d}\|_{2}\quad\text{and}\quad\|v\|_{L^{2}}=\|R\bm{d}\|_{2}\leq(a_{\bm{0}}/\rho)\|\bm{d}\|_{2},$
(56)

where $R=\text{diag}((r_{\bm{\nu}})_{\bm{\nu}\in\mathbb{Z}^{d}})$ with
$r_{\bm{\nu}}=1/(4\pi^{2}\|\bm{\nu}\|_{2}^{2}+\rho/a_{\bm{0}})$ being the
rescaling factors, and where we used the fact that $r_{\bm{\nu}}\leq
a_{\bm{0}}/\rho$ for every $\bm{\nu}\in\mathbb{Z}^{d}$. Using the triangle
inequality, we see that

$\displaystyle\|(\Delta-\rho)(u-\hat{u})\|_{L^{2}}$
$\displaystyle\leq\|(\Delta-\rho)(u-u_{\Lambda})\|_{L^{2}}+\|(\Delta-\rho)(u_{\Lambda}-\hat{u})\|_{L^{2}}$
$\displaystyle\|u-\hat{u}\|_{L^{2}}$
$\displaystyle\leq\|u-u_{\Lambda}\|_{L^{2}}+\|u_{\Lambda}-\hat{u}\|_{L^{2}}.$

These, combined with (5.1.2) and (56) (with $v=\hat{u}-u_{\Lambda}$ and
$\bm{d}=\hat{\bm{c}}-\bm{c}_{\Lambda}$), imply the following error bounds:

$\displaystyle\|(\Delta-\rho)(u-\hat{u})\|_{L^{2}}$
$\displaystyle\leq\|(\Delta-\rho)(u-u_{\Lambda})\|_{L^{2}}+\frac{9}{\sqrt{s}}\sigma_{s}(\bm{c}_{\Lambda})_{1}+\frac{1}{2}\left(\frac{9}{\sqrt{B_{\Phi}}d_{1}}+\frac{28\sqrt{B_{\Phi}}}{b_{\Phi}}\right)\|\bm{e}\|_{2}$
$\displaystyle\|u-\hat{u}\|_{L^{2}}$
$\displaystyle\leq\|u-u_{\Lambda}\|_{L^{2}}+\frac{a_{\bm{0}}}{\rho}\left(\frac{9}{\sqrt{s}}\sigma_{s}(\bm{c}_{\Lambda})_{1}+\frac{1}{2}\left(\frac{9}{\sqrt{B_{\Phi}}d_{1}}+\frac{28\sqrt{B_{\Phi}}}{b_{\Phi}}\right)\|\bm{e}\|_{2}\right).$

Summing these two inequalities and observing that, thanks to Lemma 7, for any
$v\in H^{2}(\mathbb{T}^{d})$,

$\|(\Delta-\rho)v\|_{L^{2}}+\|v\|_{L^{2}}\leq\|\Delta
u\|_{L^{2}}+(1+\rho)\|v\|_{L^{2}}\leq\sqrt{2}(1+\rho){\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}\leq\sqrt{2}(1+\rho)\|v\|_{H^{2}},$

we see that

$\displaystyle\|u-\hat{u}\|_{L^{2}}+\|(\Delta-\rho)(u-\hat{u})\|_{L^{2}}$
$\displaystyle\leq\sqrt{2}(1+\rho)\|u-u_{\Lambda}\|_{H^{2}}$
$\displaystyle\quad+\left(1+\frac{a_{\bm{0}}}{\rho}\right)\left(\frac{9}{\sqrt{s}}\sigma_{s}(\bm{c}_{\Lambda})_{1}+\frac{1}{2}\left(\frac{9}{\sqrt{B_{\Phi}}d_{1}}+\frac{28\sqrt{B_{\Phi}}}{b_{\Phi}}\right)\|\bm{e}\|_{2}\right).$
(57)

##### Step 4: Bounding the truncation error $\|\bm{e}\|_{2}$ and further
simplifications.

Following the discussion of (Wang and Brugiapaglia, 2024, page 9) and the
notation of (Adcock et al., 2022b, §7.6.1), we define

$\|\bm{e}\|_{2}=\sqrt{\frac{1}{m}\sum_{i=1}^{m}\left|\mathscr{L}[u-u_{\Lambda}](y_{i})\right|^{2}}=:E_{\Lambda,\text{disc}}(\mathscr{L}[u]).$

Now, (39) implies $m\geq 2s\log(4/\varepsilon)$, giving the following bound
for $E_{\Lambda,\text{disc}}(\mathscr{L}[u])$ with probability
$1-\varepsilon/2$:

$E_{\Lambda,\text{disc}}(\mathscr{L}[u])\leq\sqrt{2}\left(\frac{E_{\Lambda,\infty}(\mathscr{L}[u])}{\sqrt{s}}+E_{\Lambda,2}(\mathscr{L}[u])\right),$
(58)

thanks to (Adcock et al., 2022b, Lemma 7.11), where
$E_{\Lambda,\infty}(\mathscr{L}[u])=\|\mathscr{L}[u-u_{\Lambda}]\|_{L^{\infty}}$
and $E_{\Lambda,2}(\mathscr{L}[u])=\|\mathscr{L}[u-u_{\Lambda}]\|_{L^{2}}$.
Hence, (39) is sufficient so that both the rNSP and (58) hold with probability
$1-\varepsilon$ and we obtain

$\|\bm{e}\|_{2}\leq\sqrt{2}\left(\frac{\|\mathscr{L}[u-u_{\Lambda}]\|_{L^{\infty}}}{\sqrt{s}}+\|\mathscr{L}[u-u_{\Lambda}]\|_{L^{2}}\right),$
(59)

Now, we estimate the two error terms on the right-hand side. Note that, for
$1\leq p\leq\infty$, and any $v\in H^{2}(\mathbb{T}^{d})$,

$\displaystyle\|\mathscr{L}[v]\|_{L^{p}}$
$\displaystyle=\|-\nabla\cdot(a\nabla v)+\rho v\|_{L^{p}}$
$\displaystyle=\|-\nabla a\cdot\nabla v-a\Delta v+\rho v\|_{L^{p}}$
$\displaystyle\leq\|\nabla a\cdot\nabla v\|_{L^{p}}+\|a\Delta
v\|_{L^{p}}+\|\rho v\|_{L^{p}}.$ (60)

Using the inequality $\|\bm{v}\|_{1}\leq d^{1-\frac{1}{p}}\|\bm{v}\|_{p}$, for
all $\bm{v}\in\mathbb{R}^{d}$, we begin by analyzing the first two terms.
First, we see that

$\displaystyle\|\nabla a\cdot\nabla v\|_{L^{p}}$
$\displaystyle=\left\|\sum_{k=1}^{d}\frac{\partial a}{\partial
x_{k}}\frac{\partial v}{\partial
x_{k}}\right\|_{L^{p}}\leq\sum_{k=1}^{d}\left\|\frac{\partial a}{\partial
x_{k}}\frac{\partial v}{\partial
x_{k}}\right\|_{L^{p}}\leq\max_{k\in\left[d\right]}\left\|\frac{\partial
a}{\partial x_{k}}\right\|_{L^{\infty}}\sum_{k=1}^{d}\left\|\frac{\partial
v}{\partial x_{k}}\right\|_{L^{p}}$
$\displaystyle=\max_{k\in\left[d\right]}\left\|\frac{\partial a}{\partial
x_{k}}\right\|_{L^{\infty}}\left\|\left(\left\|\frac{\partial v}{\partial
x_{k}}\right\|_{L^{p}}\right)_{k=1}^{d}\right\|_{1}\leq\|a\|_{W^{1,\infty}}\cdot
d^{1-\frac{1}{p}}\left\|\left(\left\|\frac{\partial v}{\partial
x_{k}}\right\|_{L^{p}}\right)_{k=1}^{d}\right\|_{p}$ $\displaystyle\leq
d^{1-\frac{1}{p}}\|a\|_{W^{1,\infty}}\|v\|_{W^{2,p}}.$

Second, we have

$\displaystyle\|a\Delta v\|_{L^{p}}$
$\displaystyle\leq\|a\|_{L^{\infty}}\left\|\sum_{k=1}^{d}\frac{\partial^{2}v}{\partial
x_{k}^{p}}\right\|_{L^{p}}\leq\|a\|_{L^{\infty}}\sum_{k=1}^{d}\left\|\frac{\partial^{2}v}{\partial
x_{k}^{2}}\right\|_{L^{p}}=\|a\|_{L^{\infty}}\left\|\left(\left\|\frac{\partial^{2}v}{\partial
x_{k}^{2}}\right\|_{L^{p}}\right)_{k=1}^{d}\right\|_{1}$
$\displaystyle\leq\|a\|_{L^{\infty}}\cdot
d^{1-\frac{1}{p}}\left\|\left(\left\|\frac{\partial^{2}v}{\partial
x_{k}^{2}}\right\|_{L^{p}}\right)_{k=1}^{d}\right\|_{p}\leq
d^{1-\frac{1}{p}}\|a\|_{L^{\infty}}\|v\|_{W^{2,p}}.$

The third term is more easily estimated as $\|\rho
v\|_{L^{p}}\leq\|\rho\|_{L^{\infty}}\|v\|_{L^{p}}\leq\|\rho\|_{L^{\infty}}\|v\|_{W^{2,p}}.$
Altogether, when $p=\infty$, we have

$\|\mathscr{L}\left[v\right]\|_{L^{\infty}}\leq\|\nabla a\cdot\nabla
v\|_{L^{\infty}}+\|a\Delta v\|_{L^{\infty}}+\|\rho
v\|_{L^{\infty}}\leq\left(2d\|a\|_{W^{1,\infty}}+\|\rho\|_{L^{\infty}}\right)\|v\|_{W^{2,\infty}},$
(61)

and, for $p=2$,

$\|\mathscr{L}\left[v\right]\|_{L^{2}}\leq\|\nabla a\cdot\nabla
v\|_{L^{2}}+\|a\Delta v\|_{L^{2}}+\|\rho
v\|_{L^{2}}\leq\left(2\sqrt{d}\|a\|_{W^{1,\infty}}+\|\rho\|_{L^{\infty}}\right)\|v\|_{H^{2}}.$
(62)

Letting $v=u-u_{\Lambda}$ and combining (58) with the above estimates yields

$\|\bm{e}\|_{2}\leq\sqrt{2}\left(\|a\|_{W^{1,\infty}}+\|\rho\|_{L^{\infty}}\right)\left(\frac{\|u-u_{\Lambda}\|_{W^{2,\infty}}}{\sqrt{s}}+\|u-u_{\Lambda}\|_{H^{2}}\right).$

Plugging the above bound into (57) leads to the recovery guarantee (40), where

$\displaystyle C^{(1)}_{a,\rho}$
$\displaystyle=9\left(1+\frac{a_{\bm{0}}}{\rho}\right)$ $\displaystyle
C^{(2)}_{a,d,\rho}$
$\displaystyle=\frac{1}{2}\left(1+\frac{a_{\bm{0}}}{\rho}\right)\left(\frac{9}{\sqrt{B_{\Phi}}d_{1}}+\frac{28\sqrt{B_{\Phi}}}{b_{\Phi}}\right)\sqrt{2}\left(\|a\|_{W^{1,\infty}}+\|\rho\|_{L^{\infty}}+1\right)+\sqrt{2}(1+\rho).$

When $\rho<1$, the bound (41) can be obtained by invoking Lemma 7 again and
observing that, for every $v\in H^{2}(\mathbb{T}^{d})$,

$\|(\Delta-\rho)v\|_{L^{2}}+\|v\|_{L^{2}}\geq\|\Delta
v\|_{L^{2}}+(1-\rho)\|v\|_{L^{2}}\geq(1-\rho){\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}\geq(1-\rho)\sqrt{\frac{2}{3}}\|v\|_{H^{2}}.$

Combining this inequality with (57), we obtain (41), where

$C^{(3)}_{a,\rho}=\frac{C^{(1)}_{a,\rho}}{\sqrt{2/3}(1-\rho)}\quad\text{and}\quad
C^{(4)}_{a,d,\rho}=\frac{C^{(2)}_{a,d,\rho}}{\sqrt{2/3}(1-\rho)}.$

This concludes the proof.

###  5.2 Proof of Theorem 1

We proceed with a proof of the main result. For a summary of the general proof
strategy, we refer to the proof sketch after the statement of Theorem 1. Most
of the technical efforts will be devoted to constructing the class of periodic
PINNs, $\mathcal{N}$ and bounding the depth and width of its neural networks.
This will allow us to leverage the CFC convergence theory (Theorem 8) to
derive the desired error bounds (21)–(22).

##### Construction of the network class $\mathcal{N}$.

As explained in the proof sketch, we will construct networks
$\psi\in\mathcal{N}$ able to produce linear combinations of the form (14).
Letting $N=|\Lambda|$, this class is of the form

$\mathcal{N}=\left\\{\psi:\mathbb{R}^{d}\to\mathbb{C}:\psi(\bm{x})=\bm{z}^{\top}\psi_{\Lambda}(\bm{x}),\bm{z}\in\mathbb{C}^{N}\right\\},$
(63)

where $\psi_{\Lambda}(\bm{x})=(\psi_{\bm{\nu}}(\bm{x}))_{\bm{\nu}\in\Lambda}$
and $\psi_{\bm{\nu}}$ is a network replicating the $\bm{\nu}$th (rescaled)
Fourier function, i.e., $\psi_{\bm{\nu}}\equiv\Psi_{\bm{\nu}}$, with
$\Psi_{\bm{\nu}}$ defined as in (13). The main task is therefore to explicitly
construct the networks $\psi_{\bm{\nu}}$. We will do so by composing different
network modules, aimed at producing more and more complex functions of the
input variables $\bm{x}=(x_{j})_{j\in[d]}$, as follows:

  1. 1.

First, using a periodic layer $\bm{v}^{(2)}\circ\bm{q}^{(1)}$ with
$\bm{v}^{(2)}(\bm{x})=\bm{x}$ (see §2.1), we will generate univariate
trigonometric functions $\sin(2\pi x_{j})$ and $\cos(2\pi x_{j})$, for
$j\in[d]$.

  2. 2.

Then, through a network module $\bm{v}_{\text{trig}\to\text{multi-freq}}$, we
will compute univariate trignometric functions corresponding to all the
frequencies we need to generate functions in the system
$\\{\Psi_{\bm{\nu}}\\}_{\bm{\nu}\in\Lambda}$. These are given by $\cos(2\pi\nu
x_{j})$ and $\sin(2\pi\nu x_{j})$, for
$\nu\in\\{-\nu_{\max},\ldots,\nu_{\max}\\}$ and $j\in[d]$, where

$\nu_{\max}:=\max\\{|\nu_{j}|:\bm{\nu}=(\nu_{j})_{j=1}^{d}\in\Lambda\\}=n-1,$

is the maximum absolute frequency of multi-indices in the hyperbolc cross
$\Lambda=\Lambda^{\text{HC}}_{d,n}$.

  3. 3.

Finally, using a module $\bm{v}_{\text{multi-freq}\to\Lambda}$, we will
generate the rescaled Fourier functions $\Psi_{\bm{\nu}}(\bm{x})$ through
tensorization.

In summary, we have $\psi_{\Lambda}=\bm{v}_{\text{multi-
freq}\to\Lambda}\circ\bm{v}_{\text{trig}\to\text{multi-
freq}}\circ\bm{q}^{(1)}.$ After constructing $\psi_{\Lambda}$, we add a last
linear layer with weights $\bm{z}^{\top}\in\mathbb{R}^{1\times N}$ to generate
linear combinations of the form (14) (with $c_{\bm{\nu}}=z_{\bm{\nu}}$). Fig.
10 summarizes our network construction.

Figure 10: Summary of a generic network $\psi\in\mathcal{N}$.

A key ingredient of the proof is the fact that neural networks with RePU
activation can exactly replicate products of real numbers. Specifically, for
any $\ell,k\in\mathbb{N}$ with $\ell\geq 2$, (Adcock et al., 2023, Lemma 7.3),
which, in turn relies on techniques from (Opschoor et al., 2022; Schwab and
Zech, 2019), establishes the existence (along with an explicit construction)
of a feedforward neural network $P^{(k)}:\mathbb{R}^{k}\to\mathbb{R}$ with
$\mathrm{RePU}_{\ell}$ activation able to exactly reproduce the product of $k$
numbers, i.e.,

$P^{(k)}(x_{1},\ldots,x_{k})=\prod_{j=1}^{k}x_{j},\quad\forall\bm{x}=(x_{j})_{j=1}^{k}\in\mathbb{R}^{k}.$

Moreover, the width and depth of this network are such that

$\text{width}(P^{(k)})\leq p_{\ell}^{(1)}\cdot
k\quad\text{and}\quad\text{depth}(P^{(k)})\leq p^{(2)}\cdot\log_{2}(k),$ (64)

where $p_{\ell}^{(1)}>0$ depends on $\ell$ only and $p^{(2)}>0$ is a universal
constant. Note that we will not track the dependence of $P^{(k)}$ on $\ell$
since this parameter is assumed to be fixed throughout the proof. Before
showing in detail how to construct the various modules composing
$\psi_{\Lambda}$, we establish a connection between the training of periodic
PINNs in $\mathcal{N}$ and CFC approximation, which will yield the desired
recovery guarantees.

##### CFC convergence theory $\Rightarrow$ error bounds (21)–(22).

To apply Theorem 8, we leverage two important facts: (i) the spectral basis
(13) is exactly replicated by the network $\psi_{\Lambda}$ and (ii) the
weights of the last layer $\bm{z}^{\top}\in\mathbb{R}^{1\times N}$,
corresponding to the coefficients of linear combinations of the form (14), are
the only trainable parameters of $\psi$. In this setting, training the network
$\psi$ by minimizing the regularized RMSE loss in (20) with regularization
term

$\mathcal{R}(\psi)=\|\bm{z}\|_{1},$

is equivalent to solving the SR-LASSO problem (36) with $A$ and $\bm{b}$ as in
(24). Hence, we can simply apply Theorem 8 to obtain the desired conclusion.

The rest of the proof is devoted to illustrating the construction of the
periodic layer $\bm{q}^{(1)}\circ\bm{v}^{(2)}$ and the network modules
$\bm{v}_{\text{trig}\to\text{multi-freq}}$ and $\bm{v}_{\text{multi-
freq}\to\Lambda}$ in full detail.

##### Construction of $\psi_{\Lambda}$ (Step 1): the periodic layer.

Recalling the notation introduced in (4), we let $l=2$, with $\phi_{i1}=0$ and
$\phi_{i1}=-\pi/2$, for $i\in[d]$. This leads to

$\bm{x}\in\mathbb{R}^{d}\mapsto\bm{q}^{(1)}(\bm{x})=\begin{bmatrix}\cos(2\pi
x_{1})\\\ \cos(2\pi x_{1}-\pi/2)\\\ \vdots\\\ \cos(2\pi x_{d})\\\ \cos(2\pi
x_{d}-\pi/2)\end{bmatrix}=\begin{bmatrix}\cos(2\pi x_{1})\\\ \sin(2\pi
x_{1})\\\ \vdots\\\ \cos(2\pi x_{d})\\\ \sin(2\pi
x_{d})\end{bmatrix}=:\begin{bmatrix}c_{1}\\\ s_{1}\\\ \vdots\\\ c_{d}\\\
s_{d}\end{bmatrix}\in\mathbb{R}^{2d},$

where we have used the short-hand notation $c_{j}=\cos(2\pi x_{j})$ and
$s_{j}=\sin(2\pi x_{j})$. The second layer within the periodic layer is a
linear (identity) layer, i.e. $\bm{v}^{(2)}(\bm{x})=\bm{x}$. Hence, in
summary,

$\text{depth}(\bm{v}^{(2)}\circ\bm{q}^{(1)})=3\quad\text{and}\quad\text{width}(\bm{v}^{(2)}\circ\bm{q}^{(1)})=2d.$
(65)

##### Construction of $\psi_{\Lambda}$ (Step 2): the module
$\bm{v}_{\textnormal{trig}\to\textnormal{multi-freq}}$.

To generate sine and cosine functions at multiple frequencies given $c_{i}$,
$s_{i}$, we first generate powers of $c_{i}$ via an inner network submodule
module $\bm{v}_{\text{trig}\to\text{pow}_{i}}$ and then apply multi-angle
trigonometric formulas (see, e.g., Weisstein (2024)).

We start by presenting the inner module
$\bm{v}_{\text{trig}\rightarrow\text{pow}_{i}}$, which constructs the powers
of cosine for a fixed dimension $i$. This is given by

$\begin{bmatrix}c_{1}\\\ s_{1}\\\ \vdots\\\ c_{d}\\\
s_{d}\end{bmatrix}\xrightarrow[\text{restriction}]{\text{linear}}\begin{bmatrix}c_{i}\\\
s_{i}\end{bmatrix}\xrightarrow[\text{$\bm{x}\mapsto
W\bm{x}+\bm{b}$}]{\text{affine layer}}\begin{bmatrix}s_{i}\\\ 1\\\ c_{i}\\\
c_{i}\\\ \vdots\\\
c_{i}\end{bmatrix}\xrightarrow[]{\text{powers}}\begin{bmatrix}s_{i}\\\ 1\\\
c_{i}\\\ c_{i}^{2}\\\ \vdots\\\
c_{i}^{\nu_{\text{max}}}\end{bmatrix}:=\text{pow}_{i},$ (66)

where the weights and biases in $W\in\mathbb{R}^{(\nu_{\max}+2)\times 2}$ and
$\bm{b}\in\mathbb{R}^{\nu_{\max}+2}$ and the network corresponding to
$\xrightarrow{\text{powers}}$ are illustrated in Fig. 11. Note that this
network contains both RePU and linear activations, in accordance with (6).

Figure 11:  Part of the network submodule
$\bm{v}_{\text{trig}\to\text{pow}_{i}}$ defined in (66), used to generate
powers of cosine at dimension $i$.

Recalling (64) and observing that consecutive affine layers can be combined
together without altering depth (since the composition of affine maps is an
affine map), the depth and width of this module are bounded by

$\displaystyle\text{width}(\bm{v}_{\text{trig}\to\text{pow}_{i}})$
$\displaystyle\leq
3+\sum_{k=2}^{\nu_{\text{max}}}\text{width}\left(P^{(k)}\right)\leq
C_{1}\cdot\nu_{\text{max}}^{2}\cdot p_{\ell}^{(1)},$
$\displaystyle\text{depth}(\bm{v}_{\text{trig}\to\text{pow}_{i}})$
$\displaystyle\leq\text{depth}\left(P^{(\nu_{\text{max}})}\right)\leq
p^{(2)}\cdot\log_{2}(\nu_{\text{max}}),$

where $C_{1}>0$ is a universal constant.

Now, the output $\text{pow}_{i}$ of
$\bm{v}_{\text{trig}\to\text{pow}_{i}}(\bm{x})$ is transformed using multi-
angle formulae (see, e.g., Weisstein (2024)) to obtain univariate sine and
cosine functions at multiple frequencies. This corresponds to a linear layer

$\text{pow}_{i}=\begin{bmatrix}s_{i}\\\ 1\\\ c_{i}\\\ \vdots\\\
c_{i}^{\nu_{\text{max}}}\end{bmatrix}\xrightarrow[\text{formulae}]{\text{multi-
angle}}\begin{bmatrix}\cos(-\nu_{\text{max}}x_{i})\\\
\sin(-\nu_{\text{max}}x_{i})\\\ \vdots\\\ \cos(\nu_{\text{max}}x_{i})\\\
\sin(\nu_{\text{max}}x_{i})\end{bmatrix}.$

For the sake of completeness, we quickly review how to apply the multi-angle
formulae in this context. These formulae allow us to generate functions of the
form $\cos(2\pi\nu x)$ and $\sin(2\pi\nu x)$ from linear combinations of
$\sin(2\pi x)$ and powers of $\cos(2\pi x)$. In fact, for any
$\nu\in\\{-\nu_{\max},\ldots,\nu_{\max}\\}$ and $x\in\mathbb{R}$, we have

$\cos(2\pi\nu
x)=\sum_{i=0}^{\lfloor\nu/2\rfloor}\sum_{j=0}^{i}(-1)^{i-j}\binom{\nu}{2i}\binom{i}{j}[\cos(2\pi
x)]^{\nu-2(i-j)}=\sum_{k=\nu-\lfloor\nu/2\rfloor}^{\nu}a_{k}^{(\nu)}\cos^{k}(2\pi
x),$ (67)

for suitable values of $a_{k}^{(\nu)}$ and where we have used the reindexing
$k=\nu-2(i-j)$ in the second inequality and $\lfloor\cdot\rfloor$ denotes the
floor function. We can use the following analogous formula for the
construction of sine functions at multiple frequencies:

$\sin(2\nu\pi x)=\sin(2\pi
x)\cdot\sum_{i=0}^{\lfloor(\nu+1)/2\rfloor}\sum_{j=0}^{i}(-1)^{i-j}\binom{\nu}{2i+1}\binom{i}{j}[\cos(2\pi
x)]^{\nu-2(i-j)-1}.$

Note that the multi-angle formula for sines does not need powers of sines but
only powers of cosines. The above formula can be replicated by combining a
linear layer with a network $P^{(2)}$ performing the product.

Taking into account the linear layer corresponding to the multi-angle formulae
and the fact that the generation of functions $\cos(2\pi\nu x_{i})$ and
$\sin(2\pi\nu x_{i})$ must be repeated for each dimension $i\in[d]$, we obtain

$\displaystyle\text{width}(\bm{v}_{\text{trig}\to\text{multi-freq}})$
$\displaystyle=d\cdot\text{width}(\bm{v}_{\text{trig}\to\text{pow}_{i}})\leq
C_{1}\cdot p_{\ell}^{(1)}\cdot d\cdot\nu_{\text{max}}^{2},$
$\displaystyle\text{depth}(\bm{v}_{\text{trig}\to\text{multi-freq}})$
$\displaystyle=\text{depth}(\bm{v}_{\text{trig}\to\text{pow}_{i}})+\text{depth}(P^{(2)})\leq
C_{2}\cdot p^{(2)}\cdot\log_{2}(\nu_{\text{max}})$

for a universal constant $C_{2}>0$.

##### Construction of $\psi_{\Lambda}$ (Step 3): the module
$\bm{v}_{\textnormal{multi-freq}\to\Lambda}$.

The final module’s objective is to generate the rescaled Fourier functions
$\Psi_{\bm{\nu}}(\bm{x})$ in (13) given the output of
$\bm{v}_{\text{trig}\to\text{multi-freq}}(\bm{x})$. Recall that each Fourier
function is a product of complex numbers:

$F_{\bm{\nu}}(\bm{x})=\exp(2\pi
i\bm{\nu}\cdot\bm{x})=\prod_{j\in\text{supp}(\bm{\nu})}\exp(2\pi
i\nu_{j}x_{j})=\prod_{j\in\text{supp}(\bm{\nu})}(\cos(2\pi
i\nu_{j}x_{j})+i\sin(2\pi i\nu_{j}x_{j})).$

In particular, the generation of $F_{\bm{\nu}}(\bm{x})$ involves the product
of $\|\bm{\nu}\|_{0}$ complex numbers whose real and imaginary parts are given
by the outputs of $\bm{v}_{\text{trig}\to\text{multi-freq}}(\bm{x})$.
Moreover, since $\bm{\nu}\in\Lambda=\Lambda^{\text{HC}}_{d,n}$, we have

$n\geq\prod_{j=1}^{d}(|\nu_{j}|+1)=\prod_{j\in\text{supp}(\bm{\nu})}(|\nu_{j}|+1)\geq
2^{\|\bm{\nu}\|_{0}}\quad\Longrightarrow\quad\|\bm{\nu}\|_{0}\leq\log_{2}n.$
(68)

Therefore, computing each $F_{\bm{\nu}}(\bm{x})$ requires at most
$\min\\{d,\log_{2}n\\}$ products.

Now, we describe how to compute products of complex numbers using RePU
networks. We would like to compute $\prod_{j=1}^{k}z_{j}$ for
$z_{j}=x_{j}+iy_{j}\in\mathbb{C}$ for generic values
$x_{j},y_{j}\in\mathbb{R}$. To do so, we first find explicit formulas for the
real and imaginary parts of this product:

$\displaystyle\prod_{j=1}^{k}z_{j}$
$\displaystyle=\prod_{j=1}^{k}(x_{j}+iy_{j})=\sum_{\bm{j}\in\\{0,1\\}^{k}}\prod_{t=1}^{k}x_{t}^{j_{t}}(iy_{t})^{1-j_{t}}=\sum_{\bm{j}\in\\{0,1\\}^{k}}i^{k-\|\bm{j}\|_{1}}\prod_{t=1}^{k}x_{t}^{j_{t}}y_{t}^{1-j_{t}}$
$\displaystyle=\underbrace{\left(\sum_{\begin{subarray}{c}\bm{j}\in\\{0,1\\}^{k}\\\
k-\|\bm{j}\|_{1}\text{
even}\end{subarray}}(-1)^{\frac{k-\|\bm{j}\|_{1}}{2}}\prod_{t=1}^{k}x_{t}^{j_{t}}y_{t}^{1-j_{t}}\right)}_{=\operatorname*{Re}\left(\prod_{j=1}^{k}z_{j}\right)}+i\cdot\underbrace{\left(\sum_{\begin{subarray}{c}\bm{j}\in\\{0,1\\}^{k}\\\
k-\|\bm{j}\|_{1}\text{
odd}\end{subarray}}(-1)^{\frac{k-\|\bm{j}\|_{1}-1}{2}}\prod_{t=1}^{k}x_{t}^{j_{t}}y_{t}^{1-j_{t}}\right)}_{=\operatorname*{Im}\left(\prod_{j=1}^{k}z_{j}\right)}.$

This shows that $\operatorname*{Re}(\prod_{j=1}^{k}z_{j})$ and
$\operatorname*{Im}(\prod_{j=1}^{k}z_{j})$ can be computed as linear
combinations of $2^{k}$ products of $k$ real numbers. These correspond to the
terms $\prod_{t=1}^{k}x_{t}^{j_{t}}y_{t}^{1-j_{t}}$, where we observe that

$x_{t}^{j_{t}}y_{t}^{1-j_{t}}=\begin{cases}x_{t}&\text{if }j_{t}=1,\\\
y_{t}&\text{if }j_{t}=0.\end{cases}$

Hence, we can compute products of $k$ complex numbers $z_{1},\ldots,z_{k}$ via
a RePU network $P^{(k)}_{\mathbb{C}}:\mathbb{R}^{2k}\to\mathbb{C}$ that takes
the real and imaginary parts of the complex numbers $z_{j}$ as inputs and
constructed by stacking $2^{k}$ copies of $P^{(k)}$ atop, adding a linear
layer realizing the linear combinations above, and a final linear layer with
complex weights $\begin{bmatrix}1&i\end{bmatrix}$:

$\begin{bmatrix}x_{1}\\\ y_{1}\\\ \vdots\\\ x_{k}\\\
y_{k}\end{bmatrix}\xrightarrow[\text{of $P^{(k)}$}]{\text{$2^{k}$
copies}}\left[\prod_{t=1}^{k}x_{t}^{j_{t}}y_{t}^{1-j_{t}}\right]_{\bm{j}\in\\{0,1\\}^{k}}\xrightarrow[\text{layer}]{\text{linear}}\begin{bmatrix}\operatorname*{Re}\left(\prod_{j=1}^{k}z_{j}\right)\\\
\operatorname*{Im}\left(\prod_{j=1}^{k}z_{j}\right)\end{bmatrix}\xrightarrow[]{\begin{bmatrix}1&i\end{bmatrix}}\operatorname*{Re}\left(\prod_{j=1}^{k}z_{j}\right)+i\operatorname*{Im}\left(\prod_{j=1}^{k}z_{j}\right).$

Therefore, recalling (64), we obtain

$\displaystyle\text{width}(P^{(k)}_{\mathbb{C}})$
$\displaystyle=2^{k}\cdot\text{width}(P^{(k)})\leq C_{3}\cdot
p^{(1)}_{\ell}\cdot 2^{k}\cdot k$
$\displaystyle\text{depth}(P^{(k)}_{\mathbb{C}})$
$\displaystyle=\text{depth}(P^{(k)})+1\leq C_{4}\cdot
p^{(2)}\cdot\log_{2}{k},$

for come universal constants $C_{3},C_{4}>0$. To conclude, the module
$\bm{v}_{\textnormal{multi-freq}\to\Lambda}$ is such that

$\displaystyle\bm{v}_{\textnormal{multi-
freq}\to\Lambda}\left(\begin{bmatrix}\cos(2\pi\nu x_{j})\\\ \sin(2\pi\nu
x_{j})\end{bmatrix}_{j\in[d],\nu\in\\{-\nu_{\max},\ldots,\nu_{\max}\\}}\right)$
$\displaystyle=\left[r_{\bm{\nu}}\cdot
P^{(\|\bm{\nu}\|_{0})}_{\mathbb{C}}\left(\begin{bmatrix}\cos(2\pi\nu_{j}x_{j})\\\
\sin(2\pi\nu_{j}x_{j})\end{bmatrix}_{j\in\text{supp}(\bm{\nu})}\right)\right]_{\bm{\nu}\in\Lambda}$
$\displaystyle=r_{\bm{\nu}}F_{\bm{\nu}}(\bm{x}),$

where $r_{\bm{\nu}}=1/(4\pi^{2}\|\bm{\nu}\|_{2}^{2}+\rho/a_{\bm{0}})$ as
desired. Note that the networks $P^{(\|\bm{\nu}\|_{0})}_{\mathbb{C}}$, with
$\bm{\nu}\in\Lambda$, can be stacked using a construction analogous to that of
Fig. 11. Also observe that for $\bm{\nu}=\bm{0}$ we simply let
$P^{(0)}_{\mathbb{C}}\equiv 1$. In addition, the multiplication by the
rescaling factor $r_{\bm{\nu}}$ can be realized by modifying weights in the
last layer (we don’t need to add a $P^{(2)}$ module because $r_{\bm{\nu}}$ is
independent of $\bm{x}$). In summary, recalling that $N=|\Lambda|$ and using
(68) and the depth and width bounds for $P^{(n)}_{\mathbb{C}}$, the
architecture bounds for the whole network module are

$\displaystyle\text{width}(\bm{v}_{\textnormal{multi-freq}\to\Lambda})$
$\displaystyle=\sum_{\bm{\nu}\in\Lambda}\text{width}(P^{(\|\bm{\nu}\|_{0})}_{\mathbb{C}})\leq
N\cdot\text{width}(P^{(\min\\{d,\log_{2}n\\})}_{\mathbb{C}})$
$\displaystyle\leq p^{(1)}_{\ell}\cdot C_{3}\cdot
N\cdot\min\\{2^{d},n\\}\cdot\min\\{d,\log_{2}n\\}$
$\displaystyle\text{depth}(\bm{v}_{\textnormal{multi-freq}\to\Lambda})$
$\displaystyle=\max_{\bm{\nu}\in\Lambda}\text{depth}(P^{(\|\bm{\nu}\|_{0})}_{\mathbb{C}})=\text{depth}(P^{(\min\\{d,2^{n}\\})}_{\mathbb{C}})\leq
C_{4}\cdot p^{(2)}\cdot\log_{2}(\min\\{d,2^{n}\\}).$

##### Final architecture bounds for $\psi$.

Combining the architecture bounds obtained for each network module in Steps 1,
2 and 3, and recalling that
$\psi(\bm{x})=\bm{z}^{\top}\psi_{\Lambda}(\bm{x})$, we obtain

$\displaystyle\text{width}\left(\psi\right)$
$\displaystyle=\max\left\\{\text{width}(\bm{v}^{(2)}\circ\bm{q}^{(1)}),\text{width}\left(\bm{v}_{\text{trig}\rightarrow\text{multi-
freq}}\right),\text{width}\left(\bm{v}_{\text{multi-
freq}\rightarrow\Lambda}\right)\right\\}$
$\displaystyle\leq\max\left\\{2d,\;C_{1}\cdot p_{\ell}^{(1)}\cdot
d\cdot\nu_{\text{max}}^{2},\;p^{(1)}_{\ell}\cdot C_{3}\cdot
N\cdot\min\\{2^{d},n\\}\cdot\min\\{d,\log_{2}n\\}\right\\}$
$\displaystyle\leq\min\left\\{4n^{5}16^{d},e^{2}n^{2+\log_{2}d}\right\\}\cdot
c_{(\ell)}^{(1)}\cdot d.$

Using the hyperbolic cross cardinality bound (52), we see that

$\text{width}\left(\psi\right)\leq C_{5}\cdot
p_{\ell}^{(1)}\cdot\min\left\\{4n^{5}16^{d},e^{2}n^{2+\log_{2}d}\right\\}\cdot
d\cdot\min\\{2^{d},n\\},$

for some universal constant $C_{5}>0$. Moreover,

$\displaystyle\text{depth}\left(\psi\right)$
$\displaystyle=\text{depth}(\bm{v}^{(2)}\circ\bm{q}^{(1)})+\text{depth}\left(\bm{v}_{\text{trig}\rightarrow\text{multi-
freq}}\right)+\text{depth}\left(\bm{v}_{\text{multi-
freq}\rightarrow\Lambda}\right)+1$ $\displaystyle\leq C_{6}\cdot
p^{(2)}\cdot\left(\log_{2}(n)+\min\\{\log_{2}d,n\\}\right),$

where $C_{6}$ is a universal constant. Letting $c_{\ell}^{(1)}=C_{5}\cdot
p_{\ell}^{(1)}$ and $c^{(2)}=C_{6}\cdot p^{(2)}$ yields (17) and (18) and
concludes the proof. $\blacksquare$

##  6 Conclusions and open problems

We have shown a new convergence result for PINNs (Theorem 1) in the form of a
practical existence theorem for the numerical solution of high-dimensional,
periodic diffusion-reaction problems. This result establishes the existence of
a class of periodic PINNs able to achieve the same accuracy as a sparse
approximation method (namely, CFC) and using a number of training samples that
scales only logarithmically or, at worst, linearly with the PDE domain’s
dimension $d$. The mild scaling of the sample complexity with respect to $d$
is numerically confirmed through experimentation (see Fig. 5). Our practical
existence theorem relies on a new CFC convergence result for diffusion-
reaction problems (Theorem 8) and an explicit construction of periodic PINNs
able to replicate the Fourier basis. We have also experimentally confirmed the
robustness of periodic PINNs to solve high-dimensional PDEs with respect to
the network hyperparameters. Finally, we compared periodic PINNs with CFC,
showing that the latter could achieve much higher (but, sometimes, worse)
accuracy than the former depending on the sparsity properties of the PDE
solution. On the other hand, the performance of periodic PINNs is numerically
observed to be consistent across the three examples considered (see Fig. 9).

We conclude by mentioning some gaps between theory and practice and open
problem for future research. First, Theorem 1 relies on RePU or linear
activations and not on more standard ones such as ReLU or tanh. The main
reason to work with RePU activations is that they allow to exactly replicate
products and, hence, using the construction illustrated in §5.2, Fourier
functions. The argument of Theorem 1 could be generalized to more general
activations, such as ReLU or tanh. In that case, however, products (and,
hence, Fourier basis functions) could only be approximated and not exactly
replicated. This issue could be handled, see (Adcock et al., 2024) and
references therein, but it will introduce nontrivial technical difficulties
due to the fact that one would have to deal with an approximate CFC matrix,
whose error with respect to the true CFC matrix should be carefully controlled
in the analysis.

Another important aspect is the presence of sufficient conditions (15) and
(16) on the PDE coefficients $a$ and $\rho$ in Theorem 1. These are inherited
by the CFC convergence analysis (Theorem 8) and are an artifact of its proof.
It was shown in the diffusion equation case that these conditions are
sufficient but far from being necessary, see (Wang and Brugiapaglia, 2024,
§4.4). This is likely to be the case in the diffusion-reaction case as well
since the convergence analysis is based on the same argument.

The optimal network $\hat{\psi}$ of Theorem 1 is assumed to be trained by
_exactly_ minimizing a regularized RMSE loss. Of course, this is not what
happens in practice, where the loss is usually only approximately minimized
using a stochastic gradient descent method. Taking into account the error
introduced by the training algorithm in the analysis is an important open
question. Moreover, only the last layer of $\hat{\psi}$ is trained, whereas
the previous layers are explicitly constructed. This gap between theory and
practice is an intrinsic limitation of the argument that Theorem 1 relies on.

Hence, differences exist between the theoretical setting of Theorem 1 and the
practical implementation of periodical PINNs in §4, for which we have followed
a setup closer to what is commonly employed in the literature. These
differences include the choice of activation, the use of complex-valued
weights and the presence of a regularization term in the loss function. It
would be interesting to study whether a numerical implementation that follows
the setting of Theorem 1 more closely would give any practical benefit over