several frame families. ETF clearly outperforms all other frames at its
optimal $\beta$ operation point. Figure 9.1-B shows the “practical capacity”
(a quantity related to the concept of “lattice decoding” [urbanke1998lattice,
zamir2014lattice]), which is equal to (9.3) without the “$I_{k}$” term inside
the logarithm. This quantity - similar to the Shannon-transform performance
measure in (3.8) - is more sensitive to eigenvalues of $F_{S}^{\dagger}F_{S}$
near zero; it thus demonstrates more clearly the gaps between the different
frames and the optimal $\beta$ operation point for each frame [MayaProject].

A similar behavior was observed in “irregular pulse-shape” optimization for
faster-than-Nyquist TDMA signaling with FFE-DFE equalization [DadushProject].

Figure 9.1: Average capacities per resource of a NOMA/CDMA system for a fixed
activity level $p=\gamma\beta=k/n=0.5$, as a function of the loading factor
$\beta=k/m$, for several frames: random (5.3), LPF (2.7), Steiner ETF
[fickus2012steiner], and sparse [zaidel2018sparse]. All frame dimensions vary
around $100\times 200$. (Upper) regular capacity (9.3). (Lower) practical
capacity. For $\beta^{-1}\approx 1$, the capacity gain of ETF is $\sim 0.4$
bit with respect to random frame, and almost 2 bits with respect to LPF (for
practical capacity).

###  9.2 DFT codes

Linear transformation can also be used to introduce redundancy against fading
(or erasures or impulses) plus noise. Say, $m$ information symbols ${\bf
x}=(x_{1},\ldots,x_{m})^{\dagger}$ are transformed into a larger $n$-vector
signal $F^{\dagger}{\bf x}$, and the receiver gets a random $k$ subset $S$
corrupted with noise:

${\bf y}=F_{S}^{\dagger}{\bf x}+{\bf z}$ (9.4)

where ${\bf y}=(y_{1},\ldots,y_{k})^{\dagger}$ is the erased receiver input,
and $F^{\dagger}$ plays the role of the “coding” transformation. Note that in
contrast to the partially-active multiple access channel (9.2), here
$F^{\dagger}$ expands the dimension $m:n$, and the interesting regime where
the channel can be inverted is $n\geq k\geq m$. Wolf [wolf1983redundancy]
proposed the concept of DFT codes, where $F$ is a LPF frame, as a means to
mitigate the effect of impulse noise which acts similarly to the erasures in
(9.4). The redundancy introduced by the frame $F$ allows the receiver to
estimate the information symbols ${\bf x}$ after the $n-k$ impulses were
detected and removed (erased). A later work by Marshall [marshall1984coding]
hinted to the fact that irregular selection of DFT frequencies leads to better
estimation performance. See also [tulino2007gaussian, tulino2010capacity,
zanko2014iterative].

###  9.3 Space-time codes (STC)

A space-time code can be viewed as a two-dimensional - space/time - analog
modulation scheme, tailored to a vector-input block-fading channel, that
introduces diversity against random fading in the spatial dimension. Codes for
fading Gaussian channels need to preserve a sufficient Hamming distance after
the coordinates were randomly attenuated by fading [tse2005fundamentals].
Rayleigh fading gives rise to the “product measure” [tse2005fundamentals, sec.
3.2.2], which becomes the “determinant criterion” for STC [tarokh1998space].
More specifically, STC transforms a message $A$ into an $m\times n$ codword
$F_{A}$ of $m$ channel uses in time over $n$ channel uses in space ($n$
transmit antennas). The pairwise error between two such codewords $A$ and $B$
in the presence of Rayleigh fading is bounded by some (SNR-dependent) constant
divided by the determinant $\det[(F_{A}-F_{B})\cdot(F_{A}-F_{B})^{\dagger}]$,
[tarokh1998space], [tse2005fundamentals, sec. 3.3.2]. If we take $B$ to be the
zero codeword, then the goodness of the codeword $A$ is inverse proportional
to $\det[F_{A}\cdot F_{A}^{\dagger}]$. In the presence of 0/1 spatial fading,
the equivalent codword is $(F_{A})_{S}$, where $S\subset[n]$ is the set of
surviving spatial coordinates. If a random set of $k$ coordinates survive,
where $k\geq m$, then we get that the error probability in the presence of a
combination of 0/1 and Rayleigh fading is bounded by

$P_{e}(A\rightarrow 0)\leq\frac{1}{{n\choose
k}}\sum_{S:|S|=k}\frac{SNR^{-m}}{\det[(F_{A})_{S}(F_{A})_{S}^{\dagger}]}.$
(9.5)

The RHS of (9.5) can be seen as a $k$-dependent criterion for STC design,
which reduces to the usual “determinant criterion” for $k=n$ ($p=1$). Figure
9.2 illustrates this measure as a function of $p=k/n$ for several “frames”
$F_{A}$. The clear advantage of ETFs for moderate $p$ gives rise to the
question: can we choose multiple space-time codewords, such that all pairwise
differences will be “near ETF”?

Figure 9.2: The “inverse subset determinant” bound (9.5) of space-time codes,
for a fixed $\gamma=0.5$ as a function of the not-erased probability $p=k/n$,
for several frames: random (5.3), LPF (2.7), Paley ETF [fickus2012steiner],
and sparse [zaidel2018sparse]. Frame dimensions $m\times n=$: (A) $10\times
20$. (B) $85\times 170$.

###  9.4 Over-sampling and $\Delta$-$\Sigma$ modulation

Sampling of a band-limited signal beyond its Nyquist rate allows to trade off
quantization resolution (distortion) for time resolution (rate)
[zamir1995rate]. The distortion enhancement follows since the desired signal
has a low-pass spectrum in the over-sampled domain, while the quantization
noise is approximately white and can be filtered out at the reconstruction.
$\Delta$-$\Sigma$ modulation enhances this effect using a “noise shaping”
feedback filter around the quantizer, that creates a high-pass quantization
noise spectrum, hence lower reconstruction distortion for the same
quantization resolution [candy1992oversampling, ordentlich2019performance].
Interestingly, it was recognized that the mapping between the original
(Nyquist-sampled) signal to the $\Delta$-$\Sigma$ modulated signal can be
viewed as a frame expansion [benedetto2004sigma]. This opens the door to more
general $\Delta$-$\Sigma$ modulation schemes, that are robust against various
forms of interference and erasures [goyal2001quantized,
strohmer2003grassmannian, eldar2003geometrically, boufounos2008causal].

###  9.5 Multiple-description (MD) source coding

When a source is transmitted over a lossy packet network, where ACK/NAK
feedback cannot be assumed, one may wish for graceful degradation as a
function of the number of erased/received packets. [gamal1982achievable,
goyal2001multiple]. In MD coding, the source is split into $n$ “descriptions”,
where a priori unknown subset of the descriptions arrive to the receiver; it
is convenient to assume that the receiver gets either a subset of size $k$, or
all the $n$ descriptions [pradhan2004n].

A simple solution based on $\Delta$-$\Sigma$ modulation for the case ($n=2$,
$k=1$) was proposed in [ostergaard2009multiple]. The $\Delta$-$\Sigma$-two-
description coding scheme uses 1:2 oversampling ratio, where the two
descriptions are created by splitting into even and odd samples, and where the
noise shaping loop controls the distortion ratio between the $k=1$ case (even
or odd samples) and the $k=2$ case (all samples). It is shown in
[ostergaard2009multiple] that this solution achieves the optimal MD rate-
distortion region for a Gaussian source under quadratic distortion if we use
an optimal dithered (lattice) quantizer [ozarow1980source]. Naturally, one may
wish to extend this solution to a general ($n,k$) pair using $m:n$
oversampling ratio, for some $m\leq k$, with periodic allocation of the
samples to the $n$ descriptions. This extension suffers, however, from noise
amplification whenever the received descriptions do not form a uniform
sampling pattern [mashiach2010multiple, mashiach2013sampling], which is
exactly the same as the signal amplification problem in band-limited
interpolation for “source coding with erasures” (Chapter 2).

Figure 9.3: Rate-distortion performance of $\Delta$-$\Sigma$ MD coding with
over-sampling ratio $n:m=31:6$, where $k=6$ descriptions are received out of
$n=31$ transmitted descriptions, for several interpolation filters: low pass
(LP), difference-set (DS), uniform (Unif), and random (rand) spectrum
interpolation filter, compared to the “PPR bound” [pradhan2004n].

The fact that ETFs - which are highly robust against random erasures -
correspond to “irregular spectrum” (see the difference-set DFT construction in
Section 4.1.2), inspired the idea to modify the interpolation filter in the MD
$\Delta$-$\Sigma$ modulation scheme to be an “irregular interpolator”
[goren2021diversity]. Figure 9.3, taken from [goren2021diversity],
demonstrates this idea in MD image (JPEG-based) coding. See also
[palgy2010multiple].

###  9.6 Coded distributed computation (CDC)

Analog MDS (Reed-Solomon) codes were recently proposed to protect against
straggling worker nodes in distributed computation of a linear function; e.g.,
of the gradient in training a deep neural network [lee2018speeding,
raviv2020gradient]. Specifically, to obtain the matrix multiplication $XA$
with $K$ ideal workers, the master node breaks $A$ into $K$ pieces
$[A_{1},...,A_{K}]$, and distribute the partial multiplications
$XA_{1},...,XA_{K}$. If there are more workers but some are stragglers, say,
only $K$ out of $N$ compute their task in time, the master worker can use the
generator matrix of a $K:N$ “analog erasure-correction code” to expand
$A_{1},...,A_{K}$ into $\tilde{A}_{1},...,\tilde{A}_{N}$, so $XA$ can be
computed from any $K$ subset of $XA_{1},...,XA_{N}$. For example, for
$K=2,N=3$ the code is simply $\tilde{A}_{1}=A_{1},\tilde{A}_{2}=A_{2}$ and
$\tilde{A}_{3}=A_{1}+A_{2}$, and analog versions of BCH and Reed Solomon codes
provide simple encoding and decoding schemes for more general $(K,N)$ pairs
[lee2018speeding, li2020coded, raviv2020gradient]. It is insightful to think
of the expansion as multiplying $A$ by a frame $F$ with an aspect ratio $N/K$,
i.e., $\tilde{A}=AF$, where $F$ is composed of $K\times N$ square blocks, each
being a scaled identity matrix $g_{i,j}I$, for $i\in[K]$, $j\in[N]$. For
example, in the 2:3 example above all the $g$’s are 0/1, and

$F=\left[\begin{array}[]{ccc}I&0&I\\\ 0&I&I\end{array}\right].$

Now suppose that beyond the stragglers problem, we also face noise due to
finite word-length computation. That is, the master node needs to reconstruct
$XA$ from

$Y_{S}=(XAF)_{S}+Z$

where $S$ denotes the subset of non straggling workers and $Z$ denotes the
computation noise. Our goal is to find coding frames $F$ which are resilient
to noise, i.e., numerically stable [fahim2019numerically,
ramamoorthy2019numerically]. The challenge though is to keep the simplicity of
encoding and decoding of analog codes inspired from digital erasure correction
codes [babadi2011spectral, calderbank2010construction, raviv2020gradient], yet
enjoy the better robustness to analog noise of ETFs.

Figure 9.4: Comparison of noise amplification Vs. the frame redundancy
$\gamma^{-1}$ for several polynomial codes: uniform sampled unit circle
(USPC), non-uniform sampled unit circle (NUSPC), and non-consecutive
(difference-set) powers, compared to the MANOVA curve. Figure is taken from
[RoyeeYosibashPaper].

###  9.7 Neural networks

The robustness of ETF subsets, revealed in [haikin2016analog,
haikin2017random], inspired the work of Bank and Girayes
[bank2018relationship, bank2020ETF] on training of deep neural networks with
Dropout. Their scheme regularizes the training of the parameters in each
iteration, so that they form a set of vectors with a low coherence, just like
a subset from an ETF.

###  9.8 Frames in compressed sensing: phase transition

Pose the compressed sensing measurement as $x\mapsto Ax=y$. We observe $y$ and
wish to recover $x$ for $A$ known. Denote our chosen signal reconstruction
scheme by $\mathcal{R}:t\mapsto x$. For example, $\mathcal{R}$ can be
$\mathcal{R}:y\mapsto\text{argmin}_{y=Ax}\left|\left|x\right|\right|_{1}$. The
basic question in (noiseless) compressed sensing is - for which
(matrix,signal)  pair $(A,x)$ do we have $\mathcal{R}(y)=x$? When the
underlying model includes randomness in the pair $(A,x)$, one can ask this
question either in the worst-case sense (asking to characterize those cases
$(A,x)$ when $\mathcal{R}(y)=x$ w.p 1) or average-case sense (asking to
characterize those cases cases $(A,x)$ when $\mathcal{R}(y)=x$ with high
probability).

The Restricted Isometry Property became popular as it appeared in sufficient
conditions for worst-case recovery in compressed sensing problems. However,
these sufficient conditions are known to be far from tight, and indeed explain
only a fraction of the cases $(A,x)$ in which recovery (using $\ell_{1}$
minimization, say) is known to hold.

The situation is different when one considers average-case recovery. Donoho
and Tanner [donoho2009observed] developed a theory that places cases $(A,x)$
on a sparsity-undersampling phase plane. Using large-scale experiments they
have identified - empiricially at first - an asymptotic phase transition
phenomenon, which holds approximately for finite $(A,x)$ and exactly in the
large-system limit. An exact formula for the asymptotic phase transition was
later identified for Gaussian i.i.d sensing matrices $A$ and a variety of
signal models $x$. In a definite sense, then, the average-case analysis is
satisfying and complete.

