Adversarial Meta-Learning of Gamma-Minimax Estimators
That Leverage Prior Knowledge

We first define the discretization of the model space ${\mathcal{M}}$ that we will use. Let $\{{\mathcal{M}}_{\ell}\}_{\ell=1}^{\infty}$ be an increasing sequence of finite subsets of ${\mathcal{M}}$ such that $\bigcup_{\ell=1}^{\infty}{\mathcal{M}}_{\ell}$ is dense in ${\mathcal{M}}$. That is, $\{{\mathcal{M}}_{\ell}\}_{\ell=1}^{\infty}$ is an increasingly fine grid over ${\mathcal{M}}$. Since ${\mathcal{M}}$ is separable, such an $\{{\mathcal{M}}_{\ell}\}_{\ell=1}^{\infty}$ necessarily exists. Define

for any $d\in{\mathcal{D}}$ and $\Gamma^{\prime}\subseteq\Pi$.

Algorithm 1 describes how the grids ${\mathcal{M}}_{\ell}$ are used to compute an approximately $\Gamma$-minimax estimator in our proposed algorithms. We will show that the approximation error decays to zero as $\ell$ grows to infinity. Here and in the rest of the algorithms in the paper, for any real-valued function $f$, when we assign $\operatorname*{argmin}_{x}f(x)$ or $\operatorname*{argmax}_{x}f(x)$ to a variable, we arbitrarily pick a minimizer or maximizer if there are multiple optimizers. In practice, the user may stop the iteration at some $\ell$ and use a $\Gamma_{\ell}$-minimax estimator $d^{*}_{\ell}$ as the output estimator. We discuss the stopping criterion in more detail at the end of this section.

We note that the minimax problem in Line 3 of Algorithm 1 is nontrivial to solve, and therefore we propose two algorithms that can solve this minimax problem in Sections 3.2 and 3.3.

We assume that the following conditions hold throughout the rest of the paper.

Condition 1 holds if the sequence $\{d_{\ell}^{*}\}_{\ell=1}^{\infty}$ eventually falls in a compact set. For example, if ${\mathcal{D}}$ is a space of neural networks and we take $\varrho$ to be the Euclidean norm in the coefficient space, then we expect Condition 1 to hold if all coefficients are restricted to fall in a bounded set, which is a common restriction in theoretical analyses of neural networks (see, e.g., Goel et al.,, 2016; Zhang et al.,, 2016; Eckle and Schmidt-Hieber,, 2019) and often leads to desirable generalization bounds (see, e.g., Bartlett,, 1997; Bartlett et al.,, 2017; Neyshabur et al.,, 2017). Our theoretical results hold for any limit point $d^{*}$ in Condition 1, even if there is more than one of them.

Condition 2 also often holds. For example, when parameterized using neural networks, all estimators are continuous functions of coefficients for common activation functions such as the sigmoid or the rectified linear unit (ReLU) (Glorot et al.,, 2011) function, and therefore $d\mapsto R(d,P)$ is continuous everywhere.

We next present a sufficient condition to ensure that $d^{*}$ is $\Gamma$-minimax, so that $d^{*}_{\ell}$ is approximately $\Gamma$-minimax for sufficiently large $\ell$.

We note that, in contrast to ${\mathcal{M}}_{\ell}$, $\Omega_{\ell}$ may be an infinite set. We may expect Condition 3 to hold in many cases, especially when $P\mapsto R(d,P)$ is continuous at each $d\in{\mathcal{D}}$ and the grid ${\mathcal{M}}_{\ell}$ contains a variety of distributions that are consistent with prior information represented by $\Gamma$. We illustrate this point with two counterexamples in Appendix A. We will check the plausibility of Condition 3 for Example 2 in our simulation and data analysis in Section 5.1 for exemplar prior information; an almost identical argument shows the plausibility of Condition 3 for Example 3.

We now present the theorem on $\Gamma$-minimaxity of $d^{*}$.

To prove Theorem 1, we utilize a result in Pinelis, (2016) to establish that $r_{\sup}(d,\Gamma)$ can be approximated arbitrarily well by a discrete prior in $\Gamma$ for any $d\in{\mathcal{D}}$. This is a key ingredient in the proof of Lemma 1, which states that, for any $d\in{\mathcal{D}}$, $r_{\sup}(d,\tilde{\Gamma}_{\ell})$ converges to $r_{\sup}(d,\Gamma)$. Then, we show that the sequence $\{r_{\sup}(d^{*}_{\ell},\Gamma_{\ell})\}_{\ell=1}^{\infty}$ is nondecreasing and upper bounded by $\inf_{d\in{\mathcal{D}}}r_{\sup}(d,\Gamma)$, which is less than or equal to the $\Gamma$-maximal Bayes risk $r_{\sup}(d^{*},\Gamma)$ of the limit point $d^{*}$ of $\{d_{\ell}^{*}\}_{\ell=1}^{\infty}$ in Condition 1. Therefore, $r_{\sup}(d^{*}_{\ell},\Gamma_{\ell})$ converges to a limit. We finally use a contradiction argument to prove that this limit is greater than or equal to $r_{\sup}(d^{*},\Gamma)$, which implies Theorem 1.

