text:
\section{Introduction}
\label{sect:intro}

Consider the problem of estimating an unknown parameter $x\in\mathbb{R}^n$ based on a linear measurement $y\in\mathbb{R}^m$ corrupted by additive noise $w \in \mathbb{R}^m$. This setup is formalized through the linear measurement model
\begin{equation}
\label{eq:linear:observation}
y = H x + w,
\end{equation}
where the observation matrix $H \in \mathbb{R}^{m \times n}$ is assumed to be known. We further assume that the distribution~$\mathbb{P}_w$ of~$w$ has finite second moments and that~$w$ is independent of~$x$. Thus, the conditional distribution~$\mathbb{P}_{y|x}$ of~$y$ given $x$ is obtained by shifting~$\mathbb{P}_w$ by~$Hx$.
We emphasize that none of the subsequent results rely on a particular ordering of the dimension~$n$ of the parameter~$x$ and the dimension~$m$ of the measurement~$y$.
The linear measurement model~\eqref{eq:linear:observation} is fundamental for numerous applications in engineering ({\em e.g.}, linear systems theory \cite{ref:golnaraghi2017automatic, ref:ogata2009modern}), econometrics ({\em e.g.}, linear regression \cite{ref:stock2015introduction, ref:wooldridge2010econometric}, time series analysis \cite{ref:chatfield2016analysis, ref:hamilton1994time}), machine learning and signal processing ({\em e.g.}, Kalman filtering \cite{ref:kay1993fundamentals, ref:murphy2012machine, ref:oppenheim2015signals}) or information theory ({\em e.g.}, multiple-input multiple-output systems \cite{ref:cover2006elements, ref:mackay2003information}) etc.
In addition, model~\eqref{eq:linear:observation} also emerges naturally in many applications in operations research such as traffic management and control~\cite{ref:lint2012applications}, inventory control~\cite{ref:aviv2003time}, advertising and promotion budgeting~\cite{ref:sriram2007optimal} or resource management~\cite{ref:rubel2017robust}.

An estimator of $x$ given $y$ is a measurable function $\psi: \mathbb{R}^m \rightarrow \mathbb{R}^n$ that grows at most linearly. Thus, there exists $C>0$ such that $|\psi(y)|\le C(1+\|y\|)$ for all $y\in\mathbb{R}^m$. The function value $\psi(y)$ is the prediction of $x$ based on the measurement $y$ under the estimator $\psi$. In the following we denote the family of all estimators by~${\mc F}$. The quality of an estimator is measured by a risk function $R:{\mc F}\times\mathbb{R}^n\rightarrow \mathbb{R}$, which quantifies the mismatch between the parameter $x$ and its prediction $\psi(y)$. A popular risk function is the mean square error~(MSE)
\begin{align*}
	R (\psi, x) =  \mathds{E}_{\mathbb{P}_{y|x}} \left[ \| x - \psi(y)\|^2 \right],
\end{align*}
which defines the estimation error as the expected squared Euclidean distance between $\psi(y)$ and $x$. If $x$ was known, then $R (\psi, x)$ could be minimized directly, and the constant estimator $\psi^\star(y)\equiv x$ would be optimal. In practice, however, $x$ is unobservable. Otherwise there would be no need to solve an estimation problem in the first place. With $x$ unknown, it is impossible to minimize the MSE directly. The statistics literature proposes two complementary workarounds for this problem: the Bayesian approach and the minimax approach.

The Bayesian statistician treats $x$ as a random vector governed by a {\em prior} distribution~$\mathbb{P}_x$ that captures her beliefs about $x$ before seeing~$y$ \cite[\S~11.4]{ref:kay1993fundamentals} and solves the minimum MSE (MMSE) estimation problem
\begin{align}\label{eq:Bayes}
	\mathop{\text{minimize}}_{\psi\in{\mc F}} ~ \mathds{E}_{\mathbb{P}_x} \left[ R(\psi, x) \right] .
\end{align}
If the distribution $\mathbb{P}_x$ of~$x$ has finite second moments, then~\eqref{eq:Bayes} is solvable. In this case, the optimal estimator, which is usually termed the Bayesian MMSE estimator, is of the form $\psi^\star_{\mathcal B}(y)=\mathds{E}_{\mathbb{P}_{x|y}} [x]$, where the conditional distribution~$\mathbb{P}_{x|y}$ of~$x$ given $y$ is obtained from $\mathbb{P}_x$ and $\mathbb{P}_{y|x}$ via Bayes' theorem. However, the Bayesian MMSE estimator suffers from two conceptual shortcomings. First, $\psi^\star_{\mathcal B}$ is highly sensitive to the prior distribution~$\mathbb{P}_x$, which is troubling if the statistician has little confidence in her beliefs. Second, computing $\psi^\star_{\mathcal B}$ requires precise knowledge of the noise distribution $\mathbb{P}_w$, which is typically unobservable and thus uncertain at least to some extent. Moreover, $\psi^\star_{\mathcal B}$ may generically have a complicated functional form, and evaluating $\psi^\star_{\mathcal B}(y)$ to high precision for a particular measurement $y$ ({\em e.g.}, via Monte Carlo simulation) may be computationally challenging if the dimension of $x$ is high.

These shortcomings are mitigated if we restrict the space ${\mc F}$ of all measurable estimators in~\eqref{eq:Bayes} to the~space
\begin{align} \label{eq:affine}
	\Ac = \big\{ \psi \in{\mc F}\;:\; \existsA \in \mathbb{R}^{n \times m}, \, b \in \mathbb{R}^n \text{ with } \psi(y) = A y + b ~\forall y\in\mathbb{R}^m  \big\}
\end{align}
of all {\em affine} estimators. In this case the distributions $\mathbb{P}_x$ and $\mathbb{P}_w$ need not be fully known. Instead, in order to evaluate the optimal affine estimator $\psi_{\mathcal A}^\star(y)= A^\star y+b^\star$, it is sufficient to know the mean vectors $\mu_x$ and $\mu_w$ as well as the covariance matrices $\Sigma_x$ and $\Sigma_w$ of the distributions $\mathbb{P}_x$ and $\mathbb{P}_w$, respectively. If $H \Sigma_x H^\top + \Sigma_w \succ 0$, which is the case if the noise covariance matrix has full rank, then the coefficients of the best affine estimator can be computed in closed form. Using~\eqref{eq:linear:observation} together with the independence of $x$ and $w$ one can show that
\begin{align}
	\label{eq:optimal:affine}
	A^\star = \Sigma_x H^\top (H \Sigma_x H^\top + \Sigma_w)^{-1} \quad\text{and}\quad
	b^\star = \mu_x - A^\star (H \mu_x + \mu_w).
\end{align}
If the random vector $(x,y)$ follows a normal distribution, then the best affine estimator is also optimal among all measurable estimators. In general, however, we do not know how much optimality is sacrificed by restricting attention to affine estimators. Moreover, the uncertainty about $\mathbb{P}_x$ and $\mathbb{P}_w$ transpires through to their first- and second-order moments. As the coefficients~\eqref{eq:optimal:affine} tend to be highly sensitive to these moments, their uncertainty remains worrying.

The minimax approach models the statistician's prior knowledge concerning $x$ via a convex closed uncertainty set $\mathbb{X}\subseteq \mathbb{R}^n$ as commonly used in robust optimization. The minimax MSE estimation problem is then formulated as a zero-sum game between the statistician, who selects the estimator $\psi\in{\mc F}$ with the goal to minimize the MSE, and nature, who chooses the parameter value $x\in\mathbb{X}$ with the goal to maximize the MSE.\begin{subequations}
	\begin{align}
		\label{eq:minimax}
		\mathop{\text{minimize}}_{\psi\in{\mc F}} ~ \max_{x \in \mathbb{X}} ~ R(\psi, x)
	\end{align}
	By construction, any minimizer $\psi_{\mathcal{M}}^\star$ of~\eqref{eq:minimax} incurs the smallest possible estimation error under the worst parameter realization within the uncertainty set $\mathbb{X}$. For this reason $\psi_{\mathcal{M}}^\star$ is called a minimax estimator. Note that the MSE~$R(\psi, x)$ generically displays a complicated non-concave dependence on~$x$ for any fixed~$\psi$, which implies that nature's inner maximization problem in~\eqref{eq:minimax} is usually non-convex. Thus, we should not expect the zero-sum game~\eqref{eq:minimax} between the statistician and nature to admit a Nash equilibrium. However, the inner maximization problem can be convexified by allowing nature to play mixed (randomized) strategies, that is, by reformulating~\eqref{eq:minimax} as the (equivalent) convex-concave saddle point problem
	\begin{align}
		\label{eq:minimax-saddle}
		\mathop{\text{minimize}}_{\psi\in{\mc F}} \max_{\mathbb{Q}_x \in \mathcal{M}(\mathbb{X})} \mathds{E}_{\mathbb{Q}_x} \left[ R(\psi, x) \right] ,
	\end{align}
\end{subequations}
where $\mathcal{M}(\mathbb{X})$ stands for the family of all distributions supported on $\mathbb{X}$ with finite second-order moments. As $\mathds{E}_{\mathbb{Q}_x} [ R(\psi, x)]$ is convex in $\psi$ for any fixed $\mathbb{Q}_x$ and concave (linear) in $\mathbb{Q}_x$ for any fixed $\psi$, while ${\mc F}$ and $\mathcal{M}(\mathbb{X})$ are both convex sets, the zero-sum game~\eqref{eq:minimax-saddle} admits a Nash equilibrium $(\psi_{\mathcal{M}}^\star, \mathbb{Q}_x^\star)$ under mild technical conditions. Note that $\psi_{\mathcal{M}}^\star$ is again a minimax estimator. Moreover, $\psi_{\mathcal{M}}^\star$ is the statistician's best response to nature's choice $\mathbb{Q}_x^\star$ and vice versa. Using the terminology introduced above, this means that $\psi_{\mathcal{M}}^\star$ is the Bayesian MMSE estimator corresponding to the prior~$\mathbb{Q}_x^\star$. For this reason, $\mathbb{Q}_x^\star$ is usually referred to as the {\em least favorable prior}. Even though the minimax approach exonerates the statistician from narrowing down her beliefs to a single prior distribution $\mathbb{Q}_x$, it still requires precise information about $\mathbb{P}_w$, which may not be available in practice. On the other hand, as it robustifies the estimator against {\em any} distribution on $\mathbb{X}$, the minimax approach is often regarded as overly pessimistic. Moreover, as in the case of the Bayesian MMSE estimation problem, $\psi^\star_{\mathcal{M}}$ may generically have a complicated functional form, and evaluating $\psi^\star_{\mathcal{M}}(y)$ to high precision may be computationally challenging if the dimension of $x$ is high. A simple remedy to mitigate these computational challenges would be to restrict ${\mc F}$ to the family $\mathcal A$ of affine estimators. The loss of optimality incurred by this approximation for different choices of $\mathbb{X}$ is discussed in~\cite[\S~4]{ref:juditsky2018lectures} and the references therein.

In this paper we bridge the Bayesian and the minimax approaches by leveraging tools from distributionally robust optimization. Specifically, we study distributionally robust estimation problems of the form
\begin{align}
	\label{eq:dro_estimator}
	\mathop{\text{minimize}}_{\psi\in{\mc F}} \max_{\mathbb{Q}_x \in \mathcal Q_x} \mathds{E}_{\mathbb{Q}_x} \left[ R(\psi, x) \right] ,
\end{align}
where $\mathcal Q_x\subseteq \mathcal{M}(\mathbb{R}^n)$ is an {\em ambiguity set} of multiple (possibly infinitely many) plausible prior distributions of~$x$. Note that if the ambiguity set collapses to the singleton $\mathcal Q_x = \{\mathbb{P}_{x}\}$ for some $\mathbb{P}_x\in\mathcal{M}(\mathbb{R}^n)$, then the distributionally robust estimation problem~\eqref{eq:dro_estimator} reduces to the Bayesian MMSE estimation problem~\eqref{eq:Bayes}. Similarly, under the ambiguity set $\mathcal Q_x=\mathcal{M}(\mathbb{X})$ for some convex closed uncertainty set $\mathbb{X}\subseteq \mathbb{R}^n$, problem~\eqref{eq:dro_estimator} reduces to the minimax mean square error estimation problem~\eqref{eq:minimax-saddle}. By providing considerable freedom in tailoring the ambiguity set $\mathcal Q_x$, the distributionally robust approach thus allows the statistician to reconcile the specificity of the Bayesian approach with the conservativeness of the minimax approach.

The estimation model~\eqref{eq:dro_estimator} still relies on the premise that the noise distribution~$\mathbb{P}_w$ is precisely known, and this assumption is not tenable in practice. However, nothing prevents us from further robustifying~\eqref{eq:dro_estimator} against uncertainty in $\mathbb{P}_w$. To this end, we define $\mathcal{M}(\mathbb{R}^{n+m})$ as the family of all joint distributions of $x$ and $w$ with finite second-order moments. Moreover, we define the {\em average risk} $\mathcal{R}:{\mc F} \times \mathcal{M}(\mathbb{R}^{n+m}) \rightarrow \mathbb{R}$ through
\[
\mathcal{R}(\psi, \mathbb{P}) = \mathds{E}_{\mathbb{P}} [ \| x - \psi(H x + w) \|^2 ].
\]
If $\mathbb{P}=\mathbb{P}_x\times\mathbb{P}_w$ for some marginal distributions $\mathbb{P}_x\in\mathcal{M}(\mathbb{R}^n)$ and $\mathbb{P}_w\in\mathcal{M}(\mathbb{R}^m)$, which implies that $x$ and $w$ are independent under $\mathbb{P}$, and if $\mathbb{P}_{y|x}$ is defined as $\mathbb{P}_w$ shifted by $Hx$, then $\mathcal{R}(\psi, \mathbb{P})=  \mathds{E}_{\mathbb{P}_x} [ R(\psi, x) ]$. Thus, the average risk $\mathcal{R}(\psi, \mathbb{P})$ corresponds indeed to the risk $R(\psi, x)$ averaged under the marginal distribution~$\mathbb{P}_x$. In the remainder of this paper we will study generalized distributionally robust estimation problems of the form
\begin{align}
	\label{eq:dro}
	\mathop{\text{minimize}}_{\psi\in{\mc F}} \sup\limits_{\mathbb{Q} \in \mathbb B(\wh \PP)} \mathcal{R}(\psi, \mathbb{Q}) ,
\end{align}
where the ambiguity set $\mathbb B(\wh \PP)\subseteq \mathcal{M}(\mathbb{R}^{n+m})$ captures distributional uncertainty in both $\mathbb{P}_x$ and $\mathbb{P}_w$. Specifically, we will model $\mathbb B(\wh \PP)$ as a set of factorizable distributions $\mathbb{Q}=\mathbb{Q}_x\times \mathbb{Q}_w$ close to a nominal distribution $\wh \PP = \wh \PP_x \times \wh \PP_w$ in the sense that $\mathbb{Q}_x$ and $\mathbb{Q}_w$ are close to $\wh \PP_x$ and $\wh \PP_w$ in Wasserstein distance, respectively.

\begin{definition}[Wasserstein distance]
	\label{definition:wasserstein}
	For any $d\in\mathbb N$, the type-2 Wasserstein distance between two distributions $\mathbb{Q}_1,\mathbb{Q}_2\in\mathcal{M}(\mathbb{R}^d)$ is defined as
	\begin{equation}
	\notag
	\mathds{W}(\mathbb{Q}_1, \mathbb{Q}_2) = \inf\limits_{\pi\in\Pi(\mathbb{Q}_1,\mathbb{Q}_2)}  \left( \int_{\mathbb{R}^d \times \mathbb{R}^d} \norm{\xi_1 - \xi_2}^2\, \pi({\rm d}\xi_1, {\rm d} \xi_2) \right)^{\frac{1}{2}},
	\end{equation}
	where $\Pi(\mathbb{Q}_1,\mathbb{Q}_2)$ denotes the set of all joint distributions or couplings $\pi\in \mathcal{M}(\mathbb{R}^d\times \mathbb{R}^d)$ of the random vectors $\xi_1\in\mathbb{R}^d$ and $\xi_2\in\mathbb{R}^d$ with marginal distributions $\mathbb{Q}_1$ and  $\mathbb{Q}_2$, respectively.
\end{definition}

The dependence of the Wasserstein distance on $d$ is notationally suppressed to avoid clutter. Note that $\mathds{W}(\mathbb{Q}_1, \mathbb{Q}_2)^2$ is naturally interpreted as the optimal value of a transportation problem that determines the minimum cost of moving the distribution $\mathbb{Q}_1$ to $\mathbb{Q}_2$, where the cost of moving a unit probability mass from $\xi_1$ to $\xi_2$ is given by the squared Euclidean distance $\| \xi_1 - \xi_2\|^2$. For this reason, the optimization variable $\pi$ is sometimes referred to as a transportation plan and the Wasserstein distance as the earth mover's distance.

Formally, we define the {\em Wasserstein ambiguity set} as
\begin{align}\label{eq:Ambi}
	\begin{aligned}
		\mathbb B(\wh \PP) = \left\{ \mathbb{Q}_x \times \mathbb{Q}_w:
		\def1.2{1.2}
		\begin{array}{l@{\;}ll@{\;}l}
			\mathbb{Q}_x &\in \mathcal{M}(\mathbb{R}^n),& \, \mathds{W}(\mathbb{Q}_x, \wh \PP_x) &\leq \rho_x \\
			\mathbb{Q}_w &\in \mathcal{M}(\mathbb{R}^m),& \, \mathds{W}(\mathbb{Q}_w, \wh \PP_w) &\leq \rho_w
		\end{array}
		\right\},
	\end{aligned}
\end{align}
where $\wh \PP_x$ and $\wh \PP_w$ represent prescribed nominal distributions that could be constructed via statistical analysis or expert judgement, while the Wasserstein radii $\rho_x\ge 0$ and $\rho_w\ge 0$ constitute hyperparameters that quantify the statistician's uncertainty about the nominal distributions of $x$ and $w$. We emphasize that the distributionally robust estimation model~\eqref{eq:dro} generalizes all preceding models. Indeed, if $\rho_w=0$, then~\eqref{eq:dro} reduces to the first distributionally robust model~\eqref{eq:dro_estimator}, which in turn encompasses both the MMSE estimation problem~\eqref{eq:Bayes} (for $\rho_x=0$) and the minimax estimation problem~\eqref{eq:minimax-saddle} (for $\rho_x=\infty$) as special cases.

The distributionally robust estimation model~\eqref{eq:dro} is conceptually attractive because the hyperparameters~$\rho_x$ and~$\rho_w$ allow the statistician to specify her level of trust in the nominal prior distribution~$\wh \PP_x$ and the nominal noise distribution~$\wh \PP_w$. In the remainder of the paper we will show that~\eqref{eq:dro} is also computationally attractive. This is maybe surprising because mixtures of factorizable distributions are generally not factorizable, which implies that the Wasserstein ambiguity set $\mathbb B(\wh \PP)$ is non-convex.

