Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
93,605 characters · 19 sections · 76 citation commands
Approximate Least-Favorable Distributions and Nearly Optimal Tests via Stochastic Mirror Descent
\onehalfspacing
\doublespacing
Consider the problem of testing a null hypothesis that postulates finitely many distributions for the observed data against a single alternative hypothesis. More concretely, suppose the data is modeled as a $\mathcal{Y}$-valued random variable, denoted as $Y$, and let the hypothesis testing problem take the form:
The structure of the null and alternative hypothesis in (ref) arises in nonstandard testing problems that involve nuisance parameters; we refer the reader to the work of elliott2015nearly.\footnote{See also elliott2014pre, muller2016measuring, muller2017fixed, muller2018long, guggenberger2019more, muller2020low, dou2021generalized, li2021linear, muller2025more, muralidharan2025factorial for related problems in different applications.}
We are interested in the computational aspects of finding the most powerful test of size $\alpha$ for the testing problem in (ref).\footnote{The most powerful test of size $\alpha$ correctly rejects the null hypothesis as frequently as possible, while guaranteeing that the probability of incorrectly rejecting the null is bounded above by a prespecified constant $\alpha \in (0,1)$.} It is well known that such most powerful test is based on the likelihood ratio that replaces the null hypothesis by a single mixture distribution obtained from averaging $F_1, \ldots, F_M$ according to what the literature refers to as the least-favorable distribution; see, for example, Theorem 1 in lehmann1948most. It is also known, but perhaps less so, that the least favorable distribution can be characterized as the solution of a nonlinear convex program; see Lemma 3 in krafft1967optimale. Although these results immediately suggest a computational strategy for finding the most powerful test of size $\alpha$ for (ref), the algorithms that have been thus far suggested in the econometrics literature for this purpose---for example, chiburis2008approximately, moreira2013contributions, elliott2015nearly---do not explicitly exploit the connection between finding the least-favorable distribution and solving a nonlinear convex program.
Our first result (Theorem (ref)) shows that the structure of the convex program that defines the least-favorable distribution associated with the testing problem in (ref) is amenable to the application of the methods of (stochastic) mirror descent (henceforth, S-MD); see Chapter 6.1 in bubeck2015convex and also Chapter 3 in nemirovski1983problem for textbook references. The methods of mirror descent are a family of iterative procedures designed for finding an approximate solution of convex problems in high dimensions. These methods require repeated (but finitely many) evaluations of an unbiased estimator of the subgradient of the objective function. These methods exploit the geometry of the optimization domain, but they do not exploit any higher-order smoothness information about the optimization problem.\footnote{Mirror descent is known to outperform gradient descent in some optimization domains, such as the probability simplex; see Section 4.3 of bubeck2015convex. See also Section 5.1.1 in nemirovski2009robust.} We have emphasized the term approximate solution since a common approach in theoretical computer science, operations research, and optimization is “to relax the requirement of finding an optimal solution, and instead settle for a solution that is good enough” shmoys2011design. We follow this approach and, instead of insisting on finding the exact solution of the convex program that defines the least favorable distribution, we settle for an approximate least-favorable distribution that solves the convex program of interest, but up to a prespecified additive factor $\epsilon$ (see our Definition (ref)).
Our second result (Theorem (ref)) shows that our suggested algorithm can indeed be used to provably find---after finitely many iterations and with high probability---an approximate least-favorable distribution. Our results are probabilistic since the output of the algorithm depends on the stochastic estimator of the subgradient, and not the subgradient itself. Our theoretical results yield concrete recommendations about the algorithm's implementation, including its initial condition, its step size, the number of iterations, and the number of stochastic draws per iteration that can be used to estimate the subgradient of the objective function. We think there at least two important implications from our theoretical results. First, the number of iterations used by the algorithm---$ \left( 4(1 - \alpha)^2 / \alpha^2 \epsilon^2 \right) \cdot \ln(M) $---scales logarithmically in $M$, which means there is no theoretical sense in which the number of iterations should be expected to scale poorly as function of how many elements there are in the null hypothesis in (ref). Second, and perhaps the most striking feature of our analysis, the unbiased estimator of the gradient used by our S-MD routine---defined in Equation (ref) in Theorem (ref)---can be based on a single Monte Carlo draw from each null distributions $F_{m}$, $m=1, \ldots, M$. As expected, taking a larger number of draws improves the approximation error of the S-MD routine, but using a small number of draws reduces the computational burden of the algorithm.
Interestingly, our suggested S-MD routine can be viewed as a slight modification of the algorithm suggested by elliott2015nearly for finding a least-favorable distribution for the testing problem in (ref). We think this is an important observation since, to the best of our knowledge, the theoretical guarantees for their procedure remain unknown. We note, however, that---in addition to the recommended initial condition, the step size, the number of iterations, and the number of stochastic draws per iteration---there is at least one important conceptual difference between the two algorithms. As we will explain later, the output of our recommended S-MD algorithm is based on averaging the output of the mirror descent routine over all iterations, as opposed to just trying to obtain the least-favorable distribution from the last update as suggested by elliott2015nearly.\footnote{Averaging the trajectories of a stochastic gradient-descent routine is commonly referred to as Polyak-Ruppert averaging. See ruppert1988efficient and polyak1992acceleration. See also forneron2024estimation for a discussion of Polyak-Ruppert averaging in the context of estimation and inference by stochastic optimization of nonlinear econometric models.} While it might be possible to derive theoretical guarantees for the last iterate, this will likely require considering step sizes that vary with each iteration (and in our case, the step size is fixed throughout the whole routine).
Finally, we analyze the extent to which the output of the S-MD routine allows us to construct the most powerful test for the problem in (ref). If it were possible to extract an exact least-favorable distribution, the most powerful test could be obtained directly. However, since we only obtain an approximate least-favorable distribution we need additional work to show that we can construct a nearly optimal test. We define a test $\varphi$ to be $(\epsilon,\delta)$-nearly optimal (Definition (ref)) whenever i) the size of $\varphi$ is at most $\alpha +\delta$; and ii) up to the additive constant $\epsilon$, the test $\varphi$ rejects the null hypothesis as frequently as the most powerful test of size $\alpha$. Based on this definition, we say that a test is nearly optimal if it is $(\epsilon,\delta)$-nearly optimal for some parameters $\epsilon$ and $\delta$. Our third result (Theorem (ref)) shows that it is indeed possible to use the S-MD routine to construct a randomized nearly optimal test with high probability. The nearly optimal test is shown to be the average of each of the tests of the Neyman-Pearson form associated with each vector of multipliers generated by the S-MD routine. We present explicit expressions for its size distortion and its power loss relative to the best test. We note that an important challenge in reporting such a test is that, in principle, it requires keeping track of the history of multipliers obtained from the S-MD routine. When both $M$ and the number of iterations of S-MD are large, this could come at a significant computational cost. To address this issue, we show that there is a simple strategy to implement the average test: we randomize the number of iterations uniformly between $1$ and our recommended $T$, update the multipliers, and use the resulting test to decided whether or not to reject the null hypothesis. The resulting procedure, thus, can be thought of as the Neyman-Pearson (or likelihood ratio test) associated with the last iterate of the S-MD routine that stops randomly before our suggested number of rounds.
{ Related Literature:} When the sample space $\mathcal{Y}$ is infinite, the mathematical problem that defines the most powerful test of a given size in the problem (ref) is an infinite linear programming problem: Since the testing problems analyzed in this paper posit $M$ null distributions for the observed data, the corresponding linear program has finitely many constraints but a choice variable of infinite dimension; see Section (ref) below for details. If the data were discrete-valued or if we were to discretize it---as suggested by chiburis2008approximately, moreira2013contributions, and moreira2019optimal---then one could find the most powerful test of a given size by using any algorithm for finite linear programming. krafft1967optimale is the seminal reference for using linear programming methods to characterize the most powerful test of a given size.
elliott2015nearly also show how the most powerful test for composite hypotheses can be expressed as a minimax decision problem where a false rejection of $\mathbf{H}_0$ induces a loss of 1, and a false rejection of $\mathbf{H}_1$ induces a loss of $\phi>0$ (correct choices have loss of zero). In this problem the decision maker chooses a test $\varphi$. An adversarial nature decides which element in $\{F_1, F_2, \ldots, G\}$ to use to generate the data; consequently, a mixed strategy for nature can be represented by a vector in the simplex of $\mathbb{R}^{M+1}$. One could then use an algorithm for solving the corresponding maximin problem---for example, the Hedge algorithm suggested by fernandez2024epsilon; the fictitious play algorithm suggested by guggenberger2025numerical; or a general convex optimization routine as in chamberlain2000econometric. An important limitation of this approach is that one would need to solve the maximin problem repeatedly for different values of $\phi$, until one finds a test with correct size.
Although there are no theoretical results showing that one of these algorithms is better than another for the purpose of finding a nearly optimal test, we think that our analysis illustrates how a rich literature in optimization can be leveraged to provide theoretical results about the performance of different algorithms and also to provide practical recommendations regarding their implementation.
{ Outline:} The rest of the paper is organized as follows. Section (ref) presents notation, the statement of the hypothesis testing problem of interest, and the primal and dual optimization problems that arise when searching for the most powerful test of a given size. Section (ref) presents a formal definition of a stochastic mirror descent routine (S-MD) for convex optimization and also our main results: namely, that the algorithm provably generates an approximate least favorable distribution and a nearly optimal test. Section (ref) uses an elementary testing problem that arises in the context of the univariate Gaussian location model to illustrate our main results. Section (ref) discusses some extensions and Section (ref) concludes. The proofs of our main theorems and supporting lemmas are collected in Appendix (ref). Additional results are collected in Appendix (ref).
We first present the formal statement of the hypothesis testing problem analyzed in this paper. We follow the notation and terminology used in elliott2015nearly as close as possible. We then present the convex program associated with the least-favorable distribution.
We observe a random element $Y$ that takes values in some space $\mathcal{Y}$ endowed with $\sigma$-algebra $\mathcal{F}$. Let $\nu$ denote a $\sigma$-finite measure defined over the measurable space $(\mathcal{Y},\mathcal{F})$. Let $F_1, \ldots, F_{M}$ denote $M>1$ candidate probability measures for the distribution of $Y$ under the null hypothesis. Let $G$ be the candidate distribution of $Y$ under the alternative hypothesis. Assuming all of these distributions are absolutely continuous with respect to $\nu$, Theorem 5.5.4 in dudley02 guarantees the existence of nonnegative integrable functions $f_1, \ldots, f_M, g$, which can be taken as the probability density functions of $F_1, \ldots, F_{M}, G$ relative to $\nu$.
Based on a single observation of $Y$, the testing problem of interest is
Using the typical jargon of hypothesis testing problems, the null hypothesis in (ref) is composite, since it contains more than one possible distributions for the data. The alternative hypothesis is simple, in that it contains a single distribution.\footnote{As explained in Section 2.2 of elliott2015nearly, the density $g$ can arise by appealing to the weighted average power criterion in cases where the alternative hypothesis is composite as well.}
A statistical test for (ref) (or simply a test) is a measurable function $\varphi:\mathcal{Y} \rightarrow [0,1]$, where $\varphi(y)$ is interpreted as the probability of rejecting the null hypothesis given that data $y$ were observed. A test $\varphi$ is said to be nonrandomized if $\varphi(y) \in \{0,1\}$ for $\nu$-almost every realization of $Y$; otherwise the test is said to be randomized.
The rate of Type I error under $f_{m}$ is the probability of rejecting the null hypothesis when $Y \sim f_{m}$ and it equals $\int \varphi f_{m} d \nu$. As usual, the size of a test is the largest rate of Type I error under the null hypothesis. The power of a test is the probability of rejecting the null hypothesis when $Y \sim g$ and it equals $\int \varphi g d \nu$.
We would like to find the most powerful test of size $\alpha$ for the problem (ref). By definition, such a test correctly rejects the null hypothesis as frequently as possible, but guarantees that the probability of incorrectly rejecting the null is bounded above by the prespecified constant $\alpha \in (0,1)$. Mathematically, the problem of finding the most powerful test of size $\alpha$ can be written as:
We refer to the optimization problem in (ref) as the primal problem associated with the hypothesis testing problem in (ref). We note that the primal problem is an infinite linear programming problem, in the sense of anderson1987linear. The infinite linear program in (ref) has finitely many constraints but a choice variable of infinite dimension.
Define the Lagrangian function associated with the optimization problem (ref) as
where we refer to $\kappa \equiv (\kappa_1,...,\kappa_M) \in \mathbb{R}_{+}^{M}$ as the Lagrange multipliers (or simply, multipliers) associated with each of the inequality constraints in the primal problem (ref).
Consider thus the optimization problem on $\mathbb{R}_{+}^{M}$ with variable $\kappa \equiv (\kappa_1,...,\kappa_M)$ given by
where
We refer to the problem in (ref) as the dual problem of (ref).
We have explained how the solution to the dual problem in (ref) can be used to construct the most powerful test of a given size. We now discuss the computational aspects of solving the dual problem. We start by showing that the objective function in (ref) is convex over $\mathbb{R}^{M}_{+}$ and presenting a formula for its subgradient. In particular, we show that the vector collecting the negative of the excess rate of Type I error of the test $\varphi_{\kappa}$ in (ref) is a subgradient of $f(\cdot)$ at $\kappa$.
Lemma (ref) shows that dual problem in (ref) has, as expected, a convex objective function, and we present a simple formula for a subgradient. We now show that the dual problem can be further simplified by restricting the multipliers to belong to the bounded domain:
Lemma (ref) and (ref) show that in order to find the multipliers associated with the dual program in (ref), it is sufficient to solve a convex optimization problem over a bounded domain (in particular, an $\ell_{1}$-ball around the origin with radius $1/\alpha$) instead of solving a convex optimization problem over all of $\mathbb{R}^{M}_{+}$. An intuitive explanation of this result can be given as follows. Consider the hypothesis testing problem \[ \mathbf{H_0}: Y \sim f_0 \textrm{ vs. } \mathbf{H}_1: Y \sim f_1. \] Suppose that $Y$ is real-valued and that both $f_0$ and $f_1$ are absolutely continuous with respect to Lebesgue measure. The Neyman-Pearson lemma implies that most powerful test of size $\alpha$ rejects if and only if \[ \frac{f_1}{f_0} > c_\alpha, \] where $c_{\alpha}$ is the critical value that satisfies \[ P_0 \left( \frac{f_1}{f_0} > c_{\alpha} \right) = \alpha. \] Markov's inequality trivially implies that
Thus $c_{\alpha} \leq 1/\alpha$. Lemma (ref) can be viewed as a generalization of this simple observation. In the next section we show that---due to Lemma (ref) and (ref)---the structure of the dual problem in (ref) is amenable to the application of the methods of stochastic mirror descent.
This section presents our main results. First, we present a formal definition of a stochastic mirror descent (S-MD) routine for convex optimization. The members of the mirror descent family are indexed by what is called a mirror map. We apply S-MD to the problem in (ref) by setting the mirror map to be equal to the negative entropy. This choice is motivated by our Lemma (ref) and the fact that negative entropy is a common recommendation for convex problems over $\ell_1$ balls; see Section 5.7 in nemirovski2009robust and Example 2 in srebro2012convex. Second, we define an approximate least-favorable distribution and show that S-MD provably obtains an approximate least favorable distribution (Theorem (ref)). Finally, we define a nearly optimal test and show that the S-MD routine can be used to generate such a test (Theorem (ref)).
This section follows as closely as possible the notation in Sections 4.1 and 6.1 of bubeck2015convex. Let $\mathbb{R}^{M}_{++}$ denote the set of all strictly positive vectors in $\mathbb{R}^M$. We say that a map $\Phi:\mathbb{R}_{++}^{M} \rightarrow \mathbb{R}$ is a mirror map if it satisfies the following properties
Mirror Maps are used to build iterative algorithms for constrained optimization problems when unbiased estimators of the gradient are available. More precisely, consider the optimization problem \[ \inf_{\kappa \in \mathcal{X}} f(\kappa), \] where $f:\mathcal{X} \rightarrow \mathbb{R}$ is a convex function and $\mathcal{X} \subseteq \mathbb{R}_{+}^{M}$. Suppose that $\widehat{G}(\kappa)$ is an unbiased estimator of a subgradient of $f$ at $\kappa$, in the sense that $\mathbb{E}\left[ \widehat{G}(\kappa) \right]$ is a subgradient of the function $f$ at $\kappa$.\footnote{The expectation should be understood as being conditional on $\kappa$, since $\kappa$ is stochastic. See Chapter 6 in bubeck2015convex.} Let $D_{\Phi}(\kappa,\kappa')$ denote the Bregman divergence associated with $\Phi$, that is \[ D_{\Phi}(\kappa,\kappa') \equiv \Phi(\kappa) - \Phi(\kappa') - \nabla \Phi(\kappa') (\kappa-\kappa'). \] The stochastic mirror descent algorithm (henceforth, S-MD) given the mirror map $\varphi$ is defined as follows:
The general interpretation of the S-MD update in equation (ref) is that “the method is trying to minimize the local linearization of the function while not moving too far away from the previous point, with distances measured via the Bregman divergence of the mirror map”; see p. 301 in bubeck2015convex.
An important observation regarding Algorithm (ref) is that its output is the average value of $\kappa_{t}$ over all the iterations, and not its last value. Averaging the trajectories of a stochastic optimization routine is commonly referred to as Polyak-Ruppert averaging; see ruppert1988efficient and polyak1992acceleration. This idea goes back to seminal work on mirror descent by nemirovski1983problem.
The following theorem specializes the mirror descent routine in Algorithm (ref) to the dual problem in (ref). The mirror map is set to be equal to the negative entropy.
It is useful to make an explicit connection between the updating formula in (ref) and Equation 10, p. 782 in elliott2015nearly. Following the notation in elliott2015nearly, define \[ \mu^{t+1}_{m} \equiv \ln(\kappa_{t+1,m}). \] Taking logarithms on both sides of (ref) and using the definition of $\widehat{G}_{N}(\kappa_{t})$ yields
When $c_{t}=1$, the term $\ln(c_t)=0$, and thus, (ref) essentially matches Equation 10, p. 782 in elliott2015nearly after noting that $\widehat{G}_{m,N}(\kappa_t)$ is a Monte Carlo estimate of the negative excess rate of Type I error the test $\varphi_{\kappa_{t}}$: \[ \alpha-\int \varphi_{\kappa_t} f_m d \nu. \] Other than notation, the main difference between our expressions is the presence of the additional term $c_{t}$. In the stochastic mirror descent routine, this term is used to take into account the fact that---because of our Lemma (ref)---the optimization domain in the dual problem (ref) can be restricted to values of $\kappa$ such that $\sum_{m=1}^{M}\kappa_m \leq 1/\alpha$.
In this subsection we show that if the number of epochs ($T$) and the step size ($\eta$) are chosen appropriately, then $\bar{\kappa}_{T} \equiv (1/T) \sum_{t=1}^{T} \kappa_{t}$---obtained using Equations (ref) and (ref) in Theorem (ref)---indeed can be used to generate an approximate least-favorable distribution for the problem in (ref).
In order to formalize this statement, note that we have defined $\bar{v}$ as the value of the dual problem $\inf_{\kappa \in \mathcal{X}} f(\kappa)$, where $f(\cdot)$ is defined in (ref) and the set $\mathcal{X}$ is defined in (ref). Let
denote the probability simplex in $\mathbb{R}^{M}$.
We now present a result showing that Algorithm (ref) provably generates an approximate least-favorable distribution. For any nonnegative real number $x$, let $\lceil x \rceil$ denote the “ceiling function”; that is smallest integer larger than $x$.
Theorem (ref) shows that---even after finitely many iterations---the S-MD routine for the dual problem, \[ \inf_{\kappa \in \mathcal{X}} f(\kappa),\] generates an approximate least-favorable distribution (in the sense of our Definition (ref)) with high probability. The approximate least favorable distribution in (ref) is the direction of the multipliers $\overline{\kappa}_T$ (in analogy to what we would do if we had access to the exact solution of the dual problem). The probabilistic statement in the theorem arises due to the randomness in the gradient estimator in (ref), which makes the output of the S-MD routine behave as a random variable. Note that if we fix the frequency at which we would like to obtain an approximate least-favorable distribution (over different runs of the S-MD routine), then the number of draws used to construct the gradient estimator ($N$) determines how close we get to finding a “good enough” solution for the dual problem.
For example, suppose that $\epsilon=.1$, $\alpha=10\%$ and $M=200$, and suppose we set $\Omega = \sqrt{\ln(1/\alpha)}$, so that $1-\exp(-\Omega^2) = 1-\alpha = 90\%$. If we run the S-MD routine using $N=1$ (only one draw per density $f_m$), then with probability 90% we will obtain a \[ \left(1 + 2 \sqrt{\frac{\ln(1/\alpha)}{\ln(M)(1-\alpha)^2}} \right)\epsilon \approx 2.5 \epsilon = .25 \] least-favorable distribution. If we use $N=10$ (only ten draws per density) we get a \[ \left(1 + 2 \sqrt{\frac{\ln(1/\alpha)}{\ln(M)\cdot 10 \cdot (1-\alpha)^2}} \right)\epsilon \approx 1.5 \epsilon = .15 \] least favorable distribution. If we use $N=100$, with probability 90% we get a $1.15\epsilon = .11$-least favorable distribution. More generally, Theorem (ref) shows that, for any target probability, we can always make $N$ large enough to get as close as we would like to the desired approximation error $\epsilon$.
In our view, the most surprising part of Theorem (ref) is that even if the number of draws per density used to implement the S-MD routine are as low as $N=1$, it is still possible to get an approximate least-favorable distribution that provides a non-trivial approximate solution to the dual problem in (ref). For instance, in the example above, using only one draw per density yields an approximation error of $2.5 \epsilon =.25$ with probability 90%. The approximation error of $.25$ should be interpreted as a worst-case guarantee that applies to any testing problem of the form (ref). As we show in our illustrative example, the resulting approximation error can, in practice, be considerably smaller.
Now that we have established that the S-MD routine in Algorithm (ref) (with negative entropy as a mirror map) provably generates an approximate least-favorable distribution---in the sense of our Definition (ref)---we discuss the extent to which the S-MD routine can also be used to generate a nearly optimal test.
Before presenting a formal definition of what we mean by a nearly optimal test, it is helpful to explain why it is not entirely trivial to translate the approximate least-favorable distribution in Equation (ref) into a nearly optimal test. Let $\kappa^*_{T}$ be the multipliers that we obtain after running the S-MD routine, with $T$ and $\eta$ defined as in Theorem (ref). Consider the test of the Neyman-Pearson form, $\varphi_{\kappa^*_{T}}$, defined in (ref) based on the multipliers $\kappa^*_{T}$.\footnote{Note that $\kappa^*_{T}$ can be decomposed into its direction $\lambda^*_{T}$ as in Equation (ref) and its norm $\textrm{cv}^*_{T} \equiv \sum_{m=1}^{M}\kappa^*_{T,m}$. Thus, the test in (ref) rejects the null if and only if \[ g(y) > \textrm{cv}^*_{T} \sum_{m=1}^{M} \lambda^*_{T,m} f_m(y). \] } Since the S-MD routine never explicitly tried to enforce size control, it is possible that the size of $\varphi_{\kappa^*_{T}}$ is strictly above the nominal level $\alpha$. Mathematically, this happens because an approximate optimizer to the dual problem does not necessarily imply a feasible solution to the primal problem (let alone a nearly optimal one). This suggests that, when defining a nearly optimal test, it could be helpful to take into account i) possible violations of the required size; ii) as well as potential loss in power, relative to the optimal solution.
Let $\bar{v}$ be defined as value function of the dual problem in (ref). As we show in Section (ref) of Appendix (ref), duality holds, and $\bar{v}$ equals the power of the most powerful test of size $\alpha$.\footnote{For the sake of exposition, we deliberately write $\bar{v}$ instead of $\bar{v}(\alpha)$.}
The following theorem shows that Algorithm (ref) can provably be used to generate a nearly optimal test.
Theorem (ref) shows that the S-MD routine analyzed in this paper provably generates a nearly optimal test, in the sense of our Definition (ref). There are three aspects about our result that are worth highlighting.
First, it is rather surprising that the nearly optimal test in (ref) takes the form of an “average” test. In order to get some intuition of why this construction is helpful to obtain theoretical results, it is useful to explicitly write the dual in (ref) as the minimax problem \[ \min_{\kappa \in \mathbb{R}_{+}^{M}} \max_{\varphi} L(\varphi, \kappa), \] where $L(\varphi,\kappa)$ is the Lagrangian function defined in (ref). In this problem, the “min” player is choosing a vector of (Lagrange) multipliers, and the “max” player is choosing a test. For a fixed $\kappa$, the best response of the max player is a test of the Neyman-Pearson form $\varphi_{\kappa}$ defined in (ref). A mirror descent routine for this problem initializes the choice of $\kappa$ by the “min” player, and iteratively updates its values based on the (sub)gradient of the Lagrangian with respect to $\kappa_t$, which---by results analogous to the envelope theorem---will give the rates of Type I error of $\varphi_{\kappa_{t}}$. The most powerful test of size $\alpha$ is the solution to the maximin problem \[ \max_{\varphi} \min_{\kappa \in \mathbb{R}_{+}^{M}} L(\varphi, \kappa). \] The question is how to translate the iterates, $\{\kappa_{t}\}_{t=1}^{T}$, into a solution to the maximin problem. This question is common in the application of the mirror descent algorithm to minimax problems that arise in game theory and statistical decision theory; see fernandez2024epsilon and the discussion of matrix games in arora2012multiplicative. While these papers consider different problems to the one studied in this paper, a suggestion therein is to use the best responses of the max player $\{\varphi_{\kappa_{t}}\}_{t=1}^{T}$ and randomize over them with uniform probability. In our problem, such a construction becomes the average test in (ref); and this is what motivated us to study its performance. To the best of our knowledge, our results have not been stated elsewhere.
Second, the upper bound on the rate of Type I error features a term whose value changes with each specific run of the S-MD routine. This is not ideal, but we were not able to derive a better bound. To better understand the role of this term, consider again the example we discussed after Theorem (ref). Suppose that $\epsilon=.1$, $\alpha=10\%$ and $M=200$, and suppose we set $\Omega = \sqrt{\ln(1/\alpha)}$, so that $1-\exp(-\Omega^2) = 1-\alpha = 90\%$. If we run the S-MD routine using $N=1$ (only one draw per density $f_m$), \[ 1 + \left[1 + 2\sqrt{\frac{\ln(1/\alpha)}{\ln(M)(1-\alpha)^2}} \right]\epsilon \approx 1+2.5 \epsilon = 1.25. \] If the term
were not part of (ref), then we could conclude that with probability at least 90%, the test $\bar{\varphi}_{T}$ has size of at most $12.5\%$ (that is, there is a size distortion of 2.5%) and power that is no less than $\bar{v}$ minus 25 percent points. Again, making $N$ arbitrarily large makes the size closer to $\alpha(1+\epsilon)$ and the power at least $\bar{v}-\epsilon$. The interpretation of Theorem (ref) changes slightly when we incorporate (ref). Suppose for example that in one run of the S-MD routine the term in (ref) equals .05. Then, for that run, our best hope is that (ref) is satisfied with the larger bound $1.3$ instead of $1.25$. This means that we could see a rate or Type I error of the average test as high as 13%. As we discuss in the next section, in our illustrative example the term in (ref) tends to be small (and negative), but we do not have any theoretical guarantees for this.\footnote{The terms $-\widehat{G}_{N}(\kappa_{t})^{\top} \kappa_{t}$ is a Monte-Carlo estimate of the average excess rate of Type I error of $\varphi_{\kappa_{t}}$, evaluated at the different null densities. The term in (ref) averages these Monte-Carlo estimates over all iterations.}
Third, to report the entire test (as function of all possible data), one would have to retain the history of multipliers obtained from the S-MD routine. When both $M$ and the number of iterations of S-MD are large, this could come with significant computational and data storage expense. However, a typical use case is to perform the test on a specific data set. For this purpose, data storage requirements can be much reduced because, rather than retaining the weights for every epoch, it suffices to store the implied test results, i.e. one bit per epoch and parameter value being tested (and even less if one is content with only updating the average). Furthermore, rather than reporting a rejection probability, it usually suffices to either accept or reject, although in cases where the rejection probability is interior, this decision will then be random conditionally on the data. But this can be achieved at much lower computational expense: (ii) By a simple application of the Law of Iterated Expectations, rejecting with probability corresponding to the average test is equivalent to first drawing a “realized epoch” $t^* \sim \operatorname{unif}(\{1,\ldots,T\})$ and then executing the Neyman-Pearson test for epoch $t^*$ only. (iii) As $t^*$ can be drawn before starting the iteration, it is only necessary to execute iterations up to epoch $t^*$, threreby only executing $T/2$ iterations in expectation. By using these simplifications, it should be easy to execute the test.
An alternative to the average test is to simply report the Neyman-Pearson test in (ref) associated with the multipliers $\bar{\kappa}_{T}$ obtained as the output of Algorithm (ref) based on the mirror map $\Phi(\kappa) = \sum_{m=1}^{M} \kappa_m \ln(\kappa_{m})$ and the unbiased estimator of the gradient $\widehat{G}_{N}(\cdot)$ defined in Equation (ref) of Theorem (ref). While it is challenging to analyze the size and power properties of this test for finite $T$, in Appendix (ref) we show that, under some conditions, as $T \rightarrow \infty$ the power of the test will converge to $\bar{v}$, and it will have correct size (in a sense we make precise).
It is common to construct confidence regions by inverting tests. However, the test advocated here is in general randomized, raising the question of how to invert it. Conceptual discussions of this matter go back at least to the 1950's Stevens50,Lehmann59.\footnote{A notable difference to our setting is that these discussions were motivated by the randomized nature of optimal or exact tests in highly discrete sample spaces.} To summarize some key points, for this paragraph only let the test function $\rho(\cdot)$ also depend on the parameter value to be tested, i.e. we temporarily define $\rho:\mathcal{Y} \times \Theta \to [0,1]$, where $\rho(y,\theta)$ is the probability of rejecting the instance of $H_0$ characterized by parameter value $\theta$ given data $y$; Lehmann59 calls this the critical function. Then we can define no less than four intervals that arguably invert our test:
Mirroring discussions in the previous subsection, idea 4 is the computationally most involved and its justification is asymptotic, idea 3 will be the computationally easiest, and idea 1 is intermediate; 2 is computationally equivalent to 1. Some users might find 2 hard to interpret CasellaBerger05. We leave further analysis of the issue to future research.
In order to illustrate the performance of the SM-D routine in Algorithm (ref)---along with the implications of the theoretical guarantees in Theorem (ref) and Theorem (ref)---we consider an elementary testing problem that arises in the context of the univariate Gaussian location model. More precisely, suppose we observe a realization of the random variable \[ Y \sim \mathcal{N}(\theta,1). \] The location parameter, $\theta$, is unknown to the econometrician. Let $\Theta_0 \equiv \{\theta_{0,1}, \ldots, \theta_{0,M}\}$ be an equally-spaced grid over the interval $[-5,0]$ consisting of $M=200$ points. We order the elements of $\Theta_0$ in decreasing order, so that $\theta_{0,1} = 0$. We also define the singleton set $\Theta_1 \equiv \{\theta_1\}$.
We assume that the econometrician is interested in the following hypothesis testing problem: \[ \textbf{H}_0: \theta \in \Theta_0 \quad \textrm{vs.} \quad \textbf{H}_1: \theta \in \Theta_1. \] It is well known that the most powerful test of size $\alpha$ for this problem---which we denote as $\varphi^*_{\alpha}$---rejects the null if $Y$ is large enough. More precisely, $\varphi^*_{\alpha}(Y) \equiv \mathbf{1}\{ Y \geq z_{1-\alpha} \}$, where $z_{1-\alpha}$ is the $1-\alpha$ quantile of a standard normal. The power of this test can then be expressed in terms of the normal c.d.f. as $\Pr(N(0,1) \leq \theta_1 - z_{1-\alpha})$. Since the most powerful test of size $\alpha$ for this example is known, we can use this information to analyze the theoretical guarantees we provided in Theorem (ref) and Theorem (ref).
The Stochastic Mirror Descent (S-MD) Routine: We first obtain a nearly optimal test of size $\alpha = 10\%$. We set $\theta_1 = 2$, which means that the largest power of a test of size $\alpha=10\%$ is $\Phi(2-1.28) \approx 76.38\%$. This is also the value of the dual problem in (ref). We set $\epsilon = .1$, and use the formulae in Theorem (ref) to determine the maximum number of iterations ($T$) and the learning rate ($\eta$) for the S-MD routine:
In this example $M = 200 > (\exp\left(1\right)/\alpha) = \exp(1)\cdot10$. Thus, in accordance with our theoretical derivations, the initial condition for the S-MD routine ($\kappa_0 \in \mathbb{R}^{M}$) is chosen to be: \[ \kappa_0 = \left( 1/(\alpha M) , \ldots, 1/(\alpha M) \right)^{\top} = (.05, \ldots, .05)^{\top}. \]
The main component of the S-MD routine is the stochastic mirror descent update. We implement the unbiased estimator of the gradient in Theorem (ref) using only one draw per density; that is, $N=1$. More precisely, if we let $f_{m}(\cdot)$ denote the p.d.f. of $Y$ under the null hypothesis $\theta_{0,m} \in \Theta_0$ and we let $g(\cdot)$ be the p.d.f. of $Y$ under the alternative $\Theta_1$, the mirror descent update necessitates an unbiased estimator of the rates of Type I error of the test
In our example, an unbiased estimator for the rate of Type I error at $\theta_{0,m}$ can be succinctly obtained by sampling $Z \sim N(0,1)$ and using $\varphi_{\kappa}(Z + \theta_{0,m})$ as an estimator. More precisely, in each epoch $t$ we obtain one draw $Z_t \sim \mathcal{N}(0,1)$ and compute the mirror descent update (coordinate by coordinate) as \[ \kappa_{m,t+1} \equiv \kappa_{m,t} \cdot \exp\left( \eta \left[ \varphi_{\kappa_t}(Z_t + \theta_{0,m}) - \alpha \right] \right), \quad \text{for each } m=1, \ldots , M. \]
The intuition of the update is very simple. If $\varphi_{\kappa}(Z + \theta_{0,m})$---the unbiased estimator of the rate of Type I error of $\varphi_{\kappa_{t}}$ at $\theta_{0,m}$---is larger than $\alpha$, then $\kappa_{m,t+1}$ increases (and otherwise decreases). We also know that $\| \kappa_{t} \|_{1}$ must be less than or equal than $1/\alpha$. Thus, after the update we check if $\| \kappa_{t+1} \| \leq 1/\alpha$. If this is the case, we keep $\kappa_{t+1}$ as is; but otherwise we normalize $\kappa_{t+1}$ to guarantee that $\sum_{m=1}^{M} \kappa_{m} \leq 1/\alpha$. This gives us back the update described in part 2 of Theorem (ref).
In general, updating $\kappa_{t}$ is numerically very cheap when $N=1$, as it only involves obtaining one sample from each of the null densities (along with the evaluation of the null and alternative densities at each draw) and also evaluation of the exponential function. Using Matlab R2024a on a personal ASUS Vivobook Pro 15 @ 2.5GHz Intel Core Ultra 9 185H, it took around 280 seconds (slightly less than 5 minutes) to complete all of the $T=171,666$ iterations.
Approximate Least-Favorable Distribution: As suggested by our Theorem (ref), in order to construct an approximate least-favorable distribution we standardize the average value of $\kappa_t$ to represent it as a probability distribution. More precisely, the blue bars in Figure (ref) below correspond to \[ \lambda_{T} \equiv \frac{ \bar{\kappa}_{T} }{ \| \bar{\kappa}_{T}\|_1 }, \quad \textrm{where } \bar{\kappa}_{T} \equiv \frac{1}{T} \sum_{t=1}^{T} \kappa_{t}. \] In the testing problem we are considering, it is known that the least-favorable distribution loads all of its mass on $\theta_{0,1}$. As Figure (ref) shows, the output of the S-MD routine resembles such distribution.
Our definition of approximate least-favorable distribution in Definition (ref) makes reference to the value function of the dual problem in (ref), and not in terms of its minimizer. In this example, it is easy to show that the value of the dual (which we denoted by $\bar{v}$) equals $76.38\%$ (the power of the most powerful test of size $\alpha=10\%$). We now argue that, as expected, the distribution $\lambda_{T}$ approximately solves the dual problem, in the sense of Definition (ref). To see this, we just need to evaluate the function: \[ f(\bar{\kappa}_{T}) = \underbrace{\int \varphi_{\bar{\kappa}_{T}} g(y)dy}_{\approx 76.96\%}- \underbrace{\sum_{m=1}^{M} \bar{\kappa}_{T,m} \left( \int \varphi_{\bar{\kappa}_{T}}(y) f_m(y) dy - \alpha \right)}_{-1.01\%} \approx 77.97\%, \] where $\phi_{\bar{\kappa}_{T}}$ is the test of Neyman-Pearson form defined in (ref), and where a Monte-Carlo approximation with 100,000 draws is used for the evaluation of each of the integrals. Thus, in this example: \[f(\bar{\kappa}_{T}) \approx 77.97\% < \bar{v} + \epsilon = 76.38\% + 10\% = 86.38\%. \] This means that in our run of the S-MD routine we obtained an $\epsilon=.016$-least-favorable distribution. The quality of the approximation is much better than what we expected based on our theoretical results in Theorem (ref). This is consistent with the fact that the results in Theorem (ref) apply to every possible testing problem with $M$ null densities and a single alternative. We also conducted a Monte-Carlo simulation where we implemented the S-MD routine with different draws for the estimation of the subgradient, with 10,000 draws being used for integral evaluation in each. In all of the 100 runs we obtain a $10\%$-least favorable distribution. This is consistent that the theoretical results we presented in Theorem (ref) apply to any testing problem of the form $\eqref{eqn:testing_problem_intro}$.
Nearly Optimal Test $\bar{\varphi}_T$: Figure (ref) reports the test $\bar{\varphi}_T$ (red, solid line), which is the test defined in (ref). For comparison, we also report the test $\varphi_{\bar{\kappa}_T}$ (blue, solid line). The size of $\bar{\varphi}_{T}$ is approximately $10.22\%$, and its power is approximately $76.98\%$. This means that the test $\bar{\varphi}_{T}$ is slightly over-sized, but it has competitive power.
We also conducted 100 different runs of our S-MD routine. In all of the 100 hundred runs the size was at most $\alpha(1+\epsilon)$ and the lowest power achieved was 76.82%.
Time Comparison of Using More Draws in the S-MD Routine: We also analyzed the increased computational effort of increasing the number of draws used to evaluate the subgradient during each iteration. To do this, we re-ran our illustrative example using 1, 10, 100, and 1,000 draws in each round. Table (ref) shows the runtime as a function of the number of draws. Lastly, we also calculated the expected runtime associated with 20,000 draws in each round---the number of draws recommended by elliott2015nearly.\footnote{To be clear, they recommended to have a fixed set of draws at the beginning of the iterations and to reuse those draws via importance sampling at each iteration. While this may be computationally more feasible for certain distributions, most theoretical results on S-MD require the draws in each period to be independent of the history (otherwise, it is difficult to guarantee unbiasedness of the estimator of the subgradient). The theoretical performance guarantee with their recommendations may be derived with additional assumptions, but remains unknown to the best of our knowledge.} Our expected runtime was 369,697 seconds.\footnote{This estimate was generated by taking the time difference between the first and second round, and then multiplying by $T$. Note that this estimate is smaller than extrapolation from the last two observations based on a linear trend.} We think these results illustrate the computational gains of implementing the S-MD routine with a small number of draws to evaluate the subgradient.
Consider now a variation of the testing problem in (ref) where the alternative hypothesis is also composite, but with only $I$ possible distributions for the data:
As explained in Section 2.2 of elliott2015nearly, one can reduce the problem in (ref) to the problem in (ref) by choosing weights $w \equiv (w_1, \ldots, w_{I}) \in \Delta^{I-1}$ and defining $g \equiv \sum_{i=1}^{I} w_i g_i$. Then, the test $\varphi$ that solves (ref) can be interpreted as the test that maximizes $w$-weighted average power among all tests of size at most $\alpha$.
A common criticism of tests that maximize a weighted average power criterion (henceforth, WAP) is that they can be biased: their power for some density $g_i$ can be lower than $\alpha$moreira2013contributions,andrews2016conditional. moreira2013contributions note that one could include additional constraints in the problem (ref) and consider:
Just as before, we can define the Lagrangian function associated with problem (ref) as
where we refer to $\kappa \equiv (\kappa_1,...,\kappa_M) \in \mathbb{R}_{+}^{M}$ as the Lagrange multipliers associated with each of the inequality constraints that bound the test's rate of Type I error, and we let $\mu \equiv (\mu_1,...,\mu_I) \in \mathbb{R}_{+}^{I}$ denote the Lagrange multipliers associated with each of the inequality constraints preventing the test to be biased.
We could then proceed as we did before and define the dual optimization problem:
where
It is possible to show that the function $f(\kappa,\mu)$ is convex in its arguments. Moreover, a test $\varphi$ that achieves the maximum is the test
and a subgradient of $f(\cdot)$ at $(\kappa,\mu)$ is \[\nabla f(\kappa,\mu) \equiv -\left(\int \varphi_{\kappa,\mu} f_1 d \nu - \alpha,...,\int \varphi_{\kappa,\mu} f_M \nu d \nu - \alpha, \int \varphi_{\kappa,\mu} (-g_{1})d\nu + \alpha \ldots \int \varphi_{\kappa,\mu} (-g_{I})d\nu + \alpha \right).\]
If $\nabla f(\kappa,\mu)$ were known, the mirror descent routine (with negative entropy as mirror map) for this problem would have updates
Establishing a result similar to Theorem (ref) and (ref) is more challenging because the application of standard results would need ex-ante constraints on the $\|\cdot\|_{1}$-norm of $\kappa$ and $\mu$. But, for example, if one knew that the optimal values of $\kappa$ and $\mu$ satisfied the constraint $\| \kappa + \mu \|_{1} < 1/\alpha$ then our previous theoretical results for S-MD would apply.
We showed that---in testing problems where the null hypothesis postulates $M$ distributions for the observed data---one can use a stochastic mirror descent routine to provably obtain---after finitely many iterations---both an approximate least-favorable distribution and a nearly optimal test. The convex program that arises naturally in the testing problem in (ref) is the dual of the mathematical program that defines the most powerful test of size $\alpha$.
Our theoretical results allowed us to provide concrete recommendations about the algorithm's implementation: including its initial condition, its step size, the number of iterations, and the number of stochastic draws per iteration that can be used to approximate the subgradient of the objective function. These practical recommendations have at least two important implications. First, the number of iterations used by the algorithm scales logarithmically in $M$ (which means there is no theoretical sense in which the algorithm scales poorly as a function of the elements in the null hypothesis). Second, the algorithm can be implemented with a single stochastic draw per null density in each iteration (taking a larger number of draws improves the approximation error of the S-MD routine, but using a small number of draws reduces the computational burden of the algorithm). Importantly, our suggested algorithm coincides with a slight variation of the algorithm in elliott2015nearly.
{ Proof of Part 1 of the Theorem:} Let $\kappa_{t}$ be a realization of an arbitrary $\mathcal{X}$-valued random vector. Let $\widehat{G}_{m,N}(\kappa_t)$ be the $m$-th coordinate of $\widehat{G}_{N}(\kappa_{t})$. If suffices to show that if $(Y_{m,1}, \ldots, Y_{m,N})$ are i.i.d. random variables with distribution $Y \sim f_m$ sampled independently of the realized value of $\kappa_{t}$, then $\mathbb{E}[\widehat{G}_{m,N}(\kappa_t) | \kappa_t] = -\left(\int \varphi_\kappa f_m d\nu - \alpha \right)$.
By definition of $\widehat{G}_{m,N}(\kappa_t)$,
where the last line follows from the fact that $(Y_{m,1}, \ldots, Y_{m,N})$ are i.i.d. according to $f_m$, independently of the value of $\kappa_{t}$. Since $f_m$ is the p.d.f. of $Y$ relative to the $\sigma$-finite measure $\nu$, Problem 1, Chapter 5, p. 177 of dudley02 implies \[ \mathbb{E}\left[ \varphi_{\kappa_t}(Y_{m,n}) \Big | \kappa_{t} \right] = \int \varphi_{\kappa_{t}} f_m d \nu. \] Therefore, $\widehat{G}_N(\kappa_t)$ is an unbiased estimator of the subgradient of $f$ at the realized $\kappa_t$.
{ Proof of Part 2 of the Theorem:} Let $\widehat{G}_{N}(\kappa_{t})$ be the unbiased estimator of the subgradient of $f$ at $\kappa_{t}$. We provide an explicit solution for the problem
when $\Phi(\kappa) = \sum_{m=1}^{M} \kappa_{m} \ln(\kappa_{m})$. By definition of Bregman divergence, \[D_\Phi(\kappa,\kappa_t) = \Phi(\kappa) - \Phi(\kappa_t) - \langle \nabla \Phi(\kappa_t), \kappa - \kappa_t \rangle.\] Consequently, $\kappa_{t+1}$ is the solution to the following optimization problem
subject to the constraint \[ \sum_{m=1}^{M} \kappa_{m} \leq 1/\alpha. \] Let $\mu$ denote the Lagrange multiplier associated with this constraint. Thus, the first-order conditions of the problem for each $\kappa_{t+1,m}$ become:
where $\widehat{G}_{m,N}(\kappa_{t})$ is the m-th entry of $\widehat{G}_{N}(\kappa_{t})$. The first-order condition in (ref) can be written as \[ \kappa_{t+1,m} = \kappa_{t,m} \exp \left(-\eta \widehat{G}_{m,N}(\kappa_{t}) \right) \exp \left(-\mu \right). \] Two cases to consider. First, if \[ \sum_{m=1}^{M} \kappa_{t,m} \exp \left( -\eta \widehat{G}_{m,N}(\kappa_{t}) \right) < 1/\alpha, \] then $\mu=0$ and
Second, if \[ \sum_{m=1}^{M} \kappa_{t,m} \exp \left( -\eta \widehat{G}_{m,N}(\kappa_{t}) \right) \geq 1/\alpha, \] then $\mu>0$ and $\sum_{m=1}^{M} \kappa_{t+1,m}$ must equal $1/\alpha$. Consequently,
which can be achieved by setting \[ \mu = \ln \left( \alpha \sum_{m=1}^{M} \kappa_{t,m} \exp \left(-\eta \widehat{G}_{m,N}(\kappa_{t}) \right) \right). \]
{ Proof of Part 3 of the Theorem:} The initial condition $\kappa_1$ solves:
We re-parameterize this problem by defining \[ K \equiv \| \kappa \|_{1}, \quad p_{m} \equiv \kappa_{m}/K, \quad p = (p_1, \ldots, p_{M})^{\top}.\] Since $\kappa \in \mathbb{R}^{M}_{++}$, then $K>0$ and $w_{m}>0$ for all $m=1,\ldots, M$. Moreover, if we denote by $\Delta^{M-1}$ the simplex in $\mathbb{R}^{M}$ and use int$\left(\Delta^{M-1}\right)$ to denote its interior, the optimization problem in (ref) thus becomes the nested optimization problem
Thus, we first solve the inner problem which consists of finding the distribution in the simplex with the smallest negative entropy: \[\min_{p \in \textrm{int}\left(\Delta^{M-1}\right) } \sum_{m=1}^M p_m \ln (p_m).\] It is known that the solution of this problem is to set $p_{m}=1/M$. We verify this below for the sake of exposition. The first order conditions are \[ 1+ \ln(p_{m}) + \mu = 0,\] where $\mu$ is the Lagrange multiplier associated with $\|p\|_1 = 1$. Solving for $p_m$ yields \[ p_m = \exp \left( - (1+\mu) \right) \] which implies (by summing the left side over $m=1,\ldots, M$): \[ \frac{1}{M} = \exp \left( - (1+\mu) \right). \]\ We conclude that $p^*_{m} = 1/M$ is the optimal direction of $\kappa_1$. We now find its scale by solving the outer optimization problem
Without the constraint, the objective function has a global minimum at $K^*$ satisfying \[ \ln(K^*) + 1 - \ln(M) = 0, \] or equivalently, $K^*= M / \exp(1)$. It is also decreasing for $K<K^*$ and increasing for larger values. Therefore, the solution to the problem in (ref) is \[K^* =
,\] Thus, the initial condition is \[\kappa_{1} =
\]
The approximation error of the numerical iteration consists of two parts: optimization error and estimation error. The former is intrinsic to the optimization algorithm when applying the exact (sub)gradient, while latter is induced by the estimation error of the unknown subgradient.
\global\long\setcounter{page}{1}