We have the following corollary on the uniqueness of the $\Gamma$-minimax estimator and the convergence of $\{d^{*}_{\ell}\}_{\ell=1}^{\infty}$ for certain problems.

We prove Corollary 1 by establishing that $d\mapsto r_{\sup}(d,\Gamma)$ is strictly convex.

In practice, the user also needs to specify a stopping criterion for Algorithm 1. In Noubiap and Seidel, (2001), the authors recommended computing or approximating $r_{\sup}(d^{*}_{\ell},\Gamma)$ and stop if $r_{\sup}(d^{*}_{\ell},\Gamma)$ is sufficiently close to $r_{\sup}(d^{*}_{\ell},\Gamma_{\ell})$. However, the procedure to approximate $r_{\sup}(d^{*}_{\ell},\Gamma)$ in that work relies on the compactness of ${\mathcal{M}}$, but we do not want to assume this condition because it may restrict the applicability of the method. Therefore, we propose to use the following alternative criterion: stop if $r_{\sup}(d^{*}_{\ell},\Gamma_{\ell+1})-r_{\sup}(d^{*}_{\ell},\Gamma_{\ell})\leq\epsilon$ for a prespecified tolerance level $\epsilon>0$. This criterion was proposed but not recommended in Noubiap and Seidel, (2001) because it does not guarantee that $r_{\sup}(d^{*}_{\ell},\Gamma_{\ell})$ is close to $r_{\sup}(d^{*},\Gamma)$. For example, if ${\mathcal{M}}_{\ell+1}\setminus{\mathcal{M}}_{\ell}$ is small, it is even possible that $r_{\sup}(d^{*}_{\ell},\Gamma_{\ell+1})-r_{\sup}(d^{*}_{\ell},\Gamma_{\ell})=0$, but $d^{*}_{\ell}$ is far from being $\Gamma$-minimax. In contrast, we recommend this criterion for our proposed methods because we allow more flexibility in model specification, that is, ${\mathcal{M}}$ need not be compact. We discuss this issue in more detail in Section 4.1.

We finally remark that $r_{\sup}(d,\Gamma_{\ell})$ may be difficult to evaluate exactly. Since the risk is often an expectation, we recommend approximating $r_{\sup}(d,\Gamma_{\ell})$ for any given $d$ via Monte Carlo as follows: first, estimate risks $R(d,P)$ for all $P\in{\mathcal{M}}_{\ell}$ with a large number of Monte Carlo runs; second, estimate the corresponding least favorable prior $\pi_{d,\ell}\in\operatorname*{argmax}_{\pi\in\Gamma_{\ell}}r(d,\pi)$ using the estimated risks; third, estimate the risks $R(d,P)$ ($P\in{\mathcal{M}}_{\ell}$) again with independent Monte Carlo runs, and, finally, calculate $r(d,\pi_{d,\ell})$ with the estimated risks and the estimated least favorable prior. Using two independent estimates of the risk can remove the positive bias that would otherwise arise due to using the same data to estimate the risks and the least favorable prior.

In this section, we present methods to compute a $\Gamma_{\ell}$-minimax estimator, which corresponds to Line 3 in Algorithm 1. Gradient descent with max-oracle (GDmax) and its stochastic variant (SGDmax), which were presented in Lin et al., (2020), can be used to solve general minimax problems in Euclidean spaces. We focus on SGDmax in the main text and present GDmax in Appendix B. To apply these algorithms to find a $\Gamma_{\ell}$-m inimax estimator, we need to assume that ${\mathcal{D}}$ can be parameterized by a subset of a Euclidean space, that is, that for any $d\in{\mathcal{D}}$, there exists a real vector-valued coefficient $\beta\in^{D}$ such that $d$ may be written as $d(\beta)$. For example, ${\mathcal{D}}$ may be a neural network class. More discussions on the parameterization of ${\mathcal{D}}$ can be found in Section 4.2. In this section, in a slight abuse of notation, we define $R(\beta,P):=R(d(\beta),P)$, $r(\beta,\pi):=r(d(\beta),\pi)$ and $r_{\sup}(\beta,\Gamma_{\ell}):=r_{\sup}(d(\beta),\Gamma_{\ell})$ for a coefficient $\beta\in^{D}$, a data-generating mechanism $P\in{\mathcal{M}}$ and a prior $\pi\in\Gamma$. We assume that $\beta\mapsto R(\beta,P)$ is differentiable for all $P\in{\mathcal{M}}$, and hence so is $\beta\mapsto r(\beta,\pi)$ for all $\pi\in\Gamma$.