In [monajemi2013deterministic] the authors reported a curious phenomenon: the
same empirical apparatus that had been previously used to establish the
existence of the Donoho-Tanner phase transition for the case of Gaussian
sensning matrices, has been applied to deterministic sensing matrices based on
known frame constructions including ETFs. The authors reported overwhelming
empirical evidence, based on a very large parallel computation, that most of
the deterministic frames under study, when used as the sensing matrix $A$ in a
recovery problem and when using a variety of reconstruction algorithms,
exhibit a phase transition that agrees with that of Gaussian matrices.

This observed universality, in which many different frame constructions
exhibit a single, universal compressed sensing phase transition has not, to
the best of our knowledge, proved rigorously. However in contrast with theory
based on Restricted Isometry Property it provides a complete characterization,
within a specified framework, of the success and failure of recovery of sparse
vectors using $\ell_{1}$ minimization and related algorithms. The connection
between this universal phenomenon and the new ETF-MANOVA relation is yet to be
found.

##  Chapter 10 ETF optimality conjecture

The examples of Chapter 9 (see also Figure 1.1) showed that ETFs are superior
to all other frame families that were investigated: LPF, random and code-
based, in terms of their performance as analog codes. We are tempted to
conjecture that ETFs outperform all other frames, and achieve the optimal
frame performance characteristics (3.3), (3.9) and (3.11) defined in Chapter
3. In this chapter we provide some analytic support for the “ETF superiority”
conjecture. We show that the Welch lower bound (3.1) naturally extends to
random frame subsets and to higher-order moments (3rd and 4th moments in
(7.4)), and that the extended bound is achieved by (and sometimes only by)
ETFs.

###  10.1 Erasure Welch bound

######  Theorem 10.1.1 (Erasure Welch Bound (EWB) of order $r=2,3,4$
[haikin2018frame, MarinaThesis].).

For any $m$-by-$n$ unit-norm frame $F$, selection probability $p$, and moment
orders $r=2,3,4$, the expected sub-frame $r$th moment (7.4) is lower bounded
by

$m_{r}\geq m_{r}^{\rm ETF}$ (10.1)

where $m_{r}^{\rm ETF}$ is the corresponding value for an ETF given by Theorem
7.1.1 (using the expressions in (7.7)-(7.8), even if an ETF does not exist for
the pair $(m,n)$). Equality holds for $r=2,3$ if and only if $F$ is a tight
frame, and for $r=4$ if and only if $F$ is an ETF.

The case ($r=2,p=1$) in Theorem 10.1.1 is equivalent to the (mean-squared)
Welch bound (3.1). Thus, the EWB generalizes the Welch bound in two senses.
First, it is a bound on average of subsets of $F$ rather than on the entire
set of its pairwise moments. Second, it is a bound on higher order ($r=3,4$)
moments.111Our definition is different than the generalized Welch bound on the
average $r$th power of the absolute cross-correlation [welch1974lower].

######  Proof 10.1.2.

For the case $r=2$, the equivalent form (7.14) of the mean-squared Welch bound
(3.1), and the expression (7.13b) for $m_{2}$, imply

$\begin{split}&m_{2}\geq p+p^{2}x=m^{\rm MANOVA}(\gamma,p,2)=m_{2}^{\rm
ETF},\end{split}$ (10.2)

where the last two equations follow from (7.7b) and Theorem 7.1.1. The
equality condition for the Welch bound implies that (10.2) holds with equality
for any tight frame. The proof for $r=3,4$ is based on Jensen’s inequality,
and it appears in [haikin2018frame, MarinaThesis].

###  10.2 Minimum empirical variance and kurtosis

The EWB has an interesting implication on ETF having the most compact spectral
distribution of a randomly selected sub-frame. By (7.12) and Theorem 10.1.1,
all unit-norm frames share the same expected sub-frame 1st moment $m_{1}=p$;
thus minimizing the expected sub-frame 2nd moment $m_{2}$ amounts to
minimizing the empirical spectral variance

${\rm ESV}=m_{2}-m_{1}^{2}$ (10.3)

of a typical sub-frame (typical in the sense that its moment matches the
ensemble average), among all unit-norm frames. Moreover, all unit-norm tight
frames share the same expected sub-frame 2nd moment $m_{2}=p+p^{2}x$, and 3rd
moment $m_{3}=p+3xp^{2}+(x^{2}-x)p^{3}$; thus minimizing the expected sub-
frame 4th moment $m_{4}$ amounts to minimizing the empirical spectral kurtosis

${\rm
ESK}=\frac{m_{4}-4m_{3}m_{1}+6m_{2}m_{1}^{2}-3m_{1}^{4}}{(m_{2}-m_{1}^{2})^{2}}$
(10.4)

of a typical sub-frame, among all unit-norm tight frames. We can thus use the
MANOVA $1st$-to-$4th$ moments (7.7) to lower bound the variance and kurtosis
of the limiting sub-frame ESD of any sequence of unit-norm (tight) frames.

######  Theorem 10.2.1 (Minimum sub-frame ESV and ESK).

Consider any sequence of unit-norm frames with increasing dimensions, whose
aspect ratio converges to $\gamma$, and whose (Bernoulli-$p$ selected) sub-
frame ESD $\mu_{G_{S}}^{(n)}$ converges to some limiting law $\mu$ (see
(7.30)). The asymptotic sub-frame ESV (10.3) is then lower bounded by that of
a tight frame:

${\rm ESV}\geq p+(x-1)p^{2}$ (10.5)

where $x=1/\gamma-1$. Moreover, if the sequence of frames is tight, then the
asymptotic sub-frame ESK (10.4) is lower bounded by that of an ETF:

${\rm
ESK}\geq\frac{p+p^{2}(6x-4)+p^{3}(6x^{2}-16x+6)+p^{4}(x^{3}-7x^{2}+11x-3)}{p^{2}(x^{4}+2x^{2}+1)-p^{3}(2x^{2}+2)+p^{4}}.$
(10.6)

Recalling that the variance and kurtosis are measures for the spread of a
distribution about its mean, Theorem 10.2.1 indicates that ETFs have the most
compact typical sub-frame ESD among all unit-norm frames.

Figure 10.1: Comparison between a repetition (unit-norm tight) frame (UTF);
random (i.i.d.) frame; and ETF, where all frames are with aspect ratio
$\gamma=1/2$. The $x$-axis shows the sub-frame selection probability $p$. The
$y$-axis shows the sub-frame moments 1-6 (top), the empirical spectral
variance (bottom-left), and the empirical spectral kurtosis (bottom-right).

Figure 10.1 shows a comparison of the sub-frame moments, and the corresponding
variance and kurtosis, for several tight frames, as well as for a random
(i.i.d.) frame.

##  Chapter 11 Conclusion and discussion

In this survey we tried to demonstrate the great potential of incorporating
ideas from information theory and random matrix theory into frame theory. We
hope that it will stimulate future research along two main lines of open
questions mentioned above. One line of investigation is practical, namely,
design of robust analog codes. We raised the possibility that ETFs can be used
as robust, and possibly optimally robust, analog codes in a variety of
applications. The second line of investigation points down the ETF-MANOVA
“rabbit hole”. It concerns formal mathematical analysis and explanation of the
universality of the MANOVA distribution, as a limit of the spectrum of random
subsets of frames in a class of well structured frames, and the universality
of the rate of convergence to this limit.

#### Asymptotically robust analog codes

In all the applications we considered in Chapter 9, exact mean-case
performance is given by an expression of the form

$\mathbb{E}_{S}[\Psi(\lambda(G_{S}))]$

where $\Psi$ is a functional of the empirical spectrum, which depends on the
example. Here, $S$ denotes a subset of $[n]$, $F_{S}$ is the corresponding
subset of the frame $F$, $G_{S}=F_{S}^{\dagger}F_{S}$ is the subset Gram
matrix, $\lambda(G_{S})$ is the spectrum of $G_{S}$, and $\mathbb{E}_{S}$
denotes expectation with respect to random selection of $S$. (Note that
$\lambda(G_{S})$ is equivalent to the empirical spectral distribution
$\mu_{G_{S}}$ in (5.1).)

A common property of the three performance measures $\Psi$ mentioned in
Chapter 3, is that a subframe $F_{S}$ which is “more orthogonal”, i.e., whose
spectrum is narrower and more removed from zero, gives better performance
$\Psi(\lambda(G_{S}))$. From this perspective, “good” frames have the property
that typical subframes have a compact and zero-removed spectrum support.

Questions. Generalizing from the examples we have seen, we arrive at
fundamental new questions which span frame theory and random matrix theory:

  1. 1.

How does one calculate quantities of the form
$\mathbb{E}_{S}[\Psi(\lambda(G_{S}))]$ for various functionals $\Psi$?

  2. 2.

How can we design “good frames” - namely frames whose typical subframe has a
spectrum with support as small and as zero-removed as possible?

  3. 3.

Does there exist a “best frame”, namely, a frame with “best” spectrum (with
respect to one of the performance measures, or all of them)? If it exists, can
we characterize it?

Answers. In this survey we review recent results which attempt to answer the
questions above:

  1. 1.

For a large variety of well known frames (Table 6.1), and a large variety of
functionals $\Psi$ (Chapter 3), the quantity
$\mathbb{E}_{S}[\Psi(\lambda(G_{S}))]$ can be calculated to good and known
precision (Chapter 6). Specifically,

$\mathbb{E}_{S}\left[\Psi\left(\lambda(G_{S})\right)\right]\approx\Psi\left(f_{\beta,\gamma}\right)$

where $f_{\beta,\gamma}$ is a known cumulative distribution function that
depends on the aspect ratios of the frame and subframe in question, and
$\Psi(f)$ is defined in (6.4).

  2. 2.

Certain tight frames - of both deterministic and random design - are “good
frames” in the sense that the spectrum of their typical subframe converges to
a universal limiting spectral distribution $f_{\beta,\gamma}$ (Chapter 5),
whose support is compact and zero-removed. In fact, this universal limiting
distribution function is the famous Wachter (Jacobi) astymptotic distribution
familiar from random matrix theory (Chapters 6 and 7).

  3. 3.

There is empirical and theoretic evidence that suggests that equiangular tight
frames in fact achieve the best performance, asymptotically as the frame size
grows (Chapter 10).

#### Mathematical analysis of the ETF-MANOVA relation

More generally, this survey suggests that the spectrum of a typical subframe,
taken from a frame with certain design properties, is a fascinating object -
both from a theoretical perspective and a practical perspective. (Here,
“typical” is taken to mean a subframe selected uniformly at random from the
set of all subframes.) Our results revolved around a fascinating universality
property, whereby the spectrum of (the Gramian of) a typical subframe is
indistinguishable from that of a matrix from the MANOVA (Jacobi) random matrix
ensemble, both in finite-$n$ and asymptotically.

Interestingly, we saw that this universality property holds for a certain kind
of deterministic frames, just as it holds for certain random frames. This
means that, curiously, a famous random matrix theory law arises in a
deterministic and highly symmetric geometrical object. While one can find
hints in the literature that this random-matrix phenomena can be suspected to
arise asymptotically in these highly symmetric deterministic frames, the fact
that universality holds for finite-$n$ as well is striking in our opinion. We
highlight the many open questions and open problems that arise in this
fascinating study of spectra of random subframes.

We particularly note that, to the best of our knowledge, the observed
universality of the convergence rate of the spectrum, and of functionals of
the spectrum, to the universal limiting distribution, is a fascinating object
for further investigation, whose proof could possibly require and stimulate
new mathematical tools.

###  11.1 Summary of asymptotic results for ETF subsets

We conclude this survey by concise summaries of the asymptotic results and the
open questions.

The ETF-MANOVA relation of Chapters 6 and 7 implies that we can compute the
performance measures of large dimensional ETFs and near-ETFs as integrals of
the MANOVA distribution:

$L_{\rm
MSE}(\mbox{ETF},\gamma,p)=\log\left(\int_{\lambda_{-}}^{\lambda_{+}}\frac{1}{x}\cdot
f_{\gamma,p}(x)dx\right),$ (11.1) $L_{\rm
Shannon}(\mbox{ETF},\gamma,p)=\int_{\lambda_{-}}^{\lambda_{+}}\log(x)\cdot
f_{\gamma,p}(x)dx,$ (11.2)

and (by the convergence of the extreme eigenvalues111 While the asymptotic
results in Chapter 7 are only concerned with the bulk spectrum, our empirical
results (Chapter 6, [haikin2017random]) demonstrate a “Tracy-Widom” kind of
result ([ZeitouniBook, feier2012methods]) regarding the high-probability
convergence of the extreme eigenvalues to $\lambda_{\pm}$ of the MANOVA law. )

