solutionSolutionsolutionfile solutionfile[Papazafeiropoulos_WCL2022-0862]

#  Coverage Probability of STAR-RIS assisted Massive MIMO systems with
Correlation and Phase Errors

Anastasios Papazafeiropoulos, Zaid Abdullah, Pandelis Kourtessis, Steven
Kisseleff, Ioannis Krikidis This work was co-funded by the European Regional
Development Fund and the Republic of Cyprus through the Research and
Innovation Foundation, under the projects INFRASTRUCTURES/1216/0017 (IRIDA)
and EXCEL- LENCE/0918/0377 (PRIME). Also, this work was supported by the
Luxembourg National Research Fund (FNR) through the CORE Project under Grant
RISOTTI C20/IS/14773976.A. Papazafeiropoulos is with the Communications and
Intelligent Systems Research Group, University of Hertfordshire, Hatfield AL10
9AB, U. K., and with SnT at the University of Luxembourg, Luxembourg. P.
Kourtessis is with the Communications and Intelligent Systems Research Group,
University of Hertfordshire, Hatfield AL10 9AB, U. K. Zaid Abdullah and Steven
Kisseleff are with the SnT at the University of Luxembourg, Luxembourg. I.
Krikidis is with the IRIDA Research Centre for Communication Technologies,
Department of Electrical and Computer Engineering, University of Cyprus,
Cyprus. E-mails<EMAIL_ADDRESS>{zaid.abdullah,
<EMAIL_ADDRESS>p.kourtessis<EMAIL_ADDRESS>
###### Abstract

In this paper, we investigate a simultaneous transmitting and reflecting
reconfigurable intelligent surface (STAR-RIS) assisting a massive multiple-
input multiple-output (mMIMO) system. In particular, we derive a closed-form
expression for the coverage probability of a STAR-RIS assisted mMIMO system
while accounting for correlated fading and phase-shift errors. Notably, the
phase configuration takes place at every several coherence intervals by
optimizing the coverage probability since the latter depends on statistical
channel state information (CSI) in terms of large-scale statistics. As a
result, we achieve a reduced complexity and overhead for the optimization of
passive beamforming, which are increased in the case of STAR-RIS networks with
instantaneous CSI. Numerical results corroborate our analysis, shed light on
interesting properties such as the impact of the number of RIS elements and
the effect of phase errors, along with affirming the superiority of STAR-RIS
against reflective-only RIS.

###### Index Terms:

Reconfigurable intelligent surface (RIS), simultaneous transmission and
reflection, spatial correlation, coverage probability, 6G networks.

##  I Introduction

Reconfigurable intelligent surface (RIS) is a promising technology, enabling a
smart radio environment (SRE) to meet the challenges of sixth-generation (6G)
wireless networks [1]. It consists of a number of low-cost nearly passive
elements that can modify the phase shifts and even the amplitude of the
impinging signals to realize an SRE through a backhaul controller. RIS
favorable characteristics have led recently the research interest in meeting
certain significant targets such as maximization of the spectral and energy
efficiencies [2, 3], and maximization of the coverage probability [4, 5].
Notably, given that 6G will be built upon existing technologies in 5G such as
massive multiple-input multiple-output (mMIMO) systems, the latter suggests a
promising architecture for amalgamation with RIS [6].

Thus far, there has been a considerable amount of research on RISs. However,
most of the existing works in the literature assume independent Rayleigh
fading such as [7, 3]. In contrast, it was shown in [8], it was shown that RIS
correlation should be taken into account in practice. On this ground, several
recent research efforts have considered RIS correlation in their analysis [5,
9]. Furthermore, aiming at realistic modeling, many works have taken into
account phase-shift errors coming from the finite precision of phase-shifts
configuration [10, 2].

Recent advancements on metasurfaces have brought to the forefront the concept
of simultaneous transmitting and reflecting RISs (STAR-RISs). In particular,
the STAR-RIS can provide a full space coverage by not only reflecting in the
half space, but also refracting to space behind the RIS [11, 12]. In
particular, in [11], a general hardware model and two-channel models for the
near-field region and the far-field region of STAR-RIS have been presented
while showing that their diversity gain and coverage are greater than
conventional RIS (i.e. reflective-only RIS) assisted systems. Moreover, in
[12], three operating protocols for STAR-RIS were suggested, which are known
as energy splitting (ES), mode switching (MS), and time switching (TS).
However, none of the existing works on STAR-RIS have considered correlated
fading, phase errors, and mMIMO systems.

Against the above background, we provide in closed form the only work
obtaining the coverage probability for STAR-RIS assisted mMIMO systems with
practical channel and phase models. In particular, contrary to existing works
on STAR-RIS such as[11, 12], we formulate a model that embodies correlated
fading and phase errors to identify the realistic prospects of STAR-RIS before
its final implementation. In this respect, we derive the coverage probability
for mMIMO systems for both ES and MS protocols. Especially, compared to [11],
which assumed a single-antenna transmitter, and required no special phase-
shifts optimization, we also consider a large number of antennas at the base
station (BS). For this reason, we follow the methodology in [2, 5] based on
statistical channel state information (CSI) to optimize the passive
(reflecting and refracting) beamforming matrix (PBM) of each user equipment
(UE) at every several coherence intervals, which brings lower overhead
compared to optimizations relying on instantaneous CSI. To the best of our
knowledge, this is the only work on STAR-RIS, where both beamforming matrices
for reflection and refraction are optimized based on statistical CSI, and
thereby reducing the amount of required overhead, which can be considered as
one of the main challenges for STAR-RIS.

Notation: Vectors and matrices are described by boldface lower and upper case
symbols, respectively. The notations $(\cdot)^{\scriptscriptstyle\mathsf{T}}$,
$(\cdot)^{\scriptscriptstyle\mathsf{H}}$, and $\tr\\!\left({\cdot}\right)$
denote the transpose, Hermitian transpose, and trace operators, respectively.
Moreover, the notations $\arg\left(\cdot\right)$,
$\mathbb{E}\left[\cdot\right]$, and $\mathrm{Var}(\cdot)$ denote the argument
function, the expectation, and variance operators, respectively. The notation
$\text{diag}\left(\right)\left({\mathbf{A}}\right)$ denotes a vector with
elements equal to the diagonal elements of ${\mathbf{A}}$, while
${\mathbf{b}}\sim{\cal C}{\cal N}{({\mathbf{0}},\mathbf{\Sigma})}$ denotes a
circularly symmetric complex Gaussian vector with zero mean and covariance
matrix $\mathbf{\Sigma}$.

##  II Signal and System Models

This section presents the signal and system models of the STAR-RIS assisted
system.

###  II-A Signal Model

Regarding the description of the signal model of the STAR-RIS, let $s_{n}$
describe the incident signal on element $n\in\mathcal{N}$, where
$\mathcal{N}=\\{1,\ldots,N\\}$ is the set of RIS elements. Two independent
coefficients, denoted as the transmission and the reflection coefficients,
configure the transmitted and reflected signals in respective modes.
Especially, the transmitted ($t$) and reflected ($r$) signals by the $n$th
element can be modelled as
$t_{n}=(\sqrt{\beta_{n}^{t}}e^{j\phi_{n}^{t}})s_{n}$ and
$r_{n}=(\sqrt{\beta_{n}^{r}}e^{j\phi_{n}^{r}})s_{n}$, respectively, where
${\beta_{n}^{k}}\in[0,1]$ and $\phi_{n}^{k}\in[0,2\pi)$ express the
independent amplitude and phase-shift response of the $n$th element, and
$k\in\\{t,r\\}$ corresponds to the UE found in the transmission ($t$) or
reflection ($r$) region [11].111In practical applications, the amplitude and
phase-shifts are correlated and result in a performance loss. The study of
this correlation is left for future work. Note that the choice of
$\phi_{n}^{t}$ and $\phi_{n}^{r}$ is independent from each other, but the
amplitude adjustments are correlated based on the law of energy conservation
as

$\displaystyle\beta_{n}^{t}+\beta_{n}^{r}=1,\forall n\in\mathcal{N}.$ (1)

Operation Protocols: Herein, we present briefly the main points of the ES/MS
protocols [12]. The study of the TS protocol is left for future work.

ES protocol: All RIS elements serve simultaneously $t$ and $r$ UEs, and the
corresponding PBM for $k\in\\{t,r\\}$ is expressed as
${\bm{\Phi}}_{k}^{\mathrm{ES}}=\text{diag}\left(\right)(\sqrt{\beta_{1}^{k}}e^{j\phi_{1}^{k}},\ldots,\sqrt{\beta_{N}^{k}}e^{j\phi_{N}^{k}})\in\mathbb{C}^{N\times
N}$, where $\beta_{n}^{t},\beta_{n}^{r}\in[0,1]$,
$\beta_{n}^{t}+\beta_{n}^{r}=1$, and
$\phi_{n}^{t},\phi_{n}^{r}\in[0,2\pi),\forall n\in\mathcal{N}$.

MS protocol: The elements are divided into two groups of $N_{t}$ and $N_{r}$
elements serving UE $t$ or $r$, respectively, i.e., $N_{t}+N_{r}=N$. In such
case, the PBM for $k\in\\{t,r\\}$ is
${\bm{\Phi}}_{k}^{\mathrm{MS}}=\text{diag}\left(\right)(\sqrt{\beta_{1}^{k}}e^{j\phi_{1}^{k}},\ldots,\sqrt{\beta_{N}^{t}}e^{j\phi_{N}^{k}})\in\mathbb{C}^{N\times
N}$, where $\beta_{n}^{t},\beta_{n}^{r}\in\\{0,1\\}$,
$\beta_{n}^{t}+\beta_{n}^{r}=1$, and
$\phi_{n}^{t},\phi_{n}^{r}\in[0,2\pi),\forall n\in\mathcal{N}$. Basically,
this protocol can be considered as a special case of ES by restricting the
amplitude coefficients for transmission and reflection to binary values. Thus,
the MS protocol cannot provide the full-dimension transmission and reflection
beamforming gain as ES, however, it requires less computational complexity in
terms of PBM design.

###  II-B System Model

We consider a STAR-RIS-assisted mMIMO communication system, where a BS with a
large number of antennas, denoted by $M$, communicates with $2$ single-antenna
UEs, as shown in Fig 1.222This model can be easily extended to the multi-user
scenario as in typical mMIMO, where all UEs are divided into multiple groups
of two UEs located at the opposite sides of RIS. Especially, the incident
signals on the RIS are divided simultaneously into transmitted and reflected
signals to create full-space coverage that are received by UE $t$ located
behind the RIS, and UE $r$ located in front of the RIS, i.e., the BS and UE
$r$ are found at the same side. Notably, we assume that the UEs operate at
orthogonal frequency bands of equal sizes similar to [11]. Also, we account
for no direct links between the BS and the users due to blockages, which is
one of the main challenging practical scenarios that justifies the use of a
RIS. The phase shifts of both transmitted and reflected signals are adjusted
by a controller that exchanges information with the BS through a backhaul
link. We rely on perfect CSI, which means that the results play the role of
upper bounds of practical scenarios with imperfect CSI. This is a common
assumption on works studying coverage to allow more direct mathematical
manipulations and to focus on the transmission and reflection beamforming gain
of STAR-RIS,

Figure 1: A STAR-RIS assisted mMIMO model.

We assume narrow-band quasi-static fading channels. Specifically,
${\mathbf{h}}_{k}\in\mathbb{C}^{N\times 1}$ describes the channel between the
RIS and UE $k\in\\{t,r\\}$ while
${\mathbf{G}}=[{\mathbf{g}}_{1}\ldots,{\mathbf{g}}_{N}]\in\mathbb{C}^{M\times
N}$, expresses the LoS channel between the BS and RIS with ${\mathbf{g}}_{i},$
$i=1\in\mathcal{N}$ corresponds to the $i$th column vector of ${\mathbf{G}}$.
In particular, the $(m,n)$th entry of ${\mathbf{G}}$ can be expressed as

$\displaystyle[{\mathbf{G}}]_{m,n}$ $\displaystyle=\sqrt{\beta_{g}}$.

Note that ${\mathbf{G}}$ could be obtained by several techniques, e.g., see
[13]. The parameters $\beta_{g}$ and $\lambda$ are the path-loss of the BS-RIS
link and the carrier wavelength. Also, $d_{\mathrm{BS}}$ and
$d_{\mathrm{RIS}}$ are the inter-antenna separation at the BS and inter-
element separation at the RIS, respectively [2]. In addition, $\theta_{1,n}$,
$\psi_{1,n}$ and $\theta_{2,m}$, $\psi_{2,m}$ describe the elevation and
azimuth LoS angles of departure at the BS and arrival at the RIS with respect
to RIS element $n$, respectively.

####  II-B1 RIS correlation

Instead of assuming independent Rayleigh fading as in existing STAR-RIS works
[11, 12], we consider spatial correlation, appearing in practice [8].333The
consideration of Rician fading, which includes an LoS component, is the topic
of future work. Hence, the channel vector between the RIS and UE $k$ with
$k\in\\{t,r\\}$ can be written as

$\displaystyle{\mathbf{h}}_{k}$
$\displaystyle=\sqrt{\beta_{k}}{\mathbf{R}}_{\mathrm{RIS}}^{1/2}{\mathbf{z}}_{k},$
(2)

where $\beta_{k}$ expresses the path-loss, and
${\mathbf{R}}_{\mathrm{RIS}}\in\mathbb{C}^{N\times N}$ with
$\tr\left({\mathbf{R}}_{\mathrm{RIS}}\right)=N$ expresses the deterministic
Hermitian-symmetric positive semi-definite correlation matrix at the RIS. In
particular, if $N=N_{\mathrm{H}}N_{\mathrm{V}}$ elements with $N_{\mathrm{H}}$
the number of elements per row and $N_{\mathrm{V}}$ the number of elements per
column, while $d_{\mathrm{H}}$ and $d_{\mathrm{V}}$ are horizontal and
vertical dimensions of each elements, the $(i,j)$th element of the RIS
correlation is given by

$\displaystyle
r_{ij}=d_{\mathrm{H}}d_{\mathrm{V}}\mathrm{sinc}\left(2\|\mathbf{u}_{i}-\mathbf{u}_{j}\|/\lambda\right),$
(3)

where
$\mathbf{u}_{\epsilon}=[0,\mod(\epsilon-1,N_{\mathrm{H}})d_{\mathrm{H}},\lfloor(\epsilon-1)/N_{\mathrm{H}}\rfloor
d_{\mathrm{V}}]^{\scriptscriptstyle\mathsf{T}}$, $\epsilon\in\\{i,j\\}$.

Note that the path-losses and the correlation matrices are assumed known after
being obtained by practical methods, e.g., see [14]. Moreover,
${\mathbf{z}}_{k}\sim\mathcal{CN}\left({\mathbf{0}},{\bm{\mathrm{I}}}_{N}\right)$
describes the independent and identically distributed (i.i.d.) fast-fading
vector in the $k\mathrm{th}$ link.

####  II-B2 Phase-shift errors

The configuration of the RIS elements with infinite precision is not possible
in practice. Hence, phase errors appear, which result from imperfections in
phase-estimation and/or phase quantization that cannot be avoided [10].
Specifically, we denote $\tilde{\phi}_{n}^{k}$, with $n\in\mathcal{N}$,
$k\in\\{t,r\\}$ the i.i.d. randomly distributed phase error in $[-\pi,\pi)$ of
the $n$th phase-shift aiming at the $k$th UE. Hence, the random phase error
matrix corresponding to the PBM ${\bm{\Phi}}_{k}$ is diagonal and described by
$\widetilde{{\bm{\Phi}}}_{k}=\text{diag}\left(\right)\left(e^{j\tilde{\phi}_{1}},\ldots,e^{j\tilde{\phi}_{N}}\right)\in\mathbb{C}^{N\times
N}$. The probability density function (PDF) of $\tilde{\phi}_{n}^{k}$ is
assumed symmetric with its mean direction equal to zero, i.e.,
$\arg\left(\mathbb{E}[\mathrm{e}^{j\tilde{\phi}_{n}^{k}}]\right)=0$ [10].

In general, the Uniform and the Von Mises distributions are the main PDFs that
can describe the phase noise affecting the RIS [10]. Notably, the Uniform PDF
of phase-noise has a characteristic function (CF) equal to zero, and thus it
cannot provide any knowledge regarding the phase-estimation accuracy. In
contrast, the Von Mises PDF has a zero-mean and concentration parameter
$\kappa_{\tilde{\phi}}$, which captures the accuracy of the estimation. In
particular, the CF of the Von Mises noise CF is
$m=\frac{\mathrm{I}_{1}\\!\left(\kappa_{\tilde{\phi}}\right)}{\mathrm{I}_{0}\\!\left(\kappa_{\tilde{\phi}}\right)}$,
where $\mathrm{I}_{v}\\!\left(\kappa_{\tilde{\phi}}\right)$ is the modified
Bessel function of the first kind and order $v$.

The received signal by UE $k\in\\{t,r\\}$ through the STAR-RIS assisted
network is described by

$\displaystyle
y_{k}=\sqrt{\rho_{\mathrm{dl}}}{\mathbf{h}}_{k}^{{\scriptscriptstyle\mathsf{H}}}\widetilde{{\bm{\Phi}}}_{k}^{{\scriptscriptstyle\mathsf{H}},{\scriptscriptstyle\mathcal{X}}}{\bm{\Phi}}_{k}^{{\scriptscriptstyle\mathsf{H}},{\scriptscriptstyle\mathcal{X}}}{\mathbf{G}}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{f}}_{k}q_{k}+z_{k},$
(4)