It is often the case that $R(\beta,P)$ is expressed as an expectation. In this case, $R(\beta,P)$ may instead be approximated using Monte Carlo techniques. With $\xi$ being an exogenous source of randomness according to law $\Xi$, let $\hat{R}(\beta,P,\xi)$ be an unbiased approximation of $R(\beta,P)$ with ${\mathbb{E}}[\|\nabla_{\beta}\{\hat{R}(\beta,P,\xi)-R(\beta,P)\}\|^{2}]\leq\sigma^{2}<\infty$, where $\|\cdot\|$ denotes the $\ell_{2}$-norm in Euclidean spaces. Let $\hat{r}(\beta,\pi,\xi):=\int\hat{R}(\beta,P,\xi)\,\pi({\mathrm{d}}P)$ for $\pi\in\Gamma_{\ell}$. In this case, SGDmax (Algorithm 2) may be used to find a (locally) $\Gamma_{\ell}$-minimax estimator. Note that Algorithm 2 represents a generalization of the nested minimax AMC strategy in Luedtke et al., 2020a to $\Gamma_{\ell}$-minimax problems.

We next present two conditions needed for the validity of Algorithm 2.

Note that Condition 4 differs from Condition 2 in that the former relies on the parameterization of ${\mathcal{D}}$ in a Euclidean space equipped with the Euclidean norm, while the latter may rely on a different metric on ${\mathcal{D}}$ such as an $L^{2}$-distance.

Under these conditions, using the results in Lin et al., (2020), we can show that SGDmax yields an approximation to a local minimum of $\beta\mapsto r_{\sup}(\beta,\Gamma_{\ell})$ when the algorithms’ hyperparameters are suitably chosen. Before we formally present the theorem, we introduce some definitions related to the local optimality of a potentially nondifferentiable and nonconvex function. A real-valued function $f$ is called $q$-weakly convex if $x\mapsto f(x)+(q/2)\|x\|^{2}$ is convex ($q>0$). The Moreau envelope of a real-valued function $f$ with parameter $q>0$ is $f_{q}:x\mapsto\min_{x^{\prime}}f(x^{\prime})+\|x^{\prime}-x\|^{2}/(2q)$. A point $x$ is an $\epsilon$-stationary point ($\epsilon\geq 0$) of a $q$-weakly convex function $f$ if $\|\nabla f_{1/(2q)}(x)\|\leq\epsilon$. Similarly, a random point $x$ is an $\epsilon$-stationary point ($\epsilon\geq 0$) of a $q$-weakly convex function $f$ in expectation if ${\mathbb{E}}[\|\nabla f_{1/(2q)}(x)\|]\leq\epsilon$. If $x$ is an $\epsilon$-stationary point in expectation, we may conclude that it is an $\epsilon$-stationary point with high probability by Markov’s inequality. Lemma 3.8 in Lin et al., (2020) shows that an $\epsilon$-stationary point of $f$ is close to a point $x^{\prime}$ at which $f$ has at least one small subgradient for small $\epsilon$, so that $f(x^{\prime})$ is close to a local minimum. In other words, if an algorithm outputs an estimator $\hat{d}=d(\hat{\beta})$ such that $\hat{\beta}$ is an $\epsilon$-stationary point of $\beta\mapsto r_{\sup}(\beta,\Gamma_{\ell})$, then we know that $r_{\sup}(\hat{\beta},\Gamma_{\ell})$ is close to a local minimum of $\beta\mapsto r_{\sup}(\beta,\Gamma_{\ell})$.

We next present the validity result for Algorithm 2.

The assumption that the batch size $J=1$ is purely for convenience since increasing $J$ corresponds to decreasing variance $\sigma^{2}$. To run Algorithm 2 in practice, the user only needs to specify tuning parameters in Line 1 and all other constants in Theorem 2 need not be known. In general, a small learning rate $\eta$, a stringent accuracy $\zeta$, and a large batch size $J$ make Algorithm 2 likely to eventually reach an approximation of a local minimum of $\beta\mapsto r_{\sup}(\beta,\Gamma_{\ell})$, but computation time might increase. Similar to most numeric optimization algorithms, fine-tuning is needed to achieve a balance between convergence guarantee and computation time, but a conservative choice of tuning parameters would typically result in convergence at the cost of computation time.