$L_{\rm
StRIP}(\mbox{ETF},\gamma,p,\delta)=\left\\{\begin{array}[]{cc}1,&\mbox{if
$\delta<\min\\{1-\lambda_{-},\lambda_{+}-1\\}$}\\\
0,&\mbox{otherwise,}\end{array}\right.$ (11.3)

where $L_{\rm MSE}$, $L_{\rm Shannon}$ and $L_{\rm StRIP}$ are the performance
measures (3.2), (3.7b) and (3.11), respectively, and where the asymptotic
MANOVA density $f_{\gamma,p}(x)$ and its support edges
$(\lambda_{-},\lambda_{+})$ are given in (5.5) (with $p=\beta\gamma$). The ETF
superiority conjecture of Chapter 10 implies that the optimal measures (3.4),
(3.9) and (3.11) are given by:

$L_{\rm MSE}^{*}(\gamma,p)=L_{\rm MSE}(\mbox{ETF},\gamma,p),$ (11.4) $L_{\rm
Shannon}^{*}(\gamma,p)=L_{\rm Shannon}(\mbox{ETF},\gamma,p),$ (11.5)

and

$L_{\rm StRIP}^{*}(\gamma,p,\delta)=L_{\rm
StRIP}(\mbox{ETF},\gamma,p,\delta).$ (11.6)

The central conjecture of our work and main open question is whether ETFs are
indeed the best possible frames under random subset selection.

###  11.2 Summary of open questions and conjectures

Finally, let us list some open questions and directions for future research.

The ETF-MANOVA relation

  1. 1.

Asymptotic ETF-MANOVA theorem: find a closed-form expression for the $r$th
moment recursion formula (7.23), in order to prove the conjectured identity
(7.28).

  2. 2.

Method of proof: can we find alternative proofs for the ETF-MANOVA relation
via transform or free-probability methods? See [feier2012methods] and the
result of [raich2016eigenvalue].

  3. 3.

Sub-frame selection: compare the combinatoric “$n$ choose $k$” selection
(assumed in the empirical experiments of Chapter 6) with Bernoulli($p$)
selection (assumed in the asymptotic analysis of Chapter 7). Is (and how much)
the non-asymptotic behavior different?

  4. 4.

Symmetries of ETF: in the course of analyzing the $r$th-moment (7.4), we
identified some invariant properties of small chains of products, which seem
to hold for various ETF constructions, and may be useful to characterize the
non-asymptotic behavior of the moments [MarinaThesis].

  5. 5.

Near-ETF: how far can we relax the equi-angular and tightness conditions,
without changing the asymptotic MANOVA behavior?

  6. 6.

Non-asymptotic regime: is the expected $r$th moment of a finite frame
identical for all ETFs, or dependent on the specific ETF construction?
(Section 7.2).

  7. 7.

Universality of convergence rate: can we analyze the non-asymptotic behavior
of ETF (and near-ETF) subsets? can we prove the universal convergence rate?
(Section 6.4)

  8. 8.

Emergence of MANOVA ensemble:  what makes a typical ETF sub-frame similar to
the ratio of i.i.d. matrices (5.4) in the MANOVA definition?

  9. 9.

Beyond the limiting spectrum: find other RMT-like properties of a typical ETF
sub-frame; e.g., isotropic behavior of eigen-vectors; convergence of the
extreme eigenvalues (see footnote 1 in this chapter); diagonal of the Gram
matrix, etc.

General frame codes: asymptotic theory and bounds

  1. 10.

Low-pass frames: do randomly-selected sub-frames of LPF have a limiting ESD?

  2. 11.

Spherical codes: find similarities and distinction between maximal distance
spherical codes and ETFs.

  3. 12.

Sub-frame rectangularity: do the performance inequalities in Chapter 8
continue to hold if we fix $p$ and let both $\beta$ and $\gamma$ vary
accordingly?

  4. 13.

RMT-like behavior of frames: are there other families of deterministic frames
(for example, non-ETF harmonic frames) that exhibit a (possibly different)
limiting ESD for the subframe eigenvalues? If so, what is the set of
attainable distributions for fixed $\gamma$ and $\beta$? Is this set convex
(closed under mixtures of frame families)?

  5. 14.

Erasure Welch bound: how does the EWB (10.1) connect to Welch’s generalized
bound for higher-order moments of the cross-correlations [welch1974lower]; can
we derive one from the other?

  6. 15.

ETF superiority: find a sharper definition (beyond variance, kurtosis and
functional (“$\Psi$”) measures).

  7. 16.

Analog versus digital: are analog codes inherently inferior to digital
(Shannon) codes?

Applications for analog coding

  1. 17.

Coded distributed computing (CDC): modify the noise amplification performance
measure (3.2) to the case where one needs to recover only the input sum (or
one projection or subspace), rather than full inversion [RoyeeYosibashPaper].

  2. 18.

Source with erasures: can we obtain the full MMSE gain, and eliminate the
“$+1$” penalty term in (2.11) for a general $\beta$?

  3. 19.

Binary analog codes: can we approximate the superior ETF performance using
binary frame codes (similar to [babadi2011spectral])?

It is our hope that these questions, and others implicitly raised in this
survey, will stimulate research in this fascinating area at the intersection
of information theory, random matrices, geometry and combinatorics.

###### Acknowledgements.

Matan Gavish was partially supported by the Israel Science Foundation, grant
1523/16. Dustin G. Mixon was partially supported by AFOSR FA9550-18-1-0107 and
NSF DMS 1829955. Ram Zamir was partially supported by the Israel Science
Foundation, grants 676/15 and 2623/20. We thank Ofer Zeitouni for suggesting
the moment approach to prove the ETF-MANOVA theorem.

##  Appendix A Entropy-coded dithered quantization

To asses the rate-distortion performance (2.11) of the analog coding scheme,
we shall adopt the additive-noise “test channel” of entropy-coded
(subtractive) dithered quantization (ECDQ) [zamir1992universal]. In this
model, $Q^{({\rm dither})}(\tilde{{\bf x}})=\tilde{{\bf x}}+{\bf Z}$, where
${\bf Z}$ is an independent uniform or Gaussian noise, and the quantizer
entropy is given by the mutual information $H(Q^{({\rm dither})}(\tilde{{\bf
X}}))=I(\tilde{{\bf X}};\tilde{{\bf X}}+{\bf Z})$.111To be more precise, the
quantizer operation is given by $Q(\tilde{{\bf x}}+{\rm dither})-{\rm
dither}$, where ${\rm dither}$ is uniform over the fundamental quantizer cell.
For good high-dimensional lattice quantizers, this uniform dither distribution
tends to be white Gaussian. See [zamir2014lattice].  The equivalent error in
the important samples, ${\bf E}_{s}=\hat{{\bf X}}_{s}-{\bf X}_{s}$ is thus, by
(2.8), given by ${\bf E}_{s}=F_{S}^{\dagger}\cdot(\tilde{{\bf X}}+{\bf
Z})-{\bf X}_{s}=F_{S}^{\dagger}\cdot{\bf Z}$, implying that the mean-squared
distortion per important sample (2.2) is

$D=\frac{1}{k}\mathbb{E}\left\\{\|{\bf
E}_{s}\|^{2}\right\\}=\frac{\sigma_{z}^{2}}{k}\cdot\mbox{trace}\\{F_{S}^{\dagger}F_{S}\\}=\sigma_{z}^{2},$
(A.1)

where $\sigma_{z}^{2}$ is the quantizer mean-squared error, and the last
equation follows since $F_{S}$ has $k$ unit-norm columns.222Further
improvement can be obtained using minimum mean-squared error “Wiener”
estimation, [haikin2016analog], but this is negligible when $\sigma_{x}^{2}\gg
D$. Since the decoder is blind to the location of important samples (that
affect the correlation of the vector $\tilde{{\bf X}}$), the coding rate is
equal to that of a white Gaussian vector $\tilde{{\bf W}}$ with the same power
as $\tilde{{\bf X}}$, [lapidoth1997role]. Assuming the LS estimator (2.9), we
have $E\|\tilde{{\bf
X}}\|^{2}=\sigma_{x}^{2}\cdot\mbox{trace}\\{(F_{S}^{\dagger}F_{S})^{-1}\\}$,
so the ECDQ mutual information formula (with a Gaussian dither) becomes:

$\displaystyle R_{\rm analog}$ $\displaystyle=\frac{1}{n}I(\tilde{{\bf
W}};\tilde{{\bf W}}+{\bf Z})$ (A.2)
$\displaystyle=\frac{m}{2n}\cdot\log\left(1+\frac{\frac{1}{m}\mathbb{E}\|\tilde{{\bf
X}}\|^{2}}{\sigma_{z}^{2}}\right)$ (A.3)
$\displaystyle=\beta\cdot\frac{p}{2}\cdot\log\left(1+\frac{\sigma_{x}^{2}}{D}\cdot\frac{\mbox{trace}\\{(F_{S}^{\dagger}F_{S})^{-1}\\}}{m}\right)$
(A.4)

where $\beta\stackrel{{\scriptstyle\Delta}}{{=}}k/m$, and we used (2.2) with
$k=pn$.

##  Appendix B Information theoretic proofs for sub-frame inequalities

In proving Theorems 8.0.1 and 8.0.2 we use the following two lemmas.

######  Lemma B.0.1.

For any $m\times n$ frame $F$, and $k<n$,

$\frac{1}{{n\choose k}}\sum_{S:|S|=k}F_{S}\cdot
F_{S}^{\dagger}=\frac{k}{n}\cdot F\cdot F^{\dagger},$ (B.1)

where the average in the left hand side is over all $k$-subsets of
$\\{1,\ldots,n\\}$ as in (3.8)-(3.2).

Remark: For Bernoulli($p$) selection, i.e., when each index in
$\\{1,\ldots,n\\}$ belongs to $S$ with probability $p$ independently of the
other indices, the lemma becomes $\mathbb{E}_{S}\\{F_{S}\cdot
F_{S}^{\dagger}\\}=p\cdot F\cdot F^{\dagger}$, where
$\mathbb{E}_{S}\\{\cdot\\}$ denotes expectation with respect to the selection
of $S$.

Proof: The case $k=1$ of (B.0.1) is the standard expansion $\sum_{i=1}^{n}{\bf
f}_{i}\cdot{\bf f}_{i}^{\dagger}=F\cdot F^{\dagger}$ (multiplied by $1/n$),
where the ${\bf f}_{i}$’s are the columns of $F$. For a general $k$, let $P$
denote a random diagonal matrix, whose $\\{0,1\\}$ diagonal elements
correspond to the selection of the subset $S$. That is, $F_{S}\cdot
F_{S}^{\dagger}=F\cdot P\cdot F^{\dagger}$, where the diagonal of $P$ is
uniform over all $n$ choose $k$ binary vectors with $k$ ones and $n-k$ zeroes.
Thus $\mathbb{E}_{S}\\{F_{S}\cdot F_{S}^{\dagger}\\}=\mathbb{E}_{S}\\{F\cdot
P\cdot F^{\dagger}\\}=F\cdot\mathbb{E}_{S}\\{P\\}\cdot
F^{\dagger}=F\cdot(k/n\cdot I_{n})\cdot F^{\dagger}$, which is the right-hand
side of (B.1). $\Box$

######  Lemma B.0.2.

For $k_{1}<k_{2}<n$, if $S_{2}$ is uniformly drawn from all $k_{2}$ subsets of
$\\{1,\ldots,n\\}$, and $S_{1}$ is uniformly drawn from all $k_{1}$ subsets of
$S_{2}$, then $S_{1}$ is uniform on all $k_{1}$ subsets of $[n]$.

Proof: Since the distribution of $S_{1}$ is invariant under the action of the
symmetric group of $[n]$, it is necessarily uniform. $\Box$

Proof of Theorem 8.0.1: We first prove the inequality with respect to the
“edge”: $L(F,k)\leq L(F,n)$, and then extend to any $k_{1}<k_{2}$. The Ky Fan
inequality [CoverBook], says that if $K_{1}$ and $K_{2}$ are $m\times m$ PSD
matrices, then

$\det[a\cdot K_{1}+(1-a)\cdot K_{2}]\geq\det[K_{1}]^{a}\cdot\det[K_{2}]^{1-a}$
(B.2)

for any $0<a<1$, which by taking logarithm implies that $\log\det[K]$ is
concave [CoverBook]. Thus, starting from (3.8),

$\displaystyle L_{\rm Shannon}(F,k)$
$\displaystyle=\mathbb{E}_{S}\\{\log[\sqrt[m]{\det(F_{S}\cdot
F_{S}^{\dagger})}/(k/m)]\\}$ (B.3)
$\displaystyle\leq\log[\sqrt[m]{\det(\mathbb{E}_{S}\\{F_{S}\cdot
F_{S}^{\dagger}\\})}/(k/m)]$ (B.4) $\displaystyle=\log[\sqrt[m]{\det(k/n\cdot
F\cdot F^{\dagger})}/(k/m)]$ (B.5) $\displaystyle=\log[\sqrt[m]{\det(F\cdot
F^{\dagger})}/(n/m)]$ (B.6) $\displaystyle=L_{\rm Shannon}(F,n)$ (B.7)

where the second line is by Jensen’s inequality and the log concavity of the
determinant (B.2), and the third line is by the identity in Lemma B.0.1.
Turning to the general case, for $k_{1}<k_{2}<n$, let $S_{1}\subset S_{2}$
denote a $k_{1}$-subset of a $k_{2}$-subset $S_{2}$ of $\\{1,\ldots,n\\}$. By
the definition of the average subset Shannon transform of a frame (3.8), we
have

$\displaystyle L_{\rm Shannon}(F,k_{1})$
$\displaystyle=\mathbb{E}_{S_{1}}\\{\Psi_{\rm Shannon}(F_{S_{1}})\\}$ (B.8)
$\displaystyle=\mathbb{E}_{S_{2}}\\{E_{S_{1}|S_{2}}\\{\Psi_{\rm
Shannon}(F_{S_{1}})\\}$ (B.9) $\displaystyle=\mathbb{E}_{S_{2}}\\{L_{\rm
Shannon}(F_{S_{2}},k_{1})\\}$ (B.10)
$\displaystyle\leq\mathbb{E}_{S_{2}}\\{L_{\rm Shannon}(F_{S_{2}},k_{2})\\}$
(B.11) $\displaystyle=\mathbb{E}_{S_{2}}\\{\Psi_{\rm Shannon}(F_{S_{2}})\\}$
(B.12) $\displaystyle=L_{\rm Shannon}(F,k_{2})$ (B.13)

where $\mathbb{E}_{S_{1}|S_{2}}\\{\cdot\\}$ denotes expectation over a uniform
distribution on all $k_{1}$-subsets of a given $S_{2}$, and
$\mathbb{E}_{S_{2}}\\{\cdot\\}$ denotes expectation over a uniform
distribution on all $k_{2}$-subsets of $\\{1,\ldots,n\\}$. The second line
follows from Lemma B.0.2 by iterated expectation; the third line follows by
viewing $F_{S_{2}}$ as the full frame in the inner expectation; the inequality
follows from the first part of the proof (B.3)-(B.7) setting $k=k_{1}$ and
$n=k_{2}$; and the last two lines are again by the definition (3.8). $\Box$

Proof of Theorem 8.0.2: For an $n\times n$ covariance matrix $K$ (i.e., $K$ is
a non-negative matrix), and a $k$-subset $S\subset\\{1,\ldots,n\\}$, let
$K_{S}$ denote the corresponding $k\times k$ sub-matrix of $K$. Define the
subset trace-inverse of the covariance $K$ as the average of
$\text{trace}\left(K_{S}^{-1}\right)/k$ over all $k$-subsets:

$M_{k}^{(n)}=\frac{1}{{n\choose
k}}\sum_{S:|S|=k}\frac{1}{k}\text{trace}\left(K_{S}^{-1}\right),$ (B.14)

for $k=1,\ldots,n$. For example, the two extremes are
$M_{1}^{(n)}=1/n\sum_{i=1}^{n}(K_{ii})^{-1}$, and
$M_{n}^{(n)}=1/n\sum_{i=1}^{n}(K^{-1})_{ii}$. It is shown in
[khina2020monotonicity] that the sequence of subset trace inverses is
monotonically non decreasing, i.e.,

$M_{1}^{(n)}\leq M_{2}^{(n)}\leq\ldots\leq M_{n}^{(n)}$ (B.15)

with equality if and only if the covariance matrix $K$ is proportional to
identity.111 If the matrix $K$ is Toeplitz (corresponding to stationary
vector), then the inequality between the two edges, $M_{1}^{(n)}\leq
M_{n}^{(n)}$, follows from the arithmetic-to-harmonic means inequality applied
to the eigenvalues $\lambda_{1},\ldots,\lambda_{n}$ of $K$, and equality holds
if and only if the eigenvalues are constant, i.e., $K=\lambda I$. (This is
because by the Toeplitz property $K_{ii}=\sigma^{2}$ for all $i$, and
therefore $M_{1}^{(n)}=1/\sigma^{2}$; furthermore,
$\sigma^{2}=\text{trace}\left(K\right)/n=(\lambda_{1}+\ldots+\lambda_{n})/n$,
while the eigenvalues of $K^{-1}$ are the reciprocals of the eigenvalues of
$K$, so $M_{n}^{(n)}=(1/\lambda_{1}+\ldots+1/\lambda_{n})/n$.)  This
inequality can be thought of as the “MSE counterpart” of the subset
determinant inequality in [CoverBook] and [dembo1991information]. Using
similar tools as in [CoverBook], we can show the same inequality also for
subset log trace-inverse. Namely, letting

$L_{k}^{(n)}=\frac{1}{{n\choose
k}}\sum_{S:|S|=k}\log\left(\frac{1}{k}\text{trace}\left(K_{S}^{-1}\right)\right),$
(B.16)

we have

$L_{1}^{(n)}\leq L_{2}^{(n)}\leq\ldots\leq L_{n}^{(n)}.$ (B.17)

Proof: Starting with the right edge, noting that ${n\choose n}=1$ and
${n\choose n-1}=n$, we have

$\displaystyle L_{n}^{(n)}$
$\displaystyle=\log\left(\frac{1}{n}\text{trace}\left(K^{-1}\right)\right)$
$\displaystyle=\log\left(M_{n}^{(n)}\right)$
$\displaystyle\geq\log\left(M_{n-1}^{(n)}\right)$
$\displaystyle=\log\left(\frac{1}{n}\sum_{S:|S|=n-1}\frac{1}{n-1}\text{trace}\left(K_{S}^{-1}\right)\right)$
$\displaystyle\geq\frac{1}{n}\sum_{S:|S|=n-1}\log\left(\frac{1}{n-1}\text{trace}\left(K_{S}^{-1}\right)\right)$
$\displaystyle=L_{n-1}^{(n)},$ (B.18)

where the equality lines are by definition; the first inequality follows from
(B.15); and the second inequality follows by Jensen. We then continue
similarly to the second part of the proof of Theorem 8.0.1 to prove that
$L_{k}^{(n)}\leq L_{k+1}^{(n)}$ for any $k$. Specifically, by Lemma B.0.2,
when computing $L_{k}^{(n)}$ we first condition on a specific $(k+1)$-subset
$S_{2}$ and average over all its $k$-subsets, and then average over all
$(k+1)$-subsets $S_{2}$ of $\\{1,\ldots,n\\}$. By (B), the inner average is
upper bounded by the log trace-inverse of $K_{S_{2}}$, which becomes
$L_{k+1}^{(n)}$ after taking the outer average. $\Box$

Theorem 8.0.2 now follows from (B.17) by viewing the $F^{\dagger}\cdot F$ as
the covariance matrix $K$, so $L_{\rm MSE}(F,k)=L_{k}^{(n)}$. (Note that for
$k>m$ the matrix $F_{S}^{\dagger}F_{S}$ is singular, so $L_{k}^{(n)}$ is
infinite.)

##  Appendix C Variance of $1st$ and $2nd$ moments (proof of Theorem 7.1.3)

Using the moment-variance formula (7.15),

$\displaystyle V_{r}$ $\displaystyle\triangleq{\rm
Var}\left[\frac{1}{n}\mbox{trace}\left((FPF^{\dagger})^{r}\right)\right]$
$\displaystyle=\mathbb{E}\left[\left(\frac{1}{n}\mbox{trace}\left((FPF^{\dagger})^{r}\right)\right)^{2}\right]-m_{r}^{2}$

we obtain for $r=1$,

$\displaystyle V_{1}^{\rm ETF}$
$\displaystyle=\mathbb{E}\left[\left(\frac{1}{n}\mbox{trace}\left(FPF^{\dagger}\right)\right)^{2}\right]-\mathbb{E}\left[\frac{1}{n}\mbox{trace}\left(FPF^{\dagger}\right)\right]^{2}$
$\displaystyle=\mathbb{E}\left[\left(\frac{1}{n}\mbox{trace}\left(FPF^{\dagger}\right)\right)^{2}\right]-m_{1}^{2}$
$\displaystyle=\frac{1}{n^{2}}\mathbb{E}\left[\left(\sum_{i=1}^{n}\langle{\bf
f}_{i},{\bf f}_{i}\rangle P_{i}\right)\left(\sum_{j=1}^{n}\langle{\bf
f}_{j},{\bf f}_{j}\rangle P_{j}\right)\right]-p^{2}$
$\displaystyle=\frac{1}{n^{2}}\mathbb{E}\left[\left(\sum_{i=1}^{n}P_{i}\right)\left(\sum_{j=1}^{n}P_{j}\right)\right]-p^{2}$
$\displaystyle=\frac{1}{n^{2}}\left[np+n(n-1)p^{2}\right]-p^{2}$
$\displaystyle=\frac{1}{n}\left[p-p^{2}\right].$ (C.1)

For $r=2$, we obtain

$\displaystyle V_{2}^{\rm ETF}$
$\displaystyle=\mathbb{E}\left[\left(\frac{1}{n}\mbox{trace}\left((FPF^{\dagger})^{2}\right)\right)^{2}\right]-\mathbb{E}\left[\frac{1}{n}\mbox{trace}\left((FPF^{\dagger})^{2}\right)\right]^{2}$
$\displaystyle=\mathbb{E}\left[\left(\frac{1}{n}\mbox{trace}\left((FPF^{\dagger})^{2}\right)\right)^{2}\right]-m_{2}^{2}$
$\displaystyle=\frac{1}{n^{2}}\mathbb{E}\left[\left(\sum_{i,j=1}^{n}|\langle{\bf
f}_{i},{\bf
f}_{j}\rangle|^{2}P_{i}P_{j}\right)\left(\sum_{k,m=1}^{n}|\langle{\bf
f}_{k},{\bf f}_{m}\rangle|^{2}P_{k}P_{m}\right)\right]-(p+p^{2}x)^{2}.$ (C.2)

If $F$ is an ETF, then according to the Welch bound $|\langle{\bf f}_{i},{\bf
f}_{j}\rangle|^{2}=\frac{x}{n-1}$ for $i\neq j$. Similarly to the computation
of the moments, we split to cases by number of distinct values of
$\\{i,j,k,m\\}$.

$\displaystyle\mathbb{E}\left[\left(\sum_{i,j=1}^{n}|\langle{\bf f}_{i},{\bf
f}_{j}\rangle|^{2}P_{i}P_{j}\right)\left(\sum_{k,m=1}^{n}|\langle{\bf
f}_{k},{\bf f}_{m}\rangle|^{2}P_{k}P_{m}\right)\right]=np$
$\displaystyle+\left[n(n-1)+4n(n-1)\frac{x}{n-1}+2n(n-1)\frac{x^{2}}{(n-1)^{2}}\right]p^{2}$
$\displaystyle+\left[2n(n-1)(n-2)\frac{x}{n-1}+4n(n-1)(n-2)\frac{x^{2}}{(n-1)^{2}}\right]p^{3}$
$\displaystyle+n(n-1)(n-2)(n-3)\frac{x^{2}}{(n-1)^{2}}p^{4}.$ (C.3)

It follows that

$\displaystyle V_{2}=\frac{1}{n}\Bigl{[}p$
$\displaystyle+\left(-1+4x+\frac{2x^{2}}{n-1}\right)p^{2}+\left(-4x+\frac{4(n-2)x^{2}}{n-1}\right)p^{3}$
$\displaystyle+\left(\frac{(6-4n)x^{2}}{n-1}\right)p^{4}\Bigr{]}.$ (C.4)


# An Augmented Lagrangian Deep Learning Method for Variational Problems with
Essential Boundary Conditions

Jianguo Huang<EMAIL_ADDRESS>
School of Mathematical Sciences, and MOE-LSC, Shanghai Jiao Tong University,
Shanghai, China
Haoqin Wang<EMAIL_ADDRESS>
School of Mathematical Sciences, and MOE-LSC, Shanghai Jiao Tong University,
Shanghai, China
Tao Zhou 111Corresponding author<EMAIL_ADDRESS>
LSEC, Institute of Computational Mathematics and Scientific/Engineering
Computing,
Academy of Mathematics and Systems Science, Chinese Academy of Sciences,
Beijing, China

###### Abstract

This paper is concerned with a novel deep learning method for variational
problems with essential boundary conditions. To this end, we first reformulate
the original problem into a minimax problem corresponding to a feasible
augmented Lagrangian, which can be solved by the augmented Lagrangian method
in an infinite dimensional setting. Based on this, by expressing the primal
and dual variables with two individual deep neural network functions, we
present an augmented Lagrangian deep learning method for which the parameters
are trained by the stochastic optimization method together with a projection
technique. Compared to the traditional penalty method, the new method admits
two main advantages: i) the choice of the penalty parameter is flexible and
robust, and ii) the numerical solution is more accurate in the same magnitude
of computational cost. As typical applications, we apply the new approach to
solve elliptic problems and (nonlinear) eigenvalue problems with essential
boundary conditions, and numerical experiments are presented to show the
effectiveness of the new method.