We remark that one could also work with an alternative ambiguity set of the form
\begin{align}
	\label{eq:Ambi'}
	\mathbb B'(\wh \PP) = \left\{ \mathbb{Q}_x \times \mathbb{Q}_w: \mathbb{Q}_x \in \mathcal{M}(\mathbb{R}^n),~ \mathbb{Q}_w \in \mathcal{M}(\mathbb{R}^m), ~\mathds{W}(\mathbb{Q}_x\times \mathbb{Q}_w, \wh \PP_x\times \wh \PP_w) \leq \rho \right\},
\end{align}
which involves only a single hyperparameter $\rho\ge 0$ and is therefore less expressive but maybe easier to calibrate than~$\mathbb B(\wh \PP)$. The following lemma is instrumental to understanding the relation between~$\mathbb B(\wh \PP)$ and~$\mathbb B'(\wh \PP)$. The proof of this result is relegated to the appendix.

\begin{lemma}[Pythagoras' theorem for Wasserstein distances]
	\label{lem:wasserstein-decomposition}
	For any $\mathbb{Q}_x^1,\mathbb{Q}_x^2\in\mathcal{M}(\mathbb{R}^n)$ and $\mathbb{Q}_w^1,\mathbb{Q}_w^2\in\mathcal{M}(\mathbb{R}^m)$ we have $\mathds{W}(\mathbb{Q}_x^1\times \mathbb{Q}_w^1, \mathbb{Q}_x^2\times \mathbb{Q}_w^2)^2 = \mathds{W}(\mathbb{Q}_x^1, \mathbb{Q}_x^2)^2 + \mathds{W}(\mathbb{Q}_w^1, \mathbb{Q}_w^2)^2$.
\end{lemma}

If we denote the ambiguity sets \eqref{eq:Ambi} and \eqref{eq:Ambi'} temporarily by $\mathbb B_{\rho_x,\rho_w} (\wh \PP)$ and $\mathbb B'_\rho(\wh \PP)$ in order to make their dependence on the hyperparameters explicit, then Lemma~\ref{lem:wasserstein-decomposition} implies that
\begin{align*}
	\mathbb B'_\rho(\wh \PP) =\bigcup_{\rho_x^2+\rho_w^2\le \rho^2} \mathbb B_{\rho_x,\rho_w} (\wh \PP).
\end{align*}
This relation suggests that $\mathbb B'_\rho(\wh \PP)$ could be substantially larger than $\mathbb B_{\rho_x,\rho_w} (\wh \PP)$ for any fixed $\rho,\rho_x,\rho_w\ge 0$ with $\rho_x^2+\rho_w^2= \rho^2$ and thus lead to substantially more conservative estimators.

In the following we summarize the key contributions of this paper.
\begin{enumerate}
	\item \label{contribution1} We construct a safe approximation for the distributionally robust MMSE estimation problem~\eqref{eq:dro} by restricting attention to {\em affine} estimators and by maximizing the average risk over an {\em outer} approximation of the Wasserstein ambiguity set, which is described through first- and second-order moment conditions. We then prove that this safe approximation is equivalent to a tractable convex program.
	\item \label{contribution2} We also study a {\em dual} estimation problem, which is obtained by interchanging the minimization and maximization operations in the {\em primal} problem~\eqref{eq:dro} and thus lower bounds the optimal value of~\eqref{eq:dro}. We then construct a safe approximation for this dual problem by restricting the Wasserstein ambiguity set to contain only {\em normal} distributions. Assuming that the nominal distribution is normal, we prove that this safe approximation is again equivalent to a tractable convex program.
	\item By construction, the primal and dual estimation problems are upper and lower bounded by their respective safe approximations. We prove, however, that the optimal values of the safe approximations collapse if the nominal distribution is normal. This result has three important implications.
	\begin{enumerate}
	    \item The primal and dual estimation problems and their safe approximations are all equivalent. This implies via contributions~(\ref{contribution1}) and~(\ref{contribution2}) that both original estimation problems are tractable.
	    \item The primal estimation problem is solved by an {\em affine} estimator, and the dual estimation problem is solved by a {\em normal} distribution. In other words, we have discovered a new class of adaptive distributionally robust optimization problems for which affine decision rules are optimal.
	    \item The affine estimator and the normal distribution that solve the primal and dual estimation problems, respectively, form a {\em Nash equilibrium} for the zero-sum game between the statistician and nature. Thus, the optimal normal distribution constitutes a {\em least favorable prior}, and the optimal affine estimator represents the corresponding {\em Bayesian MMSE estimator}.
	\end{enumerate}
	We leverage these insights to prove that the optimal affine estimator can be constructed easily from the least favorable prior without the need to solve another optimization problem.
	\item We argue that our main results remain valid if the nominal distribution is any elliptical distribution.
	\item We develop a tailor-made Frank-Wolfe algorithm that can solve the dual estimation problem orders of magnitude faster than state-of-the-art general purpose solvers. We show that this algorithm enjoys a linear convergence rate. Moreover, we prove that the direction-finding subproblems can be solved in quasi-closed form, which means that the algorithm offers a favorable iteration complexity.
\end{enumerate}

We highlight that the Wasserstein ambiguity set~\eqref{eq:Ambi} is non-convex as it contains only distributions under which the signal and the noise are independent. To our best knowledge, we describe the first distributionally robust optimization model with independence conditions that admits a tractable reformulation.
We also emphasize that some of our results hold only if the nominal distribution~$\wh \PP$ is normal or elliptical. While this is restrictive, there is strong evidence that normal distributions are natural candidates for~$\wh \PP$. One reason for this is that the normal distribution has maximum entropy among all distributions with prescribed first- and second-order moments~\cite[\S~12]{ref:cover2006elements}. Therefore, it has appeal as the least prejudiced baseline model. Similarly, if the parameter~$x$ in~\eqref{eq:linear:observation} is normally distributed, then a normal distribution minimizes the mutual information between~$x$ and the observation~$y$ among all noise distributions with bounded variance~\cite[Lemma~II.2]{ref:diggavi2001worst}. In this sense, normally distributed noise renders the observations least informative. Conversely, if the noise in~\eqref{eq:linear:observation} is normally distributed, then a normal distribution maximizes the MMSE across all distributions of~$x$ with bounded variance~\cite[Proposition~15]{ref:guo2011estimation}. In this sense, normally distributed parameters are the hardest to estimate. Using normal nominal distributions thus amounts to adopting a worst-case perspective.

In the following we briefly survey the landscape of existing MMSE estimation models that have a robustness flavor. Several authors have addressed the {\em minimax} MMSE estimation problem~\eqref{eq:minimax} from the perspective of classical robust optimization \cite{ref:beck2007regularization, ref:beck2007mean, ref:eldar2008minimax, ref:eldar2004linear, ref:eldar2004robust, ref:juditsky2018nearoptimality}. To guarantee computational tractability, in all of these papers the estimators are restricted to be affine functions of the measurements. In this case, the minimax MMSE estimation problem can be reformulated as a tractable semidefinite program (SDP) if the uncertainty set $\mathbb{X}$ is an ellipsoid~\cite{ref:eldar2004linear, ref:eldar2004robust} or an intersection of two ellipsoids~\cite{ref:beck2007regularization}. Similar SDP reformulations are available if the observation matrix $H$ is also subject to uncertainty and ranges over a spectral norm ball~\cite{ref:eldar2004robust} or displays a block circulant structure, with each block ranging over a Frobenius norm ball \cite{ref:beck2007mean}. If the uncertainty set is described by an intersection of several ellipsoids, then the minimax MMSE estimation problem admits an (inexact) SDP relaxation~\cite{ref:eldar2008minimax}. Even though the restriction to affine estimators may incur a loss of optimality, affine estimators are known to be near-optimal in all of the above minimax estimation models~\cite{ref:juditsky2018nearoptimality}.

Another stream of literature investigates the {\em distributionally robust} estimation model~\eqref{eq:dro_estimator} with an ambiguous signal distribution and a crisp noise distribution. When focusing on affine estimators only, this model can be reformulated as a tractable SDP if the uncertainty in the signal distribution is characterized through spectral constraints on its covariance matrix~\cite{ref:eldar2004competitive}. This tractability result also extends to uncertain observation matrices. Similar SDP reformulations are available for the distributionally robust estimation model~\eqref{eq:dro} when both the signal and the noise distribution are ambiguous and their covariance matrices are subject to spectral constraints~\cite{ref:eldar2006robust}. Extensions to uncertain block circulant observation matrices are discussed in~\cite{ref:beck2006robust}.

Some authors have studied less structured distributionally robust estimation problems where the signal~$x$ and the measurement~$y$ are governed by a distribution that may {\em not} obey the linear measurement model~\eqref{eq:linear:observation}. In this case, the zero-sum game between the statistician and nature admits a Nash equilibrium if nature may choose any distribution that has a bounded Kullback-Leibler divergence with respect to a {\em normal} nominal distribution~\cite{ref:levy2004robust}. Intriguingly, the (affine) Bayesian MMSE estimator for the nominal distribution is optimal in this model and thus enjoys strong robustness properties. On the downside, there is no hope to improve this estimator's performance by tuning the size of the Kullback-Leibler ambiguity set. The underlying distributionally robust estimation model also serves as a fundamental building block for a robust Kalman filter~\cite{ref:levy2012robust}. Extensions to general $\tau$-divergence ambiguity sets that contain only normal distributions are reported in~\cite{ref:zorzi2017robustness}. We emphasize that all papers surveyed so far merely derive SDP reformulations or SDP relaxations that can be addressed with general purpose solvers, but none of them develops a customized solution algorithm.

The present paper extends the distributionally robust MMSE estimation model introduced in~\cite{ref:shafieezadeh2018wasserstein}, which accommodates a simple Wasserstein ambiguity set for the distribution of the signal-measurement pairs and makes no structural assumptions about the measurement noise. Note, however, that the linear measurement model~\eqref{eq:linear:observation} abounds in the literature on control theory, signal processing and information theory, implying that there are numerous applications where the measurement noise is {\em known} to be additive and independent of the signal. As we will see in Section~\ref{sect:numerical}, ignoring this structural information may result in weak estimators that sacrifice predictive performance. In Sections~\ref{sect:approx}--\ref{sect:nash} we will further see that constructing an explicit Nash equilibrium is considerably more difficult in the presence of structural information. Finally, we describe here an accelerated Frank-Wolfe algorithm that improves the sublinear convergence rate established in~\cite{ref:shafieezadeh2018wasserstein} to a linear rate. We emphasize that, in contrast to the robust MMSE estimators derived in~\cite{ref:levy2004robust, ref:zorzi2017robustness} that are insensitive to the radii of the underlying divergence-based ambiguity sets, the estimators constructed here change with the Wasserstein radii~$\rho_x$ and~$\rho_w$. Thus, using a Wasserstein ambiguity set to robustify the nominal MMSE estimation problem has a regularizing effect and leads to a parametric family of estimators that can be tuned to attain maximum prediction accuracy. Similar connections between robustification and regularization have previously been discovered in the context of statistical learning~\cite{ref:shafieezadeh2019mass-transportation} and covariance estimation~\cite{ref:nguyen2018distributionally}.

The paper is structured as follows. Sections~\ref{sect:approx} and~\ref{sect:dual} develop conservative approximations for the primal and dual distributionally robust MMSE estimation problems, respectively, both of which are equivalent to tractable convex programs.  Section~\ref{sect:nash} shows that if the nominal distribution is normal, then both approximations are exact and can be used to find a Nash equilibrium for the zero-sum game between the statistician and nature. Extensions to non-normal nominal distributions are discussed in Section~\ref{sect:elliptical}. Section~\ref{sect:algorithm} develops an efficient Frank-Wolfe algorithm for the dual MMSE estimation problem, and Section~\ref{sect:numerical} reports on numerical results.

\textbf{Notation.} For any $A \in \mathbb{R}^{d \times d}$ we use $\Tr{A}$ to denote the trace and $\|A\|=\sqrt{\Tr{A^\top A}}$ to denote the Frobenius norm of $A$. By slight abuse of notation, the Euclidean norm of $v\in\mathbb{R}^d$ is also denoted by $\|v\|$. Moreover, $I_d$ stands for the identity matrix in $\mathbb{R}^{d \times d}$. For any $A,B\in\mathbb{R}^{d\times d}$, we use $\inner{A}{B} = \Tr{A^\top B}$ to denote the inner product and $A\otimes B$ to denote the Kronecker product of $A$ and $B$. The space of all symmetric matrices in $\mathbb{R}^{d\times d}$ is denoted by $\mathbb S^d$. We use $\mathbb{S}_{+}^d$ ($\mathbb{S}_{++}^d$) to represent the cone of symmetric positive semidefinite (positive definite) matrices in $\mathbb S^d$. For any $A,B\in\mathbb S^d$, the relation $A\succeq B$ ($A\succ B$) means that $A-B\in\mathbb{S}_{+}^d$ ($A-B\in \mathbb{S}_{++}^d$). The unique positive semidefinite square root of a matrix $A \in \mathbb{S}_{+}^d$ is denoted by $A^\frac{1}{2}$. For any $A \in \mathbb S^d$, $\lambda_{\min}(A)$ and $\lambda_{\max}(A)$ denote the minimum and maximum eigenvalues of $A$, respectively.

\section{The Gelbrich MMSE Estimation Problem}
\label{sect:approx}

The distributionally robust estimation problem~\eqref{eq:dro} poses two fundamental challenges. First, checking feasibility of the inner maximization problem in~\eqref{eq:dro} requires computing the Wasserstein distances $\mathds{W}(\wh \PP_x,\mathbb{Q}_x)$ and $\mathds{W}(\wh \PP_w,\mathbb{Q}_w)$, which is \#P-hard even if~$\wh \PP_x$ and~$\wh \PP_w$ are simple two-point distributions, while~$\mathbb{Q}_x$ and~$\mathbb{Q}_w$ are uniform distributions on hypercubes \cite{ref:taskesen2019complexity}. Efficient algorithms for computing Wasserstein distances are available only if both involved distributions are discrete~\cite{ref:cuturi2013sinkhorn, ref:peyre2019computational, ref:solomon2015convolutional}, and analytical formulas are only known in exceptional cases ({\em e.g.}, if both distributions are Gaussian~\cite{ref:givens1984class} or belong to the same family of elliptical distributions~\cite{ref:gelbrich1990formula}). The second challenge is that the outer minimization problem in~\eqref{eq:dro} constitutes an infinite-dimensional functional optimization problem. In order to bypass these computational challenges, we first seek a conservative approximation for~\eqref{eq:dro} by relaxing the ambiguity set~$\mathbb B(\wh \PP)$ and restricting the feasible set~${\mc F}$. We begin by constructing an outer approximation for the ambiguity set. To this end, we introduce a new distance measure on the space of mean vectors and covariance matrices.
\begin{definition}[Gelbrich distance]
	\label{def:induced_distance}
	For any $d\in\mathbb N$, the Gelbrich distance between two tuples of mean vectors and covariance matrices $(\mu_1, \Sigma_1),(\mu_2, \Sigma_2) \in \mathbb{R}^d \times \mathbb{S}_{+}^d$ is defined as
	\[
	\mathds{G} \big( (\mu_1, \Sigma_1), (\mu_2, \Sigma_2) \big) = \sqrt{\norm{\mu_1 - \mu_2}^2 + \Tr{\Sigma_1 + \Sigma_2 - 2 \left( \Sigma_2^\frac{1}{2} \Sigma_1 \Sigma_2^\frac{1}{2} \right)^\frac{1}{2}}}.
	\]

\end{definition}

The dependence of the Gelbrich distance on $d$ is notationally suppressed in order to avoid clutter. One can show that $\mathds{G}$ constitutes a metric on $\mathbb{R}^d \times \mathbb{S}_{+}^d$, that is, $\mathds{G}$ is symmetric, non-negative, vanishes if and only if $(\mu_1,\Sigma_1)=(\mu_2,\Sigma_2)$ and satisfies the triangle inequality~\cite[pp.~239]{ref:givens1984class}.

\begin{proposition}[{Commuting covariance matrices \cite[p.~239]{ref:givens1984class}}]
	\label{prop:commuting}
	If $\mu_1,\mu_2\in\mathbb{R}^d$ are identical and $\Sigma_1, \Sigma_2\in\mathbb{S}_{+}^d$ commute ($\Sigma_1 \Sigma_2= \Sigma_2 \Sigma_1$), then the Gelbrich distance simplifies to $\mathds{G}\big( (\mu_1, \Sigma_1), (\mu_2, \Sigma_2)\big) = \norm{\sqrt{\Sigma_1} - \sqrt{\Sigma_2}}$.
\end{proposition}

While the Gelbrich distance is non-convex, the squared Gelbrich distance is convex in all of its arguments.

\begin{proposition}[Convexity and continuity of the squared Gelbrich distance]
    \label{prop:gelbrich-convexity}
    The squared Gelbrich distance $\mathds{G}\big( (\mu_1, \Sigma_1), (\mu_2, \Sigma_2)\big)^2$ is jointly convex and continuous in $\mu_1 , \mu_2 \in \mathbb{R}^d$ and $\Sigma_1, \Sigma_2 \in \mathbb{S}_{+}^d$.
\end{proposition}
\begin{proof}[Proof of Proposition~\ref{prop:gelbrich-convexity}]
    By \cite[Proposition~2]{ref:malago2018wasserstein}, the squared Gelbrich distance $\mathds{G}\big( (\mu_1, \Sigma_1), (\mu_2, \Sigma_2)\big)^2$ coincides with the optimal value of the semidefinite program
    \begin{equation*}
        \begin{array}{cl}
             \min & \| \mu_1 - \mu_2 \|_2^2 + \Tr{\Sigma_1 + \Sigma_2 - 2 C} \\[1ex]
             \st & C \in \mbb R^{d \times d} \\[1ex]
             & \begin{bmatrix} \Sigma_1 & C \\ C^\top & \Sigma_2 \end{bmatrix} \succeq 0,
        \end{array}
    \end{equation*}
    see also \cite[Section~3]{ref:dowson1982frechet}. Less general results that hold when one of the matrices $\Sigma_1$ or $\Sigma_2$ is positive definite are reported in~\cite{ref:givens1984class, ref:knott1984optimal, ref:olkin1982distance}.
    The convexity of the squared Gelbrich distance then follows from~\cite[Proposition~3.3.1]{ref:bertsekas2009convex}, which guarantees that convexity is preserved under partial minimization. Moreover, the continuity of the squared Gelbrich distance follows from the continuity of the matrix square root established in Lemma~\ref{lemma:hoelder}.

\end{proof}

Our interest in the Gelbrich distance stems mainly from the next proposition, which lower bounds the Wasserstein distance between two distributions in terms of their first- and second-order moments. We will later see that this bound becomes tight when $\mathbb{Q}_1$ and $\mathbb{Q}_2$ are normal or---more generally---elliptical distributions of the same type.

\begin{proposition}[Moment bound on the Wasserstein distance {\cite[Theorem~2.1]{ref:gelbrich1990formula}}]
	\label{prop:wasserstein}
	For any distributions $\mathbb{Q}_1$,~$\mathbb{Q}_2 \in \mathcal{M}(\mathbb{R}^d)$ with mean vectors $\mu_1$, $\mu_2 \in \mathbb{R}^d$ and covariance matrices $\Sigma_1$, $\Sigma_2 \in \mathbb{S}_{+}^d$, respectively, we have
	\begin{align} \notag
		\mathds{W}(\mathbb{Q}_1, \mathbb{Q}_2) \geq \mathds{G} \big( (\mu_1, \Sigma_1), (\mu_2, \Sigma_2) \big).
	\end{align}
\end{proposition}

Proposition~\ref{prop:wasserstein} prompts us to construct an outer approximation for the Wasserstein ambiguity set $\mathbb B(\wh \PP)$ by using the Gelbrich distance. Specifically, we define the {\em Gelbrich ambiguity set} centered at $\wh \PP=\wh \PP_x\times \wh \PP_w$ as
\begin{align*}
	\begin{aligned}
		{\mbb G}(\wh \PP) = \left\{ \mathbb{Q}_x \times \mathbb{Q}_w:
		\def1.2{1.2}
		\begin{array}{l}
			\begin{array}{l@{\;}@{\;}l@{\;}l}
				\mathbb{Q}_x \in \mathcal{M}(\mathbb{R}^n), & \mu_x\,= \mathds{E}_{\mathbb{Q}_x}[x], & \Sigma_x\,= \mathds{E}_{\mathbb{Q}_x} [x x^\top]  -\mu_x \mu_x^\top \\
				\mathbb{Q}_w \in \mathcal{M}(\mathbb{R}^m), & \mu_w=\mathds{E}_{\mathbb{Q}_w}[w], & \Sigma_w = \mathds{E}_{\mathbb{Q}_w} [w w^\top] - \mu_w \mu_w^\top
			\end{array}\\
			\begin{array}{ll}
				\mathds{G} \big( (\mu_x, \Sigma_x), (\wh{\m}_x, \covsa_x) \big)\leq \rho_x, & \mathds{G} \big( (\mu_w, \Sigma_w), (\wh{\m}_w, \covsa_w) \big) \leq \rho_w
			\end{array}
		\end{array}
		\right\},
	\end{aligned}
\end{align*}
where $\wh{\m}_x$ and $\wh{\m}_w$ denote the mean vectors and $\covsa_x$ and $\covsa_w$ the covariance matrices of $\wh \PP_x$ and $\wh \PP_w$, respectively.

\begin{corollary}[Relation between Gelbrich and Wasserstein ambiguity sets]
	\label{cor:wasserstein-in-gelbrich}
	For any $\wh \PP =\wh \PP_x\times \wh \PP_w$ with $\wh \PP_x\in  \mathcal{M}(\mathbb{R}^n)$ and $\wh \PP_w\in  \mathcal{M}(\mathbb{R}^m)$ we have $\mathbb B(\wh \PP) \subseteq {\mbb G}(\wh \PP)$.
\end{corollary}
\begin{proof}[Proof of Corollary~\ref{cor:wasserstein-in-gelbrich}]
	Select any $\mathbb{Q}=\mathbb{Q}_x\times\mathbb{Q}_w\in\mathbb B(\wh \PP)$ and define $\mu_x$ and $\mu_w$ as the mean vectors and $\Sigma_x$ and $\Sigma_w$ as the covariance matrices of $\mathbb{Q}_x$ and $\mathbb{Q}_w$, respectively. By Proposition~\ref{prop:wasserstein} we then have
	\[
	\mathds{G} \big( (\mu_x, \Sigma_x), (\wh{\m}_x, \covsa_x) \big) \le \mathds{W}(\mathbb{Q}_x, \wh \PP_x)\le \rho_x\quad\text{and}\quad
	\mathds{G} \big( (\mu_w, \Sigma_w), (\wh{\m}_w, \covsa_w) \big) \le \mathds{W}(\mathbb{Q}_w, \wh \PP_w)\le \rho_w,
	\]
	which in turn implies that $\mathbb{Q}\in{\mbb G}(\wh \PP)$. We may thus conclude that $\mathbb B(\wh \PP) \subseteq {\mbb G}(\wh \PP)$.
\end{proof}
By restricting ${\mc F}$ to the set $\Ac$ of all affine estimators while relaxing $\mathbb B(\wh \PP)$ to the Gelbrich ambiguity set ${\mbb G}(\wh \PP)$, we obtain the following conservative approximation of the distributionally robust estimation problem~\eqref{eq:dro}.
\begin{align}
	\label{eq:dro:approx}
	\mathop{\text{minimize}}_{\psi\in\Ac} \sup\limits_{\mathbb{Q} \in {\mbb G}(\wh \PP)} \mathcal{R}(\psi, \mathbb{Q})
\end{align}
From now on we will call~\eqref{eq:dro} and~\eqref{eq:dro:approx} the Wasserstein and Gelbrich MMSE estimation problems, and we will refer to their minimizers as Wasserstein and Gelbrich MMSE estimators, respectively. As the average risk $\mathcal{R}(\psi, \mathbb{Q})$ of a fixed affine estimator $\psi\in\Ac$ is convex and quadratic in the mean vector $\mu$ and affine in the covariance matrix $\Sigma$ of the distribution~$\mathbb{Q}$, the inner maximization problem in~\eqref{eq:dro:approx} is non-convex. Thus, one might suspect that the Gelbrich MMSE estimation problem is intractable. Below we will show, however, that~\eqref{eq:dro:approx} is equivalent to a finite convex program that can be solved in polynomial time. To this end, we first show that, under mild conditions, problem~\eqref{eq:dro:approx} is stable with respect to changes of its input parameters.

\begin{proposition}[Regularity of the Gelbrich MMSE estimation problem]
	\label{prop:Gelbrich}
	The Gelbrich MMSE estimation problem~\eqref{eq:dro:approx} enjoys the following regularity properties.
	\begin{enumerate}[label = $(\roman*)$]
		\item {\bf Conservativeness:} \label{prop:Gelbrich:conservative}
		Problem~\eqref{eq:dro:approx} upper bounds the Wasserstein MMSE estimation problem~\eqref{eq:dro}.

		\item {\bf Solvability:} \label{prop:Gelbrich:solvability}
		If $\rho_x>0$ and $\rho_w>0$, then the minimum of~\eqref{eq:dro:approx} is attained.

		\item {\bf Stability:} \label{prop:Gelbrich:cont}
		If $\rho_x>0$ and $\rho_w>0$, then the minimum of~\eqref{eq:dro:approx} is continuous in $(\wh{\m}_x,\wh{\m}_w,\covsa_x,\covsa_w)$.

	\end{enumerate}
\end{proposition}

The proof of Proposition~\ref{prop:Gelbrich} is lengthy and technical and is therefore relegated to the appendix. We are now ready to prove the main result of this section.

\begin{theorem}[Gelbrich MMSE estimation problem]
	\label{thm:conservative}
	The Gelbrich MMSE estimation problem~\eqref{eq:dro:approx} is equivalent to the finite convex optimization problem
	\begin{equation}
	\label{eq:VA}
	\begin{array}{cl}
		\inf
		& \dualvar_x \big( \rho_x^2 - \Tr{\covsa_x} \big) + \dualvar_x^2 \inner{[\dualvar_x I_n - (I_n- A H)^\top (I_n - AH)]^{-1}}{\covsa_x} \\
		& \hspace{1cm} + \dualvar_w \big( \rho_w^2 - \Tr{\covsa_w} \big) + \dualvar_w^2 \inner{(\dualvar_w I_m - A^\top A)^{-1}}{\covsa_w} \\[1ex]
		\st & A \in \mathbb{R}^{n \times m},\quad \dualvar_x, \dualvar_w \in \mathbb{R}_+ \\
		& \dualvar_x I_n - (I_n-A H)^\top (I_n-A H) \succ 0, \quad \dualvar_w I_m - A^\top A \succ 0.
	\end{array}
	\end{equation}
	Moreover, if $\rho_x> 0$ and $\rho_w > 0$, then~\eqref{eq:VA} admits an optimal solution\footnote{We say that $A^\star$ solves~\eqref{eq:VA} if adding the constraint $A=A^\star$ does not change the infimum of~\eqref{eq:VA}. Note that the infimum of the resulting problem over $(\dualvar_x,\dualvar_w)$ may not be attained, {\em i.e.}, the existence of a solution $A^\star$ does not imply that~\eqref{eq:VA} is solvable.}~$A^\star$, and the infimum of~\eqref{eq:dro:approx} is attained by the affine estimator $\psi^\star(y) =  A^\star y +b^\star$, where $b^\star = \wh{\m}_x - A^\star (H \wh{\m}_x + \wh{\m}_w)$.
\end{theorem}

{The \textit{strict} semidefinite inequalities in~\eqref{eq:VA} ensure that the inverse matrices in the objective function exist.}

\begin{proof}[Proof of Theorem~\ref{thm:conservative}]
	Throughout this proof we denote by $\psi_{A,b}\in\Ac$ the affine estimator $\psi_{A,b}(y)=Ay+b$ corresponding to the sensitivity matrix $A \in \mathbb{R}^{n \times m}$ and the vector $b\in\mathbb{R}^n$ of intercepts. In the following we fix some $A \in \mathbb{R}^{n \times m}$ and define $K = I_n-A H$ in order to simplify the notation. By the definitions of the average risk $\mathcal{R}(\psi, \mathbb{Q})$ and the Gelbrich ambiguity set ${\mbb G}(\widehat \mathbb{P})$, we then have
	\begin{equation}
	\label{eq:robust}
	\inf_{\substack{b}} \sup\limits_{\mathbb{Q} \in {\mbb G}(\widehat \mathbb{P})} \mathcal{R}(\psi_{A,b}, \mathbb{Q}) = \left\{
	\begin{array}{ccl}
		\inf\limits_{b}& \sup\limits_{\substack{\mu_x, \mu_w \\ \Sigma_x, \Sigma_w\succeq 0}} & \inner{K^\top K}{\Sigma_x + \mu_x \mu_x^\top} + \inner{A^\top A}{\Sigma_w + \mu_w \mu_w^\top} + b^\top b \\[-2ex]
		&& \hspace{1cm} - 2 \mu_x^\top K^\top A \mu_w - 2 b^\top ( K \mu_x - A \mu_w ) \\[1ex]
		&\st & \mathds{G} \big( (\mu_x, \Sigma_x), (\wh{\m}_x, \covsa_x) \big)^2 \leq \rho_x^2 \\
		&& \mathds{G} \big( (\mu_w, \Sigma_w), (\wh{\m}_w, \covsa_w) \big)^2 \leq \rho_w^2.
	\end{array}	\right.
	\end{equation}
	The outer minimization problem in~\eqref{eq:robust} is convex because the objective function of the minimax problem is convex in $b$ for any fixed $(\mu_x, \mu_w, \Sigma_x, \Sigma_w)$ and because convexity is preserved under maximization. Moreover, the inner maximization problem in~\eqref{eq:robust} is non-convex because its objective function is convex in~$(\mu_x, \mu_w)$. This observation prompts us to maximize over $(\mu_x,\mu_w)$ and $(\Sigma_x,\Sigma_w)$ sequentially and to reformulate~\eqref{eq:robust}~as
	\begin{equation}
	\label{eq:robust2}
	\begin{array}{cccl}
		\inf\limits_{b}& \sup\limits_{\substack{\mu_x, \mu_w\\ \| \mu_x - \wh{\m}_x \| \leq \rho_x \\\| \mu_w - \wh{\m}_w \| \leq \rho_w}}& \sup\limits_{\substack{ \Sigma_x, \Sigma_w\succeq 0}} & \inner{K^\top K}{\Sigma_x + \mu_x \mu_x^\top} + \inner{A^\top A}{\Sigma_w + \mu_w \mu_w^\top} + b^\top b \\[-4ex]
		&&& \hspace{0.2cm} - 2 \mu_x^\top K^\top A \mu_w - 2 b^\top ( K \mu_x - A \mu_w ) \\[1ex]
		&&\st & \mathds{G} \big( (\mu_x, \Sigma_x), (\wh{\m}_x, \covsa_x) \big)^2 \leq \rho_x^2 \\
		&&& \mathds{G} \big( (\mu_w, \Sigma_w), (\wh{\m}_w, \covsa_w) \big)^2 \leq \rho_w^2.
	\end{array}
	\end{equation}
	As $\| \mu_x - \wh{\m}_x \| \leq  \mathds{G} ( (\mu_x, \Sigma_x), (\wh{\m}_x, \covsa_x) )$ and as this inequality is tight for $\Sigma_x=\covsa_x$, the extra constraint $\| \mu_x - \wh{\m}_x \| \leq \rho_x$ is actually redundant and merely ensures that the maximization problem over $\Sigma_x$ remains feasible for any admissible choice of $\mu_x$. An analogous statement holds for $\mu_w$ and $\Sigma_w$.
	By the definition of the Gelbrich distance, the innermost maximization problem over $(\Sigma_x,\Sigma_w)$ in~\eqref{eq:robust2} admits the Lagrangian~dual
	\begin{equation}
	\label{eq:robust-inner}
	\begin{array}{c@{\;}c@{}l}
		\inf\limits_{\dualvar_x, \dualvar_w \geq 0} & \sup\limits_{\Sigma_x, \Sigma_w\succeq 0} & \inner{K^\top K}{\Sigma_x + \mu_x \mu_x^\top} + \inner{A^\top A}{\Sigma_w + \mu_w \mu_w^\top} + b^\top b - 2 \mu_x^\top K^\top A \mu_w - 2 b^\top ( K \mu_x - A \mu_w ) \\[-0.5ex]
		&& \hspace{1cm} + \dualvar_x \big( \rho_x^2 - \| \mu_x - \wh{\m}_x \|^2 - \Tr{\Sigma_x + \covsa_x - 2 \big( \covsa_x^\frac{1}{2} \Sigma_x \covsa_x^\frac{1}{2} \big)^\frac{1}{2}} \big) \\[2ex]
		&& \hspace{1cm} + \dualvar_w\! \big( \rho_w^2 - \| \mu_w - \wh{\m}_w \|^2 - \Tr{\Sigma_w + \covsa_w - 2 \big( \covsa_w^\frac{1}{2} \Sigma_w \covsa_w^\frac{1}{2} \big)^\frac{1}{2}} \big).
	\end{array}
	\end{equation}
	Strong duality holds by \cite[Proposition~5.5.4]{ref:bertsekas2009convex}, which applies because the primal problem has a non-empty compact feasible set. Next, we observe that the inner maximization problem in~\eqref{eq:robust-inner} can be solved analytically by using Proposition~\ref{prop:lin_obj:Lag} in the appendix, and thus the dual problem~\eqref{eq:robust-inner} is equivalent to
	\begin{align}
		\label{eq:robust-inner2}
		\begin{array}{cl}
			\inf\limits_{\substack{\dualvar_x,\dualvar_w \\ \dualvar_x I_n \succ K^\top K \\ \dualvar_w I_m \succ A^\top A}} & \inner{K^\top K}{\mu_x \mu_x^\top} + \inner{A^\top A}{\mu_w \mu_w^\top} + b^\top b
			- 2 \mu_x^\top K^\top A \mu_w - 2 b^\top ( K \mu_x - A \mu_w ) \\[-4ex]
			& \hspace{1cm} + \dualvar_x \big(\rho_x^2 - \| \mu_x - \wh{\m}_x \|^2 -\Tr{ \covsa_x} +  \dualvar_x  \inner{(\dualvar_x I_n - K^\top K )^{-1}}{\covsa_x} \big) \\[2ex]
			& \hspace{1cm} + \dualvar_w \big(\rho_w^2 - \| \mu_w - \wh{\m}_w \|^2 -\Tr{ \covsa_w} + \dualvar_w \inner{(\dualvar_w I_m - A^\top A)^{-1}}{\covsa_w} \big).
		\end{array}
	\end{align}
	Substituting~\eqref{eq:robust-inner2} back into~\eqref{eq:robust2} then allows us to reformulate the Gelbrich MMSE estimation problem~\eqref{prop:Gelbrich}~as
	\begin{align}
		\label{eq:robust3}
		\begin{array}{c@{\,}l}
			\inf\limits_{b } \, \sup\limits_{\substack{\mu_x, \mu_w \\  \| \mu_x - \wh{\m}_x \| \leq \rho_x \\ \| \mu_w - \wh{\m}_w \| \leq \rho_w}} \,
			\inf\limits_{\substack{\dualvar_x,\dualvar_w \\ \dualvar_x I_n \succ K^\top K \\ \dualvar_w I_m \succ A^\top A }} & \inner{K^\top K}{\mu_x \mu_x^\top}
			+ \inner{A^\top A}{\mu_w \mu_w^\top} + b^\top b - 2 \mu_x^\top K^\top A \mu_w - 2 b^\top ( K \mu_x - A \mu_w ) \\[-4ex]
			& \hspace{1cm} + \dualvar_x \big(\rho_x^2 - \| \mu_x - \wh{\m}_x \|^2 -\Tr{ \covsa_x} +  \dualvar_x  \inner{(\dualvar_x I_n - K^\top K )^{-1}}{\covsa_x} \big) \\[2ex]
			& \hspace{1cm} + \dualvar_w \big(\rho_w^2 - \| \mu_w - \wh{\m}_w \|^2 -\Tr{ \covsa_w} + \dualvar_w \inner{(\dualvar_w I_m - A^\top A)^{-1}}{\covsa_w} \big).
		\end{array}
	\end{align}
	The infimum of the inner minimization problem over $(\dualvar_x, \dualvar_w)$ in~\eqref{eq:robust3} is convex quadratic in~$b$. Moreover, it is concave in~$(\mu_x,\mu_w)$ because $K^\top K - \dualvar_x I_n\prec 0$ and $A^\top A - \dualvar_w I_m\prec 0$ for any feasible choice of $(\dualvar_x,\dualvar_w)$ and because concavity is preserved under minimization. Finally, the feasible set for~$(\mu_x,\mu_w)$ is convex and compact. By Sion's classical minimax theorem, we may therefore interchange the infimum over $b$ with the supremum over $(\mu_x, \mu_w)$. The minimization problem over $b$ thus reduces to an unconstrained (strictly) convex quadratic program that has the unique optimal solution $ b = K \mu_x - A \mu_w$. Substituting this expression back into~\eqref{eq:robust3} then yields
	\begin{align}
		\label{eq:robust4}
		\begin{array}{c@{\,}l}
			\sup\limits_{\substack{\mu_x, \mu_w \\  \| \mu_x - \wh{\m}_x \| \leq \rho_x \\ \| \mu_w - \wh{\m}_w \| \leq \rho_w}} \,
			\inf\limits_{\substack{\dualvar_x,\dualvar_w \\ \dualvar_x I_n \succ K^\top K \\ \dualvar_w I_m \succ A^\top A }} &
			\dualvar_x \big( \rho_x^2 - \| \mu_x - \wh{\m}_x \|^2 - \Tr{\covsa_x} \big) + \dualvar_x^2 \inner{(\dualvar_x I_n - K^\top K)^{-1}}{\covsa_x}  \\[-4ex]
			& \hspace{0.5cm} + \dualvar_w \big( \rho_w^2 - \| \mu_w - \wh{\m}_w \|^2 - \Tr{\covsa_w} \big) + \dualvar_w^2 \inner{(\dualvar_w I_m - A^\top A)^{-1}}{\covsa_w}.
		\end{array}
	\end{align}
	It is easy to verify that the resulting maximization problem over $(\mu_x, \mu_w)$ is solved by $\mu_x=\wh{\m}_x$ and $\mu_w= \wh{\m}_w$. Substituting the corresponding optimal value into~\eqref{eq:robust} finally yields
	\[
	\inf_{\substack{b}} \sup\limits_{\mathbb{Q} \in {\mbb G}(\widehat \mathbb{P})} \mathcal{R}(\psi_{A,b}, \mathbb{Q}) = \left\{
	\begin{array}{c@{\,}l}
	\inf\limits_{\substack{\dualvar_x,\dualvar_w \\ \dualvar_x I_n \succ K^\top K \\ \dualvar_w I_m \succ A^\top A }} &
	\dualvar_x \big( \rho_x^2 - \Tr{\covsa_x} \big) + \dualvar_x^2 \inner{(\dualvar_x I_n - K^\top K)^{-1}}{\covsa_x}  \\[-3ex]
	& \hspace{0.5cm} + \dualvar_w \big( \rho_w^2 - \Tr{\covsa_w} \big) + \dualvar_w^2 \inner{(\dualvar_w I_m - A^\top A)^{-1}}{\covsa_w}.
	\end{array}
	\right.
	\]
	From the above equation and the definition of $K$ it is evident that the Gelbrich MMSE estimation problem
	\begin{equation} \label{eq:dro:approx:2}
	\inf_{\psi\in\Ac} \sup\limits_{\mathbb{Q} \in {\mbb G}(\wh \PP)} \mathcal{R}(\psi, \mathbb{Q}) = \inf_{A, b} \sup\limits_{\mathbb{Q} \in {\mbb G}(\wh \PP)} \mathcal{R}(\psi_{A,b}, \mathbb{Q})
	\end{equation}
	is indeed equivalent to the finite convex optimization problem~\eqref{eq:VA}.

	Assume now that $\rho_x>0$ and $\rho_w > 0$. In this case we know from Proposition~\ref{prop:Gelbrich}~\ref{prop:Gelbrich:solvability} that the Gelbrich MMSE estimation problem~\eqref{eq:dro:approx:2} admits an optimal affine estimator $\psi^\star(y)=A^\star y+b^\star$ for some $A^\star\in\mathbb{R}^{n\times m}$ and $b^\star \in\mathbb{R}^m$. The reasoning in the first part of the proof then implies that $A^\star$ solves~\eqref{eq:VA}. Moreover, it implies that $b^\star$ is optimal in~\eqref{eq:robust} when we fix $A=A^\star$. As~\eqref{eq:robust} is equivalent to~\eqref{eq:robust3} and as the unique optimal solution of~\eqref{eq:robust3} for $A=A^\star$ is given by $b = \wh{\m}_x - A^\star (H \wh{\m}_x + \wh{\m}_w)$, we may finally conclude that
	\[
	b^\star = \wh{\m}_x - A^\star (H \wh{\m}_x + \wh{\m}_w).
	\]
	By reversing these arguments, one can further show that if $A^\star$ solves~\eqref{eq:VA} and $b^\star$ is defined as above, then the affine estimator $\psi^\star(y)=A^\star y+b^\star$ is optimal in~\eqref{eq:dro:approx:2}. This observation completes the proof.
\end{proof}

Using Schur complement arguments, the convex program~\eqref{eq:VA} can be further simplified to a standard semidefinite program (SDP), which can be addressed with off-the-shelf solvers.

\begin{corollary}[SDP reformulation] \label{cor:primal:refor}
	The Gelbrich MMSE estimation problem~\eqref{eq:dro:approx} is equivalent to the SDP
	\begin{equation}
	\label{eq:program:SDP:linear}
	\begin{array}{cl}
		\inf & \dualvar_x \big( \rho_x^2 - \Tr{\covsa_x} \big) + \Tr{U_x} + \dualvar_w \big( \rho_w^2 - \Tr{\covsa_w} \big) + \Tr{U_w} \vspace{1mm}\\
		\st & A \in \mathbb{R}^{n \times m}, \quad \dualvar_x, \dualvar_w \in \mathbb{R}_+ \\
		& U_x \in \mathbb{S}_{+}^{n}, \quad V_x \in \mathbb{S}_{+}^{n} , \quad U_w \in \mathbb{S}_{+}^{m}, \quad V_w \in \mathbb{S}_{+}^{m} \\
		& \begin{bmatrix}
			U_x ~ & \dualvar_x \covsa_x^\frac{1}{2} \\
			\dualvar_x \covsa_x^\frac{1}{2} ~ & V_x
		\end{bmatrix} \succeq 0, \quad\,\,
		\begin{bmatrix}
			\dualvar_x I_n - V_x ~ & ~ I_n - H^\top A^\top \\
			I_n - A H ~& ~ I_n
		\end{bmatrix} \succeq 0 \\
		& \begin{bmatrix}
			U_w ~ & \dualvar_w \covsa_w^\frac{1}{2} \\
			\dualvar_w \covsa_w^\frac{1}{2} ~ & V_w
		\end{bmatrix} \succeq 0, \quad
		\begin{bmatrix}
			\dualvar_w I_m - V_w ~ & A^\top \\
			A ~ & I_n
		\end{bmatrix} \succeq 0.
	\end{array}
	\end{equation}
\end{corollary}
\begin{proof}[Proof of Corollary~\ref{cor:primal:refor}]
	Define the extended real-valued function $h_w:\mathbb{R}^{n\times m}\times \mathbb{R}_+\rightarrow (-\infty,\infty]$ through
	\[
	h_w(A,\dualvar_w) = \left\{ \begin{array}{ll} \dualvar_w^2 \inner{(\dualvar_w I_m - A^\top A)^{-1}}{\covsa_w} &\text{if } \dualvar_w I_m - A^\top A \succ 0, \\
	\infty & \text{otherwise.} \end{array}\right.
	\]
	If $\dualvar_w I_m - A^\top A \succ 0$, then, we have
	\begin{align}
		h_w(A,\dualvar_w) =&\inf_{U_w\succeq 0} \left\{ \Tr{U_w}\; :\; U_w \succeq \dualvar_w^2 \covsa_w^\frac{1}{2} (\dualvar_w I_m - A^\top A)^{-1} \covsa_w^\frac{1}{2} \right\} \notag \\
		=&\inf_{U_w\succeq 0, V_w\succ 0} \left\{ \Tr{U_w}\;:\; U_w \succeq \dualvar_w^2 \covsa_w^\frac{1}{2} V_w^{-1} \covsa_w^\frac{1}{2}, ~\dualvar_w I_m - A^\top A \succeq V_w
		\right\} \notag \\
		=&\inf_{U_w\succeq 0, V_w\succ 0} \left\{ \Tr{U_w}\;:\;
		\begin{bmatrix}
			\dualvar_w I_m - V_w & A^\top \\
			A & I_n
		\end{bmatrix} \succeq 0, ~
		\begin{bmatrix}
			U_w & \dualvar_w \covsa_w^\frac{1}{2} \\
			\dualvar_w \covsa_w^\frac{1}{2} & V_w
		\end{bmatrix} \succeq 0
		\right\},
		\label{eq:SDP-reformulation}
	\end{align}
	where the first equality holds due to the cyclicity of the trace operator and because $U_w \succeq \bar U_w$ implies $\Tr{U_w} \geq \Tr{\bar U_w}$ for all $U_w,\bar U_w\succeq 0$, the second equality holds because $V_w\succeq \bar V_w$ is equivalent to $V_w^{-1}\preceq \bar V_w^{-1}$ for all $V_w,\bar V_w\succ 0$, and the last equality follows from standard Schur complement arguments; see, {\text{\emph{e.g.}}}, \cite[\S~A.5.5]{ref:boyd2004convex}. If $\dualvar_w I_m - A^\top A \nsucc 0$, on the other hand, then the first matrix inequality in~\eqref{eq:SDP-reformulation} implies that $V_w$ must have at least one non-positive eigenvalue, which contradicts the constraint $V_w\succ 0$. The SDP~\eqref{eq:SDP-reformulation} is therefore infeasible, and its infimum evaluates to $\infty$. Thus, $h_w(A,\dualvar_w)$ coincides with the optimal value of the SDP~\eqref{eq:SDP-reformulation} for all $A\in\mathbb{R}^{n\times m}$ and $\dualvar_w\in\mathbb{R}_+$.

	A similar SDP reformulation can be derived for the function $h_x:\mathbb{R}^{n\times m}\times \mathbb{R}_+\rightarrow (-\infty,\infty]$ defined through
	\[
	h_x(A,\dualvar_x)= \left\{ \begin{array}{ll} \dualvar_x^2 \inner{[\dualvar_x I_n - (I_n - A H)^\top (I_n - A H)]^{-1}}{\covsa_x} & \text{if } \dualvar_x I_n - (I_n - A H)^\top (I_n - A H) \succ 0,\\
	\infty & \text{otherwise.}\end{array} \right.
	\]
	The claim now follows by substituting the SDP reformulations for $h_w(A,\dualvar_w)$ and $h_x(A,\dualvar_x)$ into~\eqref{eq:VA}. In doing so, we may relax the strict semidefinite inequalities $V_w\succ 0$ and $V_x\succ 0$ to weak inequalities $V_w\succeq 0$ and $V_x\succeq 0$, which amounts to taking the closure of the (non-empty) feasible set and does not change the infimum of problem~\eqref{eq:VA}. This observation completes the proof.
\end{proof}
\begin{remark}[Numerical stability]
	The SDP~\eqref{eq:program:SDP:linear} requires the square roots of the nominal covariance matrices as inputs. Unfortunately, iterative methods for computing matrix square roots often suffer from numerical instability in high dimensions. As a remedy, one may replace those matrix inequalities in~\eqref{eq:program:SDP:linear} that involve $\covsa^\frac{1}{2}_x$ and $\covsa^\frac{1}{2}_w$ with
	\[
	\begin{bmatrix}
	U_x ~ & \dualvar_x \Lambda_x^\top \\
	\dualvar_x \Lambda_x ~ & V_x
	\end{bmatrix} \succeq 0 \qquad \text{and} \qquad \begin{bmatrix}
	U_w ~ & \dualvar_w \Lambda_w^\top \\
	\dualvar_w \Lambda_w ~ & V_w
	\end{bmatrix} \succeq 0,
	\]
	where $\Lambda_x$ and $\Lambda_w$ represent the lower triangular Cholesky factors of $\covsa_x$ and $\covsa_w$, respectively. Thus, we have $\covsa_x = \Lambda_x \Lambda_x^\top$ and $\covsa_w = \Lambda_w \Lambda_w^\top$. We emphasize that $\Lambda_x$ and $\Lambda_w$ can be computed reliably in high dimensions.
\end{remark}

\section{The Dual Wasserstein MMSE Estimation Problem over Normal
	Priors}
\label{sect:dual}

We now examine the dual Wasserstein MMSE estimation
problem
\begin{align}
	\label{eq:dual-dro}
	\mathop{\text{maximize}}_{\mathbb{Q} \in \mathbb B(\wh \PP)}\inf\limits_{\psi\in{\mc F}}  \mathcal{R}(\psi, \mathbb{Q}) ,
\end{align}
which is obtained from~\eqref{eq:dro} by interchanging the order of minimization and maximization. Any maximizer $\mathbb{Q}^\star$ of this dual estimation problem, if it exists, will henceforth be called a {\em least favorable prior}. Unfortunately, problem~\eqref{eq:dual-dro} is
generically intractable. Below we will demonstrate, however,
that~\eqref{eq:dual-dro} becomes tractable if the nominal
distribution~$\wh \PP$ is normal.

\begin{definition}[Normal distributions]
	\label{def:normal-dist}
	We say that $\mathbb P$ is a normal distribution on $\mathbb{R}^d$ with mean $\mu\in\mathbb{R}^d$ and covariance matrix $\Sigma\in\mathbb{S}_{+}^d$, that is, $\mathbb P=\mathcal N(\mu,\Sigma)$, if $\mathbb P$ is supported on $\text{\rm supp}(\mathbb P)=\{\mu+Ev: v\in\mathbb{R}^k\}$, and if the density function of $\mathbb P$ with respect to the Lebesgue measure on $\text{\rm supp}(\mathbb P)$ is given by
	\[
	\varrho_{\mathbb P}(\xi) = \frac{1}{\sqrt{(2\pi)^k\det(D)}} e^{-(\xi-\mu)^\top ED^{-1}E^\top(\xi-\mu)},
	\]
	where $k=\text{\rm rank}(\Sigma)$, $D\in \mathbb{S}_{++}^k$ is the diagonal matrix of the positive eigenvalues of $\Sigma$, and $E\in\mathbb{R}^{d\times k}$ is the matrix whose columns correspond to the orthonormal eigenvectors of the positive eigenvalues of $\Sigma$.
\end{definition}

Definition~\ref{def:normal-dist} also accounts for degenerate normal distributions with singular covariance matrices. We now recall some basic properties of normal distributions that are
crucial for the results of this paper.

\begin{proposition}[{Affine transformations
		\cite[Theorem~2.16]{ref:fang1990symmetric}}]
	\label{prop:normal-affine}
	If $\xi \in\mathbb{R}^d$ follows the normal distribution
	$\mathcal N(\mu, \Sigma)$, while $A \in \mathbb{R}^{k \times d}$ and
	$b \in \mathbb{R}^k$, then $A \xi + b\in\mathbb{R}^k$ follows the
	normal distribution $\mathcal N(A \mu + b, A
	\Sigma A^\top)$.
\end{proposition}

\begin{proposition}[{Affine conditional expectations
		\cite[Corollary~5]{ref:cambanis1981theory}}]
	\label{prop:normal-cond-exp}
	Assume that $\xi \in\mathbb{R}^d$ follows the normal distribution $\mathbb{P}=
	\mathcal N(\mu, \Sigma)$ and that
	\begin{align*}
		\xi = \begin{bmatrix} \xi_1 \\ \xi_2 \end{bmatrix}, \qquad
		\mu = \begin{bmatrix} \mu_1 \\ \mu_2 \end{bmatrix} \qquad
		\text{and} \qquad \Sigma = \begin{bmatrix}
			\Sigma_{11} & \Sigma_{12} \\ \Sigma_{21} & \Sigma_{22}
		\end{bmatrix},
	\end{align*}
	where $\xi_1,\mu_1\in\mathbb{R}^{d_1}$, $\xi_2,\mu_2\in\mathbb{R}^{d_2}$,
	$\Sigma_{11}\in\mathbb{R}^{d_1\times d_1}$, $\Sigma_{22}\in\mathbb{R}^{d_2\times d_2}$ and
	$\Sigma_{12}=\Sigma_{21}^\top \in\mathbb{R}^{d_1\times d_2}$ for some
	$d_1,d_2\in\mathbb N$ with $d_1+d_2= d$. Then, there exist
	$A\in\mathbb{R}^{d_1\times d_2}$ and $b\in\mathbb{R}^{d_1}$ such that
	$\mathds{E}_\mathbb{P}[\xi_1| \xi_2]=A \xi_2 + b $ $\mathbb{P}$-almost surely.
\end{proposition}

Another useful but lesser known property of normal distributions is
that their Wasserstein distances can be expressed analytically in terms of the distributions' first- and second-order moments.

\begin{proposition}[Wasserstein distance between normal distributions {\cite[Proposition~7]{ref:givens1984class}}]
	\label{prop:normal:distance}
	The Wasserstein distance between two normal distributions
	$\mathbb{Q}_1 = \mathcal N(\mu_1, \Sigma_1)$ and $\mathbb{Q}_2 = \mathcal N(\mu_2, \Sigma_2)$ equals the Gelbrich distance between their mean vectors and covariance matrices, that is, $\mathds{W}(\mathbb{Q}_1, \mathbb{Q}_2) = \mathds{G}( (\mu_1, \Sigma_1), (\mu_2,
	\Sigma_2) )$.
\end{proposition}

Assume now that the nominal distributions of the parameter $x\in\mathbb{R}^n$ and the noise $w\in\mathbb{R}^m$ are normal, that is, assume that $\wh \PP_x = \mathcal N(\widehat
\mu_x, \covsa_x)$ and $\wh \PP_w = \mathcal N(\widehat \mu_w,
\covsa_w)$.
Thus, the joint nominal distribution $\wh \PP = \wh \PP_{x} \times
\wh \PP_{w}$ is also normal, that is,
\begin{align}
	\label{eq:nominal:normal}
	\wh \PP=\mathcal N(\widehat \mu, \covsa) \qquad \text{where}
	\qquad \widehat \mu = \begin{bmatrix} \widehat \mu_x \\ \widehat \mu_w
	\end{bmatrix} \qquad \text{and} \qquad \covsa = \begin{bmatrix}
		\covsa_{x} & 0 \\ 0 & \covsa_{w} \end{bmatrix}.
\end{align}

Armed with the fundamental results on normal distributions summarized above, we are now ready to address the dual Wasserstein MMSE estimation problem~\eqref{eq:dual-dro} with a normal nominal distribution. In analogy to Section~\ref{sect:approx}, where we proposed the Gelbrich MMSE estimation problem as an easier conservative approximation for the original {\em primal} estimation problem~\eqref{eq:dro}, we will now construct an easier conservative approximation for the original {\em dual} estimation problem~\eqref{eq:dual-dro}. To this end, we define the restricted ambiguity set
\begin{align}
	\notag
	\mathbb B_{\mathcal N}(\wh \PP) = \left\{ \mathbb{Q}_x \times \mathbb{Q}_w \in \mathcal{M}(\mathbb{R}^{n})\times
	\mathcal{M}(\mathbb{R}^{m}) :
	\begin{array}{l}
		\exists \Sigma_x\in\mathbb{S}_{+}^n,~ \Sigma_w\in\mathbb{S}_{++}^m \text{ with }\\
		\mathbb{Q}_x = \mathcal N(\wh{\m}_x, \Sigma_x), ~\mathbb{Q}_w =
		\mathcal N(\wh{\m}_w, \Sigma_w), \\
		\mathds{W} (\mathbb{Q}_x, \wh \PP_x ) \leq \rho_x,~\mathds{W} (\mathbb{Q}_w, \wh \PP_w) \leq \rho_w
	\end{array}
	\right\}.
\end{align}

By construction, $\mathbb B_{\mathcal N}(\wh \PP)$ contains all {\em normal}
distributions $\mathbb{Q}=\mathbb{Q}_x \times \mathbb{Q}_w$ from within the original
Wasserstein ambiguity set $\mathbb B(\wh \PP)$ that have the same mean vector
$(\wh{\m}_x, \wh{\m}_w)$ as the nominal
distribution~$\wh \PP=\wh \PP_x\times\wh \PP_w$, and where the covariance matrix of $\mathbb{Q}_w$ is strictly positive definite.
Thus, we have $\mathbb B_{\mathcal N}(\wh \PP) \subseteq \mathbb B(\wh \PP)$. Note also that $\mathbb B_{\mathcal N}(\wh \PP)$ is non-convex because mixtures of normal distributions usually fail to be normal.

By restricting the original Wasserstein ambiguity set $\mathbb B(\wh \PP)$ to
its subset $\mathbb B_{\mathcal N}(\wh \PP)$, we obtain the following conservative
approximation for the dual Wasserstein MMSE estimation
problem~\eqref{eq:dual-dro}.
\begin{align}
	\label{eq:dual-dro:conservative}
	\mathop{\text{maximize}}_{\mathbb{Q} \in \mathbb
		B_{\mathcal N}(\wh \PP)}\inf\limits_{\psi\in{\mc F}}  \mathcal{R}(\psi, \mathbb{Q})
\end{align}
We will henceforth refer to~\eqref{eq:dual-dro:conservative} as the dual Wasserstein MMSE estimation problem {\em over normal priors}. The
following main theorem shows that~\eqref{eq:dual-dro:conservative} is
equivalent to a finite convex optimization problem.
\begin{theorem}[Dual Wasserstein MMSE estimation problem over
	normal priors]
	\label{thm:least-favorable-prior}
	Assume that the Wasserstein ambiguity set $\mathbb B_{\mathcal N}(\wh \PP)$ is centered at a normal distribution $\wh \PP$ of the form~\eqref{eq:nominal:normal}. Then, the dual  Wasserstein MMSE estimation problem over normal priors~\eqref{eq:dual-dro:conservative} is equivalent to the finite convex optimization problem
	\begin{equation}
	\label{eq:program:dual}
	\begin{array}{cl}
		\sup & \Tr{\Sigma_x - \Sigma_x H^\top \left( H \Sigma_x H^\top +
			\Sigma_w \right)^{-1} H \Sigma_x} \\
		\st &  \Sigma_x \in \mathbb{S}_{+}^{n}, \quad \Sigma_w \in \mathbb{S}_{++}^{m} \\[0.5ex]
		& \Tr{\Sigma_x + \covsa_x - 2 \big( \covsa_x^\frac{1}{2}
			\Sigma_x \covsa_x^\frac{1}{2} \big)^\frac{1}{2}} \leq \rho_x^2 \\[0.5ex]
		& \Tr{\Sigma_w + \covsa_w - 2 \big( \covsa_w^\frac{1}{2}
			\Sigma_w \covsa_w^\frac{1}{2} \big)^\frac{1}{2}} \leq \rho_w^2.
	\end{array}
	\end{equation}
	If $\covsa_w \succ 0$, then~\eqref{eq:program:dual} is solvable, and the maximizer denoted by~$(\Sigma_x^\star, \Sigma_w^\star)$ satisfies $\Sigma_x^\star \succeq \lambda_{\min}(\covsa_x) I_n$ and $\Sigma_w^\star \succeq \lambda_{\min}(\covsa_w) I_m$. Moreover, the supremum of~\eqref{eq:dual-dro:conservative} is attained by the normal distribution $\mathbb{Q}^\star=\mathbb{Q}^\star_x\times\mathbb{Q}^\star_w$ defined through $\mathbb{Q}^\star_x = \mathcal N(\wh{\m}_x, \Sigma_x^\star)$ and $\mathbb{Q}^\star_w= \mathcal N(\wh{\m}_w, \Sigma_w^\star)$.
\end{theorem}

\begin{proof}[Proof of Theorem~\ref{thm:least-favorable-prior}]
	If $(x,w)$ is governed by a normal distribution $\mathbb{Q} \in
	\mathbb B_{\mathcal N}(\wh \PP)$, then the linear transformation $(x, y) = (x, Hx +
	w)$ is also normally distributed by virtue of
	Proposition~\ref{prop:normal-affine}, and the average risk $\mathcal{R}(\psi,
	\mathbb{Q})$ is minimized by the Bayesian MMSE estimator $\psi^\star_{\mathcal
		B}(y)=\mathds{E}_{\mathbb{P}_{x|y}} [x]$, which is affine due to
	Proposition~\ref{prop:normal-cond-exp}. Thus, in the dual Wasserstein
	MMSE estimation problem with normal priors, the set~${\mc F}$ of {\em
		all} estimators may be restricted to the set~$\Ac$ of all {\em affine}
	estimators without sacrificing optimality, that is,
	\begin{equation}
	\label{eq:dual-simplification}
	\sup\limits_{\mathbb{Q} \in \mathbb B_{\mathcal N}(\wh \PP)} \inf\limits_{\psi\in{\mc F}} \mathcal{R}(\psi,
	\mathbb{Q}) = \sup\limits_{\mathbb{Q} \in \mathbb B_{\mathcal N}(\wh \PP)} \inf\limits_{\psi\in\Ac}  \mathcal{R}(\psi,
	\mathbb{Q}).
	\end{equation}
	As the average risk $\mathcal{R}(\psi, \mathbb{Q})$ of an affine estimator
	$\psi\in\Ac$ simply evaluates the expectation of a quadratic function in
	$(x,w)$, it depends on $\mathbb{Q}$ only through its first and second moments.
	Moreover, as $\mathbb{Q}$ and $\wh \PP$ are normal distributions, their Wasserstein distance coincides with
	the Gelbrich distance between their mean vectors and covariance
	matrices; see Proposition~\ref{prop:normal:distance}. Thus, the
	maximization problem over $\mathbb{Q}$ on the right hand side
	of~\eqref{eq:dual-simplification} can be recast as an equivalent
	maximization problem over the first and second moments of~$\mathbb{Q}$.
	Specifically, by the definitions of $\mathcal{R}(\psi, \mathbb{Q})$ and $ \mathbb
	B_\phi(\wh \PP)$ we~find
	\begin{align*}
		\sup\limits_{\mathbb{Q} \in \mathbb B_{\mathcal N}(\wh \PP)} \inf\limits_{\psi\in\Ac} \mathcal{R}(\psi, \mathbb{Q})
		= & \left\{
		\begin{array}{cl}
			\sup\limits_{\Sigma_x, \Sigma_w} & \inf\limits_{\substack{A, K \\ K =
					I_n-A H}} \, \inf\limits_{b} ~\; \inner{K^\top K}{\Sigma_x + \wh{\m}_x
				\wh{\m}_x^\top} + \inner{A^\top A}{\Sigma_w + \wh{\m}_w \wh{\m}_w^\top} +
			b^\top b \\[-2ex]
			& \hspace{2.5cm} - 2 \wh{\m}_x^\top K^\top A \wh{\m}_w -
			2b^\top (K\wh{\m}_x - A \wh{\m}_w)\\[1ex]
			\st & \mathds{G}\big((\wh{\m}_x, \Sigma_x), (\wh{\m}_x,
			\covsa_x)\big)\leq \rho_x,~\mathds{G}\big( (\wh{\m}_w, \Sigma_w), (\wh{\m}_w,
			\covsa_w) \big) \leq \rho_w \\
			& \Sigma_x \succeq 0,~\Sigma_w \succ
			0,        \end{array}
		\right.
	\end{align*}
	where the auxiliary decision variable $K = I_n-A H$ has
	been introduced to simplify the objective function. The innermost minimization problem over $b$ constitutes an unconstrained (strictly) convex quadratic program that has the unique optimal solution $b= K\wh{\m}_x - A \wh{\m}_w$. Substituting this minimizer back into the objective function of the above problem and recalling the definition of the Gelbrich distance then yields
	\begin{align}
		\sup\limits_{\mathbb{Q} \in \mathbb B_{\mathcal N}(\wh \PP)} \inf\limits_{\psi\in\Ac} \mathcal{R}(\psi,
		\mathbb{Q}) =& \left\{
		\begin{array}{cl}
			\sup\limits_{\Sigma_x, \Sigma_w} & \inf\limits_{\substack{A, K \\ K =
					I_n-A H}} \; \inner{K^\top K}{\Sigma_x} + \inner{A^\top
				A}{\Sigma_w} \\
			\st & \Tr{\Sigma_x + \covsa_x - 2 \big( \covsa_x^\frac{1}{2}
				\Sigma_x \covsa_x^\frac{1}{2} \big)^\frac{1}{2}} \leq \rho_x^2 \\[2ex]
			& \Tr{\Sigma_w + \covsa_w - 2 \big( \covsa_w^\frac{1}{2}
				\Sigma_w \covsa_w^\frac{1}{2} \big)^\frac{1}{2}} \leq \rho_w^2 \\[2ex]
			& \Sigma_x \succeq 0,~\Sigma_w \succ 0.
		\end{array}
		\right.    \label{eq:least1}
	\end{align}
	By using the equality $K = I_n-A H$ to eliminate $K$, the inner minimization problem in~\eqref{eq:least1} can be reformulated as an unconstrained quadratic program in~$A$. As $\Sigma_w \succ 0$, this quadratic program is strictly convex, and an elementary calculation reveals that its unique optimal solution is given by
	\[
	A^\star =  \Sigma_x H^\top \left( H {\Sigma_x} H^\top +
	\Sigma_w\right)^{-1}.
	\]
	Substituting $A^\star$ as well as the corresponding auxiliary decision variable $K^\star=I_n-A^\star H$ into the objective function of~\eqref{eq:least1} finally yields the postulated convex program~\eqref{eq:program:dual}.

	Assume now that $\covsa_w \succ 0$, and define
	\[
	\mathcal S_x = \left\{ \Sigma_x \in \mathbb{S}_{+}^n: \mathds{G}\big((\wh{\m}_x, \Sigma_x), (\wh{\m}_x,
	\covsa_x)\big)\leq \rho_x \right\}\quad \text{and}\quad
	\mathcal S_w = \left\{ \Sigma_w \in \mathbb{S}_{+}^m: \mathds{G}\big( (\wh{\m}_w, \Sigma_w), (\wh{\m}_w,
	\covsa_w) \big) \leq \rho_w \right\}.
	\]
	Equations~\eqref{eq:dual-simplification} and~\eqref{eq:least1} imply that
	\begin{align}
		\sup\limits_{\mathbb{Q} \in \mathbb B_{\mathcal N}(\wh \PP)} \inf\limits_{\psi\in{\mc F}} \mathcal{R}(\psi,
		\mathbb{Q}) \le &  \sup\limits_{\Sigma_x \in \mathcal S_x} \sup\limits_{\Sigma_w \in \mathcal S_w}  \inf\limits_{\substack{A, K \\ K =
				I_n-A H}} \; \inner{K^\top K}{\Sigma_x} + \inner{A^\top
			A}{\Sigma_w} \notag \\
		=& \sup\limits_{\substack{\Sigma_x \in \mathcal S_x \\ \Sigma_x \succeq \lambda_{\min}(\covsa_x) I_n}} \sup\limits_{\substack{\Sigma_w \in \mathcal S_w \\ \Sigma_w \succeq \lambda_{\min}(\covsa_w) I_m}} \inf\limits_{\substack{A, K \\ K =
				I_n-A H}} \; \inner{K^\top K}{\Sigma_x} + \inner{A^\top
			A}{\Sigma_w}, \label{eq:least2}
	\end{align}
	where the inequality holds because we relax the requirement that $\Sigma_w$ be strictly positive definite, and the equality follows from applying Lemma~\ref{lemma:monotone loss} consecutively to each of the two maximization problems. If $\covsa_w \succ 0$, then problem~\eqref{eq:least2} constitutes a restriction of~\eqref{eq:least1} and therefore provides also a lower bound on the dual Wasserstein MMSE estimation problem. In summary, we thus have
	\begin{align}
		\sup\limits_{\mathbb{Q} \in \mathbb B_{\mathcal N}(\wh \PP)} \inf\limits_{\psi\in{\mc F}} \mathcal{R}(\psi,
		\mathbb{Q}) =& \left\{
		\begin{array}{cl}
			\sup\limits_{\Sigma_x, \Sigma_w} & \inf\limits_{\substack{A, K \\ K =
					I_n-A H}} \; \inner{K^\top K}{\Sigma_x} + \inner{A^\top
				A}{\Sigma_w} \\
			\st & \Tr{\Sigma_x + \covsa_x - 2 \big( \covsa_x^\frac{1}{2}
				\Sigma_x \covsa_x^\frac{1}{2} \big)^\frac{1}{2}} \leq \rho_x^2 \\[2ex]
			& \Tr{\Sigma_w + \covsa_w - 2 \big( \covsa_w^\frac{1}{2}
				\Sigma_w \covsa_w^\frac{1}{2} \big)^\frac{1}{2}} \leq \rho_w^2 \\[2ex]
			& \Sigma_x \succeq\lambda_{\min}(\covsa_x) I_n, ~\Sigma_w \succeq \lambda_{\min}(\covsa_w) I_m.
		\end{array}
		\right. \label{eq:least3}
	\end{align}
	This reasoning implies that if $\covsa_w \succ 0$, then the constraints $\Sigma_x \succeq \lambda_{\min}(\covsa_x) I_n$ and $\Sigma_w \succeq \lambda_{\min}(\covsa_w) I_m$ can be appended to problem~\eqref{eq:least1} and, consequently, to problem~\eqref{eq:program:dual} without altering their common optimal value. Problem~\eqref{eq:program:dual} with the additional constraints $\Sigma_x \succeq \lambda_{\min}(\covsa_x) I_n$ and $\Sigma_w \succeq \lambda_{\min}(\covsa_w) I_m$ has a continuous objective function over a compact feasible set and is thus solvable. Any of its optimal solutions is also optimal in problem~\eqref{eq:program:dual}, which has no redundant constraints. Thus, problem~\eqref{eq:program:dual} is solvable.

	It remains to show that $\mathbb{Q}^\star$ as constructed in the theorem statement is optimal in~\eqref{eq:dual-dro:conservative}. The feasibility of $(\Sigma_x^\star, \Sigma_w^\star)$ in~\eqref{eq:program:dual} implies that $\mathbb{Q}^\star \in \mathbb B_{\mathcal N}(\wh \PP)$, and thus $\mathbb{Q}^\star$ is feasible in~\eqref{eq:dual-dro:conservative}. Moreover, we have
	\begin{align}
		\label{eq:q-optimality}
		\sup\limits_{\mathbb{Q} \in  \mathbb B_{\mathcal N} (\wh \PP)} \, \inf\limits_{\psi \in {\mc F}} \,
		\mathcal{R}(\psi, \mathbb{Q}) &\ge \inf\limits_{\psi \in {\mc F}} \, \mathcal{R}(\psi, \mathbb{Q}^\star) =
		\Tr{\Sigma_x^\star - \Sigma_x^\star H^\top \left( H \Sigma_x^\star H^\top +
			\Sigma_w^\star \right)^{-1} H \Sigma_x^\star},
	\end{align}
	where the equality follows from elementary algebra, recalling that the affine estimator $\psi(y) = A^\star y+b^\star$ with
	\[
	A^\star =  \Sigma_x^\star  H^\top \left( H {\Sigma_x^\star}
	H^\top + \Sigma_w^\star\right)^{-1}\quad \text{and} \quad b^\star=
	\mu_x -  A^\star (H \wh{\m}_x + \wh{\m}_w)
	\]
	is the Bayesian MMSE estimator for the normal distribution
	$\mathbb{Q}^\star$. As the right hand side of~\eqref{eq:q-optimality} coincides
	with the maximum of~\eqref{eq:program:dual} and as
	problem~\eqref{eq:program:dual} is equivalent to the dual Wasserstein
	MMSE estimation problem~\eqref{eq:dual-dro:conservative} over normal
	priors, we may thus conclude that the
	inequality in~\eqref{eq:q-optimality} is tight. Thus, we find
	\[
	\sup\limits_{\mathbb{Q} \in  \mathbb B_{\mathcal N}(\wh \PP)} \, \inf\limits_{\psi \in {\mc F}} \,
	\mathcal{R}(\psi, \mathbb{Q}) = \inf\limits_{\psi \in {\mc F}} \, \mathcal{R}(\psi, \mathbb{Q}^\star),
	\]
	which in turn implies that $\mathbb{Q}^\star$ is optimal
	in~\eqref{eq:dual-dro:conservative}. This observation completes the proof.
\end{proof}

\begin{remark}[Singular covariance matrices]
	A nonlinear SDP akin to~\eqref{eq:program:dual} has been derived in~\cite{ref:shafieezadeh2018wasserstein} under the stronger assumption that the covariance matrix of the nominal distribution $\wh \PP$ is non-degenerate, which implies that $\covsa_x \succ 0$ and $\covsa_w \succ 0$. However, the weaker condition $\covsa_w \succ 0$ is sufficient to ensure that the matrix inversion in the objective function of problem~\eqref{eq:program:dual} is well-defined. Therefore, Theorem~\ref{thm:least-favorable-prior} remains valid if the nominal covariance matrix $\covsa_x$ is singular, which occurs in many applications.
	On the other hand, it is common to require that $\covsa_w = \sigma^2 I_m$ for some $\sigma>0$, see, e.g.,~\cite{ref:chang2000adaptive}.
\end{remark}

Corollary~\ref{corol:dual:refor} below asserts that the convex program~\eqref{eq:program:dual} admits a canonical linear SDP reformulation. The proof is omitted as it relies on standard Schur complement arguments familiar from the proof of Corollary~\ref{cor:primal:refor}.
\begin{corollary}[SDP reformulation] \label{corol:dual:refor}
	Assume that the Wasserstein ambiguity set $\mathbb B_{\mathcal N}(\wh \PP)$ is centered at a normal distribution $\wh \PP$ of the form~\eqref{eq:nominal:normal} with noise covariance matrix $\covsa_w \succ 0$. Then, the dual Wasserstein MMSE estimation problem~\eqref{eq:dual-dro:conservative} over normal priors is equivalent to the SDP
	\begin{equation} \label{eq:SDP:dual}
	\begin{array}{cl}
		\max & \Tr{\Sigma_x} - \Tr{U}  \\[1ex]
		\st & \Sigma_x \in \mathbb{S}_{+}^n, \, \Sigma_w \in \mathbb{S}_{+}^m, \, V_x \in
		\mathbb{S}_{+}^n, \, V_w \in \mathbb{S}_{+}^m, \, U \in \mathbb{S}_{+}^n \\[1ex]
		& \begin{bmatrix} \covsa_x^\frac{1}{2} \Sigma_x \covsa_x^\frac{1}{2} &
			V_x \\ V_x & I_n \end{bmatrix} \succeq 0, \quad
		\begin{bmatrix} \covsa_w^\frac{1}{2} \Sigma_w \covsa_w^\frac{1}{2} & V_w
			\\ V_w & I_m \end{bmatrix} \succeq 0\\[3ex]
		& \Tr{\Sigma_x + \covsa_x - 2V_x} \leq \rho_x^2, \quad
		\Tr{\Sigma_w + \covsa_w - 2V_w} \leq \rho_w^2 \\[2ex]
		& \begin{bmatrix} U & \Sigma_x H^\top \\ H \Sigma_x & H \Sigma_x
			H^\top + \Sigma_w \end{bmatrix} \succeq 0, \quad \Sigma_x \succeq
		\lambda_{\min}(\covsa_x) I_n, \quad \Sigma_w \succeq
		\lambda_{\min}(\covsa_w) I_m.
	\end{array}
	\end{equation}
\end{corollary}
We emphasize that the lower bounds on $\Sigma_x $ and $\Sigma_w$ are redundant but have been made explicit in~\eqref{eq:SDP:dual}.
\section{Nash Equilibrium and Optimality of Affine Estimators}
\label{sect:nash}
If $\wh \PP$ is a normal distribution of the form~\eqref{eq:nominal:normal}, then we have
\begin{equation} \label{eq:strong:duality}
\inf\limits_{\psi\in\Ac} \sup\limits_{\mathbb{Q} \in \mathbb G(\wh \PP)} \mathcal{R}(\psi, \mathbb{Q}) \ge
\inf\limits_{\psi\in{\mc F}} \sup\limits_{\mathbb{Q} \in \mathbb B(\wh \PP)} \mathcal{R}(\psi, \mathbb{Q}) \ge
\sup\limits_{\mathbb{Q} \in \mathbb B(\wh \PP)}\inf\limits_{\psi\in{\mc F}}  \mathcal{R}(\psi, \mathbb{Q}) \ge
\sup\limits_{\mathbb{Q} \in \mathbb B_{\mathcal N}(\wh \PP)}\inf\limits_{\psi\in{\mc F}}  \mathcal{R}(\psi, \mathbb{Q}),
\end{equation}
where the first inequality follows from the inclusions $\Ac \subseteq {\mc F}$ and $\mathbb B(\wh \PP) \subseteq \mathbb G(\wh \PP)$, the second inequality exploits weak duality, and the last inequality holds due to the inclusion $\mathbb B_{\mathcal N}(\wh \PP) \subseteq \mathbb B(\wh \PP)$. Note that the leftmost minimax problem is the Gelbrich MMSE estimation problem~\eqref{eq:dro:approx} studied in Section~\ref{sect:approx}, and the rightmost maximin problem is the dual Wasserstein MMSE estimation problem~\eqref{eq:dual-dro:conservative} over normal priors studied in Section~\ref{sect:dual}. We also highlight that these restricted primal and dual estimation problems sandwich the original Wasserstein estimation problems~\eqref{eq:dro} and~\eqref{eq:dual-dro}, which coincide with the second and third problems in~\eqref{eq:strong:duality}, respectively. The following theorem asserts that all inequalities in~\eqref{eq:strong:duality} actually collapse to equalities.

\begin{theorem}[Sandwich theorem]
	\label{thm:sandwich}
	If $\wh \PP$ is a normal distribution of the form~\eqref{eq:nominal:normal}, then the optimal values of the restricted primal and dual estimation problems~\eqref{eq:dro:approx} and~\eqref{eq:dual-dro:conservative} coincide, i.e.,
	\[
	\inf\limits_{\psi\in\Ac} \sup\limits_{\mathbb{Q} \in \mathbb G(\wh \PP)} \mathcal{R}(\psi, \mathbb{Q}) = \sup\limits_{\mathbb{Q} \in \mathbb B_{\mathcal N}(\wh \PP)}\inf\limits_{\psi\in{\mc F}}  \mathcal{R}(\psi, \mathbb{Q}).
	\]
\end{theorem}
\begin{proof}[Proof of Theorem~\ref{thm:sandwich}]
	By Theorem~\ref{thm:conservative}, the Gelbrich MMSE estimation problem~\eqref{eq:dro:approx} can be expressed as
	\begin{align*}
		\inf\limits_{\psi\in\Ac} \sup\limits_{\mathbb{Q} \in {\mbb G}(\wh \PP)} \mathcal{R}(\psi, \mathbb{Q}) =
		\left\{
		\begin{array}{ccl}
			\inf\limits_{\substack{A, K \\ K = I_n - A H}} & 	\inf\limits_{\substack{\dualvar_x,\dualvar_w \\ \dualvar_x I_n \succ K^\top K \\ \dualvar_w I_m \succ A^\top A }} &
			\dualvar_x \big( \rho_x^2 - \Tr{\covsa_x} \big) + \dualvar_x^2 \inner{(\dualvar_x I_n - K^\top K)^{-1}}{\covsa_x}  \\[-4ex]
			&& \hspace{0.5cm} + \dualvar_w \big( \rho_w^2 - \Tr{\covsa_w} \big) + \dualvar_w^2 \inner{(\dualvar_w I_m - A^\top A)^{-1}}{\covsa_w},
		\end{array}
		\right.
	\end{align*}
	where the auxiliary variable $K = I_n-A H$ has been introduced to highlight the problem's symmetries.
	Next, we introduce the feasible sets
	\[
	\mathcal S_x = \left\{ \Sigma_x \in \mathbb{S}_{+}^n: \Tr{\Sigma_x + \covsa_x - 2 \big( \covsa_x^\frac{1}{2} \Sigma_x \covsa_x^\frac{1}{2} \big)^\frac{1}{2}} \leq \rho_x^2  \right\}
	\]
	and
	\[
	\mathcal S_w = \left\{ \Sigma_w \in \mathbb{S}_{+}^m: \Tr{\Sigma_w + \covsa_w - 2 \big( \covsa_w^\frac{1}{2} \Sigma_w \covsa_w^\frac{1}{2} \big)^\frac{1}{2}} \leq \rho_w^2  \right\},
	\]
	both of which are convex and compact by virtue of Lemma~\ref{lemma:compact:FS}.
	Using Proposition~\ref{prop:quadratic}\,\ref{prop:quad:dual} to reformulate the inner minimization problem over $\gamma_x$ and $\gamma_w$, we then obtain
	\begin{align}
		\inf\limits_{\psi\in\Ac} \sup\limits_{\mathbb{Q} \in {\mbb G}(\wh \PP)} \mathcal{R}(\psi, \mathbb{Q}) &=
		\inf\limits_{\substack{A, K \\ K = I_n - A H}} \sup\limits_{\substack{\Sigma_x \in \mathcal S_x \\ \Sigma_w \in \mathcal S_w}} ~\inner{K^\top K}{\Sigma_x} + \inner{A^\top A}{\Sigma_w} \notag
		\\
		&= \sup\limits_{\Sigma_x \in \mathcal S_x} \inf\limits_{\substack{A, K \\ K = I_n - A H}} \sup\limits_{ \Sigma_w \in \mathcal S_w} ~\inner{K^\top K}{\Sigma_x} + \inner{A^\top A}{\Sigma_w} \notag
	\end{align}
	where the second equality holds due to Sion's minimax theorem~\cite{ref:sion1958minimax}. Define now the auxiliary function
	\[
	    f(A) = \sup\limits_{\Sigma_w \in \mathcal S_w}~\inner{A^\top A}{\Sigma_w}.
	\]
	As $\Sigma_w \succeq 0$ for any $\Sigma_w \in \mathcal S_w$, $f$ constitutes a pointwise maximum of convex functions and is therefore itself convex. In addition, as the set $\mathcal S_w$ is compact by Lemma~\ref{lemma:compact:FS}, $f$ is everywhere finite and thus continuous thanks to \cite[Theorem~2.35]{ref:rockafellar2010variational}. This allows us to conclude that
	\begin{align}
		\inf\limits_{\psi\in\Ac} \sup\limits_{\mathbb{Q} \in {\mbb G}(\wh \PP)} \mathcal{R}(\psi, \mathbb{Q})
		&= \sup\limits_{\Sigma_x \in \mathcal S_x} \inf\limits_{\substack{A, K \\ K = I_n - A H}} ~\inner{K^\top K}{\Sigma_x} + f(A) \notag \\
		&= \sup\limits_{\substack{\Sigma_x \in \mathcal S_x \\ \Sigma_x \succeq \lambda_{\min}(\covsa_x) I_n}} \inf\limits_{\substack{A, K \\ K = I_n - A H}} ~\inner{K^\top K}{\Sigma_x} + f(A) \notag \\
		&=  \inf\limits_{\substack{A, K \\ K = I_n - A H}} \sup\limits_{\substack{\Sigma_x \in \mathcal S_x \\ \Sigma_x \succeq \lambda_{\min}(\covsa_x) I_n}} \sup\limits_{ \Sigma_w \in \mathcal S_w} ~\inner{K^\top K}{\Sigma_x} + \inner{A^\top A}{\Sigma_w}, \notag
	\end{align}
	where the second equality exploits Lemma~\ref{lemma:monotone loss}, and the last equality follows from Sion's minimax theorem~\cite{ref:sion1958minimax}. Another (trivial) application of Lemma~\ref{lemma:monotone loss} then allows us to append the constraint $\Sigma_w \succeq \lambda_{\min}(\covsa_w) I_m$ to the maximization problem over $\Sigma_w$. Sion's minimax theorem~\cite{ref:sion1958minimax} finally implies that
	\begin{align*}
		\inf\limits_{\psi\in\Ac} \sup\limits_{\mathbb{Q} \in {\mbb G}(\wh \PP)} \mathcal{R}(\psi, \mathbb{Q})&= \sup\limits_{\substack{\Sigma_x \in \mathcal S_x \\ \Sigma_x \succeq \lambda_{\min}(\covsa_x) I_n}} \sup\limits_{\substack{ \Sigma_w \in \mathcal S_w \\ \Sigma_w \succeq \lambda_{\min}(\covsa_w) I_m}} \inf\limits_{\substack{A, K \\ K = I_n - A H}} ~\inner{K^\top K}{\Sigma_x} + \inner{A^\top A}{\Sigma_w} \notag \\
		&= \sup\limits_{\mathbb{Q} \in \mathbb B_{\mathcal N}(\wh \PP)} \inf\limits_{\psi\in{\mc F}} \mathcal{R}(\psi, \mathbb{Q}), \notag
	\end{align*}
	where the last equality has already been established in the proof of Theorem~\ref{thm:least-favorable-prior}; see Equation~\eqref{eq:least3}. Thus, the claim follows.
\end{proof}

Theorem~\ref{thm:sandwich} suggests that solving any of the restricted estimation problems is tantamount to solving both original primal and dual estimation problems. This intuition is formalized in the following corollary.

\begin{corollary}[Nash equilibrium] \label{corol:nash}
	If $\wh \PP$ is a normal distribution of the form~\eqref{eq:nominal:normal} with $\covsa_w \succ 0$, then the affine estimator~$\psi^\star$ that solves~\eqref{eq:dro:approx} is optimal in the primal Wasserstein MMSE estimation problem~\eqref{eq:dro}, while the normal distribution~$\mathbb{Q}^\star$ that solves~\eqref{eq:dual-dro:conservative} is optimal in the dual Wasserstein MMSE estimation problem~\eqref{eq:dual-dro}. Moreover, $\psi^\star$ and $\mathbb{Q}^\star$ form a Nash equilibrium for the game between the statistician and nature, that is,
	\begin{equation}
	\label{eq:nash}
	\mathcal{R}(\psi^\star, \mathbb{Q}) \le \mathcal{R}(\psi^\star, \mathbb{Q}^\star) \le \mathcal{R}(\psi, {\mathbb{Q}}^\star) \quad \forall \psi \in {\mc F},~ \mathbb{Q} \in  \mathbb B(\wh \PP) \,.
	\end{equation}
\end{corollary}
\begin{proof}[Proof of Corollary~\ref{corol:nash}]
	As $\covsa_w \succ 0$, the Gelbrich MMSE estimation problem~\eqref{eq:dro:approx} is solved by the affine estimator $\psi^\star$ defined in Theorem~\ref{thm:conservative}, and the dual Wasserstein MMSE estimation problem~\eqref{eq:dual-dro} over normal priors is solved by the normal distribution $\mathbb{Q}^\star$ defined in Theorem~\ref{thm:least-favorable-prior}. Thus, we have
	\begin{align*}
		\mathcal{R} (\psi^\star, \mathbb{Q}^\star) \ge \inf\limits_{\psi \in {\mc F}} \mathcal{R} (\psi, \mathbb{Q}^\star) = \max\limits_{\mathbb{Q} \in \mathbb B_{\mathcal N}(\wh \PP)} \inf\limits_{\psi\in{\mc F}} \mathcal{R}(\psi, \mathbb{Q}) = \min\limits_{\psi\in\Ac} \sup\limits_{\mathbb{Q} \in {\mbb G}(\wh \PP)} \mathcal{R}(\psi, \mathbb{Q}) = \sup\limits_{\mathbb{Q} \in {\mbb G}(\wh \PP)} \mathcal{R}(\psi^\star, \mathbb{Q}) \ge \mathcal{R} (\psi^\star, \mathbb{Q}^\star),
	\end{align*}
	where the three equalities follow from the definition of $\mathbb{Q}^\star$,  Theorem~\ref{thm:sandwich} and the definition of $\psi^\star$, respectively. As the left and the right hand sides of the above expression coincide, we may then conclude that
	\[
	\mathcal{R}(\psi^\star, \mathbb{Q}) \le \mathcal{R}(\psi^\star, \mathbb{Q}^\star) \le \mathcal{R}(\psi, {\mathbb{Q}}^\star) \quad \forall \psi \in {\mc F},~ \mathbb{Q} \in  \mathbb G(\wh \PP).
	\]
	Moreover, as $\mathbb B(\wh \PP) \subseteq \mathbb G(\wh \PP)$, the above relation implies~\eqref{eq:nash}.

	It remains to be shown that $\psi^\star$ and $\mathbb{Q}^\star$ solve the primal and dual Wasserstein MMSE estimation problems~\eqref{eq:dro} and~\eqref{eq:dual-dro}, respectively. As for $\psi^\star$, we have
	\[
	\sup\limits_{\mathbb{Q} \in \mathbb B(\wh \PP)} \mathcal{R}(\psi^\star, \mathbb{Q}) \le \sup\limits_{\mathbb{Q} \in {\mbb G}(\wh \PP)} \mathcal{R}(\psi^\star, \mathbb{Q}) = \inf\limits_{\psi\in\Ac} \sup\limits_{\mathbb{Q} \in \mathbb G(\wh \PP)} \mathcal{R}(\psi, \mathbb{Q}) = \inf\limits_{\psi\in{\mc F}} \sup\limits_{\mathbb{Q} \in \mathbb B(\wh \PP)} \mathcal{R}(\psi, \mathbb{Q}),
	\]
	where the inequality holds because $\mathbb B(\wh \PP) \subseteq {\mbb G}(\wh \PP)$. The first equality follows from  the definition of $\psi^\star$, while the second equality exploits Theorem~\ref{thm:sandwich}, which implies that all inequalities in~\eqref{eq:strong:duality} are in fact equalities. This reasoning shows that $\psi^\star$ is optimal in~\eqref{eq:dro}. The optimality of $\mathbb{Q}^\star$ in~\eqref{eq:dual-dro} can be proved similarly.
\end{proof}

Corollary~\ref{corol:nash} implies that $\psi^\star$ can be viewed as a Bayesian estimator for the least favorable prior $\mathbb{Q}^\star$ and that $\mathbb{Q}^\star$ represents a worst-case distribution for the optimal estimator $\psi^\star$. Next, we will argue that $\psi^\star$ can not only be constructed from the solution of the convex program~\eqref{eq:VA}, which is equivalent to the Gelbrich MMSE estimation problem~\eqref{eq:dro:approx}, but also from the solution of the convex program~\eqref{eq:program:dual}, which is equivalent to the dual MMSE estimation problem~\eqref{eq:dual-dro:conservative} over normal~priors. This alternative construction is useful because problem~\eqref{eq:program:dual} is amenable to highly efficient first-order methods to be derived in Section~\ref{sect:algorithm}.

\begin{corollary}[Dual construction of the optimal estimator]
	\label{corol:alternative}
	If $\wh \PP$ is a normal distribution of the form~\eqref{eq:nominal:normal} with $\covsa_w\succ 0$, and $(\Sigma_x^\star, \Sigma_w^\star)$ is a maximizer of~\eqref{eq:program:dual}, then the affine estimator $\psi^\star(y) = A^\star y + b^\star$ with
	\begin{equation}
	\label{eq:alternative:def}
	A^\star =  \Sigma_x^\star H^\top \left( H {\Sigma_x^\star} H^\top + \Sigma_w^\star\right)^{-1} \quad \text{and} \quad
	b^\star= \wh{\m}_x - A^\star (H \wh{\m}_x + \wh{\m}_w)
	\end{equation}
	solves the Wasserstein MMSE estimation problem~\eqref{eq:dro}.
\end{corollary}

\begin{proof}[Proof of Corollary~\ref{corol:alternative}]
	Define $\psi^\star$ as the affine estimator that solves~\eqref{eq:dro:approx} and $\mathbb{Q}^\star$ as the normal distribution that solves~\eqref{eq:dual-dro:conservative}. By Corollary~\ref{corol:nash}, the second inequality in~\eqref{eq:nash} holds for all admissible estimators $\psi\in{\mc F}$, which implies that $\psi^\star \in \arg\min_{\psi\in{\mc F}} \mathcal{R}(\psi,\mathbb{Q}^\star)$, that is, $\psi^\star$ solves the Bayesian MMSE estimation problem corresponding to $\mathbb{Q}^\star$. As any Bayesian MMSE estimator satisfies $\psi^\star(y)=\mathds{E}_{\mathbb{Q}^\star_{x|y}} [x]$ for $\mathbb{Q}^\star$-almost all $y$ and as $\Sigma^\star_w\succ 0$, we may use the known formulas for conditional normal distributions to conclude that the unique affine Bayesian MMSE estimator for $\mathbb{Q}^\star$ is of the form $\psi^\star(y) = A^\star y + b^\star$ with parameters defined as in~\eqref{eq:alternative:def}.
\end{proof}

\section{Non-normal Nominal Distributions}
\label{sect:elliptical}

We will first show that the results of Sections~\ref{sect:approx}--\ref{sect:nash} remain valid if $\wh \PP$ is an arbitrary elliptical (but maybe non-normal) distribution. To this end, we first review some basic results on elliptical distributions.

\begin{definition}[Elliptical distributions]
	\label{definition:elliptical}
	The distribution $\mathbb{P}$ of $\xi\in\mathbb{R}^d$ is called elliptical if the
	characteristic function $\Phi_\mathbb{P}(t)= \mathds{E}_\mathbb{P}[ \exp(i t^\top \xi)]$
	of $\mathbb{P}$ is given by $\Phi_\mathbb{P}(t)=\exp(i t^\top \mu) \phi(t^\top S
	t)$ for some location parameter ${\mu \in \mathbb{R}^d}$, dispersion matrix
	$S \in \mathbb{S}_{+}^{d}$ and characteristic generator $\phi:\mathbb{R}_+\! \rightarrow
	\mathbb{R}$. In this case we write~${\mathbb{P}=\mathcal{E}^d_\phi(\mu, S)}$. The class of
	all $d$-dimensional elliptical distributions with characteristic
	generator $\phi$ is denoted by~$\mathcal{E}^{d}_\phi$.
\end{definition}

The class of elliptical distributions was introduced
in~\cite{ref:Kelker-1970} with the aim to generalize the family of
normal distributions, which are obtained by setting the characteristic
generator to~$\phi(u)=e^{-u/2}$. We emphasize that, unlike the
moment-generating function $M_\mathbb{P}(t)= \mathds{E}_\mathbb{P}[ \exp(t^\top \xi)]$,
the characteristic function $\Phi_\mathbb{P}(t)$ is always finite for all
$t\in\mathbb{R}^d$ even if some moments of $\mathbb{P}$ do not exist. Thus,
Definition~\ref{definition:elliptical} is general enough to cover also
heavy-tailed distributions with non-zero tail dependence coefficients
\cite{ref:hult2002advances}. Examples of elliptical distributions include the Laplace, logistic and $t$-distribution etc. Useful theoretical properties
of elliptical distributions are discussed in~\cite{ref:cambanis1981theory,
	ref:fang1990symmetric}. We also highlight that elliptical distributions are
central to a wide spectrum of diverse applications ranging from genomics
\cite{ref:posekany2011biological} and medical imaging \cite{ref:ruttimann1998statistical}
to finance \cite[\S~6.2.1]{ref:jondeau2007financial}, to name a few.

If the dispersion matrix $S\in\mathbb{S}_{+}^d$ has rank $r$, then there
exists $\Lambda\in\mathbb{R}^{d\times r}$ with $S=\Lambda\Lambda^\top$, and
there exists a generalized inverse $\Lambda^{-1}\in\mathbb{R}^{r\times d}$ with
$\Lambda^{-1} \Lambda=I_r=\Lambda^\top(\Lambda^{-1})^\top$. One easily
verifies that if $\xi\in\mathbb{R}^d$ follows an elliptical distribution
$\mathbb{P}=\mathcal{E}^d_\phi(\mu, S)$, then $\tilde \xi=
\Lambda^{-1}(\xi-\mu)\in\mathbb{R}^r$ follows the spherically symmetric
distribution $\tilde\mathbb{P}=\mathcal{E}^d_\phi(0, I_r)$ with characteristic function
$\Phi_{\tilde\mathbb{P}}(t)=\phi(\|t\|^2)$. Thus, the choice of the
characteristic generator $\phi$ is constrained by the implicit condition
that $\phi(\|t\|^2)$ must be an admissible characteristic function. For
instance, the normalization of probability distributions necessitates
that $\phi(0) = 1$, while the dominated convergence theorem implies that
$\phi$ must be continuous etc. As any distribution is uniquely
determined by its characteristic function, and as $\phi(\|t\|^2)$
depends only on the norm of~$t$, the spherical distribution $\tilde \mathbb{P}$
is indeed invariant under rotations. This implies that $\mathds{E}_{\tilde
	\mathbb{P}}[\tilde \xi]=0$ and, via the linearity of the expectation, that
$\mathds{E}_{\mathbb{P}}[\xi]=\mu$ provided that $\tilde\xi$ and $\xi$ are
integrable, respectively. Thus, the location parameter $\mu$ of an
elliptical distribution coincides with its mean vector whenever the mean
exists. By the definition of the characteristic function, the covariance
matrix of $\tilde \mathbb{P}$, if it exists, can be expressed as
\[
\tilde \Sigma=-\left.\nabla_t^2\Phi_{\tilde \mathbb{P}}(t)\right|_{t=0} =
-\left.\nabla_t^2\phi(\|t\|^2)\right|_{t=0}= -2 \phi'(0) I_r,
\]
where $\phi'(0)$ denotes the right derivative of $\phi(u)$ at $u=0$.
Hence, $\tilde \Sigma$ exists if and only if $\phi'(0)$ exists and is
finite. Similarly, the covariance matrix of $\mathbb{P}$ is given by $\Sigma=-2
\phi'(0)S$, if it exists \cite[Theorem~4]{ref:cambanis1981theory}. Below
we will focus on elliptical distributions with finite first- and
second-order moments ({\em i.e.}, we will only consider characteristic
generators with $|\phi'(0)|<\infty$), and we will assume that
$\phi'(0)=-\frac{1}{2}$, which ensures that the dispersion matrix $S$
equals the covariance matrix $\Sigma$. The latter assumption does not
restrict generality. In fact, changing the characteristic generator to
$\phi(\frac{-u}{2\phi'(0)})$ and the dispersion matrix to $-2
\phi'(0)S$ has no impact on the elliptical distribution $\mathbb{P}$ but
matches the dispersion matrix $S$ with the covariance matrix $\Sigma$.

The elliptical distributions inherit many desirable properties from the normal distributions but are substantially more expressive as they include also heavy- and light-tailed distributions. For example, any class of elliptical distributions with a common characteristic generator is closed under affine transformations and affine conditional expectations; see {\em e.g.}, \cite[Theorem~1 and Corollary~5]{ref:cambanis1981theory}. Moreover, the Wasserstein distance between two elliptical distributions with the same characteristic generator equals the Gelbrich distance between their mean vectors and covariance matrices \cite[Theorem~2.4]{ref:gelbrich1990formula}. Thus, the Propositions~\ref{prop:normal-affine}, \ref{prop:normal-cond-exp} and~\ref{prop:normal:distance} extend verbatim from the class of normal distributions to {\em any} class of elliptical distributions that share the same characteristic generator. For the sake of brevity, we do not restate these results for elliptical distributions.

The above discussion suggests that the results of Sections~\ref{sect:approx}--\ref{sect:nash} carry over almost verbatim to MMSE estimation problems involving elliptical nominal distributions. In the following we will therefore assume that
\begin{align}
	\label{eq:nominal:elliptical}
	\wh \PP=\mathcal{E}^{n+m}_\phi(\wh{\m}, \covsa) \qquad \text{with}
	\qquad \widehat \mu = \begin{bmatrix} \widehat \mu_x \\ \widehat \mu_w
	\end{bmatrix} \qquad \text{and} \qquad \covsa = \begin{bmatrix}
		\covsa_{x} & 0 \\ 0 & \covsa_{w} \end{bmatrix},
\end{align}
where $\phi$ denotes a prescribed characteristic generator. As the class of all elliptical distributions with characteristic generator $\phi$ is closed under affine transformations, the marginal distributions $\wh \PP_x$ and $\wh \PP_w$ of $x$ and $w$ under $\wh \PP$ are also elliptical distributions with the same characteristic generator $\phi$.

Note that while the signal~$x$ and the noise~$w$ are uncorrelated under $\wh \PP$ irrespective of~$\phi$, they fail to be independent unless~$\wh \PP$ is a normal distribution. When working with generic elliptical nominal distributions, we must therefore abandon any independence assumptions.
Otherwise, the ambiguity set would be empty for small radii $\rho_x$ and $\rho_w$. This insight prompts us to redefine the Wasserstein ambiguity set as
\begin{equation}
\label{eq:Ambi-elliptical}
\mathbb B(\wh \PP) = \left\{ \mathbb{Q}\in \mathcal{M}(\mathbb{R}^{n+m}):
\mathds{E}_\mathbb{Q}[xw^\top]= \mathds{E}_\mathbb{Q}[x]\cdot\mathds{E}_\mathbb{Q}[w]^\top,~\mathds{W}(\mathbb{Q}_x, \wh \PP_x) \leq \rho_x,~ \mathds{W}(\mathbb{Q}_w, \wh \PP_w) \leq \rho_w \right\},
\end{equation}
which relaxes the independence condition in~\eqref{eq:Ambi} and merely requires $x$ and $w$ to be uncorrelated. When using the new ambiguity set~\eqref{eq:Ambi-elliptical} to model the distributional uncertainty, we can again compute a Nash equilibrium between the statistician and nature by solving a tractable convex optimization problem.

\begin{theorem}[Elliptical distributions] \label{thm:nash-elliptical}
	Assume that $\wh \PP$ is an elliptical distribution of the form~\eqref{eq:nominal:elliptical} with characterisic generator $\phi$ and noise covariance matrix $\covsa_w \succ 0$, and define the ambiguity set $\mathbb B(\wh \PP)$ as in~\eqref{eq:Ambi-elliptical}. If $(\Sigma_x^\star, \Sigma_w^\star)$ solves the finite convex program~\eqref{eq:dual-dro:conservative}, then the affine estimator $\psi^\star(y) = A^\star y + b^\star$ with
	\[
	A^\star =  \Sigma_x^\star H^\top \left( H {\Sigma_x^\star} H^\top + \Sigma_w^\star\right)^{-1} \quad \text{and} \quad
	b^\star= \wh{\m}_x - A^\star (H \wh{\m}_x + \wh{\m}_w)
	\]
	solves the Wasserstein MMSE estimation problem~\eqref{eq:dro}, while the elliptical distribution
	\[
	\mathbb{Q}^\star=\mathcal E_\phi^{n+m}(\wh{\m},\Sigma^\star)\quad \text{with} \quad \Sigma^\star= \begin{bmatrix}
	\Sigma_{x}^\star & 0 \\ 0 & \Sigma_{w^\star} \end{bmatrix}
	\]
	solves the dual Wasserstein MMSE estimation problem~\eqref{eq:dual-dro}. Moreover, $\psi^\star$ and $\mathbb{Q}^\star$ form a Nash equilibrium for the game between the statistician and nature, that is,
	\[
	\mathcal{R}(\psi^\star, \mathbb{Q}) \le \mathcal{R}(\psi^\star, \mathbb{Q}^\star) \le \mathcal{R}(\psi, {\mathbb{Q}}^\star) \quad \forall \psi \in {\mc F},~ \mathbb{Q} \in  \mathbb B(\wh \PP) \,.
	\]
\end{theorem}
\begin{proof}[Proof of Theorem~\ref{thm:nash-elliptical}]
	The proof replicates the arguments used to establish Theorems~\ref{thm:conservative}, \ref{thm:least-favorable-prior} and \ref{thm:sandwich} as well as Corollary~\ref{corol:nash} with obvious minor modifications. Details are omitted for brevity.
\end{proof}

Theorem~\ref{thm:nash-elliptical} asserts that the optimal estimator depends only on the first and second moments of the nominal elliptical distribution~$\wh \PP$ but {\em not} on its characteristic generator. Whether~$\wh \PP$ displays heavier or lighter tails than a normal distribution has therefore no impact on the prediction of the signal. Note, however, that the characteristic generator of~$\wh \PP$ determines the shape of the least favorable prior.

If the nominal distribution fails to be elliptical, the minimum of the Gelbrich MMSE estimation problem~\eqref{eq:dro:approx} may strictly exceed the maximum of the dual Wasserstein MMSE estimation problem~\eqref{eq:dual-dro:conservative} over normal priors. Note that in this case the ambiguity set $\mathbb B_{\mathcal N}(\wh \PP)$ may even be empty. Moreover, while typically suboptimal for the original Wasserstein MMSE estimation problem~\eqref{eq:dro}, the usual affine estimator constructed from a solution of the nonlinear SDP~\eqref{eq:program:dual} still solves the Gelbrich MMSE estimation problem~\eqref{eq:dro:approx}.
\begin{proposition}[Non-elliptical nominal distributions] \label{prop:equivalence}
    	Suppose that $\wh \PP = \wh \PP_{x} \times \wh \PP_{w}$, where $\wh \PP_x$ and $\wh \PP_w$ are arbitrary signal and noise distributions with mean vectors $\wh{\m}_x$ and $\wh{\m}_w$ and covariance matrices $\covsa_x \succeq 0$ and $\covsa_w \succ 0$, respectively. Then, the nonlinear SDP~\eqref{eq:program:dual} is solvable, and for any optimal solution~$(\Sigma_x^\star, \Sigma_w^\star)$ of~\eqref{eq:program:dual} the affine estimator $\psi^\star(y) = A^\star y + b^\star$ with
	\[
	A^\star =  \Sigma_x^\star H^\top \left( H {\Sigma_x^\star} H^\top + \Sigma_w^\star\right)^{-1} \quad \text{and} \quad
	b^\star= \wh{\m}_x - A^\star (H \wh{\m}_x + \wh{\m}_w)
	\]
	solves the Gelbrich MMSE estimation problem~\eqref{eq:dro:approx}.
\end{proposition}
\begin{proof}[Proof of Proposition~\ref{prop:equivalence}]
Denote by~$\wh \PP' = \wh \PP'_{x} \times \wh \PP'_{w}$ the normal distribution with the same first and second moments as~$\wh \PP$. As $\covsa_w\succ 0$, the nonlinear SDP~\eqref{eq:program:dual} is then solvable by virtue of Theorem~\ref{thm:least-favorable-prior}. Theorem~\ref{thm:sandwich} further implies that the first inequality in~\eqref{eq:strong:duality} with $\wh \PP'$ instead of $\wh \PP$ collapses to the equality
    \begin{equation}
        \label{eq:gebrich-meets-wasserstein}
        \inf\limits_{\psi\in\Ac} \sup\limits_{\mathbb{Q} \in \mathbb G(\wh \PP')} \mathcal{R}(\psi, \mathbb{Q})= \inf\limits_{\psi\in{\mc F}} \sup\limits_{\mathbb{Q} \in \mathbb B(\wh \PP')} \mathcal{R}(\psi, \mathbb{Q}).
    \end{equation}
    In addition, Corollary~\ref{corol:alternative} ensures that the affine estimator $\psi^\star$ defined in the proposition statement solves the modified Wasserstein MMSE estimation problem with normal nominal distribution $\wh \PP'$ on the right hand side of~\eqref{eq:gebrich-meets-wasserstein}. Because $\psi^\star$ is affine, it is also feasible in the modified Gelbrich MMSE estimation problem on the left hand side. In addition, the average risk of any affine estimator depends only on the mean vectors and covariance matrices of $x$ and $w$. If we denote by $T$ the mean-covariance projection that maps any distribution $\mathbb{Q} = \mathbb{Q}_x \times \mathbb{Q}_w$ of $(x,w)$ to the mean vectors and covariance matrices of $x$ and $w$ under $\mathbb{Q}$, then the images of the ambiguity sets $\mathbb G(\wh \PP')$ and $\mathbb B(\wh \PP')$ under $T$ coincide by Proposition~\ref{prop:normal:distance}.
    These observations imply that the affine estimator $\psi^\star$ also solves the estimation problem on the left hand side of~\eqref{eq:gebrich-meets-wasserstein}.
    As $\wh \PP$ and $\wh \PP'$ share the same first and second moments, the Gebrich ball ${\mbb G}(\wh \PP)$ around the generic distribution~$\wh \PP$ coincides with the Gelbrich ambiguity set ${\mbb G}(\wh \PP')$ around the normal distribution~$\wh \PP'$. Thus, we find
    \begin{align*}
          \sup\limits_{\mathbb{Q} \in {\mbb G}(\wh \PP)} \mathcal{R}(\psi^\star, \mathbb{Q}) = \sup\limits_{\mathbb{Q} \in {\mbb G}(\wh \PP')} \mathcal{R}(\psi^\star, \mathbb{Q}) =  \inf\limits_{\psi\in\Ac} \sup\limits_{\mathbb{Q} \in \mathbb G(\wh \PP')} \mathcal{R}(\psi, \mathbb{Q}) =  \inf\limits_{\psi\in\Ac} \sup\limits_{\mathbb{Q} \in \mathbb G(\wh \PP)} \mathcal{R}(\psi, \mathbb{Q}),
    \end{align*}
    where the second equality holds because $\psi^\star$ solves the Gelbrich MMSE estimation problem with the normal nominal distribution $\wh \PP'$ on the left hand side of~\eqref{eq:gebrich-meets-wasserstein}. Hence, $\psi^\star$ solves the Gelbrich MMSE estimation problem~\eqref{eq:dro:approx} with the generic nominal distribution $\wh \PP$.
\end{proof}

\section{Numerical Solution of Wasserstein MMSE Estimation Problems}
\label{sect:algorithm}

By Corollaries~\ref{cor:primal:refor} and~\ref{corol:dual:refor}, the primal and dual Wasserstein MMSE estimation problems~\eqref{eq:dro} and~\eqref{eq:dual-dro} can be addressed with off-the-shelf SDP solvers. Unfortunately, however, general-purpose interior-point methods quickly run out of memory when the signal dimension~$n$ and the noise dimension~$m$ grow. It is therefore expedient to look for customized first-order algorithms that can handle larger problem instances.

In this section we develop a Frank-Wolfe method for the nonlinear SDP~\eqref{eq:program:dual}, which is equivalent to the dual Wasserstein MMSE estimation problem~\eqref{eq:dual-dro}. This approach is meaningful because any solution to~\eqref{eq:program:dual} allows us to construct both an optimal estimator as well as a least favorable prior that form a Nash equilibrium in the sense of Corollary~\ref{corol:nash}; see also Corollary~\ref{corol:alternative}. Addressing the nonlinear SDP~\eqref{eq:program:dual} directly with a Frank-Wolfe method has great promise because the subproblems that identify the local search directions can be shown to admit quasi-closed form solutions and can therefore be solved very quickly.

In Section~\ref{sec:fully-adaptive-FW} we first review three variants of the  Frank-Wolfe algorithm corresponding to a static, an adaptive and a more flexible {\em fully} adaptive stepsize rule, and we prove that the fully adaptive rule offers a linear convergence guarantee under standard regularity conditions. In Section~\ref{sect:FW-MMSE} we then show that the nonlinear SDP~\eqref{eq:program:dual} is amenable to the fully adaptive Frank-Wolfe algorithm and can thus be solved efficiently.

\subsection{Frank-Wolfe Algorithm for Generic Convex Optimization Problems}
\label{sec:fully-adaptive-FW}
Consider a generic convex minimization problem of the form
\begin{equation}
\label{eq:generic-convex}
f^\star = \min_{s \in \mathcal S}~f(s)
\end{equation}
with a convex compact feasible set $\mathcal S \subseteq \mathbb{R}^d$ and a convex differentiable objective function $f:\mathcal S\rightarrow \mathbb R$. We assume that for each precision $\delta \in [0, 1]$ we have access to an inexact oracle $F: \mathcal S \rightarrow \mathcal S$ that maps any~$s\in\mathcal S$ to a $\delta$-approximate solution of an auxiliary problem linearized around~$s$. More precisely, we assume that
\begin{align}\label{oracle}
	\left( F(s) - s \right)^\top \nabla f(s) \leq \delta \min\limits_{z \in \mathcal S}~ \left( z - s \right)^\top \nabla f(s).
\end{align}
By the standard optimality condition for convex optimization problems, the minimum on the right hand side of~\eqref{oracle} vanishes if and only if~$s$ solves the original problem~\eqref{eq:generic-convex}. Otherwise, the minimum is strictly negative. If $\delta = 1$, then the oracle returns an exact mininizer of the linearized problem. If $\delta=0$, on the other hand, then the oracle returns any solution that is weakly preferred to $s$ in the linearized problem. Given an oracle satisfying~\eqref{oracle}, one can design a Frank-Wolfe algorithm whose iterates obey the recursion
\begin{align}\label{frank-wolfe}
	s_{t+1} = s_t + \eta_t (F(s_t) - s_t) \quad \forall t\in\mathbb N\cup\{0\},
\end{align}
where $s_0\in\mathcal S$ is an arbitrary initial feasible solution, $\delta$ is a prescribed precision, and $\eta_t \in[0,1]$ is a stepsize that may depend on the current iterate~$s_t$.
The Frank-Wolfe algorithm was originally developed for quadratic programs~\cite{ref:frank1956algorithm} and later extended to general convex programs with differentiable objective functions and compact convex feasible sets \cite{ref:levitin1966constrained, ref:demyanov1970approximate, ref:dunn1978conditional, ref:dunn1979rates, ref:dunn1980convergence}. Convergence guarantees for the Frank-Wolfe algorithm typically rely on the assumption that the gradient of $f$ is Lipschitz continuous \cite{ref:levitin1966constrained, ref:dunn1979rates, ref:dunn1980convergence, ref:garber2015faster, ref:freund2016new}, that $f$ has a bounded curvature constant \cite{ref:clarkson2010coresets, ref:jaggi2013revisiting}, or that the gradient of $f$ is H\"older continuous \cite{ref:nesterov2018complexity}.

Throughout this section we will assume that the decision variable can be represented as $s=(s^{[1]},\ldots, s^{[K]})$, where $s^{[k]}\in\mathbb R^{d_k}$ and $\sum_{k=1}^{K} d_k=d$. Moreover, we will assume that the feasible set~$\mathcal S = \times_{k=1}^K \mathcal S^{[k]}$ constitutes a $K$-fold Cartesian product, where the marginal feasible set~$\mathcal S^{[k]}\subseteq \mathbb R^{d_k}$ is convex and compact for each $k=1,\ldots, K$. This assumption is unrestrictive because we are free to set $K=1$ and $\mathcal S^{[1]}=\mathcal S$. For ease of notation, we use from now on $\nabla_{[k]}$ to denote the partial gradient with respect to the subvector~$s^{[k]}\in\mathcal S^{[k]}$,~$k=1,\ldots,K$.

The subsequent convergence analysis will rely on the following regularity conditions.

\begin{assumption} [Regularity conditions]
	\label{a:FW} ~
	\begin{enumerate}[label = $(\roman*)$]
		\item \label{a:FW:smooth}
		The objective function is $\beta$-smooth for some $\beta>0$, i.e.,
		\[
		\|\nabla f(s) - \nabla f(\bar s)\| \leq \beta \|s - \bar s\|\quad \forall s,\bar s\in\mathcal S.
		\]

		\item \label{a:FW:set}
		The marginal feasible sets are $\alpha$-strongly convex with respect to $f$ for some $\alpha>0$, i.e.,
		\[
		\theta s^{[k]} + (1-\theta) \bar s^{[k]} - \theta(1-\theta)\frac{\alpha}{2}\left\|s^{[k]} - \bar s^{[k]} \right\|^2 \frac{\nabla_{[k]} f(s)}{\|\nabla_{[k]} f(s)\|} \in \mathcal S^{[k]} \quad \forall s, \bar s \in \mathcal S,~\theta \in [0,1],~k=1,\ldots, K.
		\]

		\item \label{a:FW:lower}
		The objective function is $\varepsilon$-steep for some $\varepsilon > 0$, i.e.,
		\[
		\| \nabla_{[k]} f(s) \| \ge \varepsilon\quad \forall s \in \mathcal S,~k=1,\ldots, K.
		\]
	\end{enumerate}
\end{assumption}

Assumption~\ref{a:FW}\,\ref{a:FW:set} relaxes the standard strong convexity condition prevailing in the literature, which is obtained by setting $K=1$ and requiring that the condition stated here remains valid when the normalized gradient~${\nabla_{[1]} f(s) / \|\nabla_{[1]} f(s)\|}$ is replaced with any other vector in the Euclidean unit ball, see, {\em e.g.},~\cite[Equation~(25)]{ref:journee2010generalized}. We emphasize that our weaker condition is sufficient for the standard convergence proofs of the Frank-Wolfe algorithm but is necessary for our purposes because the feasible set of problem~\eqref{eq:program:dual} fails to be strongly convex in the traditional sense. Similarly, Assumption~\ref{a:FW}\,\ref{a:FW:lower} generalizes the usual $\varepsilon$-steepness condition from the literature, which is recovered by setting~$K=1$, see, {\em e.g.}, \cite[Assumption~1]{ref:journee2010generalized}. Under this assumption the gradient never vanishes on $\mathcal S$, and the minimum of~\eqref{eq:generic-convex} is attained on the boundary of~$\mathcal S$.

In the following we will distinguish three variants of the Frank-Wolfe algorithm with different stepsize rules. The {\em vanilla Frank-Wolfe} algorithm employs the harmonically decaying static stepsize
\begin{equation*}
	\eta_t = \frac{2}{2 + t},
\end{equation*}
which results in a sublinear $\mathcal O(1/t)$ convergence whenever Assumption~\ref{a:FW}\,\ref{a:FW:smooth} holds \cite{ref:frank1956algorithm, ref:dunn1978conditional}. The {\em adaptive Frank-Wolfe} algorithm uses the stepsize
\begin{align}\label{adp_step}
	\eta_t = \min \left\{1, \frac{(s_t - F (s_t))^\top \nabla f(s_t)}{\beta \| s_t - F(s_t) \|^2 } \right\},
\end{align}
which adapts to the iterate $s_t$. If all of the Assumptions~\ref{a:FW}\,\ref{a:FW:smooth}--\ref{a:FW:lower} hold, then the adaptive Frank-Wolfe algorithm enjoys a linear~$\mathcal O (c^t)$ convergence guarantee, where $c\in(0,1)$ is an explicit function of the oracle precision~$\delta$, the smoothness parameter~$\beta$, the strong convexity parameter~$\alpha$ and the steepness parameter~$\varepsilon$~\cite{ref:levitin1966constrained,ref:garber2015faster}. Note that the stepsize~\eqref{adp_step} is constructed as the unique solution of the univariate quadratic program
\begin{equation*}
	\min_{\eta \in [0, 1]} ~ f(s_{t}) - \eta \big(s_{t} - F(s_{t})\big)^\top \nabla f(s_{t}) + \frac{1}{2}{\beta \eta^2} \left\| s_{t} - F(s_{t}) \right\|^2,
\end{equation*}
which minimizes a quadratic majorant of the objective function $f$ along the line segment from $s_t$ to $F(s_t)$.

The adaptive stepsize rule~\eqref{adp_step} has undergone further scrutiny in \cite{ref:pedregosa2018stepsize}, where it was discovered that one may improve the algorithm's convergence behavior by replacing the global smoothness parameter~$\beta$ in~\eqref{adp_step} with an adaptive smoothness parameter~$\beta_t$ that captures the smoothness of $f$ along the line segment from~$s_t$ to~$F(s_t)$. This extra flexibility is useful because~$\beta_t$ can be chosen smaller than the unnecessarily conservative global smoothness parameter~$\beta$ and because $\beta_t$ is easier to estimate than~$\beta$, which may not even be accessible.

Following~\cite{ref:pedregosa2018stepsize}, we will henceforth only require that $\beta_t>0$ satisfies the inequality
\begin{align} \label{eq:stepsize_full}
	f\Big( s_t - \eta_t(\beta_t) \big(s_t - F(s_t)\big) \Big) \le f(s_t) - \eta_t(\beta_t) \big(s_t - F(s_t)\big)^\top \nabla f(s_t) + \frac{1}{2}{\beta_t \eta_t(\beta_t)^2} \big\|s_t - F(s_t)\big\|^2,
\end{align}
where $\eta_t(\beta_t)$ is defined as the adaptive stepsize~\eqref{adp_step} with $\beta$ replaced by $\beta_t$. As it adapts both to $s_t$ and $\beta_t$, we will from now on refer to $\eta_t=\eta_t(\beta_t)$ as the {\em fully} adaptive stepsize. The above discussion implies that~\eqref{eq:stepsize_full} is always satisfiable if Assumption~\ref{a:FW}\,\ref{a:FW:smooth} holds, in which case one may simply set~$\beta_t$ to the global smoothness parameter~$\beta$. In practice, however, the inequality~\eqref{eq:stepsize_full} is often satisfiable for much smaller values~$\beta_t \ll \beta$ that may not even be related to the smoothness properties of the objective function. A close upper bound on the smallest~$\beta_t>0$ that satisfies~\eqref{eq:stepsize_full} can be found efficiently via backtracking line search. Specifically, the {\em fully adaptive Frank-Wolfe} algorithm sets $\beta_t$ to the smallest element of the discrete search space $\frac{\beta_{t-1}}{\zeta}\cdot\{1, \tau, \tau^2, \tau^3,\ldots\}$ that satisfies~\eqref{eq:stepsize_full}, where  $\tau>1$ and $\zeta>1$ are prescribed line search parameters. A detailed description of the fully adaptive Frank-Wolfe algorithm in pseudocode is provided in Algorithm~\ref{algorithm:FAFW}.

It has been shown in~\cite{ref:pedregosa2018stepsize} that Algorithm~\ref{algorithm:FAFW} enjoys the same sublinear~$\mathcal O(1/t)$ convergence guarantee as the vanilla Frank-Wolfe algorithm when Assumption~\ref{a:FW}\,\ref{a:FW:smooth} holds. Below we will leverage techniques from~\cite{ref:levitin1966constrained, ref:garber2015faster} to show that Algorithm~\ref{algorithm:FAFW} offers indeed a linear convergence rate if all of the Assumptions~\ref{a:FW}\,\ref{a:FW:smooth}--\ref{a:FW:lower} hold.

\begin{table}[th]
	\begin{minipage}{0.71\columnwidth}
		\begin{algorithm}[H]
			\caption{Fully adaptive Frank-Wolfe algorithm}
			\label{algorithm:FAFW}
			\begin{algorithmic}
				\REQUIRE initial feasible point $s_0 \in \mathcal S$, initial smoothness parameter $\beta_{-1} > 0$ \\
				\hspace{2.2em}
				line search parameters $\tau > 1$, $\zeta > 1$, initial iteration counter $t=0$ \\[0.5ex]
				\WHILE{stopping criterion is not met} \vspace{0.25em}
				\STATE solve the oracle subproblem to find $\tilde s_t = F(s_t)$
				\STATE set $d_t \leftarrow \tilde s_{t} - s_t$ and $ g_t \leftarrow -d_t^\top{\nabla f(s_{t})}$ \vspace{0.1em}
				\STATE set $\beta_t \leftarrow \beta_{t-1} / \zeta$ and $\eta \leftarrow \min \{1, g_t/(\beta_t \| d_t \|^2) \} $
				\WHILE{$ \displaystyle f(s_t + \eta d_t) > f(s_t) - \eta g_t + \frac{\eta^2 \beta_t}{2} \| d_t \|^2 $}
				\STATE $\beta_t \leftarrow \tau \beta_t$ and $\eta \leftarrow \min \{1, g_t/(\beta_t \| d_t \|^2) \}$
				\ENDWHILE
				\STATE set $\eta_t \leftarrow \eta$ and $s_{t+1} \leftarrow s_t + \eta_t d_t$
				\STATE set $t \leftarrow t + 1$
				\ENDWHILE \vspace{0.5ex}
				\ENSURE $s_t$
			\end{algorithmic}
		\end{algorithm}
	\end{minipage}
\end{table}

\begin{theorem}[Linear convergence of the fully adaptive Frank-Wolfe algorithm]\label{theorem:algorithm}
	If Assumption~\ref{a:FW} holds and $\overline{\beta} = \max \{ \tau \beta, \beta_{-1} \}$, then Algorithm~\ref{algorithm:FAFW} enjoys the linear convergence guarantee
	\[f(s_t) - f^\star \le \max \left\{ 1 - \frac{\delta}{2} , 1 - \frac{(1-\sqrt{1 - \delta})\alpha \varepsilon}{4 \overline \beta} \right\}^t (f(s_{0}) - f^\star) \quad \forall t\in\mathbb N.\]
\end{theorem}
The proof of Theorem~\ref{theorem:algorithm} relies on the following preparatory lemma.

\begin{lemma}[Bounds on the surrogate duality gap]
	\label{lem:FW}
	The surrogate duality gap~$g_t = -d_t^\top \nabla f(s_t)$ corresponding to the search direction~$d_t = F(s_t) - s_t $ admits the following lower bounds.
	\begin{enumerate}[label = $(\roman*)$,itemsep = 2mm]
		\item \label{lem:FW:convex}
		If the objective function $f$ is convex, then $g_t \ge \delta(f(s_t) - f^\star)$.
		\item \label{lem:FW:set}
		If the marginal feasible sets are $\alpha$-strongly convex with respect to $f$ for some $\alpha>0$ in the sense of Assumption~\ref{a:FW}\,\ref{a:FW:set}, then
		\[
		g_t \geq \min_{k \in \{ 1, \hdots, K \}} \frac{(1-\sqrt{1 - \delta}) \alpha} {2 \delta} \| d_t \|^2 \|\nabla_{[k]} f(s_t) \|.
		\]
	\end{enumerate}
\end{lemma}

\begin{proof}[Proof of Lemma~\ref{lem:FW}]
	By the definition of~$g_t$ we have
	\begin{align}\label{eq:lem:gk}
		g_t = \big(s_t - F(s_t)\big)^\top \nabla f(s_t) \ge \delta \big(s_t - s\big)^\top \nabla f(s_t) \qquad \forall s \in \mathcal S,
	\end{align}
	where the inequality follows from the defining property~\eqref{oracle} of the inexact oracle with precision~$\delta$. Setting~$s$ in~\eqref{eq:lem:gk} to a global minimizer~$s^\star$ of~\eqref{eq:generic-convex} then implies via the first-order convexity condition for~$f$ that
	$$g_t \ge \delta \big(s_t - s^\star\big)^\top \nabla f(s_t) \ge \delta\big(f(s_t) - f^\star\big)  \,.$$
	This observation establishes assertion~\ref{lem:FW:convex}. To prove assertion~\ref{lem:FW:set}, we first rewrite the estimate~\eqref{eq:lem:gk} as
	\begin{equation}
	\label{eq:cartesian}
	g_t \ge \delta (s_t - s )^\top \nabla f(s) = \delta \sum_{t=1}^T \big(s_t^{[k]} - s^{[k]} \big)^\top \nabla_{[k]} f(s_t) \quad \forall s \in \mathcal S \,.
	\end{equation}
	In the following, we denote by $F^{[k]}:\mathcal S\to \mathcal S^{[k]}$ the $k$-th suboracle for $k=1,\ldots, K$, which is defined through the identity $F=(F^{[1]},\ldots, F^{[K]})$. Similarly, for any $\theta\in[0,1]$, we define $s(\theta)=(s^{[1]}(\theta), \ldots,s^{[K]}(\theta))$ through
	\begin{align*}
		s^{[k]}(\theta) = \theta F^{[k]}(s_t) + (1-\theta)s_t^{[k]} - \frac{\alpha}{2}\theta(1-\theta)\left\|F^{[k]}(s_t) - s_t^{[k]}\right\|^2 \frac{\nabla_{[k]} f(s_t)}{\|\nabla_{[k]} f(s_t)\|}.
	\end{align*}
	By Assumption~\ref{a:FW}\,\ref{a:FW:set}, we have $s^{[k]}(\theta)\in\mathcal S^{[k]}$ for every $k=1,\ldots, K$. Thanks to the rectangularity of the feasible set this implies that $s(\theta)\in\mathcal S$. Setting $s$ in~\eqref{eq:cartesian} to $s(\theta)$, we thus find
	\begin{align*}
		g_t &\ge \delta \sum_{t=1}^T \Big(\theta \big(s_t^{[k]} - F^{[k]}(s_t)\big) + \frac{\alpha}{2} \theta(1-\theta)\left\|F^{[k]}(s_t) - s_t^{[k]}\right\|^2\frac{\nabla_{[k]} f(s_t)}{\|\nabla_{[k]} f(s_t)\|} \Big)^\top \nabla_{[k]} f(s_t)  \\
		& = \delta \left( \theta g_t + \frac{\alpha}{2} \theta (1-\theta) \left[ \sum_{k=1}^K \left\|F^{[k]}(s_t) - s_t^{[k]}\right\|^2 \|\nabla_{[k]} f(s_t)\| \right] \right) \\
		&\geq \delta \left( \theta g_t + \frac{\alpha}{2} \theta (1-\theta) \|F(s_t) - s_t \|^2 \Big(\min_{k \in \{ 1, \hdots, K \} } \|\nabla_{[k]} f(s_t)\| \Big) \right) \quad \forall \theta \in [0,1] \,,
	\end{align*}
	where the equality follows from the definition of~$g_t$, and the last inequality exploits the Pythagorean theorem. Reordering the above inequality to bring $g_t$ to the left hand side yields
	\begin{align}\label{eq:lem:gk2}
		g_t \ge \min_{k \in \{ 1, \hdots, K \} } \frac{\alpha}{2} \|F(s_t) - s_t\|^2 \|\nabla_{[k]} f(s_t)\| \frac{\delta \theta(1-\theta)}{1 - \delta \theta}  \qquad \forall \theta \in [0,1].
	\end{align}
	A tedious but straightforward calculation shows that the lower bound on the right hand side of~\eqref{eq:lem:gk2} is maximized by $\theta^\star=(1 - \sqrt{1 - \delta}) / \delta$. Assertion~\ref{lem:FW:set} then follows by substituting $\theta^\star$ into~\eqref{eq:lem:gk2}.
\end{proof}

\begin{proof}[Proof of Theorem~\ref{theorem:algorithm}]
	By Assumption~\ref{a:FW}\,\ref{a:FW:smooth} the function $f$ is $\beta$-smooth, and thus one can show that
	\begin{align}
		\label{eq:surrogate-objective}
		f(s_t + \eta d_t) \leq f(s_t) - \eta g_t + \frac{\eta^2 \beta}{2} \| d_t \|^2 \quad \forall \eta \in [0,1] \,,
	\end{align}
	where the surrogate duality gap~$g_t\ge 0$ and the search direction~$d_t\in\mathbb R^d$ are defined as in Lemma~\ref{lem:FW}. We emphasize that~\eqref{eq:surrogate-objective} holds in fact for all $\eta\in\mathbb R$. However, the next iterate $s_{t+1}=s_t + \eta d_t$ may be infeasible unless $\eta \in [0,1]$.  The inequality~\eqref{eq:surrogate-objective} implies that any $\beta_t\ge \beta$ satisfies the condition of the inner while loop of Algorithm~\ref{algorithm:FAFW}, and thus the loop must terminate at the latest after $\lceil \log(\zeta\beta/\beta_{-1})/\log(\tau)\rceil$ iterations, outputting a smoothness parameter~$\beta_t$ and a stepsize~$\eta_t$ that satisfy the inequality~\eqref{eq:stepsize_full}.
	We henceforth denote by $h_t = f(s_t) - f^\star$ the suboptimality of the $k$-th iterate and note that
	\begin{equation}
	\label{eq:h_t-recursion}
	h_{t+1} = f(s_t+\eta_t d_t)-f(s_t)+ h_t \le -g_t + \frac{1}{2}\beta_t\eta_t^2 \|d_t\|^2+h_t\,,
	\end{equation}
	where the inequality exploits~\eqref{eq:stepsize_full} and the definitions of $g_t$ and $d_t$. In order to show that $h_t$ decays geometrically, we distinguish the cases $(i)$ $g_t / (\beta_t \|d_t\|^2) \ge 1$ and $(ii)$ $g_t / (\beta_t \|d_t\|^2) < 1$. In case~$(i)$, the stepsize $\eta_t$ defined in~\eqref{eq:stepsize_full} satisfies $\eta_t = \min \{1, g_t / (\beta_t \|d_t\|^2) \} = 1$, and thus we have
	\begin{equation}
	\label{eq:h_t-recursion_i}
	h_{t+1} \leq \left( \frac{\beta_t \| d_t \|^2}{2 g_t} - 1 \right) g_t +h_t \leq -\frac{g_t}{2}+h_t \leq  \left( 1 - \frac{\delta}{2} \right) h_t,
	\end{equation}
	where the first inequality follows from~\eqref{eq:h_t-recursion}, while the third inequality holds due to Lemma~\ref{lem:FW}\,\ref{lem:FW:convex}.

	In case~$(ii)$, the stepsize satisfies~$\eta_t = g_t / \beta_t \|d_t\|^2 < 1$, and thus we find
	\begin{align}
		\nonumber
		h_{t+1} &\leq - g_t + \frac{g_t^2}{2\beta_t \| d_t \|^2} + h_t \leq - \frac{g_t^2}{2 \beta_t \| d_t \|^2} + h_t \leq \left(1 - \frac{\delta g_t}{2 \beta_t \| d_t \|^2}\right) h_t\\
		\label{eq:h_t-recursion_ii}
		&\le \left( 1 - \min_{k \in \{ 1, \hdots, K \} } \frac{(1-\sqrt{1 - \delta}) \alpha}{4 \beta_t} \|\nabla_{[k]} f(s_t)\|\right) h_t
		\le \left( 1 - \frac{(1-\sqrt{1 - \delta}) \alpha\varepsilon}{4 \overline\beta} \right) h_t\,,
	\end{align}
	where the first and the second inequalities follow from~\eqref{eq:h_t-recursion} and from multiplying $-g_t$ with $\eta_t<1$, respectively, while the third and the fourth inequalities exploit Lemmas~\ref{lem:FW}\,\ref{lem:FW:convex} and~\ref{lem:FW}\,\ref{lem:FW:set}, respectively. The last inequality in~\eqref{eq:h_t-recursion_ii} holds because of Assumption~\ref{a:FW}\,\ref{a:FW:lower} and because $\beta_t \leq \overline{\beta}$ for all $t\in\mathbb N$; see~\cite[Proposition~2]{ref:pedregosa2018stepsize}. By the estimates~\eqref{eq:h_t-recursion_i} and~\eqref{eq:h_t-recursion_ii}, the suboptimality of the current iterate decays at least~by
	\[
	\max \left\{ 1 - \frac{\delta}{2} , 1 - \frac{(1-\sqrt{1 - \delta})\alpha \varepsilon}{4 \overline \beta} \right\} <1
	\]
	in each iteration of the algorithm. This observation completes the proof.
\end{proof}

\subsection{Frank-Wolfe Algorithm for Wasserstein MMSE Estimation Problems}
\label{sect:FW-MMSE}

We now use the fully adaptive Frank-Wolfe algorithm of Section~\ref{sec:fully-adaptive-FW} to solve the nonlinear SDP~\eqref{eq:program:dual}, which is equivalent to the dual Wasserstein MMSE estimation problem over normal priors. Recall from Corollary~\ref{corol:alternative} that any solution of~\eqref{eq:program:dual} can be used to construct a least favorable prior and an optimal estimator that form a Nash equilibrium. Unlike the generic convex program~\eqref{eq:generic-convex}, the nonlinear SDP~\eqref{eq:program:dual} is a convex {\em maximization} problem. This prompts us to apply Algorithm~\ref{algorithm:FAFW} to the convex minimization problem obtained from problem~\eqref{eq:program:dual} by turning the objective function upside down.

Throughout this section we assume that $\covsa_x\succ 0$, $\covsa_w\succ 0$, $\rho_x>0$ and $\rho_w>0$, which implies via Theorem~\ref{thm:least-favorable-prior} that the nonlinear SDP~\eqref{eq:program:dual} is solvable and can be reformulated more concisely as
\begin{equation}
\label{eq:program:dual:concise}
\max_{\Sigma_x\in\mathcal S^+_x, \Sigma_w\in\mathcal S^+_w}~ f(\Sigma_x, \Sigma_w)\,,
\end{equation}
where the objective function $f:\mathcal S^+_x \times \mathcal S^+_w\to \mathbb R$ is defined through
\[
f(\Sigma_x, \Sigma_w) = \Tr{\Sigma_x - \Sigma_x H^\top \left( H \Sigma_x H^\top + \Sigma_w \right)^{-1} H \Sigma_x}\,,
\]
and where the separate feasible sets for $\Sigma_x$ and $\Sigma_w$ are given by
\[
\mathcal S^+_x = \left\{ \Sigma_x \in \mathbb{S}_{+}^n: \Tr{\Sigma_x + \covsa_x - 2 \big( \covsa_x^\frac{1}{2} \Sigma_x \covsa_x^\frac{1}{2} \big)^\frac{1}{2}} \leq \rho_x^2,~\Sigma_x\succeq \lambda_{\min}(\covsa_x)I_n \right\}
\]
and
\[
\mathcal S^+_w = \left\{ \Sigma_w \in \mathbb{S}_{+}^m: \Tr{\Sigma_w + \covsa_w - 2 \big( \covsa_w^\frac{1}{2} \Sigma_w \covsa_w^\frac{1}{2} \big)^\frac{1}{2}} \leq \rho_w^2,~\Sigma_w\succeq \lambda_{\min}(\covsa_w)I_m  \right\}\,,
\]
respectively. One readily verifies that $f$ is concave and differentiable. Moreover, in the terminology of Section~\ref{sec:fully-adaptive-FW}, the feasible set of the nonlinear SDP~\eqref{eq:program:dual:concise} constitutes a Cartesian product of $K=2$ marginal feasible sets~$\mathcal S^+_x$ and~$\mathcal S^+_w$, both of which are convex and compact thanks to Lemma~\ref{lemma:compact:FS}. Note that~$\mathcal S_x^+$ and~$\mathcal S_w^+$ constitute restrictions of the feasible sets~$\mathcal S_x$ and~$\mathcal S_w$, respectively, which appeared in the proofs of Theorems~\ref{thm:least-favorable-prior} and~\ref{thm:sandwich}. The oracle problem that linearizes the objective function of the nonlinear SDP~\eqref{eq:program:dual:concise} around a fixed feasible solution $\Sigma_x\in\mathcal S^+_x$ and $\Sigma^+_w\in\mathcal S_w$ can now be expressed concisely as
\begin{equation}
\label{eq:oracle}
\max_{L_x\in\mathcal S^+_x, L_w\in\mathcal S^+_w} \inner{L_x - \Sigma_x}{D_x} + \inner{L_w - \Sigma_w}{D_w}\,,
\end{equation}
where $D_x = \nabla_{\Sigma_x} f\big(\Sigma_x, \Sigma_w\big)$ and $D_w =\nabla_{\Sigma_w} f\big(\Sigma_x, \Sigma_w\big)$ represent the gradients of $f$ with respect to~$\Sigma_x$ and~$\Sigma_w$. Lemma~\ref{lem:Taylor-f} offers analytical formulas for $D_x$ and $D_w$ and shows that they are both positive semidefinite.

The oracle problem~\eqref{eq:oracle} is manifestly separable in $L_x$ and $L_w$ and can therefore be decomposed into a sum of two structurally identical marginal subproblems. The Frank-Wolfe algorithm is an ideal method to address the nonlinear SDP~\eqref{eq:program:dual:concise} because these two marginal oracle subproblems admit quasi-closed form solutions. Specifically, Proposition~\ref{prop:quadratic}\,\ref{prop:quad:cov_min} in the appendix implies that problem~\eqref{eq:oracle} is uniquely solved by
\[
L_x^\star = (\gamma_x^\star)^2 (\gamma_x^\star I_n - D_x)^{-1} \covsa_x (\gamma_x^\star I_n - D_x)^{-1} \quad \text{and}\quad
L_w^\star = (\gamma_w^\star)^2 (\gamma_w^\star I_m - D_w)^{-1} \covsa_w (\gamma_w^\star I_m - D_w)^{-1},
\]
where $\gamma_x^\star\in(\lambda_{\max}(D_x),\infty)$ and $\gamma_w^\star\in(\lambda_{\max}(D_w),\infty)$ are the unique solutions of the algebraic equations
\begin{equation}
\label{eq:algebraic-equations}
\rho_x^2 - \inner{\covsa_x}{\big( I_n - \gamma_x^\star (\gamma_x^\star I_n - D_x)^{-1} \big)^2} = 0\quad \text{and}\quad
\rho_w^2 - \inner{\covsa_w}{\big( I_m - \gamma_w^\star (\gamma_w^\star I_m - D_w)^{-1} \big)^2} = 0,
\end{equation}
respectively. In practice, these algebraic equations need to be solved numerically.
However, the numerical errors in~$\gamma_x^\star$ and~$\gamma_w^\star$ must be contained to ensure that~$L_x^\star$ and~$L_w^\star$ give rise to a $\delta$-approximate solution for~\eqref{eq:oracle} in the sense of~\eqref{oracle}. In the following we will show that $\delta$-approximate solutions for each of the two oracle subproblems in~\eqref{eq:oracle} and for each $\delta\in(0,1)$ can be computed with an efficient bisection algorithm.

\begin{table}[th]
	\begin{minipage}{0.60\columnwidth}
		\begin{algorithm}[H]
			\caption{Bisection algorithm for the oracle subproblem}
			\label{algorithm:bisection}
			\begin{algorithmic}
				\REQUIRE nominal covariance matrix $\covsa \in \mathbb{S}_{++}^d$, radius $\rho \in\mathbb R_{++}$, \\	\hspace{2.3em}  reference covariance matrix $\Sigma \in \mathbb{S}_{+}^d$ feasible in~\eqref{eq:oracle}, \\
				\hspace{2.3em} gradient matrix $D \in \mathbb{S}_{+}^d$, $D\neq 0$, precision $\delta \in (0, 1)$,\\
				\hspace{2.3em} dual objective function $\varphi(\gamma)$ defined in Theorem~\ref{thm:oracle}
				\vspace{0.21em}
				\STATE set $\lambda_1\leftarrow \lambda_{\max}(D)$, and let $v_1\in\mathbb R^d$ be an eigenvector for $\lambda_1$
				\STATE set $\underline\gamma \leftarrow \lambda_1 ( 1 + (v_1^\top \covsa v_1)^\frac{1}{2} / \rho ) $ and $\overline\gamma \leftarrow \lambda_1 ( 1 + \Tr{\covsa}^{\frac{1}{2}} / \rho ) $
				\REPEAT
				\STATE Set $\tilde \gamma \leftarrow (\overline\gamma + \underline\gamma) / 2$ and $\tilde L \leftarrow (\tilde \gamma)^2 (\tilde \gamma I_d - D)^{-1} \covsa (\tilde \gamma I_d - D)^{-1}$
				\STATE \textbf{if } $\frac{{\rm d} \varphi}{{\rm d}\gamma}(\tilde \gamma)<0$ \textbf{ then } Set $\underline\gamma \leftarrow \tilde \gamma$ \textbf{ else } Set $\overline\gamma \leftarrow \tilde \gamma$ \textbf{ endif}
				\UNTIL{$\frac{{\rm d} \varphi}{{\rm d}\gamma}(\tilde \gamma)>0$ and $\inner{\tilde L - \Sigma}{D}  \geq \delta\,\varphi(\tilde \gamma) }$ \vspace{0.21em}
				\ENSURE $\tilde L$
			\end{algorithmic}
		\end{algorithm}
	\end{minipage}
\end{table}

\begin{theorem}[Approximate oracle] \label{thm:oracle}
	For any fixed $\rho\in \mathbb R_{++}$,  $\covsa\in\mathbb{S}_{++}^d$ and $D \in\mathbb{S}_{+}^d$, $D\neq 0$, consider the generic oracle subproblem
	\begin{equation}
	\label{eq:oracle-primal}
	\begin{array}{cl}
		\displaystyle \max_{L \in \mathbb{S}_{+}^d} & \inner{L - \Sigma}{D}\\[-1.5ex]
		\st & \Tr{L + \covsa - 2 \big( \covsa^\frac{1}{2} L \covsa^\frac{1}{2} \big)^\frac{1}{2}} \leq \rho^2,~L\succeq \lambda_{\min}(\covsa)I_d\,,
	\end{array}
	\end{equation}
	where $\Sigma\in\mathbb{S}_{+}^d$ represents a feasible reference solution. Moreover, denote the feasible set of problem~\eqref{eq:oracle-primal} by~$\mathcal S^+$, let $\delta\in(0,1)$ be the desired oracle precision, and define $\varphi(\gamma)= \gamma(\rho^2 + \inner{\gamma(\gamma I_d - D)^{-1} - I_d}{\covsa})-\inner{\Sigma}{D}$ for any $\gamma>\lambda_{\max}(D)$. Then, Algorithm~\ref{algorithm:bisection} returns in finite time a matrix $\tilde L\in\mathbb{S}_{+}^d$ with the following~properties.
	\begin{enumerate}[label = $(\roman*)$,itemsep = 2mm]
		\item Feasibility: $\tilde L\in\mathcal S^+$
		\item $\delta$-Suboptimality: $\inner{\tilde L - \Sigma}{D} \geq \delta \max_{L \in \mathcal S^+} ~ \inner{L - \Sigma}{D}$
	\end{enumerate}
\end{theorem}
\begin{proof}[Proof of Theorem~\ref{thm:oracle}]
	Proposition~\ref{prop:quadratic}\,\ref{prop:quad:cov_min} in the appendix guarantees that the lower bound on $L$ in~\eqref{eq:oracle-primal} is redundant and can be omitted without affecting the problem's optimal value. By Proposition~\ref{prop:quadratic}\,\ref{prop:quad:dual}, the oracle subproblem~\eqref{eq:oracle-primal}
	thus admits the strong Lagrangian dual
	\[
	\min_{\gamma > \lambda_{\max}(D)}~\varphi(\gamma)\,,
	\]
	where the convex and differentiable function $\varphi(\gamma)$ is defined as in the theorem statement. In the following we denote by~$\lambda_1>0$ the largest eigenvalue of $D$ and let $v_1\in\mathbb R^d$ be a corresponding eigenvalue. By Proposition~\ref{prop:quadratic}\,\ref{prop:quad:cov_min}, the dual oracle subproblem admits a minimizer $\gamma^\star\in[\underline \gamma,\overline\gamma]$ that is uniquely determined by the first-order optimality condition $\frac{{\rm d} \varphi}{{\rm d}\gamma}(\gamma^\star)=0$, where
	\[
	\underline\gamma= \lambda_1\left( 1+\sqrt{v_1^\top\covsa v_1}/\rho \right)\quad \text{and}\quad\overline\gamma= \lambda_1\left( 1+\sqrt{\Tr{\covsa}}/\rho\right),
	\]
	while the primal problem~\eqref{eq:oracle-primal} admits a unique maximizer $L^\star=L(\gamma^\star)$, where
	\[
	L(\gamma) = \gamma^2 (\gamma I_d - D)^{-1} \covsa (\gamma I_d - D)^{-1}.
	\]
	From the proof of Proposition~\ref{prop:quadratic}\,\ref{prop:quad:cov_min} it is evident that $L(\gamma)\succ \lambda_{\min}(\covsa)I_d$ for every $\gamma>0$.
	A direct calculation further shows that
	\begin{align*}
		\frac{{\rm d}\varphi}{{\rm d}\gamma}(\gamma)=\rho^2 - \inner{\covsa}{\big( I_d - \gamma (\gamma I_d - D)^{-1} \big)^2}
		=\rho^2- \Tr{L(\gamma) + \covsa - 2 \big( \covsa^\frac{1}{2} L(\gamma) \covsa^\frac{1}{2} \big)^\frac{1}{2}}.
	\end{align*}
	Recalling that $\varphi(\gamma)$ is convex and that its derivative is non-negative for all~$\gamma\ge \gamma^\star$, the above reasoning implies that $L(\gamma)\in\mathcal S^+$ for all $\gamma\ge \gamma^\star$. Note also that the optimal value of the primal problem~\eqref{eq:oracle-primal} is non-negative because $\Sigma \in \mathcal S^+$. The continuity of $\inner{L(\gamma)-\Sigma}{D}$ at $\gamma=\gamma^\star$ thus ensures that there exists $\delta'>0$ with
	\[
	\inner{L(\gamma)- \Sigma}{D} \leq \delta \inner{L(\gamma^\star)- \Sigma}{D} = \delta \max_{L \in \mathcal S^+} ~ \inner{L - \Sigma}{D}\quad
	\forall \gamma\in[\gamma^\star,\gamma^\star+\delta'].
	\]
	In summary, computing a feasible and $\delta$-suboptimal matrix $\tilde L\in\mathbb{S}_{+}^d$ is tantamount to finding $\tilde \gamma\in[\gamma^\star,\gamma^\star+\delta']$. Algorithm~\ref{algorithm:bisection} uses bisection over the interval $[\underline \gamma, \overline\gamma]$ to find a $\tilde\gamma$ with these properties.
\end{proof}

Theorem~\ref{thm:oracle} complements~\cite[Theorem~3.2]{ref:shafieezadeh2018wasserstein}, which constructs an approximate oracle for a nonlinear SDP similar to~\eqref{eq:program:dual:concise} that offers an {\em additive} error guarantee. The {\em multiplicative} error guarantee of the oracle constructed here is needed to ensure the linear convergence of the fully adaptive Frank-Wolfe algorithm. Next, we prove that the nonlinear SDP~\eqref{eq:program:dual:concise} satisfies all regularity conditions listed in Assumption~\ref{a:FW}.

\begin{proposition}[Regularity conditions of the nonlinear SDP~\eqref{eq:program:dual:concise}] \label{prop:regularity}
	If $\rho_x \in \mbb R_{++}$, $\rho_w \in \mbb R_{++}$, $\covsa_x \in \mathbb{S}_{+}^n$ and $\covsa_w \in \mathbb{S}_{+}^m$, then the nonlinear SDP~\eqref{eq:program:dual:concise} obeys the following regularity conditions.

	\begin{enumerate}[label = $(\roman*)$,itemsep = 2mm]
		\item \label{prop:regularity-i} The objective function of problem~\eqref{eq:program:dual:concise} is $\beta$-smooth in the sense of Assumption~\ref{a:FW}\,\ref{a:FW:smooth}, where
		\[
		\beta = 2\lambda_{\min}^{-1}(\covsa_w) \left( C + C \; \lambda_{\max}^2(H^\top H) + \lambda_{\max}(H^\top H) \right),
		\]
		which depends on the auxiliary constant $C = \lambda_{\max} ( H^\top H ) \cdot \lambda_{\min}^{-2}(\covsa_w) \cdot ( \rho_x +  \Tr{\covsa_x}^\frac{1}{2} )^4$.

		\item \label{prop:regularity-ii} The marginal feasible sets $\mathcal S_x^+$ and $\mathcal S^+_w$ of problem~\eqref{eq:program:dual:concise} are $\alpha$-strongly convex with respect to $-f$ in the sense of Assumption~\ref{a:FW}\,\ref{a:FW:set}, where $\alpha = \min \, \{\alpha_x, \alpha_w \}$, which depends on the auxiliary constants
		\[
		\alpha_x = \frac{\lambda_{\min}^{\frac{5}{4}}(\covsa_x)}{2 \rho_x \big( \rho_x + \Tr{\covsa_x}^\frac{1}{2} \big)^{\frac{7}{2}} } \quad\text{and}\quad \alpha_w = 		\frac{\lambda_{\min}^{\frac{5}{4}}(\covsa_w)}{2 \rho_w \big( \rho_w + \Tr{\covsa_w}^\frac{1}{2} \big)^{\frac{7}{2}} }\,.
		\]
		\item \label{prop:regularity-iii}
		The objective function of problem~\eqref{eq:program:dual:concise} is $\varepsilon$-steep in the sense of Assumption~\ref{a:FW}\,\ref{a:FW:lower}, where $\varepsilon = \min \, \{\varepsilon_x, \varepsilon_w \}$, which depends on the auxiliary constants
		\begin{align*}
			\varepsilon_x &= \left( \frac{\lambda_{\min}(\covsa_w)}{\big( \rho_x + \Tr{\covsa_x}^\frac{1}{2} \big)^2 \lambda_{\max}(H^\top H)+ \big(\rho_w + \Tr{\covsa_w}^{\frac{1}{2}} \big)^2} \right)^2
		\end{align*}
		and
		\begin{align*}
			\varepsilon_w &= \lambda_{\max}(H^\top H) \left( \frac{\lambda_{\min}(\covsa_x)}{\big( \rho_w + \Tr{\covsa_w}^\frac{1}{2} \big)^2 + \lambda_{\min}(\covsa_x) \lambda_{\max}(H^\top H)} \right)^2.
		\end{align*}
	\end{enumerate}
\end{proposition}

\begin{proof}[Proof of Proposition~\ref{prop:regularity}]
	The proof repeatedly uses the fact that, for any $A \in \mbb R^{d_1 \times d_2}$ and $ B \in \mathbb{S}_{+}^{d_2}$, we have
	\begin{equation}\label{eq:eig}
	\lambda_{\max}(A B A^\top)
	= \lambda_{\max}(A^\top A B)
	\leq \lambda_{\max}(A^\top A) \, \lambda_{\max}(B).
	\end{equation}
	The equality in~\eqref{eq:eig} holds because all eigenvalues of $A B A^\top$ are non-negative and because the non-zero spectrum of $A B A^\top$ is identical to that of $A^\top A B$ due to \cite[Proposition~4.4.10]{ref:bernstein2009matrix}. The inequality follows from the observation that $\lambda_{\max}(A^\top A)$ and $\lambda_{\max}(B)$ coincide with the operator norms of the positive semidefinite matrices $A^\top A$ and $B$, respectively.

	As for assertion~\ref{prop:regularity-i}, recall first that the objective function$~f$ of the nonlinear SDP~\eqref{eq:program:dual} is concave. In order to show that $f$ is $\beta$-smooth for some $\beta>0$, it thus suffices to prove that the largest eigenvalue of the positive semidefinite Hessian matrix of $-f$ admits an upper bound uniformly across $\mathcal S^+_x \times \mathcal S^+_w$. By Lemma~\ref{lem:Taylor-f}, the partial gradients of~$f$ evaluated at~$\Sigma_x\succ 0$ and~$\Sigma_w\succ 0$ are given by
    \begin{align*}
	D_x &= \nabla_{\Sigma_x} f\big(\Sigma_x, \Sigma_w\big) = (I_n - \Sigma_x H^\top G^{-1}H )^\top (I_n - \Sigma_x H^\top G^{-1}H)\\
	D_w &= \nabla_{\Sigma_w} f\big(\Sigma_x, \Sigma_w\big) = G^{-1} H \Sigma_x^2 H^\top G^{-1}\,,
    \end{align*}
    where $G= H\Sigma_x H^\top +\Sigma_w$. Moreover, the Hessian matrix
    \[
	\mathcal H =
	\begin{bmatrix}
	\mathcal H_{xx} & \mathcal H_{xw} \\
	\mathcal H_{xw}^\top & \mathcal H_{ww}
	\end{bmatrix} \succeq 0
	\]
	of the convex function~$-f$ evaluated at $\Sigma_x\succ 0$ and~$\Sigma_w\succ 0$ consists of the submatrices
    \begin{align*}
    \begin{array}{l@{\;}l@{\;}l@{\;\;}l}
	\mathcal H_{xx} & = -\nabla^2_{xx} f(\Sigma_x, \Sigma_w) &=& 2 D_x \otimes H^\top G^{-1} H \\[0.5ex]
	\mathcal H_{xw} &= -\nabla^2_{xw} f(\Sigma_x, \Sigma_w) &=& H^\top G^{-1} \otimes (H^\top D_w - \Sigma_x H^\top G^{-1}) +  (H^\top D_w - \Sigma_x H^\top G^{-1}) \otimes H^\top G^{-1} \\[0.5ex]
	\mathcal H_{ww} & = -\nabla^2_{ww} f(\Sigma_x, \Sigma_w) &=& 2 D_w \otimes G^{-1},
	\end{array}
	\end{align*}
	where $\nabla_x$ and $\nabla_w$ are used as shorthands for the nabla operators with respect to $\vect (\Sigma_x)$ and $\vect (\Sigma_w)$, respectively. To construct an upper bound on $\lambda_{\max}(\mathcal H)$ uniformly across $\mathcal S^+_x \times \mathcal S^+_w$, we note first that
	\begin{equation}
	\label{eq:hessian-bound}
	\lambda_{\max}(\mathcal H) \leq \lambda_{\max}(\mathcal H_{xx}) + \lambda_{\max}(\mathcal H_{ww}) = 2 \left( \lambda_{\max} \left( D_x \right) \lambda_{\max} \left( H^\top G^{-1} H \right) + \lambda_{\max} \left( D_w \right) \lambda_{\max} \left( G^{-1} \right) \right),
	\end{equation}
	where the inequality follows from \cite[Fact~5.12.20]{ref:bernstein2009matrix} and the subadditivity of the maximum eigenvalue, whereas the equality exploits the trace rule of the Kronecker product  \cite[Proposition~7.1.10]{ref:bernstein2009matrix}. In the remainder of the proof, we derive an upper bound for each term on the right hand side of the above expression.

	By the definition of $G$ and because $\Sigma_w\in\mathcal S^+_w$, we have $G \succeq \lambda_{\min}(\Sigma_w) I_m \succeq \lambda_{\min}(\covsa_w) I_m$. As $\covsa_w\succ 0$ by assumption, we may thus conclude that $\lambda_{\max} ( G^{-1} ) \leq \lambda_{\min}^{-1}(\covsa_w)$, which in turn implies via~\eqref{eq:eig} that
	\[
	\lambda_{\max} ( H^\top G^{-1} H ) \leq \lambda_{\max}(G^{-1}) \, \lambda_{\max}(H H^\top) \leq \lambda_{\min}^{-1}(\covsa_w) \, \lambda_{\max}(H^\top H).
	\]
	Similarly, by the definition of~$D_w$ we find
	\begin{align*}
		\lambda_{\max} ( D_w )
		&= \lambda_{\max} (G^{-1} H \Sigma_x^2 H^\top G^{-1}) \\
		&\leq \lambda_{\max}^2 ( \Sigma_x ) \, \lambda_{\max}^2 ( G^{-1} ) \, \lambda_{\max} ( H^\top H ) \leq \big( \rho_x + \Tr{\covsa_x}^{\frac{1}{2}} \big)^4 \, \lambda_{\min}^{-2} ( \covsa_w ) \, \lambda_{\max} ( H^\top H ),
	\end{align*}
	where the first inequality follows from applying the estimate~\eqref{eq:eig} twice, while the last inequality reuses the bound on~$\lambda_{\max} ( G^{-1} )$ derived above and exploits Lemma~\ref{lemma:compact:FS}. Finally, by the definition of~$D_w$ we have
	\begin{align*}
		\lambda_{\max} ( D_x ) &\leq 1+\lambda_{\max} ( -H^\top G^{-1} H \Sigma_x )+\lambda_{\max} ( -\Sigma_x H^\top G^{-1} H) + \lambda_{\max} ( H^\top G^{-1} H \Sigma_x^2 H^\top G^{-1} H ) \\&\leq 1 + \lambda_{\max} ( H^\top G^{-1} H \Sigma_x^2 H^\top G^{-1} H ) \\
		&\leq 1 + \lambda_{\max}^2 ( \Sigma_x ) \, \lambda_{\max}^2 ( G^{-1} ) \, \lambda_{\max}^2 ( H^\top H ) \\
		&\leq 1 + \big( \rho_x + \Tr{\covsa_x}^{\frac{1}{2}} \big)^4 \, \lambda_{\min}^{-2}(\covsa_w) \, \lambda_{\max}^2 ( H^\top H ).
	\end{align*}
	where the first inequality holds due to the subadditivity of the maximum eigenvalue and~\cite[Proposition~4.4.10]{ref:bernstein2009matrix}, which implies that the nonzero spectra of $-\Sigma_x H^\top G^{-1} H$ and $-H^\top G^{-1} H \Sigma_x$ are both real and coincide with the nonzero spectrum of the negative semidefinite matrix $-\Sigma_x^\frac{1}{2} H^\top G^{-1} H \Sigma_x^\frac{1}{2}$. The third inequality follows from applying the estimate~\eqref{eq:eig} three times, and the fourth inequality reuses the bound on~$\lambda_{\max} ( G^{-1} )$ and exploits Lemma~\ref{lemma:compact:FS}. Substituting all the above bounds into~\eqref{eq:hessian-bound} completes the proof of assertion~\ref{prop:regularity-i}.

	As for assertion~\ref{prop:regularity-ii}, we first show that the feasible set $\mathcal S^+_x$ is $\alpha_x$-strongly convex with respect to $-f$ in the sense of Assumption~\ref{a:FW}\,\ref{a:FW:set}. To see this, fix any $\Sigma_x, \Sigma_x' \in \mathcal S^+_x$ and $\theta \in [0,1]$, and set
	\begin{align}\label{cov_theta}
		\Sigma_{\theta} = \theta \Sigma_x + (1-\theta) \Sigma_x' + \theta(1-\theta)\frac{\alpha_x}{2} \| \Sigma_x - \Sigma_x' \|^2 \frac{D_x}{\| D_x \|},
	\end{align}
	where $\alpha_x>0$ is defined as in the proposition statement, and $D_x$ denotes again the partial gradient of $f$ with respect to $\Sigma_x$. To prove strong convexity of $\mathcal S^+_x$ with respect to $-f$, we will show that $\Sigma_{\theta}\in \mathcal S^+_x$. Note first that $\Sigma_\theta \succeq \lambda_{\min}(\covsa_x) I_n$ because $\Sigma_x, \Sigma_x' \in \mathcal S^+_x$ and because $D_x$ is positive semidefinite. Next, define
	\[
	\overline {\mathcal S}_x = \left\{ \Sigma_x \in \mathbb{S}_{+}^n: \lambda_{\min}(\covsa_x) I_n \preceq \Sigma_x \preceq (\rho_x + \Trace\big[\covsa_x\big]^{\frac{1}{2}})^2 I_n  \right\},
	\]
	and note that $\mathcal S^+_x \subseteq \overline{\mathcal S}_x$ by Lemma~\ref{lemma:compact:FS}. Moreover, \cite[Theorem~1]{ref:bhatia2018strong} implies that the function $g_x:\mathbb{S}_{+}^n\to\mathbb R$ defined through $g_x(\Sigma_x)=\Tr{\Sigma_x + \covsa_x - 2(\covsa_x^{\frac{1}{2}} \Sigma_x \covsa_x^{\frac{1}{2}})^{\frac{1}{2}}}$ is $\kappa_1$-strongly convex and $\kappa_2$-smooth over~$\overline {\mathcal S}_x$, where
	\begin{equation*}
		\kappa_1 = \frac{\lambda_{\min}^{\frac{1}{2}}(\covsa_x)}{2 (\rho_x + \Trace\big[ \covsa_x \big]^\frac{1}{2})^3} \quad\text{and}\quad
		\kappa_2 = \frac{\rho_x + \Trace\big[\covsa_x\big]^\frac{1}{2}}{2 \lambda_{\min}^{\frac{3}{2}}(\covsa_x)}.
	\end{equation*}
	By \cite[Theorem~12]{ref:journee2010generalized}, the sublevel set $\mathcal S_x^+=\{\Sigma_x\in\overline{\mathcal S}_x:g_x(\Sigma_x)\le \rho_x^2\}$ is thus strongly convex in the canonical sense---relative to $\overline{\mathcal S}_x$---with convexity parameter $\alpha_x=\kappa_1 / (\sqrt{2 \kappa_2} \rho_x)$. This insight implies that $g_x(\overline \Sigma_x)\le \rho_x^2$, which in turn ensures that $\overline\Sigma_x\in\mathcal S_x^+$. As $\theta\in[0,1]$ was chosen arbitrarily, we may conclude that~$\mathcal S_x^+$ is $\alpha_x$-strongly convex with respect to~$f$. Using an analogous argument, one can show that~$\mathcal S^+_w$ is $\alpha_w$-strongly convex with respect to $f$, where $\alpha_w>0$ is defined as in the proposition statement. In summary, $\mathcal S^+_x \times \mathcal S^+_w$ is therefore $\alpha$-strongly convex with respect to $f$ in the sense of Assumption~\ref{a:FW}\,\ref{a:FW:set}, where $\alpha=\min \{ \alpha_x, \alpha_w \}$.

	In order to prove assertion~\ref{prop:regularity-iii}, we will establish lower bounds on $\lambda_{\max}(D_x)$ and $\lambda_{\max}(D_w)$ uniformly across $\mathcal S^+_x \times \mathcal S^+_w$. The claim then follows from the observation that $\lambda_{\max}(D_x)\leq \|D_x\|$ and $\lambda_{\max}(D_w)\leq\|D_w\|$. We first derive a lower bound on~$\lambda_{\max}(D_x)$. To this end, set $T_1 = I_n - \Sigma_x H^\top G^{-1}H$ and $T_2 = H \Sigma_x H^\top G^{-1}$, and note that both $T_1$ and $T_2$ have real spectra thanks to~\cite[Proposition~4.4.10]{ref:bernstein2009matrix}. As $D_x = T_1^\top T_1$, we find
	\begin{align} \label{eq:Dx-bound-1}
		\lambda_{\max}(D_x)
		= \lambda_{\max} (T_1 T_1^\top)
		\geq \max_{\lambda \in \spec(T_1)} | \lambda |^2 = \max_{\lambda \in \spec(T_2)} | 1 - \lambda |^2 = \lambda_{\max}^2(I - T_2) = \lambda_{\max}^2(\Sigma_w G^{-1}),
	\end{align}
	where $\spec(T)$ denotes the eigenvalue spectrum of any square matrix $T$. The inequality in~\eqref{eq:Dx-bound-1} follows from Browne's theorem \cite[Fact~5.11.21]{ref:bernstein2009matrix}, the second equality holds because the non-zero spectrum of $\Sigma_x H^\top G^{-1}H$ matches that of $ H \Sigma_x H^\top G^{-1}$ thanks to~\cite[Proposition~4.4.10]{ref:bernstein2009matrix}, and the last equality follows from the identity
	$T_2= I_m - \Sigma_w G^{-1}$.
	Notice that all eigenvalues of $\Sigma_w G^{-1}$ are real because $T_2$ has a real spectrum.

	The estimate~\eqref{eq:Dx-bound-1} implies that a uniform lower bound on the largest eigenvalue of $D_x$ can be obtained by maximizing the largest eigenvalue of $\Sigma_w G^{-1}$ over $\mathcal S^+_x \times \mathcal S^+_w$. By the definition of $G$, we have
	\begin{equation}
	\label{eq:Dx-bound-2}
	\begin{aligned}
		\inf_{\substack{\Sigma_x \in \mathcal S^+_x \\ \Sigma_w \in \mathcal S^+_w}} \lambda_{\max}(\Sigma_w G^{-1})
		&\ge \inf_{\substack{\Sigma_x \in \mathcal S^+_x \\ \Sigma_w \in \mathcal S^+_w}} \lambda_{\min} \big(\Sigma_w (H \Sigma_x H^\top + \Sigma_w)^{-1} \big) \\
		&\ge \inf_{\substack{\Sigma_x \in \mathcal S^+_x \\ \Sigma_w \in \mathcal S^+_w}} \lambda_{\min}(\Sigma_w) \,\lambda_{\min} \left( (H \Sigma_x H^\top + \Sigma_w)^{-1} \right) \\
		&\ge \inf_{\substack{\Sigma_x \in \mathcal S^+_x \\ \Sigma_w \in \mathcal S^+_w}} \frac{\lambda_{\min}(\Sigma_w)}{\lambda_{\max}(\Sigma_x) \lambda_{\max}(H^\top H)+ \lambda_{\max}(\Sigma_w)} \\
		&\ge \frac{\lambda_{\min}(\covsa_w)}{\big( \rho_x + \Trace\big[\covsa_x\big]^\frac{1}{2} \big)^2 \lambda_{\max}(H^\top H)+ \big(\rho_w + \Trace\big[\covsa_w\big]^{\frac{1}{2}} \big)^2}, \end{aligned}
	\end{equation}
	where the second inequality holds because $\lambda_{\min}(T)=\lambda_{\max}^{-1}(T^{-1})$ for any $T\succ 0$ and because the maximum eigenvalue of a positive definite matrix coincides with its operator norm. The third inequality exploits~\eqref{eq:eig} and the subadditivity of the maximum eigenvalue, and the last inequality follows from Lemma~\ref{lemma:compact:FS}.
	Combining~\eqref{eq:Dx-bound-1} and~\eqref{eq:Dx-bound-2} shows that
	$ \| D_x \| \geq \lambda_{\max}(D_x) \geq \varepsilon_x, $
	where $\varepsilon_x$ is defined as in the proposition statement.

	Using similar arguments, we can also derive a uniform lower bound on~$\lambda_{\max}(D_w)$. Specifically, we have
	\begin{equation}
	\label{eq:lb1}
	\begin{aligned}
		\lambda_{\max}(D_w) \lambda_{\max}(H^\top H)
		&\geq \lambda_{\max}(H^\top D_w H)
		= \lambda_{\max} \big(H^\top G^{-1} H \Sigma_x \, (H^\top G^{-1} H \Sigma_x)^\top \big)\\
		&\geq \lambda_{\max}^2(H \Sigma_x H^\top G^{-1})
		= \left( 1 - \lambda_{\min}(\Sigma_w G^{-1}) \right)^2\\ & = \left( 1 - \frac{1}{1 + \lambda_{\max}(H \Sigma_x H^\top \Sigma_w^{-1})} \right)^2\\ & = \left( 1 - \frac{1}{1 + \lambda_{\max}\left(\Sigma_w^{-\frac{1}{2}} H \Sigma_x H^\top \Sigma_w^{-\frac{1}{2}}\right)} \right)^2,
	\end{aligned}
	\end{equation}
	where the first inequality follows from~\eqref{eq:eig}, and the first equality holds due to the definition of $D_w$. Moreover, the second inequality exploits Brown's theorem~\cite[Fact~5.11.21]{ref:bernstein2009matrix}, and the second equality uses the definition of~$G$. Finally, the third equality follows from the relation $\lambda_{\min}(\Sigma_w G^{-1}) = \lambda_{\max}^{-1}(G \Sigma_w)$, and the fourth equality holds due to~\cite[Proposition~4.4.10]{ref:bernstein2009matrix}.
	A uniform lower bound on $D_w$ can thus be obtained from the estimate
	\[
	\Sigma_w^{-\frac{1}{2}} H \Sigma_x H^\top \Sigma_w^{-\frac{1}{2}} \succeq \lambda_{\min}(\Sigma_x) \Sigma_w^{-\frac{1}{2}} H H^\top \Sigma_w^{-\frac{1}{2}} \succeq \frac{\lambda_{\min}(\Sigma_x)}{\lambda_{\max}(\Sigma_w)} H H^\top,
	\]
	which implies via Lemma~\ref{lemma:compact:FS} that
	\begin{align}
		\label{eq:lb2}
		\inf_{\substack{\Sigma_x \in \mathcal S^+_x \\ \Sigma_w \in \mathcal S^+_w}} \lambda_{\max}(\Sigma_w^{-\frac{1}{2}} H \Sigma_x H^\top \Sigma_w^{-\frac{1}{2}})
		&\geq \inf_{\substack{\Sigma_x \in \mathcal S^+_x \\ \Sigma_w \in \mathcal S^+_w}} \frac{\lambda_{\min}(\Sigma_x)}{\lambda_{\max}(\Sigma_w)} \lambda_{\max}(H^\top H) \geq \frac{\lambda_{\min}(\covsa_x)}{\big(\rho_w + \Trace\big[\covsa_w\big]^{\frac{1}{2}} \big)^2} \lambda_{\max}(H^\top H).
	\end{align}
	Combining \eqref{eq:lb1} and \eqref{eq:lb2} shows that $\| D_w \| \geq \lambda_{\max}(D_w) \geq \varepsilon_w$, where $\varepsilon_w$ is defined as in the proposition statement. We thus conclude that $f$ is $\varepsilon$-steep in the sense of Assumption~\ref{a:FW}\,\ref{a:FW:lower} with $\varepsilon = \min \{\varepsilon_x, \varepsilon_w \}$.
\end{proof}

By Theorem~\ref{theorem:algorithm}, which is applicable because of Proposition~\ref{prop:regularity}, the fully adaptive Frank-Wolfe algorithm (see Algorithm~\ref{algorithm:FAFW}) solves the (minimization problem equivalent to the) nonlinear SDP~\eqref{eq:program:dual:concise} at a linear convergence rate. Moreover,  Theorem~\ref{thm:oracle} ensures that the oracle problem~\eqref{eq:oracle}, which needs to be solved in each iteration of Algorithm~\ref{algorithm:FAFW}, can be solved highly efficiently via bisection (see Algorithm~\ref{algorithm:bisection}).

We emphasize that if $\covsa_x$ is singular, then the strong convexity parameter~$\alpha$ of Assumption~\ref{a:FW}\,\ref{a:FW:set} vanishes, and therefore Theorem~\ref{theorem:algorithm} is no longer applicable. In that case, however, Algorithm~\ref{algorithm:FAFW} is still guaranteed to converge, albeit at a sublinear rate; see \cite{ref:pedregosa2018stepsize} for further details.

\section{Numerical Experiments}
\label{sect:numerical}

All experiments are run on an Intel XEON CPU with 3.40GHz clock speed and 16GB of RAM. All (linear) SDPs are solved with MOSEK~8 using the YALMIP interface~\cite{ref:lofberg2004YALMIP}. In order to ensure the reproducibility of our experiments, we make all source codes available at \url{https://github.com/sorooshafiee/WMMSE}.

\subsection{Scalability of the Frank-Wolfe Algorithm}
\label{sec:scalability}
We first compare the convergence behavior of the Frank-Wolfe algorithm developed in Section~\ref{sect:algorithm} against that of MOSEK. Each experiment consists of $10$ independent simulation runs, in all of which we fix the signal and noise dimensions to $n=m=d$ and the Wasserstein radii to $\rho_x = \rho_w = \sqrt{d}$ for some~$d\in \mathbb N$. In each simulation run we randomly generate two nominal covariance matrices $\covsa_x$ and $\covsa_w$ as follows. First, we sample $Q_x$ and $Q_w$ from the standard normal distribution on $\mbb R^{d\times d}$, and we denote by $R_x$ and $R_w$ the orthogonal matrices whose columns correspond to the orthonormal eigenvectors of $Q_x + Q_x^\top$ and $Q_w+Q_w^\top$, respectively. Then, we set $\covsa_x = R_x \Lambda_x (R_x)^\top$ and $\covsa_w= R_w \Lambda_w R_w^\top$, where $\Lambda_x$ and $\Lambda_w$ are diagonal matrices whose main diagonals are sampled uniformly from $[1,5]^d$ and $[1,2]^d$, respectively. Finally, we set $\wh{\m}_x=0$ and $\wh{\m}_w=0$. The Wasserstein MMSE estimator can then be computed either by solving the nonlinear SDP~\eqref{eq:program:dual} with a Frank-Wolfe algorithm or by solving the linear SDP~\eqref{eq:SDP:dual} with MOSEK. Figures~\ref{fig:algorithm:a} and~\ref{fig:algorithm:b} show the execution time and the number of iterations needed by the vanilla, adaptive and fully adaptive versions of the Frank-Wolfe (FW) algorithm as well as by MOSEK to drive the (surrogate) duality gap below~$10^{-3}$. MOSEK runs out of memory for all dimensions~$d>100$. Figure~\ref{fig:algorithm:c} visualizes the empirical convergence behavior of the three different Frank-Wolfe algorithms. We observe that the fully adaptive Frank-Wolfe algorithm finds highly accurate solutions already after $20$ iterations for problem instances of dimension $d=1{,}000$.

\begin{figure*}
	\centering
	\subfigure[Scaling of execution time]{\label{fig:algorithm:a} \includegraphics[width=0.31\columnwidth]{time.pdf}} \hspace{0pt}
	\subfigure[Scaling of iteration count]{\label{fig:algorithm:b} \includegraphics[width=0.31\columnwidth]{iteration.pdf}} \hspace{0pt}
	\subfigure[Convergence for $d=1000$]{\label{fig:algorithm:c} \includegraphics[width=0.31\columnwidth]{convergence.pdf}}
	\caption{Scalability properties of different methods for computing the Wasserstein MMSE estimator (shown are the means (solid lines) and the ranges (shaded areas) of the respective performance measures across $10$ simulation runs)}
	\label{fig:algorithm}
\end{figure*}

\subsection{The Value of Structural Information}
In the second experiment we study the predictive power of different MMSE estimators.
The experiment consists of $100$ independent simulation runs. In each run, we use the same procedure as in Section~\ref{sec:scalability} to generate two random covariance matrices $\Sigma_x$ and $\Sigma_w$ of dimensions $n=m = 50$, and we set the true signal and noise distributions to $\mathbb{P}_x = \mathcal N(0, \Sigma_x)$ and $\mathbb{P}_w = \mathcal N(0, \Sigma_w)$, respectively. Next, we define $\covsa_x$ and $\covsa_w$ as the sample covariance matrices corresponding to $100$ independent samples from $\mathbb{P}_x$ and $\mathbb{P}_w$, respectively. Moreover, we set $H=I_n$. To assess the value of structural information, we compare the Wasserstein MMSE estimator proposed in this paper against the Bayesian MMSE estimator associated with the nominal signal and noise distributions and the unstructured Wasserstein MMSE estimator proposed in~\cite{ref:shafieezadeh2018wasserstein}. The latter uses a single Wasserstein ball to model the ambiguity of the joint distribution of~$x$ and~$y$, thereby ignoring the structural information that~$w=Hy-x$ is independent of~$x$. Both the structured and unstructured Wasserstein MMSE estimators collapse to the nominal Bayesian MMSE estimator when the underlying Wasserstein radii are set to zero. Recall also that the nominal Bayesian MMSE estimator is optimal in distributionally robust estimation problems whose ambiguity sets are defined via information divergences \cite{ref:levy2004robust, ref:levy2012robust, ref:zorzi2016robust, ref:zorzi2017robustness}. This robustness property makes the nominal Bayesian MMSE estimator an interesting benchmark.

We quantify the performance of a given estimator by its {\em regret}, which is defined as the difference between the estimator's average risk and the least possible average risk of {\em any} estimator under the unknown true distributions $\mathbb{P}_x$ and $\mathbb{P}_w$. Note that the minimum average risk is attained by the (affine) Bayesian MMSE estimator corresponding to the (normal) distributions $\mathbb{P}_x$ and $\mathbb{P}_w$.  Figure~\ref{fig:perf:a} shows the regret of the Wasserstein MMSE estimator with $\rho_x=\rho$ and $\rho_w=0$, the Wasserstein MMSE estimator with $\rho_x=0$ and $\rho_w=\rho$ and the unstructured Wasserstein MMSE estimator from~\cite{ref:shafieezadeh2018wasserstein} with Wasserstein radius $\rho$ for all $\rho\in[0.1,10]$. The solid lines represent the averages and the shaded areas the ranges of the regret across all $100$ simulation runs. The regret of the structured Wasserstein MMSE estimator with Wasserstein radii $\rho_x,\rho_w\in[0.1,10]$ averaged across all $100$ simulation runs is visualized by the surface plot in Figure~\ref{fig:perf:b}.

We observe that the average regret of the nominal Bayesian MMSE estimator amounts to~$16.7$ (the leftmost value of all curves in Figure~\ref{fig:perf:a}), while the best unstructured Wasserstein MMSE estimator attains a significantly lower average regret of~$13.1$ (the minimum of the blue curve in Figure~\ref{fig:perf:a}). The best structured Wasserstein MMSE estimator without noise ambiguity ($\rho_w=0$) displays a similar performance, attaining an average regret of~$13.2$ (the minimum of the red curve in Figure~\ref{fig:perf:a}), while the one without signal ambiguity ($\rho_x=0$) further reduces the average regret by more than $50\%$ to $6.0$ (the minimum of the yellow curve in Figure~\ref{fig:perf:a}). Finally, the best among {\em all} structured Wasserstein MMSE estimators, which is obtained by tuning both $\rho_x$ and $\rho_w$, attains an even lower average regret of~$2.2$ (the minimum of the surface plot in Figure~\ref{fig:perf:b}). This experiment confirms our hypothesis that structural (independence) information as well as information about distributional ambiguity can improve the predictive power of MMSE estimators.

Unlike in data-driven optimization, where the nominal distribution and the radii of the ambiguity sets can be tuned from data, we assume here that the nominal distribution and the radii of the ambiguity sets reflect the modeler's prior distributional information. Thus, they are reminiscent of prior distributions in Bayesian statistics.
The radii could be tuned using standard cross validation techniques, for example, if we had access to samples $(\hat x_i, \hat y_i)$, $i=1,\ldots, N$ drawn independently from the true joint distribution of $(x,y)$ under~$\mathbb{P}$. This conflicts, however, with our central assumption that only $y$ is observable. In this case, it is fundamentally impossible to assess the empirical performance of an estimator and to tune the radii via cross validation.

\begin{figure*}[t]
	\centering
	\subfigure[Regret of different Wasserstein MMSE estimators involving a single hyperparameter $\rho$.]{\label{fig:perf:a} \includegraphics[width=0.45\columnwidth]{fig2-all.pdf}} \hspace{0pt}
	\subfigure[Regret of the Wasserstein MMSE estimator involving two hyperparameters $\rho_x$ and $\rho_w$.]{\label{fig:perf:b} \includegraphics[width=0.45\columnwidth]{fig2-3d.pdf}}
	\caption{Regret of different estimators averaged across 100 simulation runs.}
	\label{fig:perf}
\end{figure*}

\noindent \textbf{Acknowledgments.} We are grateful to Erick Delage for valuable comments on an earlier version of this paper. This research was supported by the Swiss National Science Foundation grant number BSCGI0\_157733.

\subsection{\@startsection{subsection}{2}%
	\z@{.5\linespacing\@plus.7\linespacing}{.5\linespacing}%
	{\normalfont\large\bfseries}}

\def\subsubsection{\@startsection{subsubsection}{3}%
	\z@{.5\linespacing\@plus.7\linespacing}{.5\linespacing}%
	{\normalfont\itshape}}

\RequirePackage[colorlinks,citecolor=blue,urlcolor=blue]{hyperref}
\usepackage{dsfont, amssymb, amsmath, subfigure, graphicx, enumitem}
\usepackage{amsfonts,dsfont,mathtools, mathrsfs, amsthm, accents}
\usepackage[amssymb, thickqspace]{SIunits}
\usepackage{enumitem}

\usepackage{amsthm}
\numberwithin{equation}{section}
\theoremstyle{plain}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{definition}[theorem]{Definition}
\newtheorem{assumption}[theorem]{Assumption}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{remark}[theorem]{Remark}
\newtheorem{example}[theorem]{Example}

\usepackage[numbers]{natbib}
\usepackage{mathrsfs, mathtools, fancyhdr, cases}
\usepackage{float, color, bm, bbm}
\usepackage{subfigure, graphicx, epstopdf, transparent}
\usepackage{multirow, titlecaps, bigstrut}
\usepackage{algorithm, algorithmic}
\renewcommand{\algorithmicrequire}{\textbf{Input:}}
\renewcommand{\algorithmicensure}{\textbf{Output:}}
\usepackage{footnote}
\makesavenoteenv{tabular}
\makesavenoteenv{table}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{lmodern}

redpajama_set_name: RedPajamaArXiv