We note that Line 3 in Algorithm 2 may be inconvenient to implement because linear program solvers often do not use stochastic optimization. Therefore, we propose to use a convenient variant (Algorithm 6 in Appendix B), where the stochastic maximization step (Line 3 in Algorithm 2) is replaced by solving a linear program where the objective is approximated via Monte Carlo. This variant has similar validity under similar conditions. We also note that the two uniform Lipschitz continuity conditions (4 and 5) heavily rely on the fact that ${\mathcal{M}}_{\ell}$ is finite and the compactness of a set containing the coefficients. Nevertheless, the latter compactness restriction is common in theoretical analyses of neural networks (see, e.g., Goel et al.,, 2016; Zhang et al.,, 2016; Eckle and Schmidt-Hieber,, 2019). Moreover, these two conditions are sufficient conditions for the validity of the gradient-based methods, namely SGDmax, our variant of SGDmax and GDmax; a guarantee similar to these validity results might hold when two conditions are violated.

We finally remark that other algorithms similar to SGDmax can be applied, for example, (stochastic) gradient descent ascent with projection (Lin et al.,, 2020), (stochastic) mirror descent ascent, or accelerated (stochastic) mirror descent ascent (Huang et al.,, 2021). It is of future research interest to develop gradient-based methods to solve minimax problems with convergence guarantees under weaker conditions.

The algorithms in Section 3.2 may be convenient in many cases, but the requirements such as parameterization of the space ${\mathcal{D}}$ of estimators in a Euclidean space, differentiability of the risk function $R$ with respect to the coefficients $\beta$, and uniform Lipschitz continuity may be restrictive for certain problems. In this section, we propose an alternative algorithm, fictitious play, that avoids these requirements. We also present its convergence results.

Brown, (1951) introduced fictitious play as a means to find the value of a zero-sum game, that is, the optimal mixed strategy for both players and their expected gains. Robinson, (1951) then proved that fictitious play can be used to iteratively solve a two-player zero-sum game for a saddle point that is a pair of mixed strategies where both players have finitely many pure strategies. Our problem of finding a $\Gamma$-minimax estimator may also be viewed as a two-player zero-sum game where one player chooses a prior from $\Gamma$ and the other player chooses an estimator from ${\mathcal{D}}$. If we assume that, for the $\Gamma$-minimax problem at hand, the pair of both players’ optimal strategies is a saddle point, which holds in many minimax problems (e.g., v. Neumann,, 1928; Fan,, 1953; Sion,, 1958), then fictitious play may also be used to find a $\Gamma$-minimax estimator. Since $\Gamma$ may be too rich to allow for feasible implementation of fictitious play, we propose to use this algorithm to find a $\Gamma_{\ell}$-minimax estimator.

In the fictitious play algorithm in Robinson, (1951), the two players take turns to play the best pure strategy against the mixture of the opponent’s historic pure strategies, and the final output is a pair of mixtures of the two players’ historic pure strategies. Since this algorithm aims to find minimax mixed strategies, we consider stochastic estimators. That is, consider the Borel $\sigma$-field $\mathcal{F}$ over ${\mathcal{D}}$ and let $\varPi$ denote the set of all probability distributions on the measurable space $({\mathcal{D}},\mathcal{F})$. We define $\overline{{\mathcal{D}}}$ to be the space of stochastic estimators with each element taking the following form: first draw an estimator from ${\mathcal{D}}$ according to a distribution $\varpi\in\varPi$ with an exogenous random mechanism and then use the estimator to obtain an estimate based on the data. Note that we may write any $\overline{d}\in\overline{{\mathcal{D}}}$ as $\overline{d}(\varpi)$ for some $\varpi\in\varPi$. We consider estimators in $\overline{{\mathcal{D}}}$ throughout this section, with the definition of $\Gamma$-minimaxity extended in the natural way, so that $\overline{d}^{*}=\overline{d}(\varpi^{*})\in\overline{{\mathcal{D}}}$ is $\Gamma$-minimax if $r_{\sup}(\overline{d}^{*},\Gamma)=\min_{\overline{d}\in\overline{{\mathcal{D}}}}r_{\sup}(\overline{d},\Gamma)$; we similarly extend all other definitions from Section 2. We assume that there exists $\pi^{*}_{\ell}\in\Gamma_{\ell}$ ($\ell=1,2,\ldots$) such that