Keywords. The augmented Lagrangian method; Deep learning; Variational
problems; Saddle point problems; Essential Boundary Conditions.

AMS subject classifications: 65N25, 65N30, 68U99.

##  1 Introduction

Variational problems play important roles in various industrial and
engineering applications, with typical examples including partial differential
equations (PDEs) and eigenvalue problems. Many classical numerical methods
have been developed for such problems, e.g., the finite difference method, the
spectral method, and the finite element method. The first two methods are
generally used for solving problems over regular domains while the latter one
is particularly suitable for problems in irregular domains [11, 6]. In recent
years, deep learning based techniques have been widely used to solve a variety
of variational problems [32, 40, 24, 16, 35, 39, 23, 9, 28, 14, 18].
Historically, related studies can date back to the 1990s [26, 12]. We also
refer the reader to [13] and the references therein for a comprehensive review
on machine learning from the perspective of computational mathematics. For
such kind of methods, deep neural networks (DNNs) are exploited to
parameterize the PDE solutions and appropriate parameters are identified by
minimizing an optimization problem formulated from the PDEs. The most
significant feature of those methods is that they are mesh-free, and their
approximation capacity has been well studied in recent years [21, 3, 15, 19,
33].

For variational problems with natural boundary conditions, one doesn’t need to
impose these conditions on the admissible functions [16], so that the DNNs can
easily be used for approximation. However, for variational problems with
essential boundary conditions, these conditions should be imposed on the
admissible functions, and this gives rise to a significant difficulty since
one cannot enforce the boundary condition in a simple way even at the
interpolation nodes for a neural network function. It is worth noting that
even in the context of finite element methods, this is also a very tough
issue. In fact, one has to use Nitsche’s trick [30], developed further by
Stenberg [36], to handle this issue. As far as we know, there are two main
strategies to overcome the bottleneck in deep learning framework:

  * •

The first strategy is to construct neural network functions that satisfy the
essential boundary conditions exactly. For instance, if the boundary condition
is given by $u=g$ on the boundary $\Gamma$, then we construct the approximate
function by

