A Distributionally Robust Estimator that
Dominates the Empirical Average

The aim of this paper is to deploy the basic insights of DRO in order to devise an empirical estimator that outperforms the (commonly used) sample average w.r.t. the mean squared error; the latter being a standard figure of merit for measuring losses in statistics and machine learning (Girshick and Savage, 1951; James and Stein, 1961; Devroye et al., 2003; Wang and Bovik, 2009; Verhaegen and Verdult, 2007; Kuhn et al., 2019).

To set up the stage for our investigation, let $\mathrm{B}(\mathbb{R}^{d})$ be the set of Borel probability measures on $\mathbb{R}^{d}$ and $D=\{X_{1},\dots,X_{n}\}$ a set of independent identically distributed (i.i.d.) random vectors realizing $X\sim\mu\in\mathrm{B}(\mathbb{R}^{d})$. The corresponding to $D$ empirical measure is $\mu_{n}=n^{-1}\sum^{n}_{k=1}\delta_{X_{i}}$, where $\delta_{x}$ is the Delta measure of unit mass located at $x\in\mathbb{R}^{d}$.

By defining the real-valued map $f(\xi,\widehat{\mathrm{X}})\equiv\mathbb{E}_{\xi}\|\widehat{\mathrm{X}}-X\|^{2}_{2},$ $\xi\in\mathrm{B}(\mathbb{R}^{d})$, we obtain the following characterization of the empirical average

as the solution to the surrogate of the “true” problem

when $\mu$ is unknown. Plausibly then, we are interested in the performance of (1) on $\mu$:

where¹¹1We assume that the first order moment is finite $f^{\star}=f(\mu,\mathbb{E}X)$ is the optimal value of (2). Although (1) dominates among unbiased estimators for normally distributed low dimensional data, in the case of skewed or heavy tailed measures, it exhibits a rather sub-optimal performance and therefore, it leaves room for improvement. Such an improvement could be potentially achieved by an appropriate re-weighting $\nu_{n}\in\mathrm{B}(\mathbb{R}^{d})$, $\nu_{n}\ll\mu_{n}$ on $D$, such that

The most prominent approach for doing so is via distributionally robust optimization. As a paradigm, DRO (Scarf et al., 1957) has gained interest due to its connections to regularization, generalization, robustness and statistics (Staib and Jegelka, 2019), (Blanchet et al., 2021), (Bertsimas and Sim, 2004), (Blanchet et al., 2021), (Blanchet and Shapiro, 2023), (Nguyen et al., 2023),(Blanchet et al., 2024).

Let $\mathrm{F}(\beta,\rho)$ be the field of DRO problem instances $\mathrm{P}$ of the form

where $\beta_{\rho}(\mu)\subseteq\mathrm{B}(\mathbb{R}^{d})$ is the uncertainty set of radius $\rho\geq 0$ centered at $\mu$. The two main approaches for representing the uncertainty set, recently unified by (Blanchet et al., 2023), are the divergence approach, and the Wasserstein approach. In the former one, distributional shifts are measured in terms of likelihood ratios (Bertsimas and Sim, 2004; Bayraksan and Love, 2015; Namkoong and Duchi, 2016; Duchi and Namkoong, 2018; Van Parys et al., 2021), while in the latter, the uncertainty set is represented by a Wasserstein ball (Esfahani and Kuhn, 2015; Gao and Kleywegt, 2023; Lee and Raginsky, 2018; Kuhn et al., 2019).

Let $\widehat{\mathrm{X}}^{\star}$ be a solution of $\mathrm{P}\in\mathrm{F}(\beta,\rho)$. Based on the discussion so far, we are interested in an ideal sub-field of instances $\mathrm{F}^{\star}(\beta,\rho)\subset\mathrm{F}(\beta,\rho)$ formed by the following design specifications: First, if $\mathrm{P}\in\mathrm{F}^{\star}(\beta,\rho)$, then there exist a saddle point $(\widehat{\mathrm{X}}^{\star}(\xi^{\star}),\xi^{\star}(\widehat{\mathrm{X}}^{\star}))$. Then, the first design specification implies that necessarily

In some cases, robustness against distributional shifts is as safeguarding against risky events (naturally coordinated by the distribution tail). That is, solving a DRO problem is equivalent to solving a risk-averse problem. This brings us to the second specification which requires $\widehat{\mathrm{X}}^{\star}$ to be the closed form solution of a risk-averse problem (over $\mu$). Combining the aforesaid specifications would allow us to implement the empirical counterpart $\alpha(D;\xi^{\star}_{\rho,n})$ of $\widehat{\mathrm{X}}^{\star}$ from data $D$ drawn from $\mu$. That being said, in this paper we study the following motivating question:

Is there an instance $\mathrm{P}\in\mathrm{F}^{\star}(\beta,\rho)$ such that $\pi_{\alpha}(\mu,\xi^{\star}_{\rho,n}(\beta))<\pi_{\alpha}(\mu,\mu_{n})$?

The aforesaid dual relation between risk-averse and DRO suggests a place to start our design. In particular, our methodology is as follows: We are going to state a risk-averse problem (which in principle is easier to design), that is also solvable in closed form. Then, we are going to show that under certain conditions, its solution also solves a certain, precisely defined DRO problem. Then, we are going to verify our initial assertion by studying the performance of the corresponding empirical solution. We start by considering the following risk-constrained extension of (2):

The risk-averse problem (7) is a quadratically constrained quadratic program and admits the closed form solution (see Section 2)

where $\lambda^{\star}\geq{0}$. Risk-averse optimization refers to the optimization of risk measures. Unlike expectations, these are (often) convex functionals on appropariate spaces of random variables and deliver optimizers that take into account the risky events of the underlying distribution. Unlike expectations, a convex risk measure can be represented through the Legendre transform as the supremum of its affine minorants. This representation essentially characterizes the risk measure as a supremum over a certain class of affine functionals, specifically over a subset of the dual space of signed measures (Shapiro et al., 2021). For the particular class of coherent risk measures, this set contains probability measures and thus, safeguarding for risky events is as being robust against distributional shifts. In other words, solving a risk-constrained problem on the population at hand is as solving a min-max problem where expectations are taken for an appropriate re-weighting of that population.

That said, the challenge is that minimization of a risk measure does not always correspond to a DRO formulation. As an example, problem (7) is closely related to the mean-variance risk-measure (Shapiro et al., 2021, p. 275), which does not lead to a DRO formulation in general. Mean-variance optimization was initially utilized by (Markowitz, 1952), in portfolio selection. Later on, (Abeille et al., 2016) deploys mean-variance for dynamic portfolio allocation problems, and (Kalogerias et al., 2020) stressed the importance of safeguarding against risky events in mean squared error Bayesian estimation. DRO with $\chi^{2}$-divergence is almost equivalent as controlling by the variance (Gotoh et al., 2018; Lam, 2016; Duchi and Namkoong, 2019) were the worst case measure is computed exactly in $O(n\log(n))$ (Staib and Jegelka, 2019). Although minimization of the mean-variance risk-measure does not correspond to an element of $\mathrm{F}(\beta,\rho)$, the mean deviation risk measure of order $p=2$ does (see Theorem 1).

At this point, we state informally the first result of this paper (see Theorem 1), which establishes that the candidate problem (7) corresponds to an instance $\mathrm{P}\in\mathrm{F}^{\star}(\beta,\rho)$:

As discussed later, Theorem 1 states that for the corresponding levels of risk $\varepsilon$ of (7), the closed form solution of (7) solves all $\mathrm{P}\in\mathrm{F}(\chi^{2},\lambda^{\star})$. We answer the main question of the paper in the one dimensional case (see Theorem 2). We find that in the one dimensional case, the empirical counterpart of (8) dominates the empirical average for all probability distributions with negative excess kurtosis, $\kappa$.

All proofs of the results presented hereafter can be found in the supplementary material.

We begin analyzing the risk-constrained mmse problem (7) by first referring the reader to (Kalogerias et al., 2019, Lemma 1), where the well-definiteness of the statement (7) is ensured under the following regularity condition

By noticing that the constraint in (7) can be written as a perfect square, we establish a quadratic reformulation that on the one hand, allows us to interpret the solution to (7) as a stereographic projection, and on the other, sets the stage for a derivative-free solution:

Lemma 1 shows the equivalence of (7) to the convex quadratic program (11) allowing to interpret the former in geometric terms based on the standard inner product in $\mathbb{R}^{d}$. Within this setting, the solution to (7) is a stereographic projection of the $\mathrm{mmse}-$sphere onto $\mathbb{R}^{d}$. Further, we can do slightly more by applying Pythagora’s theorem to the objective in (11) to obtain