In other words, $(\overline{d}^{*},\pi^{*}_{\ell})$ is a saddle point of $r$ in $\overline{{\mathcal{D}}}\times\Gamma_{\ell}$. Under this condition and the further conditions that ${\mathcal{D}}$ is convex and $d\mapsto R(d,P)$ is convex for all $P\in{\mathcal{M}}$, it is possible to use a $\Gamma$-minimax estimator over the richer class $\overline{{\mathcal{D}}}$ of stochastic estimators to derive a $\Gamma$-minimax estimator over the original class ${\mathcal{D}}$. Indeed, for any $\overline{d}(\varpi)\in\overline{{\mathcal{D}}}$ and $P\in{\mathcal{M}}$, by Jensen’s inequality, $R(\overline{d}(\varpi),P)=\int R(d,P)\,\varpi({\mathrm{d}}d)\geq R(\underline{\overline{d}}(\varpi),P)$ where $\underline{\overline{d}}(\varpi):=\int d\,\varpi({\mathrm{d}}d)\in{\mathcal{D}}$ is the average of the stochastic estimator $\overline{d}(\varpi)$; that is, the risk of $\underline{\overline{d}}(\varpi)$ is never greater than that of $\overline{d}(\varpi)$. Therefore, we may use the fictitious play algorithm to compute $\overline{d}(\varpi^{*}_{\ell})$ for each $\ell$ and further apply Algorithm 1 to compute $\overline{d}(\varpi^{*})$. After that, we may take $\underline{\overline{d}}(\varpi^{*})$ as the final output deterministic estimator.

Algorithm 3 presents the fictitious play algorithm for finding a $\Gamma_{\ell}$-minimax estimator in $\overline{{\mathcal{D}}}$. Note that $\Gamma_{\ell}$ is convex, and hence $\pi$ always lies in $\Gamma_{\ell}$ throughout the iterations. In practice, we may initialize $\varpi$ as a point mass at an initial estimator in ${\mathcal{D}}$. In addition, similarly to Robinson, (1951), we may replace Line 5 with $d^{\dagger}_{(t)}\leftarrow\operatorname*{argmin}_{d\in{\mathcal{D}}}r(d,\pi_{(t)})$, that is, minimizing the Bayes risk with the most recently updated prior rather than with the previous prior.

We next present a convergence result for this algorithm.

Robinson, (1951) proved a similar case for two-player zero-sum games where each player has finitely many pure strategies. In contrast, in our problem, each player may have infinitely many pure strategies. A natural attempt to prove Theorem 3 would be to consider finite covers of $\bar{{\mathcal{D}}}$ and $\Gamma_{\ell}$, that is, $\bar{{\mathcal{D}}}=\bigcup_{i=1}^{I}{\mathcal{D}}_{i}$ and $\Gamma_{\ell}=\bigcup_{j=1}^{J}\Pi_{j}$, such that the range of $r(d,\pi)$ in each ${\mathcal{D}}_{i}$ and $\Pi_{j}$ is small (say less than $\epsilon$), bin pure strategies into these subsets, and then apply the argument in Robinson, (1951) to these bins. The collection of ${\mathcal{D}}_{i}$ and $\Pi_{j}$ may be viewed as finitely many approximated pure strategies to $\Gamma_{\ell}$ and $\bar{{\mathcal{D}}}$ up to accuracy $\epsilon$, respectively. Unfortunately, we found that this approach fails. The problem arises because Robinson, (1951) inducted on $I$ and $J$, and, after each induction step, the corresponding upper bound becomes twice as large. Unlike the case with finitely many pure strategies that was considered in Brown, (1951) and Robinson, (1951), as the desired approximation accuracy $\epsilon$ approaches zero, the numbers of approximated pure strategies, $I$ and $J$, may diverge to infinity, and so does the number of induction steps. Therefore, the resulting final upper bound is of order $2^{I+J}\epsilon$ and generally does not converge to zero as $\epsilon$ tends to zero. To overcome this challenge, we instead control the increase in the relevant upper bound after each induction step more carefully so that the final upper bound converges to zero as $\epsilon$ decreases to zero, despite the fact that $I$ and $J$ may diverge to infinity.

We remark that, because Line 5 of Algorithm 3 typically involves another layer of iteration in addition to that over $t$, this algorithm will often be more computationally intensive than SGDmax. Nevertheless, Algorithm 3 provides an approach to construct $\Gamma_{\ell}$-minimax estimators in cases where these other algorithms cannot be applied, for example, in settings where the risk is not differentiable in the parameters indexing the estimator or uniform Lipschitz conditions fail. In our numerical experiments, we have implemented this algorithm in the context of mean estimation (Appendix C).