$\phi(\bm{x};\bm{\theta})=\ell(\bm{x})\psi(\bm{x};\bm{\theta})+\bar{g}(\bm{x}),$
(1.1)

where $\ell(\bm{x})$ is a known function such that on $\Gamma$ it holds
$\ell(\bm{x})=0$, $\bar{g}$ is the extension of $g$ to the whole domain, and
$\psi(\bm{x};\bm{\theta})$ is another neural network function that is used to
approximate the solution in the domain. The main limitation of this approach
is that for problems with complex (non-regular) domains, it is in general not
easy to find explicit functions $\ell$ and $\bar{g}.$ For details, one can
refer to [4] and references therein. It is worth mentioning that based on the
formulation (1.1), one can also introduce an additional neural network
function on the boundary $\Gamma$ to approximate $g$ by a least squares
approach [34].

  * •

Another strategy is the penalty method, where a penalty term (with a penalty
parameter $\beta$) is included into the objective functional to enforce the
boundary condition [32, 16, 35, 39]. This method is easy to implement.
Theoretically, the penalty parameter $\beta$ should be chosen large enough,
however, this may make the optimization problem become ill-conditioned [31].
We also mention the deep Nitsche method proposed in [27], where Nitsche’s
variational formula is used for the second order elliptic problems to avoid a
large penalty parameter.

In this work, we intend to present an augmented Lagrangian deep learning
(ALDL) method to handle variational problems with essential boundary
conditions. For this purpose, we shall first rewrite the original problem as a
minimax problem associated with a feasible augmented Lagrangian, which can be
solved by the augmented Lagrangian method in an infinite dimensional setting
[17, 31]. We then express the primal and dual variables with two individual
DNN functions, respectively, and train the associated network parameters with
the stochastic optimization method (based on the augmented Lagrangian method).
It is worth noting that we require to solve a least square problem in order to
update the dual variable, and this step can be viewed as a nonlinear
projection in the corresponding parameter space. As typical applications, we
apply the ALDL method to solve elliptic problems and eigenvalue problems with
essential boundary conditions. Numerical results indicate that the ALDL method
admits two main advantages compared to the penalty method:

  * •

The choice of the penalty parameter is flexible and robust.

  * •

The numerical solution is more accurate in the same magnitude of computational
cost.

The rest of this paper is organized as follows. In Section 2, we introduce the
variational problem and its minimax formulation. As typical cases, an elliptic
PDE and eigenvalue problems are presented. In Section 3, we recall the
augmented Lagrangian method in an infinite dimensional setting, and then
propose the augmented Lagrangian deep learning method. In Section 4, numerical
examples are reported to show the performance of the proposed method. Finally,
we provide with some concluding remarks in Section 5.

##  2 The variational problem and its primal-dual formulation

To begin, we first introduce some notations for later uses. For a real Hilbert
space $V$ equipped with a norm $\|\cdot\|_{V}$, we denote by
$\langle\cdot,\cdot\rangle_{V}$ the induced inner product over $V$. We use the
standard symbols and notations for Sobolev space and their norms or semi-norms
and refer the reader to the reference [1] for details. Let
$\Omega\subset\mathbb{R}^{d}$ be a bounded domain with Lipschitz boundary,
where $d$ is a natural number. We also denote by $\Gamma$ its boundary and
$\bar{\Omega}$ the closure of $\Omega$, respectively. We let $B$ be a bounded
linear operator from $V$ to another Hilbert space $W$ and a typical example in
our variational problems is $W=L^{2}(\Gamma)$.

Throughout this paper, we consider the following variational problem

$\min_{v\in V_{g}}J(v),$ (2.1)

where $J(v)$ is a nonlinear functional over $V_{g}$ and

$V_{g}=\\{v\in V:{B}v=g\ \mbox{on}\ \Gamma\\}.$

###  2.1 Primal-dual formulation

To deal with the constraint in the admissible set $V_{g}$, the augmented
Lagrangian method [17] suggests to consider an augmented Lagrangian function
as following

$\mathcal{L}_{\beta}(v,\mu)=J(v)-\langle\mu,{B}v-g\rangle_{W}+\frac{\beta}{2}\|{B}v-g\|_{W}^{2},$
(2.2)

where $\mu\in W$ is a Lagrange multiplier function (the dual variable) and
$\beta$ is a positive constant. Then we obtain the following minimax problem:

$\min_{v\in V}\max_{\mu\in W}\mathcal{L}_{\beta}(v,\mu).$ (2.3)

By a direct manipulation, we have

$\max_{\mu\in W}\mathcal{L}_{\beta}(v,\mu)=\begin{cases}J(v),\ &v\in V_{g},\\\
+\infty,\ &v\notin V_{g}.\end{cases}$

Consequently, the variational problem (2.1) is equivalent to the minimax
problem (2.3). In other words, if $(u,\lambda)\in V\times W$ is a solution of
the problem (2.3), then $u$ is a solution of the problem (2.1). On the
contrary, if $u$ is a solution of the problem (2.1), there exists a function
$\lambda\in W$ such that $(u,\lambda)$ is a solution of the problem (2.3).
Furthermore, we assume that $(u,\lambda)$ is a saddle point of the Lagrangian
$\mathcal{L}_{\beta}(\cdot,\cdot)$, i.e.,

$\min_{v\in V}\max_{\mu\in
W}\mathcal{L}_{\beta}(v,\mu)=\mathcal{L}_{\beta}(u,\lambda)=\max_{\mu\in
W}\min_{v\in V}\mathcal{L}_{\beta}(v,\mu).$ (2.4)

Notice that it is rather difficult to show the existence of a saddle point of
a general functional and one important technique is the Ky Fan-Sion theorem
[8]. If we write

$\mathcal{F}_{\beta}(\mu)=\min_{v\in V}\mathcal{L}_{\beta}(v,\mu),$

then one can turn to solve the following dual problem

$\max_{\mu\in W}\mathcal{F}_{\beta}(\mu).$ (2.5)

###  2.2 Some applications

We now present three typical applications of the above primal dual
formulation.

####  2.2.1 Second order elliptic PDEs

The first example is the second order elliptic equation:

$\begin{cases}-\operatorname*{div}\big{(}\bm{A}(\bm{x})\nabla
u(\bm{x})\big{)}+c(\bm{x})u(\bm{x})=f(\bm{x})&\mbox{\rm in}\ \Omega,\\\
u(\bm{x})=g(\bm{x}),&\mbox{\rm on}\ \Gamma,\end{cases}$ (2.6)

where $\bm{A}(\bm{x})\in C^{1}(\bar{\Omega})$ is uniformly elliptic and
$c(\bm{x})\in L^{2}({\Omega})$ is nonnegative over $\Omega$.

The variational formula of (2.6) is

$\min_{v\in V_{g}}J(v),$ (2.7)

where

$J(v)=\frac{1}{2}\int_{\Omega}[\bm{A}(\bm{x})\nabla v(\bm{x})\cdot\nabla
v(\bm{x})+c(\bm{x})v^{2}(\bm{x})-2f(\bm{x})v(\bm{x})]\,{\rm d}x,$

and

$V_{g}=\\{v\in H^{1}(\Omega):v=g\ \mbox{a.e. on}\ \Gamma\\}.$

Based on the abstract setting given in the last subsection, we can rewrite the
variational formula (2.7) as the following problem

$\min_{v\in V}\max_{\mu\in
W}\mathcal{L}_{\beta}(v,\mu),\quad\mathcal{L}_{\beta}(v,\mu)=J(v)-\int_{\Gamma}\big{[}\mu(\bm{x})\big{(}v(\bm{x})-g(\bm{x})\big{)}\big{]}\,{\rm
d}x+\frac{\beta}{2}\int_{\Gamma}\big{[}v(\bm{x})-g(\bm{x})\big{]}^{2}\,{\rm
d}x,$ (2.8)

where $V=H^{1}(\Omega),W=L^{2}(\Gamma)$, and $\beta$ is a positive parameter.
As shown in [10], there exists a unique saddle point for the above problem.

####  2.2.2 Linear eigenvalue problems

The second example is the eigenvalue problem. Suppose we want to find the
smallest eigenvalue and its eigenfunction for a positive self-adjoint
differential operator, e.g.,

$\begin{cases}-\nabla\cdot\big{(}p(\bm{x})\nabla
u(\bm{x})\big{)}+q(\bm{x})u(\bm{x})=\rho u(\bm{x})&\mbox{\rm in}\ \Omega,\\\
u(\bm{x})=0\quad&\mbox{\rm on}\ \Gamma,\end{cases}$ (2.9)

where $p(\bm{x})\in C^{1}(\bar{\Omega})$ is uniformly elliptic and
$q(\bm{x})\in C(\bar{\Omega})$ is nonnegative over $\Omega$.

The variational formula of (2.9) (cf. [2]) is

$\min_{v\in V_{0}}J(v),\quad J(v)=\frac{\int_{\Omega}[p(\bm{x})\nabla
v(\bm{x})\cdot\nabla v(\bm{x})+q(\bm{x})v^{2}(\bm{x})]\,{\rm
d}x}{\int_{\Omega}v^{2}(\bm{x})\,{\rm d}x},$ (2.10)

where $V_{0}=H_{0}^{1}(\Omega)$.

Suppose $(\rho,u)$ is the solution of (2.9). Then it is easy to check
$(\rho,cu)$ is also a solution of (2.9), where $c$ is a non-zero real number.
This motivates us to find the normalized eigenfunction, i.e.,

$\min_{v\in
V_{0}}\tilde{J}(v),\quad\tilde{J}(v)=\int_{\Omega}[p(\bm{x})\nabla\tilde{v}(\bm{x})\cdot\nabla\tilde{v}(\bm{x})+q(\bm{x})\tilde{v}^{2}(\bm{x})]\,{\rm
d}x,$ (2.11)

where $\tilde{v}=v/\|v\|_{0,\Omega}$. As mentioned above, we can reformulate
(2.11) as the following minimax problem

$\min_{v\in V}\max_{\mu\in
W}\mathcal{L}_{\beta}(v,\mu),\quad\mathcal{L}_{\beta}(v,\mu)=\tilde{J}(v)-\int_{\Gamma}\mu(\bm{x})\tilde{v}(\bm{x})\,{\rm
d}x+\frac{\beta}{2}\int_{\Gamma}\tilde{v}^{2}(\bm{x})\,{\rm d}x,$ (2.12)

where $V=H^{1}(\Omega)$, $W=L^{2}(\Gamma)$ and $\beta$ is a positive
parameter.

####  2.2.3 Nonlinear eigenvalue problems

The third example is a nonlinear eigenvalue problem:

$\begin{cases}-\nabla\cdot\big{(}\bm{A}(\bm{x})\nabla
u(\bm{x})\big{)}+V(\bm{x})u(\bm{x})+u^{3}(\bm{x})=\rho
u(\bm{x})&\quad\mbox{in}\ \Omega,\\\ u(\bm{x})=0&\quad\mbox{on}\
\partial\Omega,\\\ \|u\|_{0,\Omega}=1,\end{cases}$ (2.13)

where $\bm{A}(\bm{x})\in(L^{\infty}(\Omega))^{d\times d}$ is symmetric and
uniformly elliptic, and $V(\bm{x})\in L^{2}(\Omega)$.

The variational formula of (2.13) (cf. [7]) is

$\min_{\begin{subarray}{c}v\in V_{0}\\\
\|v\|_{0,\Omega}=1\end{subarray}}J(v),\quad
J(v)=\frac{1}{2}{\int_{\Omega}\big{[}\bm{A}(\bm{x})\nabla v(\bm{x})\cdot\nabla
v(x)+V(\bm{x})v^{2}(\bm{x})+v^{4}(\bm{x})\big{]}}\,{\rm d}x,$ (2.14)

where $V_{0}=H_{0}^{1}(\Omega)$. According to reference [7], the ground state
non-negative solution $(\rho,u)$ of (2.14) is unique for $1\leq d\leq 3$.
Similar to linear eigenvalue problems, we use the normalization technique to
relax the constraint $\|v\|_{0,\Omega}=1$, and then reformulate the
variational problem (2.14) as

$\min_{v\in
V_{0}}\tilde{J}(v),\quad\tilde{J}(v)=\frac{1}{2}{\int_{\Omega}\big{[}\bm{A}(\bm{x})\nabla\tilde{v}(\bm{x})\cdot\nabla\tilde{v}(x)+V(\bm{x})\tilde{v}^{2}(\bm{x})+\tilde{v}^{4}(\bm{x})\big{]}}\,{\rm
d}x,$ (2.15)

where $\tilde{v}=v/\|v\|_{0,\Omega}$. As discussed in Section 2.1, we can
rewrite the variational formula (2.15) as the minimax problem

$\min_{v\in V}\max_{\mu\in
W}\mathcal{L}_{\beta}(v,\mu),\quad\mathcal{L}_{\beta}(v,\mu)=\tilde{J}(v)-\int_{\Gamma}\mu(\bm{x})\tilde{v}(\bm{x})\,{\rm
d}x+\frac{\beta}{2}\int_{\Gamma}\tilde{v}^{2}(\bm{x})\,{\rm d}x,$ (2.16)

where $V=H^{1}(\Omega)$, $W=L^{2}(\Gamma)$ and $\beta$ is a positive
parameter.

##  3 An augmented Lagrangian deep learning method

Before presenting our augmented Lagrangian deep learning method, we first
recall the augmented Lagrangian method for the minimax problem (2.3) in a
Hilbert space setting [17].