where by $\zeta^{-}$ we mean the opposite vector to $\zeta$. According to Lemma 1, the solution of (12) (and therefore of (7)) is the point lying on the ellipse of radius $\sqrt{\bar{\varepsilon}}$ and of center $\frac{1}{4}{\Sigma}^{-1}\zeta$, with the least distance from the orthogonal projection $\mathbb{E}X$. Being guided by the geometric picture, the problem is feasible for any $\bar{\varepsilon}\geq{0}$. In particular, for all $\bar{\varepsilon}\geq\bar{\varepsilon}_{c}=\|{\mathbb{E}X-\frac{1}{4}{\Sigma}^{-1}\zeta^{-}\|^{2}_{4{\Sigma}}}$, $\widehat{\mathrm{X}}^{*}=\mathbb{E}X$, while for $\bar{\varepsilon}=0$, a solution exists only when $\mathbb{E}X=\frac{1}{4}\Sigma^{-1}\zeta^{-}$. Moving forward, for all $\bar{\varepsilon}\in(0,\bar{\varepsilon}_{c})$ the solution is the unique touching point of the circle and the ellipse and it is obtained by equating the “equilibrium forces”. However, based on Lemma 1, we may proceed derivative-free (see Appendix A).

Because by construction, (7) aims to reduce the estimation error variability, it is expected to do so by shifting the optimal estimates towards the areas suffering high loss, that is, towards the tail of the underlying distribution. The next result characterizes the bias of (13): The risk constrained estimator is a shifted version of the conditional expectation by a fraction of the Fisher’s moment coefficient of skewness plus the cross third-order statistics of $X$.

In an effort to minimize the squared error on average while maintaining its variance small, (13) attributes the estimation error made by the conditional expectation to the skewness of the conditional distribution along the eigen-directions of the conditional expected error. In particular, in the one-dimensional case, cross third-order statistics vanish and therefore,

Corollary 1 reveals the “internal mechanism” with which (13) biases towards the tail of the underline distribution. Per error eigen-direction, (13) compensates for the high losses relative to the expectation, by a fraction of the Fisher’s skewness coefficient and a term $\mathrm{T}_{i}$, that captures the cross third-order statistics of the state components. Although (14) involves the cross third-order statistics of the projected state, these statistics cannot emerge artificially from a coordinate change such as $U^{\top}{X}$, but are intrinsic to the state itself. Besides, a meaningfull measure of risk should always be coordinate change invariant.

Motivated by the geometric picture, and the Lagrangean of (7), we note that there exists a unique random element such that (13) is viewed as its corresponding average (see Appendix A). Fortunately, such an element can be expressed as a Radon-Nikodym derivative w.r.t. $\mu$.

Combining Lemma 2 and (15), we obtain ³³3Where we explicitly refer to the measure w.r.t. which expectations are taken.

In words, averaging on the re-distribution $\xi$ is the same as biasing the expectation w.r.t. the initial measure $\mu$ towards the tail of $\mu$ based on a fraction of the Fisher’s moment coefficient of skewness. This fact can be considered as a strong indication of another manifestation of duality between risk-averse and distributionally robust optimization. In addition, for the one dimensional case, assuming that $\mu$ is positively skewed, we obtain the bound

That is, the expectations of neighboring populations that can be tracked with re-weighting are determined by a fraction of the skewness of the population at hand. In the case of empirical measures, (19) implies that lack of mass in the area of the tail is seen as the presence of risky events.

In this section we show that (16) (and therefore the corresponding solution of (7)) is the solution to a $\chi^{2}-$ DRO problem for all $\chi^{2}-$radii less than a computable constant. We start, by taking into account the equivalence between (7) and the mean-deviation risk measure of order $p=2$ (see e.g. (Shapiro et al., 2021)). To not overload the notation, let us declare $\|{X}-\widehat{\mathrm{X}}\|^{2}_{2}=Z_{\widehat{\mathrm{X}}}$. Then, (7) reads

or equivalently

It is a standard exercise to verify that the constraint in (21) is convex w.r.t. $\widehat{\mathrm{X}}$ which leads us to the following result.

Moving forward, strong Lagrangean duality yields

where

is the mean-deviation risk measure of order $p=2$ w.r.t. $\mu$ (see e.g. (Shapiro et al., 2021)), with sensitivity constant, the optimal Lagrange multiplier $\lambda^{\star}_{md}(\varepsilon)$.