where $\rho_{\mathrm{dl}}\geq 0$ is the transmit power at BS to UE $k$, and
$q_{k}$ is the data symbol with $\mathbb{E}[|q_{k}|^{2}]=1$, while
${\mathbf{f}}_{k}\in\mathcal{C}^{M\times 1}$ is the linear precoding vector,
and $z_{k}\sim\mathcal{CN}(0,N_{0})$ is the additive white Gaussian noise with
zero mean and variance of $N_{0}$.

Notably, the cascaded channel vector including the phase-shift errors between
the large BS and UE $k\in\\{t,r\\}$ can be expressed as
$\bar{{\mathbf{h}}}_{k}={\mathbf{G}}{\bm{\Phi}}_{k}^{{\scriptscriptstyle\mathcal{X}}}\widetilde{{\bm{\Phi}}}_{k}^{{\scriptscriptstyle\mathcal{X}}}{\mathbf{h}}_{k}\in\mathbb{C}^{M\times
1}$, where $\mathcal{X}=\mathrm{ES},\mathrm{MS}$. It is distributed as
$\mathcal{CN}({\mathbf{0}},\bar{{\mathbf{R}}}_{k})$, where
$\bar{{\mathbf{R}}}_{k}=\beta_{k}{\mathbf{G}}{\bm{\Phi}}_{k}^{{\scriptscriptstyle\mathcal{X}}}\tilde{{\mathbf{R}}}_{\mathrm{RIS}}^{{\scriptscriptstyle\mathcal{X}}}{\bm{\Phi}}_{k}^{{\scriptscriptstyle\mathsf{H}},{\scriptscriptstyle\mathcal{X}}}{\mathbf{G}}^{{\scriptscriptstyle\mathsf{H}}}$.
Note that
$\tilde{{\mathbf{R}}}_{\mathrm{RIS}}^{{\scriptscriptstyle\mathcal{X}}}=\mathbb{E}[{\bm{\Phi}}_{k}^{{\scriptscriptstyle\mathcal{X}}}{\mathbf{R}}_{\mathrm{RIS}}{\bm{\Phi}}_{k}^{{\scriptscriptstyle\mathsf{H}},{\scriptscriptstyle\mathcal{X}}}]=m^{2}{\mathbf{R}}_{\mathrm{RIS}}+\left(1-m^{2}\right){\bm{\mathrm{I}}}_{N}$,
where $m$ denotes the CF of the phase error [2, Eq. 13]. In the case of the
Uniform distribution, i.e, $m=0$, ${\mathbf{R}}_{k}$ does not depend on the
phase-shifts and RIS cannot be optimized. A similar observation is met when
RIS optimization is based on statistical CSI [2]. Hence, apart from its
practical meaning, RIS correlation should be taken into account to exploit
RIS.

##  III Main Results

This section presents the signal-to-noise ratio (SNR), the coverage
probability, and its optimization for STAR-RIS assisted systems. Generally,
the coverage probability for $k$th UE, ${P}_{\mathrm{c}}^{k}$ with
$k\in\\{t,r\\}$, is the probability that the received SNR is greater than a
threshold $T$, i.e., it is described by
${P}_{\mathrm{c}}^{k}=\mathrm{Pr}\left(\gamma_{k}>T\right)$, where
$\mathrm{Pr}(\cdot)$ denotes probability.

###  III-A SNR

First, we take advantage of channel hardening, because UEs do not have any
knowledge of the instantaneous CSI in practice, but they are aware of their
statistics [15]. Note that channel hardening appears in mMIMO systems as the
number of BS antennas increases. Hence, by resorting to the application of the
use-and-then-forget bounding technique [15], we write (4) as