The main idea here is to find the saddle point by solving the dual problem
(2.5). More precisely, one may first fix the primal variable $v_{k}$ at the
$k$th iteration, and update the dual variable to $\mu_{k+1}$. Then we fix the
approximate dual variable $\mu_{k+1}$ and find an approximate minimizer
$v_{k+1}$. The method is summarized in Algorithm 1.

Algorithm 1 The augmented Lagrangian method.

Input: The parameter $\beta_{0}>0$, tolerance $\tau_{0}>0$, increase parameter
$\alpha\geq 1$, initial guess $v_{0}$, $\mu_{0}$, max iteration number
$Epoch$.

for $k=0,1,\cdots,Epoch$ do

Update Lagrange multiplier function by
$\mu_{k+1}=\mu_{k}-\beta_{k}({B}v_{k}-g).\qquad(\star)$

Fix $\mu_{k+1}$ and find an approximate minimizer $v_{k+1}$ of
$\mathcal{L}_{\beta}(v,\mu_{k+1})$ such that
$\|\partial_{1}\mathcal{L}_{\beta}(v_{k+1},\mu_{k+1})\|\leq\tau_{k}$.

Update penalty parameter by $\beta_{k+1}=\alpha\beta_{k}$ .

Select new tolerance $\tau_{k+1}$.

end for

Output: $u=v_{Epoch+1}$, $\lambda=\mu_{Epoch+1}$.

We now present some basic concepts of the deep neural networks (DNNs). In this
paper, we shall adopt the residual neural network (ResNet) proposed in [20] to
approximate the variational problem. The ResNet can be formulated as follow:

$\bm{h}_{0}=\bm{V}\bm{x},\
\bm{h}_{\ell}=\bm{h}_{\ell-1}+\sigma(\bm{W}_{\ell}\bm{h}_{\ell-1}+\bm{b}_{\ell}),\
\ell=1,2,\dots,L,\ {\phi}(\bm{x};\bm{\theta})=\bm{a}^{T}\bm{h}_{L},$

where $\bm{V}\in\mathbb{R}^{N\times d}$, $\bm{W}_{\ell}\in\mathbb{R}^{N\times
N}$, $\bm{b}_{\ell}\in\mathbb{R}^{N}$ for $\ell=1,\dots,L$,
$\bm{a}\in\mathbb{R}^{N}$. $\sigma(x)$ is a non-linear activation function.
Here, $L$ is the depth of the ResNet, and $N$ is the width of the network.
$\bm{\theta}=\\{\bm{V},\bm{a},\bm{W}_{\ell},\bm{b}_{\ell}:1\leq\ell\leq L\\}$
denotes the set of all parameters in $\bm{\phi}$, which uniquely determines
the neural network.

To present the augmented Lagrangian deep learning method, we express the
primal and dual variables with two individual DNN functions, respectively,
i.e.,

$u(\bm{x})\approx\phi^{u}(\bm{x};\bm{\theta}_{u}),\quad\lambda(\bm{x})\approx\phi^{\lambda}(\bm{x};\bm{\theta}_{\lambda}).$

To derive the numerical solution, it suffices for us to determine parameters
$\bm{\theta}_{u}$ and $\bm{\theta}_{\lambda}$. To this end, we may closely
follow the strategies given in Algorithm 1. Notice that one can directly
follow Algorithm 1 to update the primal variable, however, for the dual
variable, Algorithm 1 is not directly applicable since we need to update the
parameters of the neural networks (not a function itself). In other words, we
cannot directly update the parameters $\bm{\theta}_{\mu}$ of the dual variable
$\phi^{\mu}$ through the equation $(\star),$ i.e.,

$\phi^{\mu}_{k+1}\,\Leftarrow\,\phi^{\mu}_{k}-\beta_{k}({B}\phi^{v}_{k}-g).$

To overcome this difficulty, we propose to solve the following least squares
problem:

$\displaystyle\bm{\theta}^{\mu}_{k+1}=\operatorname*{arg\,min}_{\bm{\theta}_{\nu}}J_{\lambda}(\phi^{\nu};\phi^{\mu}_{k},\phi^{v}_{k}),\quad\phi^{\nu}:=\phi^{\lambda}(\bm{x};\bm{\theta}_{\nu}),\
\phi^{\mu}_{k}:=\phi^{\lambda}(\bm{x};\bm{\theta}^{\mu}_{k}),\
\phi^{v}_{k}:=\phi^{u}(\bm{x};\bm{\theta}^{v}_{k});$ $\displaystyle
J_{\lambda}(\phi^{\nu};\phi^{\mu}_{k},\phi^{v}_{k})=\|\phi^{\nu}-\phi^{\mu}_{k}-\beta_{k}({B}\phi^{v}_{k}-g)\|_{L^{2}(\Gamma)}^{2}.$
(3.1)

Notice that the functionals $\mathcal{L}_{\beta}(v,\mu)$ and
$J_{\lambda}(\nu;\mu,v)$ can be expressed as integrals over $\Omega$ or
$\Gamma$ in most cases, meaning that sub-optimization problems involved can be
solved by means of the stochastic gradient descent method [5, 25].

To sum up the above discussions, we summarize our augmented Lagrangian deep
learning method in Algorithm 2.

Algorithm 2 The augmented Lagrangian deep learning method (ALDL).

Input: The parameter $\beta_{0}>0$, the increase parameter $\alpha\geq 1$, the
iteration number $Epoch$, the iteration number of model $Epoch_{u}$, the
iteration number of Lagrange multiplier $Epoch_{\lambda}$.

Initialize the network parameters $\bm{\theta}^{v}_{0}$ and
$\bm{\theta}^{\mu}_{0}$ following the default random initialization of
PyTorch.

for $k=0,1,\cdots,Epoch$ do

Fix the parameter $\bm{\theta}^{v}_{k}$ and the parameter
$\bm{\theta}^{\mu}_{k}$. With the initial guess $\bm{\theta}^{\mu}_{k}$, find
an approximate solver $\bm{\theta}^{\mu}_{k+1}$ of
$\min_{\bm{\theta}_{\nu}}J_{\lambda}\big{(}\phi^{\nu};\phi^{\mu}_{k},\phi^{v}_{k}\big{)}$
by the Adam method with $Epoch_{\lambda}$ iteration steps.

Fix the parameter $\bm{\theta}^{\mu}_{k+1}$. With the initial guess
$\bm{\theta}^{v}_{k}$, find an approximate solver $\bm{\theta}^{v}_{k+1}$ of
$\min_{\bm{\theta}_{u}}\mathcal{L}_{\beta}\big{(}\phi^{u}(\cdot\
;\bm{\theta}_{u}),\phi^{\mu}_{k+1}\big{)}$ by the Adam method with $Epoch_{u}$
iteration steps.

Update penalty parameter by $\beta_{k+1}=\alpha\beta_{k}$.

end for

Output: $u=\phi^{v}_{Epoch+1}$, $\lambda=\phi^{\mu}_{Epoch+1}$.

######  Remark 3.1.

Another strategy to update the dual variable $\phi^{\mu}$ is solving the
minimax problem (2.3) in the form of DNNs by the stochastic gradient descent
ascent method (SGDA) directly. However, how to solve this problem efficiently
is still an interesting problem and one can refer to[29, 38, 37, 22] and the
reference therein for more details. We only make some numerical comparisons of
the ALDL method and the SGDA method in Section 4.

Next, we present some details for implementing the ALDL method. For the the
elliptic problem (2.6), we have

$\displaystyle\mathcal{L}_{\beta}(v,\mu)=$
$\displaystyle\frac{1}{2}\int_{\Omega}\big{[}\bm{A}(\bm{x})\nabla
v(\bm{x})\cdot\nabla
v(\bm{x})+c(\bm{x})v^{2}(\bm{x})-2f(\bm{x})v(\bm{x})\big{]}\,{\rm d}x$
$\displaystyle-\int_{\Gamma}\big{[}\mu(\bm{x})\big{(}v(\bm{x})-g(\bm{x})\big{)}\big{]}\,{\rm
d}x+\frac{\beta}{2}\int_{\Gamma}\big{[}v(\bm{x})-g(\bm{x})\big{]}^{2}\,{\rm
d}x$ $\displaystyle=$
$\displaystyle|\Omega|\mathbb{E}_{\bm{\xi}}\left[\frac{1}{2}\bm{A}(\bm{\xi})|\nabla
v(\bm{\xi})|^{2}+\frac{1}{2}c(\bm{\xi})v^{2}(\bm{\xi})-f(\bm{\xi})v(\bm{\xi})\right]$
$\displaystyle-|\Gamma|\mathbb{E}_{\bm{\eta}}\left[\mu(\bm{\eta})\big{(}v(\bm{\eta})-g(\bm{\eta})\big{)}+\frac{\beta}{2}\big{(}v(\bm{\eta})-g(\bm{\eta})\big{)}^{2}\right],$
(3.2)

where $\bm{\xi}$ and $\bm{\eta}$ are random vectors following the uniform
distribution over $\Omega$ and $\Gamma$, respectively. $|\Omega|$ and
$|\Gamma|$ are the measure of $\Omega$ and $\Gamma$, respectively.

By the definition of $J_{\lambda}(\nu;\mu,v)$, we have

$\displaystyle J_{\lambda}(\nu;\mu_{k},v_{k})$
$\displaystyle=\int_{\Gamma}\left[\nu(\bm{x})-\mu_{k}(\bm{x})-\beta\big{(}v_{k}(\bm{x})-g(\bm{x})\big{)}\right]^{2}\,{\rm
d}x$
$\displaystyle=|\Gamma|\mathbb{E}_{\bm{\eta}}\left[\big{(}\nu(\bm{\eta})-\mu_{k}(\bm{\eta})-\beta(v_{k}(\bm{\eta})-g(\bm{\eta}))\big{)}^{2}\right].$
(3.3)

Upon substituting (3.2) and (3.3) into Algorithm 2, one can then obtain the
ALDL algorithm.

As for the linear eigenvalue problem, a direct manipulation gives

$\displaystyle\mathcal{L}_{\beta}(v,\mu)$
$\displaystyle=\frac{\mathbb{E}_{\bm{\xi}}\big{[}p(\bm{\xi})|\nabla{v}(\bm{\xi})|^{2}+q(\bm{\xi})v^{2}(\bm{\xi})\big{]}}{\mathbb{E}_{\bm{\zeta}}\big{[}v^{2}(\bm{\zeta})\big{]}}-\frac{|\Gamma|\mathbb{E}_{\bm{\eta}}\big{[}\mu(\bm{\eta}){v}(\bm{\eta})\big{]}}{\sqrt{|\Omega|\mathbb{E}_{\bm{\zeta}}\big{[}v^{2}(\bm{\zeta})\big{]}}}+\frac{\beta|\Gamma|\mathbb{E}_{\bm{\eta}}\big{[}v^{2}(\bm{\eta})\big{]}}{2|\Omega|\mathbb{E}_{\bm{\zeta}}\big{[}v^{2}(\bm{\zeta})\big{]}},$
(3.4) $\displaystyle J_{\lambda}(\nu;\mu_{k},\tilde{v}_{k})$
$\displaystyle=|\Gamma|\mathbb{E}_{\bm{\eta}}\left[\big{(}\nu(\bm{\eta})-\mu_{k}(\bm{\eta})-\beta\tilde{v}_{k}(\bm{\eta})\big{)}^{2}\right]=|\Gamma|\mathbb{E}_{\bm{\eta}}\left[\left(\nu(\bm{\eta})-\mu_{k}(\bm{\eta})-\beta\frac{{v}_{k}(\bm{\eta})}{\|v_{k}\|_{0,\Omega}}\right)^{2}\right],$
(3.5)

where $\bm{\xi}$, $\bm{\zeta}$ are i.i.d. random vectors following the uniform
distribution over $\Omega$, and $\bm{\eta}$ is the random vector following the
uniform distribution over $\Gamma$. Notice that $\tilde{v}_{k}$ is fixed in
$J_{\lambda}(\nu;\mu_{k},\tilde{v}_{k})$, so $\|v_{k}\|_{0,\Omega}$ is a known
real number in $J_{\lambda}(\nu;\mu_{k},\tilde{v}_{k}).$ Again, one may obtain
the ALDL algorithm for eigenvalue problems by substituting the above formulas
into Algorithm 2.

For the nonlinear eigenvalue problem (2.13), we have

$\displaystyle\mathcal{L}_{\beta}(v,\mu)$
$\displaystyle=\frac{\mathbb{E}_{\bm{\xi}_{1}}\big{[}\bm{A}(\bm{\xi}_{1})|\nabla{v}(\bm{\xi}_{1})|^{2}+V(\bm{\xi}_{1})v^{2}(\bm{\xi}_{1})\big{]}}{\mathbb{E}_{\bm{\xi}_{2}}\big{[}v^{2}(\bm{\xi}_{2})\big{]}}+\frac{\mathbb{E}_{\bm{\xi}_{1}}\big{[}v^{4}(\bm{\xi}_{1})\big{]}}{|\Omega|\mathbb{E}_{\bm{\xi}_{2},\bm{\xi}_{3}}\big{[}v^{2}(\bm{\xi}_{2})v^{2}(\bm{\xi}_{3})\big{]}}$
$\displaystyle-\frac{|\Gamma|\mathbb{E}_{\bm{\eta}}\big{[}\mu(\bm{\eta}){v}(\bm{\eta})\big{]}}{\sqrt{|\Omega|\mathbb{E}_{\bm{\xi}_{2}}\big{[}v^{2}(\bm{\xi}_{2})\big{]}}}+\frac{\beta|\Gamma|\mathbb{E}_{\bm{\eta}}\big{[}v^{2}(\bm{\eta})\big{]}}{2|\Omega|\mathbb{E}_{\bm{\xi}_{2}}\big{[}v^{2}(\bm{\xi}_{2})\big{]}},$
(3.6) $\displaystyle J_{\lambda}(\nu;\mu_{k},\tilde{v}_{k})$
$\displaystyle=|\Gamma|\mathbb{E}_{\bm{\eta}}\left[\big{(}\nu(\bm{\eta})-\mu_{k}(\bm{\eta})-\beta\tilde{v}_{k}(\bm{\eta})\big{)}^{2}\right]=|\Gamma|\mathbb{E}_{\bm{\eta}}\left[\left(\nu(\bm{\eta})-\mu_{k}(\bm{\eta})-\beta\frac{{v}_{k}(\bm{\eta})}{\|v_{k}\|_{0,\Omega}}\right)^{2}\right],$
(3.7)