According to the previous discussion, we observe that the shift of the expectation in (14), and (18) attributes directionality to the mass re-distribution play of (A) in Theorem 1. In this subsection we make this observation more formal through the smooth structure of the $2-$Wasserstein space $\mathrm{P}_{2,\text{ac}}(\mathbb{R}^{d})$ of probability measures absolutely continuous w.r.t. Lebesgue measure with finite second moment, equipped with the $2-$Wasserstein metric $W_{2}$. In particular, the dynamic formulation of the Wasserstein distance  (Benamou and Brenier, 2000) entails the definition of a (pseudo) Riemannian structure  (Chewi, 2023). At every point $\rho\in\mathrm{P}_{2,\text{ac}}(\mathbb{R}^{d})$, the tangent space $T_{\rho}\mathrm{P}_{2,\text{ac}}(\mathbb{R}^{d}):=\big{\{}\dot{\rho}:\mathbb{R}^{d}\rightarrow\mathbb{R}\big{|}\int\dot{\rho}=0\big{\}}$ contains the re-allocation directions and it is parameterized by functions $\omega(x)$ on $\mathbb{R}^{d}$, through the elliptic operator $\omega\mapsto-\nabla\cdot(\rho\nabla\omega)$. The metric $g:\mathrm{T}\mathrm{P}_{2,\text{ac}}(\mathbb{R}^{d})\times\mathrm{T}\mathrm{P}_{2,\text{ac}}(\mathbb{R}^{d})\rightarrow C^{\infty}(\mathrm{P}_{2,\text{ac}}(\mathbb{R}^{d}))$ ⁵⁵5 By $C^{\infty}(\mathrm{P}_{2,\text{ac}}(\mathbb{R}^{d}))$ we denote the class of smooth real-valued maps defined on $\mathrm{P}_{2,\text{ac}}(\mathbb{R}^{d})$. has value at the point $\rho\in\mathrm{P}_{2,\text{ac}}(\mathbb{R}^{d})$  (Otto, 2001):

Moving forward, let $f\in C^{\infty}(\mathrm{P}_{2,\text{ac}}(\mathbb{R}^{d}))$. Then, the gradient w.r.t. the defined metric, $\mathrm{grad}f:\mathrm{P}_{2,\text{ac}}(\mathbb{R}^{d})\rightarrow\mathrm{T}\mathrm{P}_{2,\text{ac}}(\mathbb{R}^{d})$ is the unique vector field satisfying

where $\mathrm{d}f:\mathrm{P}_{2,\text{ac}}(\mathbb{R}^{d})\rightarrow\mathrm{T}^{*}\mathrm{P}_{2,\text{ac}}(\mathbb{R}^{d})$ is the differential of $f$. To derive the gradient $(\mathrm{grad}f)(\rho)$, consider the curves $\rho_{t}:[0,1]\rightarrow\mathrm{P}_{2,\text{ac}}(\mathbb{R}^{d})$ progressing through $\rho$ and realize the tangent vector $\dot{\rho}_{t}=\nabla\cdot(\rho_{t}\nabla\omega(x))$ for some $\omega(x)$. It is

We identify $(\mathrm{grad}f)(\rho)\equiv\nabla~{}\frac{\delta{f}}{\delta\rho}(x)$, or equivalently $(\mathrm{grad}f)(\rho)=-\nabla\cdot(\rho\nabla~{}\frac{\delta{f}}{\delta\rho}(x))$, where $\frac{\delta{f}}{\delta\rho}$ denotes the first variation of $f$ w.r.t. $\rho$. For the particular case of the $\chi^{2}-$divergence from $\xi$ to $\mu$, we obtain $\mathrm{grad}\chi^{2}(\cdot\|\mu)(\xi)=2\nabla\frac{\mathrm{d}\xi}{\mathrm{d}\mu},$ which when evaluated on (32) yields

By plugging $\widehat{\mathrm{X}}=\widehat{\mathrm{X}}^{\star}_{\lambda_{mv}}$ in (34) and subsequently taking expectations, we obtain (for simplicity in the one dimensional case)

By combining (15) and (35) we may further write

Equation (35) characterizes infinitesimally the optimal strategy of (A): Given the optimal play of (S), (A) changes the measure in (31) aiming to increase the $\chi^{2}-$divergence from $\nu$ to $\mu$. For a more comprehensive review in the theory of optimal transport and applications we refer the reader to (Ambrosio et al., 2005; Villani et al., 2009; Wibisono, 2018; Figalli and Glaudo, 2021; Chewi, 2023).

Now that the risk-constrained estimator has all the desired features, we want to study its performance w.r.t. the mean squared error. We do so in two steps. First, we consider a simplified version of the empirical distributionally robust estimator and we show that it outperforms the empirical average over a rather large class of models. Second (see Appendix B), we argue that with high probability this estimator is the empirical version of the distributionally robust estimator (15). The next results refer to the one dimensional case ($d=1$), with the extension to other dimensions left open for future exploration.

By considering (15) w.r.t. $\mu_{\mathrm{n}}$, we define the following algorithm