By Theorem 1, $r_{\sup}(d_{\ell}^{*},\Gamma_{\ell})\nearrow\min_{d\in{\mathcal{D}}}r_{\sup}(d,\Gamma)$ whenever Conditions 1–3 hold and the increasing sequence $\{{\mathcal{M}}_{\ell}\}_{\ell=1}^{\infty}$ is such that $\bigcup_{\ell=1}^{\infty}{\mathcal{M}}_{\ell}$ is dense in ${\mathcal{M}}$. Though this guarantee holds for all such sequences $\{{\mathcal{M}}_{\ell}\}_{\ell=1}^{\infty}$, in practice, judiciously choosing this sequence of grids of distributions can lead to faster convergence. In particular, it is desirable that the least favorable prior $\Gamma_{\ell}$ puts mass on some of the distributions in ${\mathcal{M}}_{\ell}\backslash{\mathcal{M}}_{\ell-1}$ since, if this is not the case, then $d_{\ell}^{*}$ will be the same as $d_{\ell-1}^{*}$. While we may try to arrange for this to occur by adding many new points when enlarging ${\mathcal{M}}_{\ell-1}$ to ${\mathcal{M}}_{\ell}$, it may not be likely that any of these points will actually modify the least favorable prior unless they are carefully chosen.

To better address this issue, we propose to add grid points using a Markov chain Monte Carlo (MCMC) method. Our intuition is that, given an estimator $d$, the maximal Bayes risk is likely to significantly increase if we add distributions that (i) have a high risk for $d$, and (ii) are consistent with prior information so that there exists some prior such that these distributions lie in a high-probability region. We propose to use the MCMC algorithm to bias the selection of distributions in favor of those with the above characteristics. Let $\tau:{\mathcal{M}}\rightarrow[0,\infty)$ denote a function such that $\tau(P)>\tau(P^{\prime})$ if $P$ is more consistent with prior information than $P^{\prime}$. For example, given a prior mean $\mu$ of some real-valued summary $\Psi(P)$ of $P$ and an interval $I$ that contains $\Psi(P)$ with prior probability at least 95%, we may choose $\tau:P\mapsto\phi(\Psi(P))$, where $\phi$ is the density of a normal distribution that has mean $\mu$ and places 95% of its probability mass in $I$. We call $\tau$ a pseudo-prior. Then, with the current estimator being $d$, we wish to select distributions $P$ for which $R(d,P)\tau(P)$ is large. We may use the Metropolis-Hastings-Green algorithm (Metropolis et al.,, 1953; Hastings,, 1970; Green,, 1995) to draw samples from a density proportional to $P\mapsto R(d,P)\tau(P)$. We then let ${\mathcal{M}}_{\ell}$ be equal to the union of ${\mathcal{M}}_{\ell-1}$ and the set containing all unique distributions in this sample.

Details of the proposed scheme are provided in Algorithm 4. To use this proposed algorithm, we rely on it being possible to define a sequence of parametric models $\{\tilde{\Omega}_{\ell}\}_{\ell=1}^{\infty}$ such that $\tilde{{\mathcal{M}}}:=\cup_{\ell=1}^{\infty}\tilde{\Omega}_{\ell}$ is dense in ${\mathcal{M}}_{\ell}$—this is possible in many interesting examples (see, e.g., Chen,, 2007). When combined with separability of ${\mathcal{M}}$, this condition enables the definition of an increasing sequence of grids of distributions $\{{\mathcal{M}}_{\ell}\}_{\ell=1}^{\infty}$ such that, for each $\ell$, ${\mathcal{M}}_{\ell}\subseteq\tilde{{\mathcal{M}}}$.

The following theorem on distributional convergence follows from that for the Metropolis-Hastings-Green algorithm (see Section 3.2 and 3.3 of Green,, 1995).

Therefore, if $\mathscr{L}$ corresponds to a continuous distribution with nonzero density over the parameter space of $\tilde{{\mathcal{M}}}$, then Theorem 4 implies that $\bigcup_{\ell=1}^{\infty}{\mathcal{M}}_{\ell}$ is dense in ${\mathcal{M}}$, as required by Algorithm 1.

Implementing Algorithm 4 relies on the user making several decisions. These decisions include the choice of the pseudo-prior $\tau$ and the technique used to approximate the risk $R(d,P)$ to a reasonable accuracy. Fortunately, regardless of the decisions made, Theorem 1 suggests that $r_{\sup}(d_{\ell}^{*},\Gamma_{\ell})\nearrow\min_{d\in{\mathcal{D}}}r_{\sup}(d,\Gamma)$ for a wide range of sequences $\{{\mathcal{M}}_{\ell}\}_{\ell=1}^{\infty}$. Indeed, all that theorem requires on this sequence is that the grid ${\mathcal{M}}_{\ell}$ becomes arbitrarily fine as $\ell$ increases. Though the final decisions made are not important when $\ell$ is large, we still comment briefly on the decisions that we have made in our experiments, First, we have found it effective to approximate $R(d,P)$ via a large number of Monte Carlo draws. Second, in a variety of settings, we have also identified, via numerical experiments, candidate pseudo-priors that balance high risk and consistency with prior information (see Sections 5.1 and 5.2 for details).