$\displaystyle y_{k}$
$\displaystyle=\sqrt{\rho_{\mathrm{dl}}}\big{(}\mathbb{E}[\bar{{\mathbf{h}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{f}}_{k}]q_{k}+\bar{{\mathbf{h}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{f}}_{k}q_{k}-\mathbb{E}[\bar{{\mathbf{h}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{f}}_{k}]q_{k}\big{)}+z_{k}.$
(5)

Next, by assuming that (5) represents a single-input single-output (SISO)
system, where the BS treats the unknown terms as uncorrelated additive noise,
the achievable SNR of the link between the BS and UE $k\in\\{t,r\\}$ is
written as

$\displaystyle\gamma_{k}=\frac{\big{|}\mathbb{E}[\bar{{\mathbf{h}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{f}}_{k}]\big{|}^{2}}{\mathbb{E}[\big{|}\bar{{\mathbf{h}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{f}}_{k}\big{|}^{2}]-\big{|}\mathbb{E}[\bar{{\mathbf{h}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{f}}_{k}]\big{|}^{2}+{\sigma_{0}}{}},$
(6)

where $\sigma_{0}=N_{0}/\rho_{\mathrm{dl}}$ [15]. Regarding the selection of
precoding, we apply the simple maximum ratio transmission (MRT) precoding
because it allows the derivation of closed-from expressions and extraction of
fundamental properties together with its optimality for the single UE case,
i.e.,
${\mathbf{f}}_{k}=\frac{\bar{{\mathbf{h}}}_{k}}{\sqrt{\mathbb{E}[\|\bar{{\mathbf{h}}}_{k}\|^{2}]}}$.
444Another common choice could be the more robust regularized zero-forcing
(RZF) precoder, but that could be applicable in a multi-user scenario with
interference, and it would result in intractable expressions. Hence, the
application of RZF could be the topic of future work.

###### Proposition 1

The downlink SNR of a STAR-RIS assisted mMIMO system with correlated Rayleigh
fading and phase-shift errors at UE $k\in\\{t,r\\}$ is given by

$\displaystyle\\!\\!\\!\gamma_{k}\\!=\\!\frac{\tr^{2}(\bar{{\mathbf{R}}}_{k})}{\tr(\bar{{\mathbf{R}}}_{k}^{2})+{\sigma_{0}}\tr(\bar{{\mathbf{R}}}_{k})}.$
(7)

###### Proof:

See Appendix A. ∎

###### Remark 1

As can be seen, the SNR in (7) depends only on the PBM and the statistical CSI
by means of the correlation matrix, the path-losses, and the phase-shift
errors. Also, there is a dependence on the location of UE $k$, i.e., if it is
found behind or in front of the RIS.

###  III-B Coverage Probability

Having obtained the SNR, the following proposition provides the coverage
probability.

###### Proposition 2

The coverage probability of a mMIMO STAR-RIS assisted system, accounting for
RIS correlation and phase-shift errors of UE $k\in\\{t,r\\}$, is tightly
approximated as

$\displaystyle{P}_{\mathrm{c}}^{k}$
$\displaystyle\approx\sum^{L}_{n=1}\binom{L}{n}(-1)^{n+1}e^{-\eta\frac{T}{\gamma_{k}}},$
(8)

where $\eta=L\left(L!\right)^{-\frac{1}{L}}$ with $L$ being an approximation
parameter.

###### Proof:

See Appendix B. ∎

The coverage depends on the threshold, the number $L$ defining the tightness
of the approximation, and of course, the downlink SNR with its involved
parameters.

Conventional RIS: We assume a smart surface that consists of transmitting-only
or reflective-only elements, each with $N/2$ elements with $N$ even for
simplicity. Also, $\beta_{n}^{t}=1$ or $\beta_{n}^{r}=1,\forall n$,
respectively. In Sec. IV, we compare the performance of the STAR-RIS with the
conventional RIS and illustrate the advantage of adopting STAR-RIS in wireless
systems.

###  III-C Passive beamforming matrix optimization

A STAR-RIS assisted system serves two regions simultaneously, where UEs $t$
and $r$ are found. We focus on the optimization of the total coverage
probability given by ${P}_{\mathrm{c}}^{t}+{P}_{\mathrm{c}}^{r}$. Given the
difficulty in simultaneous optimization of both ${P}_{\mathrm{c}}^{t}$ and
${P}_{\mathrm{c}}^{r}$, we rely on alternating optimization by optimizing
first each of them with respect to its PBM while fixing the other in an
iterative manner until reaching the convergence, i.e., we perform optimization
of ${P}_{\mathrm{c}}^{k}$ with respect to
${\bm{\Phi}}_{k}^{{\scriptscriptstyle\mathcal{X}}}$. Each optimization is
achieved in terms of the projected gradient ascent until converging to a
stationary point.

Hence, based on the common assumption of infinite-resolution phase-shifters,
the maximization algorithm for ${P}_{\mathrm{c}}^{k}$, $k\in\\{t,r\\}$ UE with
respect to the PBM is formulated as

$\displaystyle\begin{split}&\\!\\!\\!(\mathcal{P}1)~{}\max_{{\bm{\Phi}}_{k}^{{\scriptscriptstyle\mathcal{X}}}}~{}~{}{P}_{\mathrm{c}}^{k}\\\
&~{}~{}~{}~{}~{}\mathrm{s.t}~{}~{}~{}~{}\beta_{n}^{t}+\beta_{n}^{r}=1,\forall
n\in\mathcal{N}\\\
&~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\beta_{n}^{t},\beta_{n}^{r}\in[0,1],\phi_{n}^{t},\phi_{n}^{r}\in[0,2\pi),~{}~{}\forall
n\in\mathcal{N}.\end{split}$ (9)

By replacing $[0,1]$ with $\\{0,1\\}$, $(\mathcal{P}1)$ can be applied to the
MS scheme.

The optimization problem $(\mathcal{P}1)$ is non-convex regarding
$\phi_{n}^{k}$. The dependence on $\phi_{n}^{k}$ is found on the covariance
matrices $\bar{{\mathbf{R}}}_{k}$. Given that it is a constrained maximization
problem, we resort to the projected gradient ascent (PGA) until convergence to
a stationary point to provide its solution [2]. The transmit power constraint
guarantees its convergence. It is worthwhile to mention that the complexity of
${P}_{\mathrm{c}}^{k}$ is $\mathcal{O}\left(G(MN^{2}+M^{3}+L)\right)$. As can
be seen, it is a function of the fundamental system parameters $M$, $N$, and
$L$ with the number of BS antennas having the higher impact. Since it has to
be performed twice for each $k\in\\{t,r\\}$, the complexity is double, while,
in the case of the conventional RIS scenario, the complexity is half of the
STAR-RIS setting.

Based on PGA, let
${\mathbf{s}}_{k,l}=[\phi_{1}^{k,l},\ldots,\phi_{N}^{k,l}]^{{\scriptscriptstyle\mathsf{T}}}$
denote the vector including the phases at step $l$, while the next iteration
point provides the increase of ${P}_{\mathrm{c}}^{k}$ upon its convergence by
projecting the solution onto the closest feasible point as specified by
$\min_{\phi_{n}^{k}\in[0,2\pi),}\|{\mathbf{s}}_{k}-\tilde{{\mathbf{s}}}_{k}\|^{2}$.
The next iteration is described by

$\displaystyle\tilde{{\mathbf{s}}}_{k,l+1}$
$\displaystyle={\mathbf{s}}_{k,l}+\mu{\mathbf{v}}_{k,l},$ (10)
$\displaystyle{\mathbf{s}}_{k,l+1}^{w}$
$\displaystyle=\beta_{n}^{t}\exp\left(j\arg\left(\tilde{{\mathbf{s}}}_{k,l+1}\right)\right),$
(11)

where the parameter $\mu$ is the step size obtained at each iteration through
the backtracking line search [16]. Moreover, ${\mathbf{v}}_{k,l}$ is the
ascent direction at step $l$, which means
${\mathbf{v}}_{k,l}=\partialderivative{{P}_{\mathrm{c}}^{k}}{{\mathbf{s}}_{k,l}^{*}}$.
The following lemma provides this derivative.

###### Lemma 1

The derivative of the coverage probability ${P}_{\mathrm{c}}^{k}$ for UE
$k\in\\{t,r\\}$ with respect to ${\mathbf{s}}_{k,l}^{*}$ is given by (15) at
the top of the next page, where $\gamma_{k}$ is given by (7), and

$\displaystyle S_{k}$ $\displaystyle=\tr^{2}(\bar{{\mathbf{R}}}_{k}),$ (12)
$\displaystyle I_{k}$
$\displaystyle=\tr(\bar{{\mathbf{R}}}_{k}^{2})+{\sigma_{0}}\tr(\bar{{\mathbf{R}}}_{k}),$
(13) $\displaystyle\partialderivative{S_{k}}{{\mathbf{s}}_{k,l}^{*}}$
$\displaystyle=2\beta_{k}\tr(\bar{{\mathbf{R}}}_{k})\text{diag}\left(\right)({\mathbf{G}}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{G}}{\bm{\Phi}}_{k}^{{\scriptscriptstyle\mathcal{X}}}\tilde{{\mathbf{R}}}_{\mathrm{RIS}}^{{\scriptscriptstyle\mathcal{X}}})$
$\displaystyle\partialderivative{I_{k}}{{\mathbf{s}}_{k,l}^{*}}$
$\displaystyle=\beta_{k}\tr({\mathbf{G}}^{{\scriptscriptstyle\mathsf{H}}}\bar{{\mathbf{R}}}_{k}{\mathbf{G}}{\bm{\Phi}}_{k}^{{\scriptscriptstyle\mathcal{X}}}\tilde{{\mathbf{R}}}_{\mathrm{RIS}}^{{\scriptscriptstyle\mathcal{X}}})$
$\displaystyle+{\sigma_{0}}\beta_{k}\tr({\mathbf{G}}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{G}}{\bm{\Phi}}_{k}^{{\scriptscriptstyle\mathcal{X}}}\tilde{{\mathbf{R}}}_{\mathrm{RIS}}^{{\scriptscriptstyle\mathcal{X}}}).$
(14)

$\displaystyle\\!\partialderivative{{P}_{\mathrm{c}}^{k}}{{\mathbf{s}}_{k,l}^{*}}\\!\\!$
$\displaystyle=\\!\frac{\partialderivative{S_{k}}{{\mathbf{s}}_{k,l}^{*}}I_{k}-S_{k}\partialderivative{I_{k}}{{\mathbf{s}}_{k,l}^{*}}}{\gamma_{k}^{2}I_{k}^{2}}\sum^{L}_{n=1}\\!\binom{L}{n}\\!\left(-1\right)^{n+1}n\eta
T\mathrm{e}^{-n\eta\frac{T}{\gamma_{k}}}.$ (15)

###### Proof:

See Appendix C. ∎

##  IV Numerical Results

In this section, we provide the numerical results for the coverage of the
STAR-RIS assisted mMIMO system. We consider a uniform planar array (UPA) of
$N=60$ elements for the RIS, while a uniform linear array (ULA) with $M=40$
antennas is assumed for the BS that serves UEs $t$ and $r$. Unless otherwise
stated, we consider the following values. The path-losses are generated
according to the NLOS version of the 3GPP Urban Micro (UMi) scenario from
TR36.814 for a carrier frequency of $2.5$ GHz, and noise level $-80$ dBm [17].
Specifically, we have $\beta_{g}=C_{g}d_{g}^{-\nu_{g}},$ and
$\beta_{k}=C_{k}d_{k}^{-\nu_{k}},~{}k\in\\{t,r\\}$ with $C_{g}=26$ dB,
$C_{t}=26$ dB, $C_{r}=28$ dB, $\nu_{g}=2.1$, $\nu_{t}=2.5$, and $\nu_{r}=2.2$.
The variables $d_{g}$ and $d_{k}$ denote the distances between the BS and RIS
and between the RIS and UE $k\in\\{t,r\\}$, respectively. We assume that the
correlation matrix for the RIS is given by [8], where the horizontal and
vertical dimension of each element is $d_{\mathrm{H}}=\lambda/8$ and
$d_{\mathrm{V}}=\lambda/8$. Also, we assume that $\beta_{n}^{t}=0.4$,
$\beta_{n}^{r}=0.6,\forall n$, $\rho_{\mathrm{dl}}=6~{}\mathrm{dB}$, and
$N_{0}=-174+10\log_{10}B_{\mathrm{c}}$ with
$B_{\mathrm{c}}=200~{}\mathrm{KHz}$. Monte-Carlo (MC) simulations are carried
out to verify our analysis.

In Fig. 2, we depict the coverage probability versus the target rate for
varying number of RIS elements $N$ while studying other key properties. We
observe that a higher number of $N$, increases ${P}_{\mathrm{c}}^{k}$ as
expected. Moreover, we study the impact of phase errors for $N=25$, and we see
that, in the case of Uniform PDF, where no knowledge on the errors is known,
the RIS cannot be optimized, since the covariance matrix does not depend on
the phase shifts. Thus, the coverage is lower. However, if the phase errors
follow the Von Mises PDF ($m=0.5$), ${P}_{\mathrm{c}}^{k}$ is higher thanks to
the possibility of performing phase-optimization under such phase-noise
distribution. Also, we make a similar observation regarding the correlation
for $N=100$, i.e., we notice that a RIS correlation allows the PBM
optimization, and ${P}_{\mathrm{c}}^{k}$ increases in the case of statistical
CSI modeling. In addition, the ES scheme achieves higher coverage than the MS
scheme because the latter is a special case of ES. Furthermore, we have added
the impact of the phase error for both STAR-RIS and conventional RIS for
$m=0.5$ and $m=0.7$, and we observe that a larger $m$ increases the estimation
accuracy, i.e., the coverage is higher in both cases, while the coverage for
the STR-RIS is generally higher compared to the conventional RIS as expected.

Fig. 3 illustrates the coverage probability versus the target rate for varying
number of BS antennas $M$, i.e., $M=30,60,90$. It is shown that
${P}_{\mathrm{c}}^{k}$ increases with an increasing number of BS antennas
because of a higher beamforming gain. Also, we notice that independent
Rayleigh fading conditions result in the worst coverage because no RIS
exploitation can be achieved in the case of statistical CSI. In this
direction, we observe that correlation should be taken into account to benefit
from the RIS, while a higher correlation leads to lower coverage. For the sake
of comparison, in the same figure, we have included the coverage in the case
of a conventional RIS as presented in Sec. III-B with $N=30$ elements.
Notably, the STAR-RIS system outperforms the conventional RIS-assisted system
because more degrees of freedom for transmission and reflection can be
harnessed.

Figure 2: Coverage probability of a STAR-RIS assisted mMIMO system versus the
target rate $T$ for varying $N$ ($M=40$,
$\mathrm{d_{\mathrm{H}}=d_{\mathrm{V}}=\lambda/8}$, $m=0.5$).  Figure 3:
Coverage probability of a STAR-RIS assisted mMIMO system versus the target
rate $T$ for varying $M$ ($N=60$,
$\mathrm{d_{\mathrm{H}}=d_{\mathrm{V}}=\lambda/8}$, $m=0.5$).

##  V Conclusion

In this paper, we obtained the coverage probability of a STAR-RIS assisted
mMIMO system under the realistic assumptions of correlated Rayleigh fading and
phase-shift errors, while previous works on STAR-RIS did not account for these
inevitable effects. Especially, we derived the coverage probability for both
$t$ and $r$ UEs in terms of large-scale statistics that change at every
several coherence intervals, and thus, reduce the increased overhead of STAR-
RIS. Among others, we depicted how the fundamental parameters affect the
coverage, and the outperformance of STAR-RIS compared to reflective-only RIS.
Future works could consider Ricean channels or even millimetre-wave
transmission.

##  Appendix A Proof of Proposition 1

Based on the property
${\mathbf{x}}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{y}}=\tr({\mathbf{y}}{\mathbf{x}}^{{\scriptscriptstyle\mathsf{H}}})$,
the term in the numerator of (7) is written as

$\displaystyle\\!\\!\\!|\mathbb{E}[\bar{{\mathbf{h}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{f}}_{k}]|^{2}$
$\displaystyle\\!=\\!\Big{|}\\!\tr(\\!\frac{\mathbb{E}[\bar{{\mathbf{h}}}_{k}\bar{{\mathbf{h}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}]}{\mathbb{E}[\|\bar{{\mathbf{h}}}_{k}\|^{2}]}\\!)\Big{|}^{2}\\!$
(16) $\displaystyle\\!=\\!\tr(\bar{{\mathbf{R}}}_{k}),$ (17)

where, in (16), we have substituted the precoder ${\mathbf{f}}_{k}$. Note that
$\bar{{\mathbf{R}}}_{k}=\beta_{k}{\mathbf{G}}{\bm{\Phi}}_{k}^{{\scriptscriptstyle\mathcal{X}}}\tilde{{\mathbf{R}}}_{\mathrm{RIS}}^{{\scriptscriptstyle\mathcal{X}}}{\bm{\Phi}}_{k}^{{\scriptscriptstyle\mathsf{H}},{\scriptscriptstyle\mathcal{X}}}{\mathbf{G}}^{{\scriptscriptstyle\mathsf{H}}}$.

Regarding the first term in the denominator, we have

$\displaystyle\\!\\!\\!\\!\mathbb{E}[\big{|}\bar{{\mathbf{h}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{f}}_{k}\big{|}^{2}]$
$\displaystyle\\!=\\!\tr(\bar{{\mathbf{R}}}_{k})+\frac{\tr(\bar{{\mathbf{R}}}_{k}^{2})}{\tr(\bar{{\mathbf{R}}}_{k})},$
(18)

where (18) is obtained by [15, Lem. B.14]. The proof is concluded by
appropriate substitutions into (6).

##  Appendix B Proof of Proposition 2

The proof follows similar lines with the technique in [18]. Specifically,
based on the definition of the coverage probability of UE $k\in\\{t,r\\}$, we
have

$\displaystyle{P}_{\mathrm{c}}^{k}$
$\displaystyle\approx\mathbb{P}(\tilde{g}>\frac{T}{\gamma_{k}})$ (19)
$\displaystyle\approx 1-(1-e^{-\eta\frac{T}{\gamma_{k}}})^{L}$ (20)
$\displaystyle=\sum_{n=1}^{L}\binom{L}{n}(-1)^{n+1}e^{-\eta\frac{T}{\gamma_{k}}}.$
(21)

In (19), we have used the approximation [18], where the constant number $1$
can be replaced by a dummy Gamma variable $\tilde{g}$, which has a unit mean
and shape parameter $L$, i.e., its pdf is
$\frac{L^{-L}}{\Gamma(L)}\tilde{g}^{L-1}e^{-{L\tilde{g}}{}}$. The tightness of
this approximation increases with increasing $L$. In (20), we have applied
Alzer’s inequality, where $\eta=L\left(L!\right)^{-\frac{1}{L}}$ [18], while
in (21), we have applied the Binomial theorem.

##  Appendix C Proof of Lemma 1

First, we apply the chain rule, which gives

$\displaystyle\partialderivative{{P}_{\mathrm{c}}^{k}}{{\mathbf{s}}_{k,l}^{*}}$
$\displaystyle=\partialderivative{{P}_{\mathrm{c}}^{k}}{\gamma_{k}}\partialderivative{\gamma_{k}}{{\mathbf{s}}_{k,l}^{*}}.$
(22)

Regarding the first term in (22), we have

$\displaystyle\partialderivative{{P}_{\mathrm{c}}^{k}}{\gamma_{k}}$
$\displaystyle=\sum^{L}_{n=1}\\!\binom{L}{n}\\!\frac{\left(-1\right)^{n+1}n\eta
T}{\gamma_{k}^{2}}\mathrm{e}^{-n\eta\frac{T}{\gamma_{k}}}.$ (23)

For the derivation of the second derivative in (22), i.e.,
$\partialderivative{\gamma_{k}}{{\mathbf{s}}_{k,l}^{*}}$, we define
$\gamma_{k}=\frac{S_{k}}{I_{k}}$, where
$S_{k}=\tr^{2}(\bar{{\mathbf{R}}}_{k})$ and
$I_{k}=\tr(\bar{{\mathbf{R}}}_{k}^{2})+{\sigma_{0}}\tr(\bar{{\mathbf{R}}}_{k})$.
Hence, the derivative can be written as

$\displaystyle\partialderivative{\gamma_{k}}{{\mathbf{s}}_{k,l}^{*}}=\frac{\partialderivative{S_{k}}{{\mathbf{s}}_{k,l}^{*}}I_{k}-S_{k}\partialderivative{I_{k}}{{\mathbf{s}}_{k,l}^{*}}}{I_{k}^{2}}.$
(24)

The computation of $\partialderivative{S_{k}}{{\mathbf{s}}_{k,l}^{*}}$
follows. In particular, we have

$\displaystyle\partialderivative{S_{k}}{{\mathbf{s}}_{k,l}^{*}}=2\tr(\bar{{\mathbf{R}}}_{k})\tr(\bar{{\mathbf{R}}}_{k}^{\prime}),$
(25)

where $(\cdot)^{\prime}$ denotes the partial derivative with respect to
${\mathbf{s}}_{k,l}^{*}$. Based on the property
$\tr\left({\mathbf{A}}\text{diag}\left(\right)({\mathbf{s}}_{k,l}^{*})\right)=\left(\text{diag}\left(\right)(A)\right)^{{\scriptscriptstyle\mathsf{T}}}{\mathbf{s}}_{k,l}^{*}$,
where ${\mathbf{A}}$ is a matrix independent of ${\mathbf{s}}_{k,l}^{*}$ we
can easily show Regarding $\partialderivative{I_{k}}{{\mathbf{s}}_{k,l}^{*}}$,
we obtain

$\displaystyle\\!\\!\\!\partialderivative{I_{k}}{{\mathbf{s}}_{k,l}^{*}}$
$\displaystyle\\!=\\!2\tr(\bar{{\mathbf{R}}}_{k}^{\prime}\bar{{\mathbf{R}}}_{k})\\!+\\!{\sigma_{0}}{}\\!\tr(\bar{{\mathbf{R}}}_{k}^{\prime})\\!.\\!$
(26)

The derivative in (14), is obtained by (26) after exploiting the property used
in (LABEL:tracederiv). Substitution of the derivatives in (24) concludes the
proof.

## References

  * [1] M. Di Renzo _et al._ , “Smart radio environments empowered by reconfigurable intelligent surfaces: How it works, state of research, and the road ahead,” vol. 38, no. 11, pp. 2450–2525.
  * [2] A. Papazafeiropoulos _et al._ , “Intelligent reflecting surface-assisted MU-MISO systems with imperfect hardware: Channel estimation, beamforming design,” _IEEE Trans. Wireless Commun._ , 2021.
  * [3] C. Huang _et al._ , “Reconfigurable intelligent surfaces for energy efficiency in wireless communication,” _IEEE Transa. Wireless Commun._ , vol. 18, no. 8, pp. 4157–4170, 2019.
  * [4] C. Guo _et al._ , “Outage probability analysis and minimization in intelligent reflecting surface-assisted MISO systems,” _IEEE Commun. Lett._ , vol. 24, no. 7, pp. 1563–1567, 2020.
  * [5] A. Papazafeiropoulos _et al._ , “Coverage probability of distributed IRS systems under spatially correlated channels,” _IEEE Wireless Commun. Lett._ , vol. 10, no. 8, pp. 1722–1726, 2021.
  * [6] K. Zhi _et al._ , “Statistical CSI-based design for reconfigurable intelligent surface-aided massive MIMO systems with direct links,” vol. 10, no. 5, pp. 1128–1132.
  * [7] Q. Wu and R. Zhang, “Intelligent reflecting surface enhanced wireless network via joint active and passive beamforming,” _IEEE Trans. Wireless Commun._ , vol. 18, no. 11, pp. 5394–5409, 2019.
  * [8] E. Björnson and L. Sanguinetti, “Rayleigh fading modeling and channel hardening for reconfigurable intelligent surfaces,” _IEEE Wireless Commun. Lett._ , vol. 10, no. 4, pp. 830–834, 2021.
  * [9] A. Papazafeiropoulos, “Ergodic capacity of IRS-assisted MIMO systems with correlation and practical phase-shift modeling,” _IEEE Wireless Commun. Lett._ , vol. 11, no. 2, pp. 421–425, 2022.
  * [10] M.-A. Badiu and J. P. Coon, “Communication through a large reflecting surface with phase errors,” vol. 9, no. 2, pp. 184–188.
  * [11] J. Xu _et al._ , “STAR-RISs: Simultaneous transmitting and reflecting reconfigurable intelligent surfaces,” _IEEE Commun. Lett._ , vol. 25, no. 9, pp. 3134–3138, 2021.
  * [12] X. Mu _et al._ , “Simultaneously transmitting and reflecting (STAR) RIS aided wireless communications,” _IEEE Trans. Wireless Commun._ , 2021.
  * [13] F. Bohagen, P. Orten, and G. E. Oien, “Design of optimal high-rank line-of-sight MIMO channels,” _IEEE Trans. Wireless Commun._ , vol. 6, no. 4, pp. 1420–1425, 2007.
  * [14] D. Neumann, M. Joham, and W. Utschick, “Covariance matrix estimation in massive MIMO,” _IEEE Signal Process. Lett._ , vol. 25, no. 6, pp. 863–867, 2018.
  * [15] E. Björnson _et al._ , “Massive MIMO networks: Spectral, energy, and hardware efficiency,” _Foundations and Trends® in Signal Processing_ , vol. 11, no. 3-4, pp. 154–655, 2017.
  * [16] S. Boyd, S. P. Boyd, and L. Vandenberghe, _Convex optimization_. Cambridge university press, 2004.
  * [17] 3GPP, “Further advancements for E-UTRA physical layer aspects (Release 9),” 3GPP TS 36.814, Tech. Rep., 2017.
  * [18] T. Bai and R. W. Heath, “Coverage and rate analysis for millimeter-wave cellular networks,” _IEEE Trans. on Wireless Commun._ , vol. 14, no. 2, pp. 1100–1114, 2015.

solutionfile


f_{G^{*}_{i}-\xi_{i}^{*},I_{d}}(Z_{i}^{\prime})}{f_{G^{*}_{i}-\xi_{i}^{*},I_{d}}(Z_{i}^{\prime})\lor\rho}\right)\right\|_{2}^{2}}$

Saha and Guntuboyina (2020b, Proof of Lemma F.3) then gives

$\displaystyle\|\hat{\theta}^{*}_{\rho}(Z)-\hat{\theta}^{*}_{\rho}(Z^{\prime})\|_{F}$
$\displaystyle\leq\bar{k}^{1/2}L^{2}(\rho)\|Z-Z^{\prime}\|_{F},$

where as before $L(\rho)\coloneqq\sqrt{-\log((2\pi)^{d}\rho^{2})}$. The same
argument applies to
$\|\hat{\theta}_{\rho}^{(j)}(Z)-\hat{\theta}_{\rho}^{(j)}(Z^{\prime})\|_{F}$.
Hence $F$ is $2\bar{k}L^{2}(\rho)$-Lipschitz. By concentration of Lipschitz
functions of Gaussians and a union bound

$\mathbb{P}\left(\zeta_{3}^{2}\geq x\right)\leq
N\exp\left(-\frac{x^{2}}{8\bar{k}L^{4}(\rho)}\right).$

Integrating the tail gives

$\displaystyle\mathbb{E}\zeta_{3}^{2}\leq 8\bar{k}L^{4}(\rho)\log(eN).$ (46)

####  C.3.4 Bounding $\mathbb{E}\zeta_{4}^{2}$

Again by the scaled Tweedie’s formula (37),

$\displaystyle\mathbb{E}\|\hat{\theta}_{\rho}^{(j)}-\hat{\theta}^{*}_{\rho}\|_{F}$
$\displaystyle\leq\sqrt{\mathbb{E}\|\hat{\theta}_{\rho}^{(j)}-\hat{\theta}^{*}_{\rho}\|_{F}^{2}}\leq\sqrt{\bar{k}\sum_{i=1}^{n}\mathbb{E}_{X_{i}\sim
f_{G^{*}_{i},I_{d}}}\left\|\frac{\nabla
f_{G^{(j)}_{i},I_{d}}(X_{i})}{f_{G^{(j)}_{i},I_{d}}(X_{i})\lor\rho}-\frac{\nabla
f_{G^{*}_{i},I_{d}}(X_{i})}{f_{G^{*}_{i},I_{d}}(X_{i})\lor\rho}\right\|_{2}^{2}}$

Saha and Guntuboyina (2020a, Lemma E.1) bounds the above expectation, yielding

$\displaystyle\left(\mathbb{E}\|\hat{\theta}_{\rho}^{(j)}-\hat{\theta}^{*}_{\rho}\|_{F}\right)^{2}$
$\displaystyle\leq
C_{d}\bar{k}\sum_{i=1}^{n}\max\left\\{\left(\frac{L^{2}(\rho)}{2}\right)^{3},\left|\log
h\left(f_{G^{*}_{i},I_{d}},f_{G^{(j)}_{i},I_{d}}\right)\right|\right\\}h^{2}\left(f_{G^{*}_{i},I_{d}},f_{G^{(j)}_{i},I_{d}}\right).$
(47)

By a change of variables,

$h^{2}\left(f_{G^{*}_{i},I_{d}},f_{G^{(j)}_{i},I_{d}}\right)=h^{2}\left(f_{G^{*},\Sigma_{i}},f_{G^{(j)},\Sigma_{i}}\right).$

Using the shorthand
$h^{2}_{i}=h^{2}\left(f_{G^{*},\Sigma_{i}},f_{G^{(j)},\Sigma_{i}}\right)$ and
using $\rho=(2\pi)^{-d/2}/n$,

$\displaystyle\left(\mathbb{E}\|\hat{\theta}_{\rho}^{(j)}-\hat{\theta}^{*}_{\rho}\|_{F}\right)^{2}$
$\displaystyle\leq C_{d}\bar{k}\sum_{i=1}^{n}\max\left\\{\left(\log
n\right)^{3},-\log h_{i}\right\\}h_{i}^{2}$ (48)
$\displaystyle=C_{d}\bar{k}\left(\sum_{i:\left(\log n\right)^{3}\geq-\log
h_{i}}\left(\log n\right)^{3}h_{i}^{2}+\sum_{i:\left(\log n\right)^{3}<-\log
h_{i}}-(\log h_{i})h_{i}^{2}\right)$ $\displaystyle\leq
C_{d}\bar{k}\left(n\left(\log n\right)^{3}\delta+\sum_{i:\left(\log
n\right)^{3}<\log h_{i}^{-1}}(\log h_{i}^{-1})h_{i}^{2}\right),$

where in the last step we used
$\frac{1}{n}\sum_{i=1}^{n}h^{2}_{i}=\bar{h}^{2}(f_{G^{*},\bullet},f_{G^{(j)},\bullet})\leq\delta$.
To bound the second term, note for $n\geq 6$, $(\log n)^{3}\geq 3\log n$,
implying $h_{i}\leq n^{-3}$ for all $i$ such that $\left(\log
n\right)^{3}<\log h_{i}^{-1}$. Since $h_{i}\log h_{i}^{-1}\leq e^{-1}$ for all
$h_{i}\in[0,1]$,

$\sum_{i:\left(\log n\right)^{3}<\log h_{i}^{-1}}(\log
h_{i}^{-1})h_{i}^{2}\leq\sum_{i:\left(\log n\right)^{3}<\log
h_{i}^{-1}}\frac{1}{en^{3}}\leq\frac{1}{en^{2}}.$

The first term dominates, so

$\displaystyle\mathbb{E}\zeta_{4}^{2}\leq C_{d}\bar{k}n(\log n)^{3}\delta.$
(49)

####  C.3.5 Bounding the metric entropy $\log N$

We will actually bound the larger covering number $\log
N(\eta^{*},{\mathcal{P}}(\mathbb{R}^{d}),m^{S})$ of the space of all
probability measures ${\mathcal{P}}(\mathbb{R}^{d})$ in the metric $m^{S}$.
For any measure $G$ we let $G_{i}$ denote the measure scaled by
$\Sigma_{i}^{-1/2}$ as in the scaled Tweedie formula. For
$G,H\in{\mathcal{P}}(\mathbb{R}^{d})$,

$\displaystyle m^{S}(G,H)$
$\displaystyle\coloneqq\max_{i\in[n]}\sup_{x:\mathfrak{d}_{S}(x)\leq
M}\left\|\frac{\Sigma_{i}\nabla
f_{G,\Sigma_{i}}(x)}{f_{G,\Sigma_{i}}(x)\lor(\rho/\sqrt{|\Sigma_{i}|})}-\frac{\Sigma_{i}\nabla
f_{H,\Sigma_{i}}(x)}{f_{H,\Sigma_{i}}(x)\lor(\rho/\sqrt{|\Sigma_{i}|})}\right\|_{2}$
$\displaystyle\leq\max_{i\in[n]}\sup_{x:\mathfrak{d}_{S}(x)\leq
M}\left\|\frac{\Sigma_{i}\nabla
f_{G,\Sigma_{i}}(x)}{f_{G,\Sigma_{i}}(x)\lor(\rho/\sqrt{|\Sigma_{i}|})}-\frac{\Sigma_{i}\nabla
f_{G,\Sigma_{i}}(x)}{f_{H,\Sigma_{i}}(x)\lor(\rho/\sqrt{|\Sigma_{i}|})}\right\|_{2}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
+\max_{i\in[n]}\sup_{x:\mathfrak{d}_{S}(x)\leq M}\left\|\frac{\Sigma_{i}\nabla
f_{G,\Sigma_{i}}(x)}{f_{H,\Sigma_{i}}(x)\lor(\rho/\sqrt{|\Sigma_{i}|})}-\frac{\Sigma_{i}\nabla
f_{H,\Sigma_{i}}(x)}{f_{H,\Sigma_{i}}(x)\lor(\rho/\sqrt{|\Sigma_{i}|})}\right\|_{2}$
$\displaystyle\leq\max_{i\in[n]}\sup_{x:\mathfrak{d}_{S}(x)\leq
M}\left\|\frac{\Sigma_{i}\nabla
f_{G,\Sigma_{i}}(x)}{f_{G,\Sigma_{i}}(x)\lor(\rho/\sqrt{|\Sigma_{i}|})}\right\|_{2}\frac{\left|f_{G,\Sigma_{i}}(x)\lor(\rho/\sqrt{|\Sigma_{i}|})-f_{H,\Sigma_{i}}(x)\lor(\rho/\sqrt{|\Sigma_{i}|})\right|}{f_{H,\Sigma_{i}}(x)\lor(\rho/\sqrt{|\Sigma_{i}|})}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
+\max_{i\in[n]}\sup_{x:\mathfrak{d}_{S}(x)\leq
M}\left\|\frac{\Sigma_{i}(\nabla f_{G,\Sigma_{i}}(x)-\nabla
f_{H,\Sigma_{i}}(x))}{f_{H,\Sigma_{i}}(x)\lor(\rho/\sqrt{|\Sigma_{i}|})}\right\|_{2}$

For the first term, by Saha and Guntuboyina (2020b, Lemma F.1),

$\left\|\frac{\Sigma_{i}\nabla
f_{G,\Sigma_{i}}(x)}{f_{G,\Sigma_{i}}(x)\lor(\rho/\sqrt{|\Sigma_{i}|})}\right\|_{2}=\left\|\Sigma_{i}^{1/2}\frac{\nabla
f_{G^{*}_{i},I_{d}}(\Sigma_{i}^{-1/2}x)}{f_{G^{*}_{i},I_{d}}(\Sigma_{i}^{-1/2}y)\lor\rho}\right\|_{2}\leq\bar{k}^{1/2}L(\rho).$

Replacing $f\lor(\rho/\sqrt{|\Sigma_{i}|})$ with $\rho/\sqrt{|\Sigma_{i}|}$ in
the denominator can only make the denominator smaller, so

$\displaystyle m^{S}(G,H)$
$\displaystyle\leq\bar{k}^{1/2}\rho^{-1}L(\rho)\max_{i\in[n]}\sup_{x:\mathfrak{d}_{S}(x)\leq
M}\sqrt{|\Sigma_{i}|}\left|f_{G,\Sigma_{i}}(x)-f_{H,\Sigma_{i}}(x)\right|$
(50) $\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
+\rho^{-1}\max_{i\in[n]}\sup_{x:\mathfrak{d}_{S}(x)\leq
M}\sqrt{|\Sigma_{i}|}\left\|\Sigma_{i}(\nabla f_{G,\Sigma_{i}}(x)-\nabla
f_{H,\Sigma_{i}}(x))\right\|_{2}$ (51)
$\displaystyle\leq\bar{k}^{d/2+1/2}\rho^{-1}L(\rho)\|f_{G,\bullet}-f_{H,\bullet}\|_{\infty,S^{M}}+\bar{k}^{d/2+1}\rho^{-1}\|f_{G,\bullet}-f_{H,\bullet}\|_{\nabla,S^{M}}$
(52)

In particular, letting

$\eta^{*}=\left(\bar{k}^{d/2+1/2}L(\rho)+\bar{k}^{d/2+1}\right)\frac{\eta}{\rho}$

we have

$\log N(\eta^{*},{\mathcal{P}}(\mathbb{R}^{d}),m^{S})\leq\log
N(\eta/2,\mathbb{F},\|\cdot\|_{\infty,S^{M}})+\log
N(\eta/2,\mathbb{F},\|\cdot\|_{\nabla,S^{M}}).$

We already have a bound on $\log
N(\eta/2,\mathbb{F},\|\cdot\|_{\infty,S^{M}})$ in Lemma 16, and we bound the
other term similarly in Lemma 17 below. Combining these bounds,

$\log N(\eta^{*},{\mathcal{P}}(\mathbb{R}^{d}),m^{S})\leq
C_{d}N(a,S^{M+a})\left(\log\frac{c_{d,\underline{k},\bar{k}}}{\eta}\right)^{d+1},$

where $a=\sqrt{-2\bar{k}\log\left(\sqrt{\underline{k}\land
1}\frac{(2\pi\underline{k})^{d/2}}{5}\eta\right)}$. Take
$\eta=\rho/n=(2\pi)^{-d/2}/n^{2}$.

$a=\sqrt{4\bar{k}\log n+2\bar{k}\log\left(\frac{5}{\sqrt{\underline{k}\land
1}\underline{k}^{d/2}}\right)}\in\left[\sqrt{2\bar{k}\log
n},\sqrt{6\bar{k}\log n}\right],$

provided $1/n\leq\frac{5}{\sqrt{\underline{k}\land 1}\underline{k}^{d/2}}\leq
n$. Hence by Saha and Guntuboyina (2020b, Lemma F.6) (see the argument on page
6) gives

$\log N\leq c_{d,\underline{k},\bar{k}}(\log n)^{1+d/2}\mathrm{Vol}(S^{2M}).$

######  Lemma 17.

For all compact $S\subset\mathbb{R}^{d}$, $M>0$ and $\eta>0$ sufficiently
small,

$\log N(\eta,\mathbb{F},\|\cdot\|_{\nabla,S^{M}})\leq
C_{d}N(a,S^{a})\left(\log\frac{c_{d,\underline{k},\bar{k}}}{\eta}\right)^{d+1}$

where
$a=\sqrt{2\bar{k}\log\frac{c_{d,\overline{k},\underline{k}}^{\prime}}{\eta}}$.

###### Proof.

Fix $G\in{\mathcal{P}}(\mathbb{R}^{d})$. By Lemma 15, there is a discrete
measure $H$ supported on $S^{a}$ with at most

$l\coloneqq\left(2\lfloor 13.5a^{2}/\bar{k}\rfloor+2\right)^{d}N(a,S^{a})+1$

atoms such that

$\|f_{G,\bullet}-f_{H,\bullet}\|_{\nabla,S^{a}}\leq
a\left(1+\frac{3}{\sqrt{2\pi}}\right)(2\pi\underline{k})^{-d/2}e^{-a^{2}/(2\bar{k})}$

Now let $\mathcal{C}$ denote a minimal $\alpha$-net of $S^{a}$. Write
$H=\sum_{j}w_{j}\delta_{a_{j}}$, and define
$H^{\prime}=\sum_{j}w_{j}\delta_{b_{j}}$ where $b_{j}\in\mathcal{C}$ is the
closest element to $a_{j}$. Then

$\displaystyle\left\|\nabla f_{H,\Sigma_{i}}(x)-\nabla
f_{H^{\prime},\Sigma_{i}}(x)\right\|_{2}$
$\displaystyle\leq\sum_{j}w_{j}\left\|\nabla\varphi_{\Sigma_{i}}(x-a_{j})-\nabla\varphi_{\Sigma_{i}}(x-b_{j})\right\|_{2}$
$\displaystyle\leq\underline{k}^{-1/2}|\Sigma_{i}|^{-1/2}\sum_{j}w_{j}\left\|\nabla\varphi\left(\Sigma_{i}^{-1/2}(x-a_{j})\right)-\nabla\varphi\left(\Sigma_{i}^{-1/2}(x-b_{j})\right)\right\|_{2}$
$\displaystyle\leq\frac{\underline{k}^{-d/2-1/2}\alpha}{(2\pi)^{d/2}}\left[1+\frac{2}{e}+\frac{\alpha}{\sqrt{\underline{k}e}}\right]$

Now let $\mathcal{D}$ denote a minimal $\beta$-net of $\Delta_{l-1}$ under
$\|\cdot\|_{1}$. Let $H^{\prime\prime}=\sum_{j}w_{j}^{\prime}\delta_{b_{j}}$
where $\|w^{\prime}-w\|_{1}\leq\beta$. Then

$\left\|\nabla f_{H^{\prime},\Sigma_{i}}(x)-\nabla
f_{H^{\prime\prime},\Sigma_{i}}(x)\right\|_{2}\leq\beta\sup_{u}\|\nabla\varphi_{\Sigma_{i}}(u)\|_{2}\leq\frac{\underline{k}^{-d/2-1/2}\beta}{(2\pi)^{d/2}\sqrt{e}}.$

By triangle inequality,

$\left\|f_{G,\bullet}-f_{H^{\prime\prime},\bullet}\right\|_{\nabla,S^{M}}\leq(2\pi\underline{k})^{-d/2}\left[\left(1+\frac{3}{\sqrt{2\pi}}\right)ae^{-a^{2}/(2\bar{k})}+\frac{\alpha}{\sqrt{\underline{k}}}\left[1+\frac{2}{e}+\frac{\alpha}{\sqrt{\underline{k}e}}\right]+\frac{\beta}{\sqrt{\underline{k}e}}\right]$

Taking $a=\sqrt{2\bar{k}\log\alpha^{-1}}\geq 1$ and $\alpha=\beta=\sqrt{k\land
1}\frac{(2\pi\underline{k})^{d/2}}{5}\eta$,

$\left\|f_{G,\bullet}-f_{H^{\prime\prime},\bullet}\right\|_{\nabla,S^{M}}\leq\frac{5a\alpha}{\sqrt{k\land
1}}(2\pi\underline{k})^{-d/2}=a\eta$

The proof is completed following same steps as the proof of Lemma 5. ∎

###  C.4 Putting together the pieces

Combining (39), (42), (45), (46), and (49) and pulling out any constants
depending on $d,\underline{k},$ or $\bar{k}$,

$\displaystyle\mathbb{E}\|\hat{\theta}-\hat{\theta}^{*}\|_{F}^{2}/n$
$\displaystyle\leq(5/n)\left[\mathbb{E}\|\hat{\theta}^{*}_{\rho}-\hat{\theta}^{*}\|_{F}^{2}+\sum_{t=1}^{4}\mathbb{E}\zeta_{t}^{2}\right]$
(53) $\displaystyle\leq
c_{d,\underline{k},\bar{k}}\bigg{(}\varepsilon_{n}^{2}(M,S,G^{*})(\sqrt{\log
n})^{d-2}$ $\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ +\frac{\log n}{n}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ +(\eta^{*})^{2}+\varepsilon_{n}^{2}(M,S,G^{*})$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ +\log(eN)\frac{(\log n)^{2}}{n}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ +\varepsilon_{n}^{2}(M,S,G^{*})(\sqrt{\log
n})^{6}\bigg{)}$ $\displaystyle\leq
c_{d,\underline{k},\bar{k}}\varepsilon_{n}^{2}(M,S,G^{*})(\sqrt{\log
n})^{(d-2)\lor 6}$

This completes the proof of Theorem 8.

##  Appendix D Proofs of Theorems 10 and 12

###### Proof of Theorem 10.

We will relate the Wasserstein distance to the average Hellinger distance, so
we rely on the tools of Nguyen (2013, proof of theorem 2). Fix a symmetric
density $K$ whose Fourier transform $\widetilde{K}$ is bounded with support on
$[-1,1]^{d}$. For any $\delta>0$ define the scaled kernel
$K_{\delta}(x)=\frac{1}{\delta^{d}}K(x/\delta)$. By the triangle inequality,

$W_{2}(G^{*},{\widehat{G}_{n}})\leq
W_{2}(G^{*},G^{*}*K_{\delta})+W_{2}(G^{*}*K_{\delta},{\widehat{G}_{n}}*K_{\delta})+W_{2}({\widehat{G}_{n}}*K_{\delta},{\widehat{G}_{n}}).$

For the first and third terms, bound the minimum over all couplings by the
strong coupling:

$W_{2}^{2}(G,G*K_{\delta})=\min_{\theta\sim G,\theta^{\prime}\sim
G,\varepsilon\sim
K}\mathbb{E}\|\theta-(\theta^{\prime}+\delta\varepsilon)\|_{2}^{2}\leq\delta^{2}\mathbb{E}_{\varepsilon\sim
K}\|\varepsilon\|_{2}^{2},$

where the inequality follows from choosing the coupling where
$\theta=\theta^{\prime}$ almost surely. Letting
$m_{2}(K)=\mathbb{E}_{\varepsilon\sim K}\|\varepsilon\|_{2}^{2}$ denote the
second moment of the (unscaled) kernel, we have

$\displaystyle W_{2}(G^{*},{\widehat{G}_{n}})\leq
2\sqrt{m_{2}(K)}\delta+W_{2}(G^{*}*K_{\delta},{\widehat{G}_{n}}*K_{\delta}).$
(54)

For the second term, Villani (2008, Theorem 6.15) yields

$W_{2}^{2}(G^{*}*K_{\delta},{\widehat{G}_{n}}*K_{\delta})\leq
2\int\|x\|_{2}^{2}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\left|G^{*}*K_{\delta}-{\widehat{G}_{n}}*K_{\delta}\right|(x)$

By Nguyen (2013, Lemma 6), for any $s>2$ such that
$m_{s}(K)=\mathbb{E}\|\varepsilon\|_{2}^{s}<\infty$,

$\displaystyle W_{2}^{2}(G^{*}*K_{\delta},{\widehat{G}_{n}}*K_{\delta})$
$\displaystyle\leq
4\left\|G^{*}*K_{\delta}-{\widehat{G}_{n}}*K_{\delta}\right\|_{L_{1}}^{(s-2)/s}R^{2/s}$
$\displaystyle\leq
4\left[2\mathrm{Vol}(B_{1})^{s/(d+2s)}R^{d/(d+2s)}\left\|G^{*}*K_{\delta}-{\widehat{G}_{n}}*K_{\delta}\right\|_{L_{2}}^{2s/(d+2s)}\right]^{(s-2)/s}R^{2/s}$
$\displaystyle=4\cdot 2^{(s-2)/s}\mathrm{Vol}(B_{1})^{(s-2)/(2s+d)}\cdot
R^{d(s-2)/(s(d+2s))+2/s}\left\|G^{*}*K_{\delta}-{\widehat{G}_{n}}*K_{\delta}\right\|_{L_{2}}^{2(s-2)/(d+2s)}$
$\displaystyle\leq 8\sqrt{\mathrm{Vol}(B_{1})}\cdot
R^{d(s-2)/(s(d+2s))+2/s}\left\|G^{*}*K_{\delta}-{\widehat{G}_{n}}*K_{\delta}\right\|_{L_{2}}^{2(s-2)/(d+2s)}$

where $R\coloneqq\mathbb{E}_{\theta^{*}\sim G^{*},\varepsilon\sim
K}\|\theta^{*}+\delta\varepsilon\|_{2}^{s}+\mathbb{E}_{\theta\sim{\widehat{G}_{n}},\varepsilon\sim
K}\|\theta+\delta\varepsilon\|_{2}^{s}$.

For moments in the term $R$, use
$\mathbb{E}\|\theta+\delta\varepsilon\|_{2}^{s}\leq
2^{s}\left(\mathbb{E}\|\theta\|_{2}^{s}+\delta^{s}m_{s}(K)\right),$ so

$R\leq 2^{s}(m_{s}(G^{*})+m_{s}({\widehat{G}_{n}})+2\delta^{s}m_{s}(K)).$

The quantity $m_{s}(K)$ is regarded as a constant depending only on $s>2$ and
$d$. By assumption, the support of ${\widehat{G}_{n}}$ is contained in the
minimum bounding box of the observations, which is further contained in
$[-U,U]^{d}$ where $U=\max_{i,j}|X_{ij}|\leq
L+\max_{i,j}|X_{ij}-\theta^{*}_{ij}|$. Since
$X_{ij}-\theta^{*}_{ij}\stackrel{{\scriptstyle\text{ind}}}{{\sim}}\mathcal{N}(0,(\Sigma_{i})_{jj})$,
we have by a standard concentration argument that

$U\leq L+4\sqrt{\overline{k}\log n}$

with probability at least $1-\frac{2d}{n^{8}}$. Hence, with the same
probability

$m_{s}({\widehat{G}_{n}})=\mathbb{E}_{{\widehat{G}_{n}}}\|\theta\|_{2}^{s}\leq
d^{s/2}\mathbb{E}_{{\widehat{G}_{n}}}\|\theta\|_{\infty}^{s}\leq
d^{s/2}\left(L+4\sqrt{\overline{k}\log n}\right)^{s}.$

This same bound holds for $m_{s}(G^{*})$.

For the $\|\cdot\|_{L_{2}}$ norm
$\left\|G^{*}*K_{\delta}-{\widehat{G}_{n}}*K_{\delta}\right\|_{L_{2}}$, let
$g_{\delta}^{(i)}$ denote the inverse Fourier transform of
$\widetilde{K_{\delta}}/\widetilde{\varphi_{\Sigma_{i}}}$, so that
$G*K_{\delta}=f_{G,\Sigma_{i}}*g_{\delta}^{(i)}$. Hence, by Proposition 8.49
of Folland (1999), we have for each $i=1,\dots,n$,

$\left\|G^{*}*K_{\delta}-{\widehat{G}_{n}}*K_{\delta}\right\|_{L_{2}}\leq
2d_{TV}(f_{G^{*},\Sigma_{i}},f_{{\widehat{G}_{n}},\Sigma_{i}})\|g_{\delta}^{(i)}\|_{L_{2}}.$

Using Plancherel’s theorem and the fact that $\widetilde{K}$ is bounded on its
support of $[-1,1]^{d}$,

$\displaystyle\|g_{\delta}^{(i)}\|_{L_{2}}^{2}$
$\displaystyle=\int_{\mathbb{R}^{d}}\frac{\widetilde{K}(\delta\omega)^{2}}{\widetilde{\varphi_{\Sigma_{i}}}(\omega)^{2}}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\omega\leq
C_{d}\int_{[-1/\delta,1/\delta]^{d}}\widetilde{\varphi_{\Sigma_{i}}}(\omega)^{-2}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\omega$
$\displaystyle=C_{d}\int_{[-1/\delta,1/\delta]^{d}}\exp(\omega^{\prime}\Sigma_{i}\omega)\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\omega=C_{d}\prod_{j=1}^{d}\int_{-1/\delta}^{1/\delta}\exp((\Sigma_{i})_{jj}\omega_{j}^{2})\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\omega_{j}$ $\displaystyle\leq
C_{d}\prod_{j=1}^{d}e^{2(\Sigma_{i})_{jj}\delta^{-2}}\int_{-1/\delta}^{1/\delta}\exp(-(\Sigma_{i})_{jj}\omega_{j}^{2})\mathop{\mathrm{\mathstrut
d}}\nolimits\\!\omega_{j}\leq
C_{d}\left(\frac{\pi}{\underline{k}}\right)^{d/2}e^{2d\overline{k}\delta^{-2}}.$

Averaging over $i=1,\dots,n$,

$\displaystyle\left\|G^{*}*K_{\delta}-{\widehat{G}_{n}}*K_{\delta}\right\|_{L_{2}}$
$\displaystyle\leq
C_{d}\underline{k}^{-d/2}e^{2d\overline{k}\delta^{-2}}\cdot\frac{1}{n}\sum_{i=1}^{n}d_{TV}(f_{G^{*},\Sigma_{i}},f_{{\widehat{G}_{n}},\Sigma_{i}})$
$\displaystyle\leq
C_{d}\underline{k}^{-d/2}e^{2d\overline{k}\delta^{-2}}\bar{h}(G^{*},{\widehat{G}_{n}}).$

Combining our calculations following (54), we have

$\displaystyle W_{2}(G^{*},{\widehat{G}_{n}})$ $\displaystyle\leq
C_{d,s}\inf_{\delta\in(0,1)}\bigg{\\{}\delta+\left(2^{s}\left(d^{s/2}\left(L+4\sqrt{\overline{k}\log
n}\right)^{s}+2\delta^{s}m_{s}(K)\right)\right)^{3d/(2d+4s)}$ (55)
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
\times\left(\underline{k}^{-d/2}e^{2d\overline{k}\delta^{-2}}\bar{h}(G^{*},{\widehat{G}_{n}})\right)^{(s-2)/(d+2s)}\bigg{\\}}.$

Assume $n$ is large enough that $4\sqrt{\overline{k}\log n}\geq L$ and
$2^{s}d^{s/2}\left(4\sqrt{\overline{k}\log n}\right)^{s}\geq 2m_{s}(K)$, so

$\displaystyle W_{2}(G^{*},{\widehat{G}_{n}})$ $\displaystyle\leq
C_{d,s}\inf_{\delta\in(0,1)}\bigg{\\{}\delta+\left(\overline{k}\log
n\right)^{3sd/(4(d+2s))}\left(\underline{k}^{-d/2}e^{2d\overline{k}\delta^{-2}}\bar{h}(G^{*},{\widehat{G}_{n}})\right)^{(s-2)/(d+2s)}\bigg{\\}}.$
(56)

Choosing
$\delta^{-2}=-\frac{1}{4\overline{k}d}\log\bar{h}(G^{*},{\widehat{G}_{n}})$
(provided $\delta<1$) and $s=d+2$ yields

$\displaystyle W_{2}^{2}(G^{*},{\widehat{G}_{n}})$ $\displaystyle\leq
C_{d}\bigg{\\{}\frac{\overline{k}d}{-\log\bar{h}(G^{*},{\widehat{G}_{n}})}+\left(\overline{k}\log
n\right)^{d/2}\left(\underline{k}^{-d}\bar{h}(G^{*},{\widehat{G}_{n}})\right)^{1/12}\bigg{\\}}.$
(57)

$\varepsilon_{n}^{2}(M,S,G^{*})$ defined in (13), with $S=[-L,L]^{d}$ and
$M=\sqrt{10\overline{k}\log n}$ gives

$\varepsilon_{n}^{2}=\frac{\left(4\sqrt{10}\left(L^{2}\lor\overline{k}\right)\right)^{d}}{n}(\log
n)^{d+1}.$

By Theorem 6,

$\displaystyle W_{2}^{2}(G^{*},{\widehat{G}_{n}})$ $\displaystyle\leq
C_{d}\bigg{\\{}\frac{\overline{k}d}{\log n-\log
C_{d,\overline{k},\underline{k},L}t^{2}(\log n)^{d+1}}+\left(\overline{k}\log
n\right)^{d/2}\left(C_{d,\overline{k},\underline{k},L}t^{2}\frac{(\log
n)^{d+1}}{n}\right)^{1/24}\bigg{\\}},$

with probability at least $1-2n^{-t^{2}}$. Take $t^{2}=8$. For $n$
sufficiently large the first term dominates, $\delta<1$, and $\log n-\log
C_{d,\overline{k},\underline{k},L}8(\log n)^{d+1}\geq(\log n)/2$. ∎

###### Proof of Theorem 12.

Take $\mu=0$ by location equivariance (see Lemma 4). Write
${\widehat{G}_{n}}=\sum_{j=1}^{\hat{k}}\hat{w}_{j}\delta_{\hat{a}_{j}}$. Since
$G^{*}=\delta_{0}$ is a point mass,

$W_{2}^{2}({\widehat{G}_{n}},G^{*})=\mathbb{E}_{\vartheta\sim{\widehat{G}_{n}}}\|\vartheta\|_{2}^{2}=\sum_{j=1}^{\hat{k}}\hat{w}_{j}\|\hat{a}_{j}\|_{2}^{2}.$

We relate this to the marginal density $f_{{\widehat{G}_{n}},\bullet}$ via

$\displaystyle\int
f_{{\widehat{G}_{n}},\Sigma_{i}}(x)\|x\|_{2}^{2}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!x$
$\displaystyle=\sum_{j=1}^{\hat{k}}\hat{w}_{j}\int\varphi_{\Sigma_{i}}(x-\hat{a}_{j})\|x\|_{2}^{2}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!x$
$\displaystyle=\sum_{j=1}^{\hat{k}}\hat{w}_{j}\int\varphi_{\Sigma_{i}}(x)(\|x\|_{2}^{2}+\|\hat{a}_{j}\|_{2}^{2}+2\langle
x,\hat{a}_{j}\rangle)\mathop{\mathrm{\mathstrut d}}\nolimits\\!x$
$\displaystyle=\int
f_{G^{*},\Sigma_{i}}(x)\|x\|_{2}^{2}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!x+W_{2}^{2}({\widehat{G}_{n}},G^{*}).$

Hence for any $i\in\\{1,\dots,n\\}$,

$\displaystyle W_{2}^{2}({\widehat{G}_{n}},G^{*})$
$\displaystyle=\int(f_{{\widehat{G}_{n}},\Sigma_{i}}(x)-f_{G^{*},\Sigma_{i}}(x))\|x\|_{2}^{2}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!x$ $\displaystyle\leq
2h(f_{{\widehat{G}_{n}},\Sigma_{i}},f_{G^{*},\Sigma_{i}})\left(\int\left(f_{{\widehat{G}_{n}},\Sigma_{i}}(x)+f_{G^{*},\Sigma_{i}}(x)\right)\|x\|_{2}^{4}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!x\right)^{1/2}.$

Averaging over $i\in\\{1,\dots,n\\}$,

$W_{2}^{2}({\widehat{G}_{n}},G^{*})\leq
2\bar{h}\left(f_{{\widehat{G}_{n}},\bullet},f_{G^{*},\bullet}\right)\max_{i=1:n}\left(\int\left(f_{{\widehat{G}_{n}},\Sigma_{i}}(x)+f_{G^{*},\Sigma_{i}}(x)\right)\|x\|_{2}^{4}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!x\right)^{1/2}.$

Applying Theorem 6 with $S=\\{0\\}$ and $M=\sqrt{10\overline{k}\log n}$,

$\bar{h}^{2}\left(f_{{\widehat{G}_{n}},\bullet},f_{G^{*},\bullet}\right)\lesssim_{d,\overline{k},\underline{k}}t^{2}\frac{(\log
n)^{d+1}}{n}$

with probability at least $1-2n^{-t^{2}}$ for all $t\geq 1$.

For the remaining terms,

$\int f_{G^{*},\Sigma_{i}}(x)\|x\|_{2}^{4}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!x=\mathbb{E}\|X_{i}\|_{2}^{4}=\mathbb{E}_{Z\sim\mathcal{N}(0,I_{d})}\|\Sigma_{i}^{1/2}Z\|_{2}^{4}\leq\overline{k}^{2}\mathbb{E}_{A\sim\chi_{d}^{2}}A^{2}=\overline{k}^{2}d(d+2)$

and

$\displaystyle\int
f_{{\widehat{G}_{n}},\Sigma_{i}}(x)\|x\|_{2}^{4}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!x$
$\displaystyle=\sum_{j}\hat{w}_{j}\int\varphi_{\Sigma_{i}}(x)\|x+\hat{a}_{j}\|_{2}^{4}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!x$ $\displaystyle\leq
8\sum_{j}\hat{w}_{j}\int\varphi_{\Sigma_{i}}(x)\Big{(}\|x\|_{2}^{4}+\|\hat{a}_{j}\|_{2}^{4}\Big{)}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!x$ $\displaystyle\leq
8\overline{k}^{2}d(d+2)+8\sum_{j}\hat{w}_{j}\|\hat{a}_{j}\|_{2}^{4}.$

By our assumption on the support that each $\hat{a}_{j}\in\mathbb{B}_{\kappa
r}(\bar{X})$, each $\hat{a}_{j}$ equivalently satisfies

$\|\hat{a}_{j}-\bar{X}\|_{2}^{4}\leq(\overline{k}/\underline{k})^{4}\max_{i}\|X_{i}-\bar{X}\|_{2}^{4}.$

Noting that
$X_{i}-\bar{X}\sim\mathcal{N}\left(0,(1-2n^{-1})\Sigma_{i}+n^{-1}\bar{\Sigma}\right)$
with $\bar{\Sigma}=n^{-1}\sum_{j=1}^{n}\Sigma_{j}$, we bound

$\displaystyle\sum_{j}\hat{w}_{j}\|\hat{a}_{j}\|_{2}^{4}$ $\displaystyle\leq
8\left(\|\bar{X}\|_{2}^{4}+\max_{j}\|\hat{a}_{j}-\bar{X}\|_{2}^{4}\right)$
$\displaystyle\leq
8\left(\|\bar{X}\|_{2}^{4}+\kappa^{4}\max_{i\in[n]}\|X_{i}-\bar{X}\|_{2}^{4}\right)$
$\displaystyle\leq_{\mathrm{st}}8\overline{k}^{2}\left(n^{-2}A_{0}^{2}+\kappa^{4}\max_{i\in[n]}A_{i}^{2}\right)$
$\displaystyle\leq
16\frac{\overline{k}^{6}}{\underline{k}^{4}}\max_{i=0:n}A_{i}^{2},$

where $A_{0},A_{1},\dots,A_{n}\sim\chi_{d}^{2}$ are possibly dependent, and
$\leq_{\mathrm{st}}$ denotes stochastic inequality.

For $t\geq 1$, we use the following tail bound (see Laurent and Massart, 2000,
Lemma 1)

${\mathbb{P}}\left(\max_{i=0:n}A_{i}^{2}\geq 60t^{2}(\log n)^{2}\right)\leq
n^{-t^{2}},$

where we have used the assumption in Theorem 6 that $n\geq(2\pi)^{d/2}$ to
eliminate the dependence on $d$. We have thus shown that

$\max_{i=1:n}\left(\int\left(f_{{\widehat{G}_{n}},\Sigma_{i}}(x)+f_{G^{*},\Sigma_{i}}(x)\right)\|x\|_{2}^{4}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!x\right)^{1/2}\leq\left(9\overline{k}^{2}d(d+2)+400\frac{\overline{k}^{6}}{\underline{k}^{4}}t^{2}(\log
n)^{2}\right)^{1/2}\lesssim\frac{\overline{k}^{3}}{\underline{k}^{2}}t(\log
n)$

with probability at least $1-n^{-t^{2}}$. Combining with a union bound over
our earlier estimate,

$W_{2}({\widehat{G}_{n}},G^{*})\lesssim_{d,\overline{k},\underline{k}}t^{3/2}\frac{(\log
n)^{(d+3)/4}}{n^{1/4}}$

with probability at least $1-3n^{-t^{2}}$ for all $t\geq 1$. ∎

## References

  * Abolfathi et al. (2018) Bela Abolfathi, DS Aguado, Gabriela Aguilar, Carlos Allende Prieto, Andres Almeida, Tonima Tasnim Ananna, Friedrich Anders, Scott F Anderson, Brett H Andrews, Borja Anguiano, et al.  The fourteenth data release of the Sloan Digital Sky Survey: first spectroscopic data from the extended Baryon Oscillation Spectroscopic Survey and from the second phase of the Apache Point Observatory Galactic Evolution Experiment.  _The Astrophysical Journal Supplement Series_ , 235(2):42, 2018.
  * Akritas and Bershady (1996) Michael G Akritas and Matthew A Bershady.  Linear regression for astronomical data with measurement errors and intrinsic scatter.  _The Astrophysical Journal_ , 470(2):706, 1996\.
  * Améndola et al. (2020) Carlos Améndola, Alexander Engström, and Christian Haase.  Maximum number of modes of Gaussian mixtures.  _Information and Inference: A Journal of the IMA_ , 9(3):587–600, 2020.
  * Anderson et al. (2018) Lauren Anderson, David W Hogg, Boris Leistedt, Adrian M Price-Whelan, and Jo Bovy.  Improving Gaia parallax precision with a data-driven model of stars.  _The Astronomical Journal_ , 156(4):145, 2018\.
  * Banerjee et al. (2021) Trambak Banerjee, Luella J Fu, Gareth M James, and Wenguang Sun.  Nonparametric empirical Bayes estimation on heterogeneous data.  2021\.  URL http://faculty.marshall.usc.edu/gareth-james/Research/Nest%20Biometrika.pdf.
  * Böhning (1985) Dankmar Böhning.  Numerical estimation of a probability measure.  _Journal of statistical planning and inference_ , 11(1):57–69, 1985.
  * Böhning (2003) Dankmar Böhning.  The EM algorithm with gradient function update for discrete mixtures with known (fixed) number of components.  _Statistics and Computing_ , 13(3):257–265, 2003\.
  * Bovy et al. (2011) Jo Bovy, David W Hogg, and Sam T Roweis.  Extreme deconvolution: Inferring complete distribution functions from noisy, heterogeneous and incomplete observations.  _Annals of Applied Statistics_ , 5(2B):1657–1677, 2011.
  * Deb et al. (2021) Nabarun Deb, Sujayam Saha, Adityanand Guntuboyina, and Bodhisattva Sen.  Two-component mixture model in the presence of covariates.  _Journal of the American Statistical Association_ , pages 1–15, 2021\.
  * Dedecker and Michel (2013) Jérôme Dedecker and Bertrand Michel.  Minimax rates of convergence for Wasserstein deconvolution with supersmooth errors in any dimension.  _Journal of Multivariate Analysis_ , 122:278–291, 2013\.
  * Delaigle and Meister (2008) Aurore Delaigle and Alexander Meister.  Density estimation with heteroscedastic error.  _Bernoulli_ , 14(2):562–579, 2008.
  * Dempster et al. (1977) Arthur P Dempster, Nan M Laird, and Donald B Rubin.  Maximum likelihood from incomplete data via the EM algorithm.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 39(1):1–22, 1977.
  * Dicker and Zhao (2016) Lee H Dicker and Sihai D Zhao.  High-dimensional classification via nonparametric empirical Bayes and maximum likelihood inference.  _Biometrika_ , 103(1):21–34, 2016.
  * Doss et al. (2020) Natalie Doss, Yihong Wu, Pengkun Yang, and Harrison H Zhou.  Optimal estimation of high-dimensional Gaussian mixtures.  _arXiv preprint arXiv:2002.05818_ , 2020.
  * Dyson (1926) Frank Dyson.  A method for correcting series of parallax observations.  _Monthly Notices of the Royal Astronomical Society_ , 86:686, 1926.
  * Dytso et al. (2019) Alex Dytso, Semih Yagli, H Vincent Poor, and Shlomo Shamai Shitz.  The capacity achieving distribution for the amplitude constrained additive Gaussian channel: An upper bound on the number of mass points.  _IEEE Transactions on Information Theory_ , 66(4):2006–2022, 2019.
  * Efron (2011) Bradley Efron.  Tweedie’s formula and selection bias.  _Journal of the American Statistical Association_ , 106(496):1602–1614, 2011.
  * Efron (2012) Bradley Efron.  _Large-scale inference: empirical Bayes methods for estimation, testing, and prediction_ , volume 1.  Cambridge University Press, 2012.
  * Efron (2014) Bradley Efron.  Two modeling strategies for empirical Bayes estimation.  _Statistical Science_ , 29(2):285, 2014.
  * Efron (2019) Bradley Efron.  Bayes, oracle Bayes and empirical Bayes.  _Statistical Science_ , 34(2):177–201, 2019.
  * Efron and Hastie (2016) Bradley Efron and Trevor Hastie.  _Computer age statistical inference_ , volume 5.  Cambridge University Press, 2016.
  * Efron and Morris (1972a) Bradley Efron and Carl Morris.  Empirical Bayes on vector observations: An extension of Stein’s method.  _Biometrika_ , 59(2):335–347, 1972a.
  * Efron and Morris (1972b) Bradley Efron and Carl Morris.  Limiting the risk of Bayes and empirical Bayes estimators—Part II: The empirical Bayes case.  _Journal of the American Statistical Association_ , 67(337):130–139, 1972b.
  * Efron and Morris (1973a) Bradley Efron and Carl Morris.  Combining possibly related estimation problems.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 35(3):379–402, 1973a.
  * Efron and Morris (1973b) Bradley Efron and Carl Morris.  Stein’s estimation rule and its competitors—an empirical Bayes approach.  _Journal of the American Statistical Association_ , 68(341):117–130, 1973b.
  * Feng and Dicker (2018) Long Feng and Lee H Dicker.  Approximate nonparametric maximum likelihood for mixture models: A convex optimization approach to fitting arbitrary multivariate mixing distributions.  _Computational Statistics & Data Analysis_, 122:80–91, 2018.
  * Folland (1999) Gerald B Folland.  _Real analysis: modern techniques and their applications_ , volume 40.  John Wiley & Sons, 1999.
  * Gu and Koenker (2016) Jiaying Gu and Roger Koenker.  On a problem of Robbins.  _International Statistical Review_ , 84(2):224–244, 2016.
  * Gu and Koenker (2017) Jiaying Gu and Roger Koenker.  Rebayes: an R package for empirical Bayes mixture methods.  Technical report, cemmap working paper, 2017.
  * Heinrich and Kahn (2018) Philippe Heinrich and Jonas Kahn.  Strong identifiability and optimal minimax rates for finite mixture estimation.  _Annals of Statistics_ , 46(6A):2844–2870, 2018\.
  * Ho and Nguyen (2016) Nhat Ho and XuanLong Nguyen.  On strong identifiability and convergence rates of parameter estimation in finite mixtures.  _Electronic Journal of Statistics_ , 10(1):271–307, 2016.
  * Hogg et al. (2010) David W Hogg, Adam D Myers, and Jo Bovy.  Inferring the eccentricity distribution.  _The Astrophysical Journal_ , 725(2):2166, 2010\.
  * James and Stein (1961) W. James and Charles Stein.  Estimation with quadratic loss.  In _Proceedings of the Fourth Berkeley Symposium on Mathematical Statistics and Probability, Volume 1: Contributions to the Theory of Statistics_ , pages 361–379, Berkeley, Calif., 1961. University of California Press.  URL https://projecteuclid.org/euclid.bsmsp/1200512173.
  * Jiang et al. (2011) Jiming Jiang, Thuan Nguyen, and J Sunil Rao.  Best predictive small area estimation.  _Journal of the American Statistical Association_ , 106(494):732–745, 2011.
  * Jiang (2020) Wenhua Jiang.  On general maximum likelihood empirical Bayes estimation of heteroscedastic IID normal means.  _Electronic Journal of Statistics_ , 14(1):2272–2297, 2020.
  * Jiang and Zhang (2009) Wenhua Jiang and Cun-Hui Zhang.  General maximum likelihood empirical Bayes estimation of normal means.  _The Annals of Statistics_ , 37(4):1647–1684, 2009.
  * Jiang and Zhang (2010) Wenhua Jiang and Cun-Hui Zhang.  Empirical Bayes in-season prediction of baseball batting averages.  In _Borrowing Strength: Theory Powering Applications–A Festschrift for Lawrence D. Brown_ , pages 263–273. Institute of Mathematical Statistics, 2010.
  * Johnstone (2019) Iain M Johnstone.  Gaussian estimation: Sequence and wavelet models.  2019\.  URL http://www-stat.stanford.edu/~imj/.
  * Kelly (2012) Brandon C Kelly.  Measurement error models in astronomy.  In _Statistical challenges in modern astronomy V_ , pages 147–162. Springer, 2012.
  * Kiefer and Wolfowitz (1956) Jack Kiefer and Jacob Wolfowitz.  Consistency of the maximum likelihood estimator in the presence of infinitely many incidental parameters.  _The Annals of Mathematical Statistics_ , pages 887–906, 1956.
  * Kim and Guntuboyina (2020) Arlene KH Kim and Adityanand Guntuboyina.  Minimax bounds for estimating multivariate Gaussian location mixtures.  _arXiv preprint arXiv:2012.00444_ , 2020.
  * Kim et al. (2019) Youngseok Kim, Peter Carbonetto, Matthew Stephens, and Mihai Anitescu.  A fast algorithm for maximum likelihood estimation of mixture proportions using sequential quadratic programming.  _Journal of Computational and Graphical Statistics_ , (just-accepted):1–34, 2019.
  * Kim et al. (2020) Youngseok Kim, Peter Carbonetto, Matthew Stephens, and Mihai Anitescu.  A fast algorithm for maximum likelihood estimation of mixture proportions using sequential quadratic programming.  _Journal of Computational and Graphical Statistics_ , 29(2):261–273, 2020.
  * Koenker and Gu (2017) Roger Koenker and Jiaying Gu.  REBayes: Empirical Bayes mixture methods in R.  _Journal of Statistical Software_ , 82(8):1–26, 2017.
  * Koenker and Mizera (2014) Roger Koenker and Ivan Mizera.  Convex optimization, shape constraints, compound decisions, and empirical Bayes rules.  _Journal of the American Statistical Association_ , 109(506):674–685, 2014.
  * Laird (1978) Nan Laird.  Nonparametric maximum likelihood estimation of a mixing distribution.  _Journal of the American Statistical Association_ , 73(364):805–811, 1978.
  * Lashkari and Golland (2008) Danial Lashkari and Polina Golland.  Convex clustering with exemplar-based models.  In _Advances in neural information processing systems_ , pages 825–832, 2008.
  * Laurent and Massart (2000) Beatrice Laurent and Pascal Massart.  Adaptive estimation of a quadratic functional by model selection.  _Annals of Statistics_ , pages 1302–1338, 2000.
  * Lesperance and Kalbfleisch (1992) Mary L Lesperance and John D Kalbfleisch.  An algorithm for computing the nonparametric MLE of a mixing distribution.  _Journal of the American Statistical Association_ , 87(417):120–126, 1992.
  * Lindsay (1995) Bruce G Lindsay.  Mixture models: theory, geometry and applications.  In _NSF-CBMS regional conference series in probability and statistics_. JSTOR, 1995.
  * Lindsay and Roeder (1993) Bruce G Lindsay and Kathryn Roeder.  Uniqueness of estimation and identifiability in mixture models.  _Canadian Journal of Statistics_ , 21(2):139–147, 1993.
  * Liu and Zhu (2007) Lei Liu and Yu Zhu.  Partially projected gradient algorithms for computing nonparametric maximum likelihood estimates of mixing distributions.  _Journal of Statistical Planning and Inference_ , 137(7):2509–2522, 2007.
  * Majewski et al. (2017) Steven R Majewski, Ricardo P Schiavon, Peter M Frinchaboy, Carlos Allende Prieto, Robert Barkhouser, Dmitry Bizyaev, Basil Blank, Sophia Brunner, Adam Burton, Ricardo Carrera, et al.  The Apache Point Observatory Galactic Evolution Experiment (APOGEE).  _The Astronomical Journal_ , 154(3):94, 2017.
  * Marriott (2002) Paul Marriott.  On the local geometry of mixture models.  _Biometrika_ , 89(1):77–93, 2002.
  * Meister (2009) Alexander Meister.  _Deconvolution problems in nonparametric statistics_ , volume 193\.  Springer Science & Business Media, 2009.
  * Morris (1983) Carl N Morris.  Parametric empirical Bayes inference: theory and applications.  _Journal of the American statistical Association_ , 78(381):47–55, 1983.
  * MOSEK ApS (2019) MOSEK ApS.  Mosek optimization suite.  2019\.  URL http://docs.mosek.com/9.0/intro.pdf.
  * Nguyen (2013) XuanLong Nguyen.  Convergence of latent mixing measures in finite and infinite mixture models.  _The Annals of Statistics_ , 41(1):370–400, 2013\.
  * Pfanzagl (1988) Johann Pfanzagl.  Consistency of maximum likelihood estimators for certain nonparametric families, in particular: mixtures.  _Journal of Statistical Planning and Inference_ , 19(2):137–158, 1988.
  * Polyanskiy and Wu (2020) Yury Polyanskiy and Yihong Wu.  Self-regularizing property of nonparametric maximum likelihood estimator in mixture models.  _arXiv preprint arXiv:2008.08244_ , 2020.
  * Ratcliffe et al. (2020) Bridget L Ratcliffe, Melissa K Ness, Kathryn V Johnston, and Bodhisattva Sen.  Tracing the assembly of the Milky Way’s disk through abundance clustering.  _The Astrophysical Journal_ , 900(2):165, 2020\.
  * Ray and Lindsay (2005) Surajit Ray and Bruce G Lindsay.  The topography of multivariate normal mixtures.  _The Annals of Statistics_ , 33(5):2042–2065, 2005.
  * Redner and Walker (1984) Richard A Redner and Homer F Walker.  Mixture densities, maximum likelihood and the EM algorithm.  _SIAM review_ , 26(2):195–239, 1984.
  * Rigollet and Weed (2018) Philippe Rigollet and Jonathan Weed.  Entropic optimal transport is maximum-likelihood deconvolution.  _Comptes Rendus Mathematique_ , 356(11-12):1228–1235, 2018.
  * Ritchie and Murray (2019) James A Ritchie and Iain Murray.  Scalable extreme deconvolution.  _arXiv preprint arXiv:1911.11663_ , 2019.
  * Robbins (1950) Herbert Robbins.  A generalization of the method of maximum likelihood-estimating a mixing distribution.  In _Annals of Mathematical Statistics_ , volume 21, pages 314–315, 1950.
  * Robbins (1956) Herbert Robbins.  An empirical Bayes approach to statistics.  _Proceedings of the Third Berkeley Symposium on Mathematical Statistics and Probability_ , 1:157–163, 1956.
  * Saha and Guntuboyina (2020a) Sujayam Saha and Adityanand Guntuboyina.  On the nonparametric maximum likelihood estimator for Gaussian location mixture densities with application to Gaussian denoising.  _Annals of Statistics_ , 48(2):738–762, 2020a.
  * Saha and Guntuboyina (2020b) Sujayam Saha and Adityanand Guntuboyina.  Supplement to “On the nonparametric maximum likelihood estimator for Gaussian location mixture densities with application to Gaussian denoising”.  2020b.  URL https://doi.org/10.1214/19-AOS1817SUPP.
  * Tan (2016) Zhiqiang Tan.  Steinized empirical Bayes estimation for heteroscedastic data.  _Statistica Sinica_ , pages 1219–1248, 2016.
  * van de Geer (2000) Sara van de Geer.  _Empirical Processes in M-estimation_ , volume 6.  Cambridge university press, 2000.
  * Villani (2008) Cédric Villani.  _Optimal transport: old and new_ , volume 338.  Springer Science & Business Media, 2008.
  * Wang (2007) Yong Wang.  On fast computation of the non-parametric maximum likelihood estimate of a mixing distribution.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 69(2):185–198, 2007.
  * Weinstein et al. (2018) Asaf Weinstein, Zhuang Ma, Lawrence D Brown, and Cun-Hui Zhang.  Group-linear empirical Bayes estimates for a heteroscedastic normal mean.  _Journal of the American Statistical Association_ , 113(522):698–710, 2018.
  * Wong and Shen (1995) Wing Hung Wong and Xiaotong Shen.  Probability inequalities for likelihood ratios and convergence rates of sieve MLEs.  _The Annals of Statistics_ , 23(2):339–362, 1995\.
  * Wu and Yang (2020) Yihong Wu and Pengkun Yang.  Optimal estimation of Gaussian mixtures via denoised method of moments.  _The Annals of Statistics_ , 48(4):1981–2007, 2020.
  * Xie et al. (2012) Xianchao Xie, SC Kou, and Lawrence D Brown.  SURE estimates for a heteroscedastic hierarchical model.  _Journal of the American Statistical Association_ , 107(500):1465–1479, 2012.
  * Zhang (2009) Cun-Hui Zhang.  Generalized maximum likelihood estimation of normal mixture densities.  _Statistica Sinica_ , pages 1297–1318, 2009.