where $\bar{\lambda}$ is a non-negative and deterministic parameter that replaces the first factor in the second term of (15). For ease of notation, let us first denote $\bar{\mathrm{X}}_{n}(D)=n^{-1}\sum_{i}X_{i}$, $\mathrm{T}_{n}(D)=n^{-1}\sum_{i}(X_{i}-\bar{\mathrm{X}}_{n})^{3}$. We call (36) the simplified empirical distributionally robust estimator, which after the above notation reads $\widehat{\mathrm{X}}^{\star}_{n,\bar{\lambda}}=\bar{\mathrm{X}}_{n}+\bar{\lambda}\mathrm{T}_{n}$, and achieves mean squared error

for any $\epsilon>0$. Strict domination follows for all those positive $\bar{\lambda}$ such that

only under reasonable conditions that render the right-hand side of (38) strictly positive. Let $Q$ be the set of all probability measures with finite second order moment and $g:Q\rightarrow\mathbb{R}$, the real-valued map with $g(\mu)\equiv 3m^{2}_{2}-m_{4}$, where $m_{2}$ and $m_{4}$ is the second and fourth central moment of $\mu$, respectively. Note that the Gaussian measures are all in $g^{-1}(0)$, and define the set $\mathrm{N}\equiv g^{-1}((0,+\infty))\subset Q$. Strict domination of (36) is demonstrated by the following

As mentioned earlier, Gaussian measures are among those with $3m^{2}_{2}=m_{4}$. In this case $\bar{\lambda}=0$, the condition (37) is satisfied with equality for any $n\geq{1}$, and the mse performance of (36) reduces to that of the empirical average. On the contrary, when the data are generated from $\mu\in\mathrm{N}$, (36) outperforms the average as soon as the parameter $\bar{\lambda}$ is chosen inside the prescribed interval. Interestingly, domination remains feasible for arbitrarily large data sets as soon as $\bar{\lambda}$ is sent to zero with a rate of at least $O(1/n)$. The distributionally robust estimator achieves better mse performance compared to the sample mean by leveraging its bias term. Further, if $\bar{\lambda}\mapsto\varphi(\bar{\lambda})\equiv\mathbb{E}|\widehat{\mathrm{X}}^{\star}_{n,\bar{\lambda}}-X|^{2}$, then the slope of $\varphi$ at zero

It is $\varphi^{\prime}(0)\leq{0}$ for all $\mu\in\mathrm{N}\cup g^{-1}(0)$, and the maximum value is attained for $\mu\in g^{-1}(0)$. In that case of course the best choice is $\bar{\lambda}=0$. On the contrary, for all those probability measures $\mu\in\mathrm{N}$ strict domination follows for (see also Appendix B):

The more negative the slope at zero is, the wider the feasible set becomes. A very negative slope at zero permits large improvement in the mse performance even for very small values of $\bar{\lambda}$ thus rendering, at the same time, the estimator as unbiased as possible. Since the slope at zero is determined by $-(3m^{2}_{2}-m_{4})$, those models with large $|g(\mu)|=|3m^{2}_{2}-m_{4}|$ trade more favorably mse performance for unbiasedness. That said, for models with large value of $g$, values of $\bar{\lambda}$ close to zero achieve a non-trivial performance difference even for large number of samples. However, on the basis of (40), large data sets tend to flatten out $\varphi$ at zero, and to maintain performance $\bar{\lambda}$ has to be increased to

thus rendering $\widehat{\mathrm{X}}^{\star}_{n,\bar{\lambda}^{\star}_{n}}$ more biased. Lastly, $\bar{\lambda}^{\star}_{n}$ approaches zero as the feasible set shrinks with a rate of at least $O(1/n)$ pushing $\widehat{\mathrm{X}}^{\star}_{n,\bar{\lambda}^{\star}_{n}}$ to the empirical average.

We deployed the main insights of DRO to come up with a distributionally robust estimator; its empirical version serving as a better estimator for the mean compared to the empirical average w.r.t. the mean squared error for all platykurtic probability distributions. The aforesaid estimator biases-up to a change of coordinates-towards the tail of the distribution with the Fisher’s moment coefficient of skewness plus additional cross third order statistics and on top of that, it can be written as an expectation w.r.t. an optimal re-weighting of the original measure. This optimal re-weighting along with the optimal estimator form a saddle point of the mse cost and the bias characterizes the directionality of the infinitesimal play through the Wasserstein geometry of the space of measures. Lastly, extending Theorem 2 in the multi-dimensional setting is left for future investigation.