It is desirable to consider a rich space ${\mathcal{D}}_{0}$ of estimators to obtain an estimator with low maximal Bayes risk, and thus good general performance. However, to make numerically constructing these estimators computationally feasible, we usually have to consider a restricted space ${\mathcal{D}}$ of estimators. This approximation is justified because, if estimators in ${\mathcal{D}}$ can approximate the $Gamma$-minimax estimator in ${\mathcal{D}}_{0}$ well, then we expect the resulting excess maximal Bayes risk is small.

Feedforward neural networks (or neural networks for short) are natural options for the space of estimators because of their universal approximation property (e.g., Hornik,, 1991; Csáji,, 2001; Hanin and Sellke,, 2017; Kidger and Lyons,, 2020). However, training commonly used neural networks can be computationally intensive. Moreover, a space of neural networks is typically nonconvex, and hence it may be difficult to find a global minimizer of the maximal Bayes risk even if the risk is convex in the estimator. Therefore, the learned estimator might not perform well.

To help overcome this challenge, we advocate for utilizing available statistical knowledge when designing the space of estimators. We call estimators that take this form statistical knowledge networks. In particular, if a simple estimator is already available, we propose to use neural networks with such an estimator as a node connected to the output node. An example of such an architecture is presented in Fig. 2. In this sample architecture, each node is an activation function such as the sigmoid or the rectified linear unit (ReLU) (Glorot et al.,, 2011) function applied to an affine transformation of the vector containing the ancestors of the node. The only exception is the output node, which is again an affine transformation of its ancestors but uses the identity activation function. When training the neural network, we may initialize the affine transformation in the output layer to only give weight to the simple estimator. Under this approach, the space of estimators is a set of perturbations of an existing simple estimator. Although we may still face the challenge of nonconvexity and local optimality, we can at least expect to improve the initial simple estimator.

In the simulation we describe in Appendix C, we compared the empirical performance of several spaces of estimators. This simulation concerns the simple problem of estimating the mean of a true distribution whose support has known bounds (Example 1), and the existing simple estimator we use in the statistical neural network is the sample mean. Fig. 3 presents the trajectory of estimated Bayes risks. As shown in subfigures (b)–(d), using the statistical knowledge network, the estimator is almost $\Gamma$-minimax after a few iterations; on the other hand, it took about 1000 iterations for the feedforward neural network to reach an approximately $\Gamma$-minimax estimator. Therefore, in this simple problem where the true $\Gamma$-minimax estimator is a shifted and scaled sample mean, statistical knowledge substantially reduced the number of iterations required to obtain an approximately $\Gamma$-minimax estimator. For more complicated problems, we expect statistical knowledge to further help improve the performance of the computed estimator.

We note that we might overcome the challenge of nonconvexity and local optimality by using an extreme learning machine (ELM) (Huang et al., 2006b, ) to parameterize the estimator. ELMs are neural networks for which the weights in hidden nodes are randomly generated and are held fixed, and only the weights in the output layer are trained. Thus, the space of ELMs with a fixed architecture and fixed hidden layer weights is convex. Like traditional neural networks, ELMs have the universal approximation property (Huang et al., 2006a, ). In addition, Corollary 1 may be applied to an ELM so that the $\Gamma_{\ell}$-minimax estimator may converge to the $\Gamma$-minimax estimator. As for traditional neural networks, we may incorporate knowledge of existing statistical estimators into an ELM.

We finally remark that, besides computational intensity when constructing (i.e., learning) a $\Gamma$-minimax estimator, another important factor to be considered when choosing ${\mathcal{D}}$ is the computational intensity to evaluate the learned estimator at the observed dataset. This is another reason for our choosing neural networks or ELMs as the space of estimators. Indeed, existing software packages (e.g., Paszke et al.,, 2019) make it easy to leverage graphics processing units to efficiently evaluate the output of neural networks for any given input. Therefore, if the existing estimator being used is not too difficult to compute, then estimators parameterized using similar architectures to that displayed in Figure 2 will be able to be computed efficiently in practice. This efficiency may be especially important in settings where the estimator will be applied to many datasets, so that the cost of learning the estimator is amortized and the main computational expense is evaluating the learned estimator.