where $\bm{\xi}_{1}$, $\bm{\xi}_{2}$, $\bm{\xi}_{3}$ are i.i.d. random vectors
produced by the uniform distribution over $\Omega$, and $\bm{\eta}$ is the
random vector following the uniform distribution over $\Gamma$.

##  4 Numerical experiments

In this section, we shall present various numerical examples to illustrate the
effectiveness of the ALDL method. We shall also perform a numerical comparison
between the penalty method (PMDL) and the ALDL method. In all our numerical
examples, we shall use the ResNet with width $N=50$ and depth $L=6$ with an
activation function $\sigma(x)=\max\\{x,0\\}^{2}$ to approximate the solution
$u$, which results in 15350 unknown parameters. To reduce the computation
cost, we use a ResNet with width $N=50$ and depth $L=2$ with an activation
function $\sigma(x)=\max\\{x,0\\}^{2}$ as an approximation to the dual
variable $\lambda$, which results in 5150 parameters.

We use the Adam optimizer [25] for training, with a learning rate
$\eta=1\text{e}-3$ for elliptic PDEs and linear eigenvalue problems, and
$\eta=5\text{e}-4$ for nonlinear eigenvalue problems. The learning rate in the
optimization is adjusted in an exponentially decaying scheme, where the
decaying rate is $0.01^{\frac{1}{50000}}$. The maximum iteration number
$Epoch$ is taken as 50000 for the PMDL method, and we set $Epoch=500$,
$Epoch_{\lambda}=100$ and $Epoch_{u}=100$ for the ALDL method. The batch size
within the computation domain is 512 for 2d problems and 2048 for 3d problems.
While the number of training points on each boundary is set to be 64 for 2d
problems and is set to be 256 for 3d problems. All numerical experiments are
implemented in Python 3.7 using Pytorch 1.3 in an NVIDIA GEFORCE RTX 2080 Ti
GPU card. The code can be shared upon request.

For elliptic PDEs, we use $u$ and $u_{dl}$ to denote the exact solution and
the deep learning solution, respectively. While for eigenvalue problems, we
let $\rho$ and $u$ be the smallest eigenvalue and the corresponding
eigenfunction with $\|u\|_{0}=1$, respectively. Also, we denote by $\rho_{dl}$
and $u_{dl}$ the DNNs solutions. For a computational domain
$\Omega=(0,1)^{d}$, we divide it into small cubes with a length $h$ uniformly
and denote the set of all vertexes as $\Omega_{h}$, i.e.,
$\Omega_{h}=\\{ih:0\leq i\leq N\\}^{d}$ with $N=1/h$. In all experiments, we
shall adopt $h=2^{-6}$ to generate the set $\Omega_{h}$ as the test locations.
In order to measure the accuracy of deep learning algorithms, we introduce the
discrete maximum norm

$\|v\|_{0,\infty,h}=\max_{\bm{x}\in E}|v(\bm{x})|,\quad E\subset\Omega_{h}.$

In addition, we define the absolute error and relative error in the domain
$\Omega$ as follows:

$\mathcal{E}_{a}^{in}=\|u(\bm{x})-u_{dl}(\bm{x})\|_{0,\infty,h},\quad\mathcal{E}_{r}^{in}=\left\|\frac{u(\bm{x})-u_{dl}(\bm{x})}{u(\bm{x})}\right\|_{0,\infty,h},\quad\bm{x}\in
E=\Omega\cap\Omega_{h}.$

The absolute error and relative error on the boundary $\Gamma$ are defined by

$\mathcal{E}_{a}^{bd}=\|u(\bm{x})-u_{dl}(\bm{x})\|_{0,\infty,h},\quad\mathcal{E}_{r}^{bd}=\left\|\frac{u(\bm{x})-u_{dl}(\bm{x})}{u(\bm{x})}\right\|_{0,\infty,h},\quad\bm{x}\in
E=\Gamma\cap\Omega_{h}.$

###  4.1 Elliptic PDEs

Consider the Poisson equation with the Dirichlet boundary

$\begin{cases}-\Delta u=f\quad&\mbox{\rm in}\ \Omega,\\\ \quad
u=g\quad&\mbox{\rm on}\ \Gamma,\end{cases}$

where $\Omega=(0,1)^{d}$ for $d=2$ and $3$. We choose appropriate $f$ and $g$
such that the exact solution yields

$u(\bm{x})=\sum_{i=1}^{d}\big{(}\sin(2\pi x_{i})+1.25\big{)}.$

####  4.1.1 The 2D case

We first solve the above problem by the ALDL method with different parameters
$\beta$. The test errors are presented in Figure 1 and in Table 2. It can be
seen from Figure 1 (left) that the ALDL method is rather robust with respect
to $\beta,$ and there is no evident difference on the accuracy (within the
domain) and convergent speed when $\beta$ is chosen from $10$ to $1000$. In
other words, one may use a relatively small $\beta$ in practice. By Figure 1
(right), we can also see the ALDL method admits a very good accuracy on the
boundary.

Figure 1: The relative error in the domain (left) and on the boundary (right)
with ALDL.

$\beta$ | $\mathcal{E}_{r}^{in}$ | $\mathcal{E}_{r}^{bd}$
---|---|---
1e+1 | 7.7843e-3 | 1.9851e-3
1e+2 | 1.0998e-2 | 1.9644e-3
1e+3 | 1.6788e-2 | 2.0246e-3
Table 1:  The final relative errors of ALDL.

$\beta$ | $\mathcal{E}_{r}^{in}$ | $\mathcal{E}_{r}^{bd}$
---|---|---
2e+2 | 4.8946e-2 | 2.2827e-2
2e+3 | 1.2054e-2 | 3.7951e-3
2e+4 | 4.0631e-1 | 1.5662e-3
Table 2: The final relative errors of PMDL.

We also solve the problem with the PMDL method with different parameters
$\beta$, and the numerical results are shown in Figure 2 and Table 2. By
Figure 2 (right) we can observe that one should use a large parameter $\beta$
to get a good accuracy on the boundary. However, by Figure 2 (left) we can see
that the large the parameter $\beta$ is, the slower the converge rate. This
may due to the ill-conditioning for a relatively large $\beta$ in the penalty
methods [31].

Figure 2: The relative error in the domain (left) and on the boundary (right)
of PMDL.

As mentioned in Remark 3.1, one can solve the minimax problem (2.8) in the
form of DNNs by the SGDA method. To compare the numerical performance between
the ALDL method and the SGDA method, we list the errors and GPU time for
different methods with the same parameters in Table 3. We learn that the ALDL
method and the SGDA method admit a similar accuracy, however, the ALDL method
saves about $25\%$ more computational cost than the SGDA method.

Method | $\beta$ | $\mathcal{E}_{r}^{in}$ | $\mathcal{E}_{r}^{bd}$ | time(s)
---|---|---|---|---
ALDL | 1e+2 | 1.0998e-2 | 1.9644e-3 | 933.82
SGDA | 1e+2 | 1.0651e-2 | 2.8226e-3 | 1256.22
ALDL | 1e+3 | 1.6788e-2 | 2.0246e-3 | 938.39
SGDA | 1e+3 | 1.7710e-2 | 2.4454e-3 | 1223.07
Table 3: The comparation of different methods.

####  4.1.2 The 3D case

We next consider the three dimensional case. Similar plots are shown in Figure
3 and Figure 4, and the final approximation errors are listed in Table 5 and
Table 5. From those pictures and tables, one can draw similar conclusions as
in the two dimensional case. When both methods provide good approximations,
the accuracy of the ALDL method is about $2$ times higher than that of the
PMDL method. However, in the case of large parameters for PMDL method, the
ALDL method outperforms the PMDL method about $50$ times.

Figure 3:  The relative error in the domain (left) and on the boundary (right)
of ALDL.

$\beta$ | $\mathcal{E}_{r}^{in}$ | $\mathcal{E}_{r}^{bd}$
---|---|---
2e+0 | 6.6582e-2 | 3.5099e-3
1e+1 | 4.0046e-2 | 6.7436e-3
1e+2 | 5.9976e-2 | 1.3568e-2
Table 4:  The final relative errors of ALDL.

$\beta$ | $\mathcal{E}_{r}^{in}$ | $\mathcal{E}_{r}^{bd}$
---|---|---
2e+2 | 5.6320e-2 | 2.6998e-2
2e+3 | 2.0876e-1 | 9.2283e-3
2e+4 | 3.5052e-1 | 1.0640e-2
Table 5:  The final relative errors of PMDL.

Figure 4:  The relative error in the domain (left) and on the boundary (right)
of PMDL.

###  4.2 Linear eigenvalue problems

We next consider the following eigenvalue problem

$\begin{cases}-\Delta u=\rho u\quad&\mbox{\rm in}\ \Omega,\\\ \quad
u=0\quad&\mbox{\rm on}\ \Gamma,\end{cases}$

where $\Omega=(0,1)^{d}$ for $d=2$ and $3$. The exact smallest eigenvalue is
$d\pi^{2}$ and the corresponding eigenfunction is

$u(\bm{x})=\prod_{i=1}^{d}\sin(\pi(x_{i}-1)).$

Due to the homogeneous Dirichlet boundary, we evaluate different approaches
with absolute errors for eigenfunctions (and relative errors for eigenvalues).

####  4.2.1 The 2D case

We present the numerical results by the ALDL method with different parameters
$\beta$ in Figure 5 and Table 7. We can see that the numerical error is
approximately $10^{-3}$ for the eigenfunctions and $10^{-4}$ for eigenvalues.
Moreover, the approximation accuracy and convergence rate are robust with
respect to different $\beta.$

Figure 5: The absolute error in the domain (left) and on the boundary (right)
of ALDL.

$\beta$ | $\mathcal{E}_{a}^{in}$ | $\mathcal{E}_{a}^{bd}$ | $|\rho_{dl}-\rho|/\rho$
---|---|---|---
2e+1 | 5.1201e-3 | 1.6031e-3 | 2.0678e-4
2e+2 | 7.7265e-3 | 1.4665e-3 | 7.3460e-4
2e+3 | 6.9291e-3 | 1.4485e-3 | 2.1878e-4
Table 6:  The final absolute errors of ALDL.

$\beta$ | $\mathcal{E}_{a}^{in}$ | $\mathcal{E}_{a}^{bd}$ | $|\rho_{dl}-\rho|/\rho$
---|---|---|---
2e+2 | 5.7031e-2 | 6.0427e-2 | 7.4183e-2
2e+3 | 9.9260e-3 | 7.0166e-3 | 7.0739e-3
2e+4 | 1.8037e-2 | 1.0109e-3 | 9.7971e-4
Table 7:  The final absolute errors of PMDL.

We also show the associated results by the PMDL method in Figure 6 and Table
7, from which we may conclude that the accuracy and convergent speed of the
PMDL method are relatively sensitive to the parameter $\beta.$

Figure 6:  The absolute error in the domain (left) and on the boundary (right)
of PMDL.

We conclude from Table 7 and Table 7 that the ALDL method has 2-20 times
advantages over the PMDL method in terms of the accuracy of eigenfunctions or
eigenvalues. Similar to the PDEs, we compare the numerical performance of the
ALDL method and the SGDA method in Table 8, from which we find that the
approximations of the former are twice as accurate as that of the latter and
the ALDL takes $30\%$ less time.

Method | $\beta$ | $\mathcal{E}_{a}^{in}$ | $\mathcal{E}_{a}^{bd}$ | $|\rho_{dl}-\rho|/\rho$ | time(s)
---|---|---|---|---|---
ALDL | 2e+2 | 7.7265e-3 | 1.4665e-3 | 7.3460e-4 | 1214.60
SGDA | 2e+2 | 7.6285e-3 | 2.1653e-3 | 1.0547e-3 | 1695.77
ALDL | 2e+3 | 6.9291e-3 | 1.4485e-3 | 2.1878e-4 | 1008.00
SGDA | 2e+3 | 9.6389e-3 | 3.0370e-3 | 1.1833e-3 | 1652.99
Table 8:  The comparation of different methods.

####  4.2.2 The 3D case

We have also performed the experiments for the three dimensional case, and the
results are shown in Figure 7, Figure 8, and Table 10 \- Table 10. Similar to
the two dimensional case, the ALDL method admits a better numerical
performance in view of the accuracy. Especially, we noticed that when the
parameter $\beta$ of the PMDL method is small, the eigenvalues obtained by the
ALDL method at least are more $100$ times accurate than those obtained by the
PMDL method. When the parameter $\beta$ of the PMDL method is large, the
eigenfunctions obtained by the ALDL method at least are more $10$ times
accurate than those obtained by the PMDL method.

Figure 7:  The absolute error in the domain (left) and on the boundary (right)
of ALDL.

$\beta$ | $\mathcal{E}_{a}^{in}$ | $\mathcal{E}_{a}^{bd}$ | $|\rho_{dl}-\rho|/\rho$
---|---|---|---
1e+1 | 1.2117e-2 | 2.1592e-2 | 1.4851e-3
1e+2 | 1.3606e-2 | 2.3021e-2 | 4.2040e-4
1e+3 | 2.3689e-2 | 1.7151e-2 | 1.7392e-3
Table 9:  The final absolute errors of ALDL.

$\beta$ | $\mathcal{E}_{a}^{in}$ | $\mathcal{E}_{a}^{bd}$ | $|\rho_{dl}-\rho|/\rho$
---|---|---|---
2e+2 | 9.1493e-2 | 8.9782e-2 | 7.6310e-2
2e+3 | 2.5106e-2 | 1.5402e-2 | 7.2753e-3
5e+4 | 1.9755e-1 | 1.1124e-2 | 9.5417e-3
Table 10:  The final absolute errors of PMDL.

Figure 8:  The absolute error in the domain (left) and on the boundary (right)
of PMDL.

###  4.3 A nonlinear eigenvalue problem

Finally, we consider the nonlinear Schrödinger eigenvalue problem

$\begin{cases}-\Delta u+Vu+u^{3}=\rho u\quad&\mbox{\rm in}\ \Omega,\\\ \quad
u=0\quad&\mbox{\rm on}\ \Gamma,\\\ \|u\|_{0,\Omega}=1,\end{cases}$

where $\Omega=(0,1)^{d}$ for $d=2$ and $3$, and $V=\sum_{i=1}^{d}x_{i}^{2}$.
Since the true solutions $\rho$ and $u$ are unavailable, we shall numerically
compute the reference solutions $\rho_{ref}$ and $u_{ref}.$ More precisely, we
construct a DNN reference solution $u_{ref}$ by:

$u_{ref}(\bm{x};\bm{\theta})=\ell(\bm{x})\psi(\bm{x};\bm{\theta}),$

where $\psi(\bm{x};\bm{\theta})$ is a ResNet function with width 50 and depth
6, and

$\ell(\bm{x})=\prod_{i=1}^{d}x_{i}(1-x_{i}).$

In this way the reference solution admits an exact match on the boundary. The
reference solution $u_{ref}(\bm{x};\bm{\theta})$ is trained by the Adam
optimizer with a learning rate $\eta=5e-4,$ and we set $Epoch=100000$ with
batch size $2048$ in the domain for $2$d problem and $8192$ for $3$d problem.

####  4.3.1 The 2D case

The numerical results for the ALDL method with different parameters $\beta$
are presented in Figure 9 and Table 12. Again, it is noticed that the method
is insensitive to the choice of $\beta.$

Figure 9:  The absolute error in the domain (left) and on the boundary (right)
of ALDL.

$\beta$ | $\mathcal{E}_{a}^{in}$ | $\mathcal{E}_{a}^{bd}$ | $|\rho_{dl}-\rho_{ref}|/\rho_{ref}$
---|---|---|---
2e+0 | 5.8958e-3 | 4.8274e-4 | 3.5429e-3
2e+1 | 7.3927e-3 | 1.8059e-3 | 1.6044e-4
2e+2 | 9.2868e-3 | 2.9953e-3 | 5.2254e-4
Table 11:  The final absolute errors of ALDL.

$\beta$ | $\mathcal{E}_{a}^{in}$ | $\mathcal{E}_{a}^{bd}$ | $|\rho_{dl}-\rho_{ref}|/\rho_{ref}$
---|---|---|---
2e+2 | 3.3200e-2 | 3.3776e-2 | 4.0211e-2
2e+3 | 1.9137e-2 | 4.1947e-3 | 3.5401e-3
2e+4 | 1.2269e-1 | 2.6753e-3 | 4.1881e-3
Table 12:  The final absolute errors of PMDL.

Figure 10:  The absolute error in the domain (left) and on the boundary
(right) of PMDL.

The associated numerical results by the PMDL method are presented in Figure 10
and Table 12, from which we can see that the larger the $\beta$ is, the slower
the converge rate seems to be. In addition, the ALDL method has $2-20$ times
advantages over the PMDL method in terms of the accuracy of eigenfunctions. As
for eigenvalues, the ALDL method has more than $100$ times the accuracy
advantage over the PMDL method in some cases.

####  4.3.2 The 3D case

We have also presented the 3D simulations, and the results are shown in Figure
11 \- Figure 12 and Table 14 \- Table 14, and one can draw similar conclusions
as in the above examples.

Figure 11:  The absolute error in the domain (left) and on the boundary
(right) of ALDL.

$\beta$ | $\mathcal{E}_{a}^{in}$ | $\mathcal{E}_{a}^{bd}$ | $|\rho_{dl}-\rho_{ref}|/\rho_{ref}$
---|---|---|---
2e+0 | 2.4886e-2 | 1.6198e-2 | 3.0750e-3
2e+1 | 3.1448e-2 | 3.1834e-2 | 4.0040e-4
1e+2 | 4.7370e-2 | 4.5080e-2 | 7.2167e-4
Table 13: The final absolute errors of ALDL.

$\beta$ | $\mathcal{E}_{a}^{in}$ | $\mathcal{E}_{a}^{bd}$ | $|\rho_{dl}-\rho_{ref}|/\rho_{ref}$
---|---|---|---
2e+2 | 4.9407e-2 | 5.0387e-2 | 4.0671e-2
2e+3 | 1.6322e-1 | 2.4755e-2 | 7.3079e-3
2e+4 | 4.3183e-1 | 2.3657e-2 | 1.0364e-2
Table 14:  The final absolute errors of PMDL.

Figure 12:  The absolute error in the domain (left) and on the boundary
(right) of PMDL.

##  5 Concluding remarks

We have proposed an augmented Lagrangian deep learning method for variational
problems with essential boundary conditions. The approach relies on first
rewriting the original problem into an equivalent minimax problem, and then
expressing the primal and dual variables with two individual DNN functions.
Then, the network parameters of the primal and dual variables are trained
using the stochastic optimization method together with a projection technique.
Applications to elliptic problems and eigenvalue problems show that the ALDL
method admits many advantages over the penalty method. In our future studies,
we shall extend our ALDL approach to time dependent problems with complex
solution structures.

## Acknowledgments

J. Huang is partially supported by the National Key Research and Development
Project (Grant No. 2020YFA0709800), NSFC (Grant No. 12071289) and Shanghai
Municipal Science and Technology Major Project (2021SHZDZX0102). T. Zhou is
supported by the National Key R&D Program of China (2020YFA0712000), NSFC
(under grant numbers 11822111, 11688101), the science challenge project (No.
TZ2018001), and youth innovation promotion association (CAS).

## References

  * [1] R.A. Adams.  Sobolev Spaces.  Academic Press, New York, 1975.
  * [2] A. Ambrosetti and A. Malchiodi.  Nonlinear Analysis and Semilinear Elliptic Problems.  Cambridge University Press, Cambridge, 2007.
  * [3] A.R. Barron.  Universal approximation bounds for superpositions of a sigmoidal function.  IEEE Trans. Inform. Theory, 39:930–945, 1993.
  * [4] J. Berg and K. Nyström.  A unified deep artificial neural network approach to partial differential equations in complex geometries.  Neurocomputing, 317:28 – 41, 2018.
  * [5] L. Bottou.  Large-scale machine learning with stochastic gradient descent.  In COMPSTAT 2010, pages 177–186, 2010.
  * [6] S.C. Brenner and L.R. Scott.  The Mathematical Theory of Finite Element Methods.  Springer, New York, 1994.
  * [7] E. Cancès, R. Chakir, and Y. Maday.  Numerical analysis of nonlinear eigenvalue problems.  J. Sci. Comput., 45:90–117, 2010.
  * [8] J. Céa.  Lectures on Optimization–Theory and Algorithms.  Tata Institute of Fundamental Research, Bombay, 1978.
  * [9] F. Chen, J. Huang, C. Wang, and H. Yang.  Friedrichs learning: Weak solutions of partial differential equations via deep learning.  arXiv e-prints, page arXiv:2012.08023, 2020.
  * [10] P.G. Ciarlet.  Linear and Nonlinear Functional Analysis with Applications.  SIAM, Philadelphia, 2013.
  * [11] P.G. Ciarlet.  The Finite Element Method for Elliptic Problems.  North-Holland, Amsterdam, 1978.
  * [12] M.W.M.G. Dissanayake and N. Phan-Thien.  Neural-network-based approximations for solving partial differential equations.  Commun. Numer. Meth. Engng., 10:195–201, 1994.
  * [13] W. E.  Machine learning and computational mathematics.  Commun. Comput. Phys., 28:1639–1670, 2020.
  * [14] W. E, J. Han, and A. Jentzen.  Deep learning-based numerical methods for high-dimensional parabolic partial differential equations and backward stochastic differential equations.  Commun. Math. Stat., 5:349–380, 2017.
  * [15] W. E, C. Ma, and L. Wu.  A priori estimates of the population risk for two-layer neural networks.  Commun. Math. Sci., 17:1407–1425, 2019.
  * [16] W. E and B. Yu.  The deep Ritz method: a deep learning-based numerical algorithm for solving variational problems.  Commun. Math. Stat., 6:1–12, 2018.
  * [17] R. Glowinski and P. Le Tallec.  Augmented Lagrangian and Operator-splitting Methods in Nonlinear Mechanics.  SIAM, Philadelphia, 1989.
  * [18] J. Han, A. Jentzen, and W. E.  Solving high-dimensional partial differential equations using deep learning.  Proc. Natl. Acad. Sci. USA, 115:8505–8510, 2018.
  * [19] J. He, L. Li, J. Xu, and C. Zheng.  Relu deep neural networks and linear finite elements.  J. Comput. Math., 38:502–527, 2020.
  * [20] K. He, X. Zhang, S. Ren, and J. FigSun.  Deep residual learning for image recognition.  In CVPR 2016, pages 770–778, 2016.
  * [21] K. Hornik, M. Stinchcombe, and H. White.  Multilayer feedforward networks are universal approximators.  Neural Networks, 2:359 – 366, 1989.
  * [22] F. Huang and X. Wu.  Efficient mirror descent ascent methods for nonsmooth minimax problems.  In NeurIPS 2021, 2021.
  * [23] J. Huang, H. Wang, and H. Yang.  Int-deep: A deep learning initialized iterative method for nonlinear problems.  J. Comput. Phys., 419:109675, 2020.
  * [24] G.E. Karniadakis, I.G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, and L. Yang.  Physics-informed machine learning.  Nat. Rev. Phys., 3:422–440, 2021.
  * [25] D.P. Kingma and J. Ba.  Adam: A method for stochastic optimization.  arXiv e-prints, page arXiv:1412.6980, 2014.
  * [26] H. Lee and I.S. Kang.  Neural algorithm for solving differential equations.  J. Comput. Phys., 91:110–131, 1990.
  * [27] Y. Liao and P. Ming.  Deep Nitsche method: Deep Ritz method with essential boundary conditions.  Commun. Comput. Phys., 29:1365–1384, 2021.
  * [28] Z. Liu, W. Cai, and Z.-Q. J. Xu.  Multi-scale deep neural network (MscaleDNN) for solving Poisson-Boltzmann equation in complex domains.  Commun. Comput. Phys., 28:1970–2001, 2020.
  * [29] A. Nemirovski, A. Juditsky, G. Lan, and A. Shapiro.  Robust stochastic approximation approach to stochastic programming.  SIAM J. Optim., 19:1574–1609, 2008.
  * [30] J. Nitsche.  Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind.  Abh. Math. Sem. Univ. Hamburg, 36:9–15, 1971.
  * [31] J. Nocedal and S. J. Wright.  Numerical Optimization (second edition).  Springer, New York, 2006.
  * [32] M. Raissi, P. Perdikaris, and G.E. Karniadakis.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  J. Comput. Phys., 378:686 – 707, 2019.
  * [33] Z. Shen, H. Yang, and S. Zhang.  Deep network approximation characterized by number of neurons.  Commun. Comput. Phys., 28:1768–1811, 2020.
  * [34] H. Sheng and C. Yang.  PFNN: A penalty-free neural network method for solving a class of second-order boundary-value problems on complex geometries.  J. Comput. Phys., 428:110085, 2021.
  * [35] J. Sirignano and K. Spiliopoulos.  DGM: A deep learning algorithm for solving partial differential equations.  J. Comput. Phys., 375:1339 – 1364, 2018.
  * [36] R. Stenberg.  On some techniques for approximating boundary conditions in the finite element method.  In International Symposium on Mathematical Modelling and Computational Methods Modelling 94, volume 63, pages 139–148. 1995.
  * [37] W. Xian, F. Huang, Y. Zhang, and H. Huang.  A faster decentralized algorithm for nonconvex minimax problems.  In NeurIPS 2021, 2021.
  * [38] T. Yoon and E. K. Ryu.  Accelerated algorithms for smooth convex-concave minimax problems with ${O}{(1/k^{2})}$ rate on squared gradient norm.  In ICML 2021, pages 12098–12109, 2021.
  * [39] Y. Zang, G. Bao, X. Ye, and H. Zhou.  Weak adversarial networks for high-dimensional partial differential equations.  J. Comput. Phys., 411:109409, 2020.
  * [40] D. Zhang, L. Guo, and G. E. Karniadakis.  Learning in modal space: solving time-dependent stochastic PDEs using physics-informed neural networks.  SIAM J. Sci. Comput., 42:A639–A665, 2020.