We apply our proposed method to Example 2. In the simulation, we set the true population to be an infinite population with the same categories and same proportions as the sample studied in Miller and Wiegert, (1989), which consists of 1088 observations in 188 categories. This setting is the same as the simulation setting in Shen et al., (2003). We set the sample size to be $n=100$ and the size of the new sample to be $m=200$. In this setting, the expected number of new categories in the new sample unconditionally on the observed sample, namely $\Phi(P_{0}):={\mathbb{E}}_{P_{0}}[\Psi(P_{0})(\mathbf{X}^{*})]$, can be analytically computed and equals $48.02$. We note that this quantity can also be computed via simulation: (i) sample $n$ and $m$ individuals with replacement from the dataset in Miller and Wiegert, (1989), (ii) count the number of new categories in the second sample, and (iii) repeat steps (i) and (ii) many times and compute the average.

It is well known that this prediction problem is difficult when $m>n$, and we run this simulation to investigate the potential gain from leveraging prior information by computing a Gamma-minimax estimator for such difficult or even ill-posed problems. We consider three sets of prior information:

1. 1.

   strongly informative: prior mean of $\Phi(P)$ in $[45,50]$, $\geq 95\%$ prior probability that $\Phi(P)$ lies in $[40,55]$;

2. 2.

   weakly informative: prior mean of $\Phi(P)$ in $[40,55]$, $\geq 95\%$ prior probability that $\Phi(P)$ lies in $[30,65]$; and

3. 3.

   almost noninformative: prior mean of $\Phi(P)$ in $[35,60]$, $\geq 95\%$ prior probability that $\Phi(P)$ lies in $[20,75]$.

We note that a traditional Bayesian approach would require specifying a prior on ${\mathcal{M}}$, including the total number of categories and the proportion of each category, which may be difficult in practice.

We check the plausibility of Condition 3 in this context. We take the strongly informative prior information as an example. Take $\Omega_{\ell}$ to be the collection of multinomial distributions with at most $\ell$ categories. It is obvious that $\bigcup_{\ell=1}^{\infty}\Omega_{\ell}={\mathcal{M}}$. Let $d\in{\mathcal{D}}$ be fixed and $\pi\in\tilde{\Gamma}_{\ell}$ be a fixed prior with finite support, that is, $\pi=\sum_{j=1}^{J}q_{j}\delta(Q_{j})$ where $\delta(\cdot)$ denotes the point mass distribution, $Q_{j}\in\Omega_{\ell}$, $q_{j}>0$ and $\sum_{j=1}^{J}q_{j}=1$. Let $\epsilon>0$ be an arbitrary small number such that $\sum_{j=1}^{J}q_{j}\Phi(Q_{j})\leq 50-\epsilon$ or $\sum_{j=1}^{J}q_{j}\Phi(Q_{j})\geq 45+\epsilon$. Since $\bigcup_{\ell=1}^{\infty}{\mathcal{M}}_{\ell}$ is dense in ${\mathcal{M}}$ and $\Phi$ is continuous, there exists a sufficiently large $i$ such that, for every distribution $Q_{j}$, there exists $P_{j}\in{\mathcal{M}}_{i}\cap\Omega_{\ell}$ satisfying the following:

• $\square$

  $|\Phi(P_{j})-\Phi(Q_{j})|\leq\epsilon$;

• $\square$

  if $\Phi(Q_{j})\in[40,55]$, then $\Phi(P_{j})\in[40,55]$;

• $\square$

  $|R(d,P_{j})-R(d,Q_{j})|\leq\epsilon$.

Take $\pi_{i}$ to be $\sum_{j=1}^{J}q_{j}\delta(P_{j})$. Then it is easy to verify that $|\sum_{j=1}^{J}q_{j}\Phi(P_{j})-\sum_{j=1}^{J}q_{j}\Phi(Q_{j})|\leq\epsilon$ and thus $\sum_{j=1}^{J}q_{j}\Phi(P_{j})\in[45,50]$; moreover, $\Phi(Q_{j})\in[40,55]$ implies that $\Phi(P_{j})\in[40,55]$ and therefore $\sum_{j=1}^{J}q_{j}{\mathbbm{1}}(\Phi(P_{j})\in[40,55])\geq\sum_{j=1}^{J}q_{j}{\mathbbm{1}}(\Phi(Q_{j})\in[40,55])\geq 95\%$. Thus, $\pi_{i}\in\tilde{\Gamma}_{i\mid\ell}$. Moreover, $|r(d,\pi)-r(d,\pi_{i})|\leq\epsilon$. Therefore, $r(d,\pi_{i})\rightarrow r(d,\pi)$ as $i\rightarrow\infty$ and Condition 3 holds.