Faster Algorithms and Constant Lower Bounds for the Worst-Case Expected Error
Year: 2022

Abstract:
The study of statistical estimation without distributional assumptions on data values, but with knowledge of data collection methods was recently introduced by Chen, Valiant and Valiant (NeurIPS 2020). In this framework, the goal is to design estimators that minimize the worst-case expected error. Here the expectation is over a known, randomized data collection process from some population, and the data values corresponding to each element of the population are assumed to be worstcase. Chen, Valiant and Valiant show that, when data values are `∞-normalized, there is a polynomial time algorithm to compute an estimator for the mean with worst-case expected error that is within a factor π2 of the optimum within the natural class of semilinear estimators. However, their algorithm is based on optimizing a somewhat complex concave objective function over a constrained set of positive semidefinite matrices, and thus does not come with explicit runtime guarantees beyond being polynomial time in the input. In this paper we design provably efficient algorithms for approximating the optimal semilinear estimator based on online convex optimization. In the setting where data values are `∞-normalized, our algorithm achieves a π2 -approximation by iteratively solving a sequence of standard SDPs. When data values are `2-normalized, our algorithm iteratively computes the top eigenvector of a sequence of matrices, and does not lose any multiplicative approximation factor. Further, using experiments in settings where sample membership is correlated with data values (e.g. "importance sampling" and "snowball sampling"), we show that our `2-normalized algorithm gives a similar advantage over standard estimators as the original `∞-normalized algorithm of Chen, Valiant and Valiant, but with much lower computational complexity. We complement these positive results by stating a simple combinatorial condition which, if satisfied by a data collection process, implies that any (not necessarily semilinear) estimator for the mean has constant worst-case expected error.

text: The study of statistical estimation without distributional assumptions on data values, but with knowledge of data collection methods was recently introduced by Chen, Valiant and Valiant (NeurIPS 2020). In this framework, the goal is to design estimators that minimize the worst-case expected error. Here the expectation is over a known, randomized data collection process from some population, and the data values corresponding to each element of the population are assumed to be worstcase. Chen, Valiant and Valiant show that, when data values are `∞-normalized, there is a polynomial time algorithm to compute an estimator for the mean with worst-case expected error that is within a factor π2 of the optimum within the natural class of semilinear estimators. However, their algorithm is based on optimizing a somewhat complex concave objective function over a constrained set of positive semidefinite matrices, and thus does not come with explicit runtime guarantees beyond being polynomial time in the input. In this paper we design provably efficient algorithms for approximating the optimal semilinear estimator based on online convex optimization. In the setting where data values are `∞-normalized, our algorithm achieves a π2 -approximation by iteratively solving a sequence of standard SDPs. When data values are `2-normalized, our algorithm iteratively computes the top eigenvector of a sequence of matrices, and does not lose any multiplicative approximation factor. Further, using experiments in settings where sample membership is correlated with data values (e.g. "importance sampling" and "snowball sampling"), we show that our `2-normalized algorithm gives a similar advantage over standard estimators as the original `∞-normalized algorithm of Chen, Valiant and Valiant, but with much lower computational complexity. We complement these positive results by stating a simple combinatorial condition which, if satisfied by a data collection process, implies that any (not necessarily semilinear) estimator for the mean has constant worst-case expected error.

heading: 1 Introduction
text: Standard methods in statistical analysis are often based on the assumption that data values are drawn independently from some underlying distribution, and may perform poorly when this assumption is violated. Methods from robust statistics often modify this assumption, for example by allowing a small fraction of data values to be arbitrary outliers while the bulk of the data values remain independent samples. However, there are natural settings in which data values may be strongly correlated to sample membership, so that distributional assumptions (e.g. independence of data values) can lead to inaccurate results. Furthermore, inaccuracy of predictions in such settings can have serious societal consequences–for example predictions of election outcomes or the spread of a contagious disease.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
Recently, Chen, Valiant, and Valiant [4] introduced a framework to capture statistical estimation in such difficult settings, by assuming that the data values are worst-case, but that the estimation algorithm is able to leverage knowledge of the randomized data collection process. The framework is simple to describe. There is a set {1, . . . n} of n indices with corresponding data values x = {x1, . . . xn} along with a probability distribution P over pairs of subsets A,B ⊆ {1, . . . n}. Here A is called the sample set and B is called the target set. A sample (A,B) is drawn from P and the values xA = {xi | i ∈ A} are revealed. The goal is then to estimate the value of some function g(xB) given only A, B and xA. The quality of an estimator in this framework is measured by its worst-case expected error. Here the data values x are worst-case and the expectation is with respect to the distribution P .
The worst-case expected error framework captures several natural settings where data values may be correlated to sample membership. Importance sampling is the setting where the target set B is always equal to the whole population {1, . . . n}, and the sample set A is chosen by independently sampling each element i ∈ {1, . . . n} with probability pi. In the worst-case expected error framework, the values xi may be arbitrarily correlated to the probability pi of being sampled. In snowball sampling [12], the elements of {1, . . . n} correspond to the vertices of a graph (e.g. a social network). A sample set A is chosen by first picking a few initial vertices. Next, each chosen vertex recruits a randomly chosen subset of its neighbors into the sample, and this process is repeated until a desired sample size is reached. The target set B can be either the whole population or a subset given by a few additional iterations of the recruitment process. Here it is natural to assume that the data values xi at neighboring vertices will be highly correlated. Finally, selective prediction[10, 18] is a temporal sampling method where the population {1, . . . n} corresponds to time steps, and the goal is to predict the average of future data values given the past (e.g. the average change in the stock market). The target set B is some time window {t, . . . t+ w} and the sample A is {1, . . . t}. If the random process for choosing the starting point t and length w of the time window is chosen appropriately, sub-constant error is attainable even when the xi are chosen to be worst-case values bounded by a constant[10, 18].
In [4], the authors design an algorithm for computing estimators for the mean of the target set B. In particular, they restrict their attention to the class of semilinear estimators i.e. estimators which are restricted to be a linear combination of the sample values xA, but where the weights of the linear combination may depend arbitrarily on P , A, and B. In the setting where the data values are bounded, the authors design an algorithm that outputs a semilinear estimator with worst-case expected error at most a π2 factor larger than that of the optimal semilinear estimator. They then demonstrate that the estimator output by their algorithm has significantly lower expected error than many standard estimators in the importance sampling, snowball sampling, and selective prediction settings where data values are correlated to sample membership. However, several open questions remain, even for the case of computing the optimal semilinear estimator for the mean.
First, the algorithm in [4] is a concave maximization over a constrained set of positive semidefinite matrices, which is solved using a general purpose convex programming solver. In particular this means that no explicit bounds on the runtime of the algorithm are given beyond being polynomial in n. This leads us to ask:
Question 1: Is there a simpler algorithm to compute a semilinear estimator with approximately optimal worst-case expected error that comes with explicit runtime guarantees?
Second, the choice of `∞-normalization (i.e. assuming the data values are bounded) is natural for some settings such as polling or testing for infectious disease. However, `2-normalization (i.e. 1 n ∑ i x 2 i ≤ 1) is also a natural choice for settings where some data values may be much larger than others such as predicting stock prices. This motivates:
Question 2: Is there an efficient algorithm to compute a semilinear estimator with approximately optimal worst-case expected error in the `2-bounded case?
Finally, the setting of the worst-case expected error is quite challenging, and so one would expect that for many randomized data collection processes it is not possible to achieve non-trivial worst-case expected error. Thus, one might ask:
Question 3: Are there simple properties of a data collection process that ensure that no estimator can achieve sub-constant worst-case expected error?

heading: 1.1 Our Results
text: We answer the first two questions by designing simple algorithms based on online gradient descent for approximating the optimal semilinear estimator in both the `∞-bounded and the `2-bounded case. To be able to describe our algorithms we first formally define the input.
Definition 1.1 (Definition 2 in [4]). A joint sample-target distribution P over {1, . . . n} is given by the uniform distribution over m pairs (Ai, Bi) where Ai, Bi ⊆ {1, . . . n}.
Because any distribution P on pairs (A,B) may be approximated to arbitrary accuracy by a uniform distribution this is a reasonable choice of parametrization. In many settings the m uniform pairs are obtained by sampling sufficiently many times from the known randomized data-collection process (see supplementary material for details). Thus, one should think of m as being much larger than n, say at least Ω(n2). Our algorithms’ runtimes also depend on a parameter ρ which we define to be n times the worst-case expected error of the optimal semilinear estimator.
Theorem 1.2 (Informal–see Theorem 2.3). There is an algorithm based on online gradient descent for approximating the optimal semilinear estimator to within a multiplicative π2 and additive error for `∞-bounded values. The algorithm runs in O(ρ2 log ρ) iterations each taking Õ(mnω−1 + n7/2) time, where ω is the current matrix multiplication exponent.
Note that when there exists a semilinear estimator achieving worst-case expected error O( 1n ), the algorithm requires O(1) iterations. Each iteration of the algorithm requires solving a standard SDP in n dimensions with n constraints (in fact the only constraints are that the diagonal entries of the psd matrix are constrained to all be equal to one).
Theorem 1.3 (Informal–see Theorem 2.4). There is an algorithm based on online gradient descent for approximating the optimal semilinear estimator to within an additive error for `2-bounded values. The algorithm runs in O(ρ2 log ρ) iterations each taking O(mn) time.
In addition to the improved runtime, each iteration of the algorithm for `2-bounded values only requires computing the eigenvector corresponding to the largest eigenvalue of an n × n matrix, making the algorithm practical and simple to implement. We provide the details of both algorithms in Section 2.
Recall that the SDP algorithm from [4] obtained lower expected error than several standard estimators on various synthetic datasets where data values are correlated to sample membership. To compare to these empirical results obtained for the `∞-bounded case, we implement our algorithm for `2-bounded values and perform experiments on the same synthetic datasets from [4]. We find that the expected error of our algorithm on these synthetic datasets matches that of the SDP algorithm from [4], while being simpler to implement and having lower computational complexity. See Section 5 for details.
Finally, in Section 4 we answer the third question by providing a simple combinatorial condition for a sample-target distribution P , which ensures that any (not necessarily semilinear) estimator has constant worst-case expected error.

heading: 1.2 Related Work
text: The most closely related work is [4] which defined the framework for the worst-case expected error and designed the first algorithms in this setting. The papers [10, 18] on selective prediction are also closely related as they design efficient estimators in the worst-case expected error setting for a specific choice of data collection distribution.
There has been extensive recent work on the theory and practice of robust estimation in high dimensions, where some fraction of data values are arbitrarily corrupted and the rest remain independent samples from the distribution [5, 16, 3, 6, 7, 8]. Most closely related to our work is the recent paper of Hopkins, Li and Zhang on robust mean estimation using online convex optimization/regret minimization [13]. While that work focuses on a quite different problem setting, the high-level idea of using methods in online convex optimization to design simple and efficient algorithms for robust estimation is closely connected to our work.

heading: 1.3 Preliminaries
text: We focus on the setting where the data values x are real numbers and the goal is to estimate the mean of the target set. Since the mean scales linearly with the data values, it is natural to assume that the values x lie in some bounded set D ⊆ Rn. For our purposes D will either be the `∞-ball or the `2-ball.
Definition 1.4. Given a sample-target distribution P and a bounded set D ⊆ Rn, an estimator f(xA, A,B) is a real-valued function which takes as input the data values xA and the sample-target index sets (A,B). The worst-case expected error of f on P is given by
max x∈D
1
m m∑ i=1 (f(xAi , Ai, Bi)−mean(xBi)) 2
We use ‖v‖ to denote the `2-norm of a vector v and ‖v‖∞ to denote the `∞-norm. For a matrix M we use ‖M‖ to denote the operator norm. For a subspace W ⊆ Rn, we use W⊥ to denote the orthogonal complement subspace, and we write ΠW for the orthogonal projection onto W .

heading: 2 Algorithms for the Worst-case Expected Error
text: In this section we design algorithms for approximating the best semilinear estimator for the mean using tools from online convex optimization. As in [4] we focus on the class of semilinear estimators which compute a linear combination of the data values xA where the weights of the linear combination can depend on the sets A and B. For convenience we introduce the notation Wi ⊆ Rn to denote the subspace of vectors that have non-zero coordinates only on indices in Ai. Then a semilinear estimator takes the form f(xAi , Ai, Bi) = 〈ai, x〉 where each vector ai ∈Wi. To further simplify notation we denote by bi ∈ Rn the vector which takes value 1|Bi| on coordinates in Bi and is zero otherwise, so that mean(xBi) = 〈bi, x〉. Using this notation we have: Definition 2.1. Given a sample-target distribution P and a set D ⊆ Rn the worst-case expected error of the optimal semilinear estimator is given by
opt(P ) = min {ai∈Wi}mi=1 max x∈D
1
m m∑ i=1 〈ai − bi, x〉2
SDP Relaxation. The first step to approximating the best semilinear estimator is to introduce an SDP relaxation for the inner maximization in Definition 2.1. To simplify notation we make the following definition.
Definition 2.2. Let a = {a1, . . . am} where ai ∈Wi for all i. Define
M(a) = 1
m ∑ i (ai − bi)(ai − bi)>.
Then the inner maximization in Definition 2.1 can be written more simply as
max x∈D
x>M(a)x
When D = {x ∈ Rn | ‖x‖∞ ≤ 1} the inner maximization can be relaxed to the semidefinite program
sdp∞(a) = max X 0,Xj,j=1
〈M(a), X〉 (1)
which has an optimal value within a factor π2 of the true worst-case expected error as shown via the rounding method of [17]. See [4] for a self-contained explanation. Therefore to approximate the best semilinear estimator it suffices to approximately solve the min-max optimization problem
min a sdp∞(a) (2)
Algorithm via Online Gradient Descent. Our first theorem states that we can approximately minimize sdp∞(a) with an application of online gradient descent. Theorem 2.3. Let P be a sample-target distribution, ρ = n · opt(P ), and > 0. There is an algorithm which computes a set of m vectors a′ = {a′1, . . . a′m} such that
sdp∞(a ′) ≤ min
a sdp∞(a) + . The algorithm runs in O ( ρ2 log ρ 2 ) iterations, and Õ(ρ 2 log ρ 2 (mn
ω−1 + n7/2 log(1/ )) time, where ω is the current matrix multiplication exponent.
Note that the parameter ρ depends on the optimal value of the worst-case expected error. In particular, if the optimal worst-case expected error is O( 1n ) then ρ = O(1), and the iteration count of the algorithm is O( 1 2 ). The algorithm is an application of online gradient descent where in each step we solve an instance of sdp∞(a) and use the solution as a convex cost function. The analysis is based on standard regret bounds for online gradient descent.
Algorithm 1: Online gradient descent for minimizing sdp∞(a). Input: A sample-target distribution P , accuracy parameter , upper bound on optimum p. Output: A sequence of m vectors a = {a1, . . . am} with ai ∈Wi for all i.
1 Let r = √
πmp 2 2 Let β = ∑m i=1‖ΠW⊥i bi‖ 2. 3 Let a(1)i = 1 |Ai|1Ai for all i. 4 for t from 1 to 36π 2n2p2
2 do 5 Let a(t) = {a(t)1 , . . . a (t) m }, and let ηt = mn√t .
6 Let X(t) be a ( 1 + 10 ) -approximate solution to sdp∞(a
(t)). 7 a
(t+1) i ← a (t) i − ηt 2mΠWiX (t)(a (t) i − bi) for all i. 8 Set λ(t) = min { 1, √ r2−β∑m
i=1‖a (t+1) i −ΠWibi‖2
} .
9 a (t+1) i ← λ(t)a (t+1) i + (1− λ(t))ΠWibi.
10 return a(t ∗) where t∗ = argmint〈M(a(t)), X(t)〉.
We now give the analysis of the algorithm using standard regret bounds for online gradient descent (see the text [11]). See the supplementary material for the referenced lemmas.
Proof. The main observation is that the algorithm is precisely online gradient descent where the convex cost function observed at each step is given by ft(a) = 〈M(a), X(t)〉. Indeed viewing a as a vector in nm dimensions (one for each ai) the component of the gradient ∇ft(a) corresponding to ai is 2mΠWiX (t)(ai − bi). Using the positive semidefinite Grothendieck inequality, one can show
that for the optimal solution a∗ lies in Br∗(b), the `2-ball of radius r∗ = √ πm opt(P ) 2 centered at b = {b1, . . . , bm} (Lemma A.1). Thus given an upper bound p ≥ opt(P ) it is sufficient to limit the search to Br(b) for r = √ πmp
2 . The last lines of the loop are just projection onto this ball (Lemma A.3). Finally, for a in Br(b), a straightforward calculation (Lemma A.2) gives ‖∇ft(a)‖ ≤ 2nrm . Thus the algorithm is online gradient descent with feasible set diameter D = 2r and gradients bounded by G = 2nrm . Thus by the textbook analysis of online gradient descent [11] we have
1
T T∑ t=1 ft(a (t))− min a′∈Br(b) 1 T T∑ t=1 ft(a ′) ≤ 3GD 2 √ T = 3πnp√ T .
Letting a∗ = argmina sdp∞(a) (which we know is contained in Br(b)) we have
1
T T∑ t=1 ft(a (t))− 1 T T∑ t=1 ft(a ∗) ≤ 3πnp√ T .
Noting that ft(a∗) = 〈M(a∗), X(t)〉 ≤ sdp∞(a∗) yields
1
T T∑ t=1 ft(a (t)) ≤ sdp∞(a∗) + 3πnp√ T ≤ sdp∞(a∗) + 2
where the last line follows from the choice of T = 36π 2n2p2
2 . For t∗ = argmint ft(a (t)) we therefore
have 〈M(a(t ∗)), X(t ∗)〉 = ft(a(t ∗)) ≤ sdp∞(a∗) +
2 .
Since the value of sdp∞(a ∗) is bounded, the fact thatX(t ∗) is a multiplicative (1+ 10 )-approximation implies that X(t
∗) is an additive 2 -approximation to sdp∞(a (t∗)) (see Lemma A.4). Thus we
conclude that
sdp∞(a (t∗)) = max
X 0,Xjj=1 〈M(a(t
∗)), X〉 ≤ 〈M(a(t ∗)), X(t ∗)〉+ 2 ≤ sdp∞(a∗) +
i.e. a(t ∗) is an -additive-approximation of mina sdp∞(a).
To analyze the runtime note that each iteration requires approximately solving an instance of sdp∞(a (t)), multiplying the solutionX(t) by each vector a(t)i −b (t) i , and then rescaling by λ
(t). First, the matrix M(a) can be constructed in time O(mnω−1) by grouping the vectors a(t)i − b (t) i into m n blocks of size n×n, and using fast matrix multiplication for each block. Next, using the fact thatM(a) is positive semi-definite, the approximate solution X(t) can be computed in time Õ(n7/2 log(1/ )) using the interior point SDP solver of [14] (see Lemma A.4 for details). Finally, we can group the vectors a(t)i − b (t) i columnwise into m n matrices of size n× n, and then use fast matrix multiplication to compute the product of each matrix with X(t). This results in a runtime of O(mnω−1) where ω is the current matrix multiplication exponent. Finally, λ(t) can be computed in O(mn) time. Putting it all together there are O(p 2n2 2 ) iterations each of which takes Õ(mn ω−1 + n7/2 log(1/ )) time, given an upper bound p ≥ opt(P ). We can find an upper bound that is at most 2 opt(P ) by starting with p = 1n and repeatedly doubling at most log(n opt(P )) = O(log ρ) times. This yields the final iteration count of (ρ 2 log ρ 2 ).
`2-norm Bounded Values. We now turn to the setting where the data values are bounded in `2-norm i.e. when D = {x ∈ Rn | ‖x‖ ≤ √ n}. In this case the inner maximization in Definition 2.1 is equal to n times the maximum eigenvalue of M(a)
sdp2(a) = max ‖x‖2≤ √ n x>M(a)x = n‖M(a)‖. (3)
The choice of normalization is such that the second moment is one i.e. 1n ∑n j=1 x 2 j = 1. A subtle difference in the `2-bounded setting is that opt(P ) may not be bounded by a constant for certain pathological examples. This is unlike the `∞-bounded setting where choosing a to be the all zeros vector immediately gives an upper bound of opt(P ) ≤ 1. However, if the target distribution is the full population mean (i.e. bi = 1n1 for all i), then setting a to all zeros and applying Cauchy-Schwarz yields
opt(P ) ≤ max ‖x‖2≤ √ n
1
m m∑ i=1 〈bi, x〉2 ≤ 1 m m∑ i=1 1 n · n = 1.
More generally, it only makes sense to compute estimators approximating the worst-case expected error when the achievable error goes to zero with n i.e. opt(P ) = o(1). Therefore, for the `2-bounded setting we make the additional assumption that opt(P ) ≤ 1. We now state our main theorem for the `2-bounded setting. Theorem 2.4. Let P be a sample-target distribution with opt(P ) ≤ 1, ρ = n · opt(P ), and > 0. There is an algorithm which computes a set of m vectors a′ = {a′1, . . . a′m} such that
sdp2(a ′) ≤ min
a sdp2(a) + . The algorithm runs in O ( ρ2 log ρ 2 ) iterations and takes O ( ρ2 log ρ 3 mn ) time.
If the optimal worst-case expected error is opt(P ) = O( 1n ), then O( 1 2 ) iterations suffices, and each iteration takes O(mn ) time. The algorithm for the `2-bounded case is also an application of online gradient descent, where in each step we solve an instance of sdp2(a) and use it as a convex cost function. However, the fact that sdp2(a) corresponds to computing the maximum eigenvalue of an n× n matrix leads to the improved runtime per iteration of O(mn ).
Algorithm 2: Online gradient descent for minimizing sdp2(a). Input: A sample-target distribution P , accuracy parameter , upper bound on optimum p. Output: A sequence of m vectors a = {a1, . . . am} with ai ∈Wi for all i.
1 Let r = √ mp 2 Let β = ∑m i=1‖ΠW⊥i bi‖ 2. 3 Let a(1)i = 1 |Ai|1Ai for all i. 4 for t from 1 to 36n 2p2
2 do 5 Let a(t) = {a(t)1 , . . . a (t) m }, and let ηt = mn√t .
6 Let X(t) = xtx>t , where xt is a ( 1 + 10 ) -approximate solution to sdp2(a
(t)). 7 a
(t+1) i ← a (t) i − ηt 2mΠWiX (t)(a (t) i − bi) for all i. 8 Set λ(t) = min { 1, √ r2−β∑m
i=1‖a (t+1) i −ΠWibi‖2
} .
9 a (t+1) i ← λ(t)a (t+1) i + (1− λ(t))ΠWibi.
10 return a(t ∗) where t∗ = argmint〈M(a(t)), X(t)〉.
To summarize briefly, the `2-bounded algorithm can be essentially obtained by modifying Algorithm 1 by letting X(t) = xtx>t where xt is an 2 -additive-approximate eigenvector of M(a
(t)). The improved running time follows from two points. First, the classical power method can compute the required eigenvector in mn time. Second, computing the gradient requires multiplyingX
(t)(a (t) i −bi)
for all i, just as in Algorithm 1. However, since in the `2-bounded case X(t) = vv> is rank one, each multiplication can be carried out in O(n) time, for a total cost of O(mn). The full details of the analysis can be found in Section B in the supplementary material.

heading: 3 The algorithm of Chen, Valiant and Valiant
text: At this point it is instructive to compare Algorithm 1 with the original SDP-based algorithm of [4]. In Appendix C of [4] the authors show how to approximate any sample-target distribution P by drawing at mostm = poly(n ) samples. Thusm should be generally thought of as a large polynomial in n. The approach taken in [4] to approximate the optimal semilinear estimator is based on first replacing the inner maximization in Definition 2.1 with sdp∞(a). This only incurs an error of at most a factor of π2 , and additionally means that the problem is convex in a = {a1, . . . am} and linear in X . Thus the min-max theorem applies and the min and max can be exchanged, resulting in an inner minimization over a which can be solved explicitly. This gives rise to a somewhat complicated semidefinite program where the objective is a concave maximization depending on the inverse of the variable matrix.
The authors note that this SDP can be converted to a more standard SDP by using the Schur complement to re-express the matrix inverse. However this conversion increases the number of constraints in the SDP to at leastm. The best known interior point solver [14] has runtime Õ( √ n(mn2 +mω+nω)) for a general SDP with an n dimensional PSD matrix variable and m constraints. The dominant term in our setting is Õ( √ nmω). For example, if m = O(n3) then the runtime for the general SDP solver
in [4] is Õ(n7.61), where we have used the fact that the current matrix multiplication exponent is approximately ω ≈ 2.37. The runtime of Algorithm 1 in this setting is O(ρ2n4.37). Depending on the achievable worst-case expected error, ρ can range from O(1) to O(n), so the runtime of Algorithm 1 ranges from O(n4.37) to O(n6.37), all of which are faster than O(n7.61) for [4].
To summarize, after making the appropriate reductions, the algorithm of [4] requires solving an SDP with O(m) constraints. In contrast, the SDP solved in each iteration of Algorithm 1 has exactly n
constraints. Therefore, given that SDP solvers have runtimes polynomial in the number of constraints, and that in our setting we should think of m n, Algorithm 1 will tend to have a faster runtime even though it solves an SDP in each iteration.

heading: 4 A constant lower bound
text: In this section we identify a simple combinatorial condition that ensures the worst-case expected error of any estimator for the mean is at least a constant. Of course it is possible to construct trivial examples of sample-target distributions P for which a constant lower bound holds. For example, consider a distribution P where the samples Ai ⊆ {1, . . . n2 } and Bi ⊆ { n 2 + 1, . . . n} for all i. In this case, given any estimator f(xA, A,B), a worst-case adversary is free to set the data values xj for j ≥ n2 arbitrarily in order to maximize (f(xAi , Ai, Bi)−mean(xBi))
2, as these values never appear as input to f . Indeed, since the estimator has zero probability of ever observing one of the data values xj for j ≥ n2 , there is no hope of estimating the mean of these values under worst-case assumptions. To go beyond the trivial case it makes sense to require that every data value has some non-negligible probability under P of being included in a sample set A. We will show that even in this case, the worst-case expected error can be constant. In fact, the condition we identify can hold for distributions P where each coordinate i is equally likely to be included in A. Definition 4.1. Let α > 0 be a constant. A sample-target distribution P is α-non-expanding if there exists a subset S ⊆ {1, . . . n} such that for an α-fraction of the pairs (Ai, Bi) exactly one of the following holds:
1. Ai ⊆ S and Bi ∩ S = ∅
2. Ai ∩ S = ∅ and Bi ⊆ S
For example, the definition is satisfied with α = 1 by the distribution P which half the time picks a uniform random subset A from the first half of indices and B from the second half, and half the time does the opposite. More subtle examples are possible where the subset S is arbitrary and α < 1, so that a constant fraction of sample and target sets can have arbitrary intersection with S. Theorem 4.2. Let P be an α-non-expanding sample-target distribution. For any estimator f(xA, A,B)
max x:‖x‖∞=1
1
m m∑ i=1 (f(xAi , Ai, Bi)−mean(xBi)) 2 ≥ α 4
Proof. Let S be the subset provided by Definition 4.1. Without loss of generality we may assume that Ai ⊆ S and Bi ∩ S = ∅ at least for at least an α2 fraction of pairs (Ai, Bi). Otherwise we could just switch the roles of S and S̄. Let c be the median of f(1Ai , Ai, Bi) on this α 2 fraction. If c is positive let E be the at least α4 fraction of indices i for which f(1Ai , Ai, Bi) ≥ c. If c is negative let E be the at least α4 fraction of indices i for which f(1Ai , Ai, Bi) < c. Set xj = 1 for all j ∈ S and xj = − sign(c) for all j /∈ S. Since Ai ⊆ S and Bi ∩ S = ∅ for all i ∈ E we have
1
m m∑ i=1 (f(xAi , Ai, Bi)−mean(xBi)) 2 ≥ 1 m ∑ i∈E (f(xAi , Ai, Bi)−mean(xBi)) 2
= 1
m ∑ i∈E (f(1Ai , Ai, Bi) + sign(c)) 2
≥ α 4 (c+ sign(c)) 2 ≥ α 4
where the final inequality follows from the fact |c+ sign(c)| ≥ 1 for all c ∈ R.

heading: 5 Experimental Results
text: In this section we empirically evaluate Algorithm 2 in settings where data values are correlated to inclusion in a sample. In [4] the authors show that solutions to argmina sdp∞(a) clearly outperform
standard choices for estimators in such settings. However, their algorithm requires using the general purpose convex programming solver MOSEK [2] via the CVXPY package [9, 1], which can be quite slow and memory intensive. We implement Algorithm 2 in Python, and run the algorithm on the same synthetic datasets as those in [4]. For the algorithm of [4] we use the publicly available code at https://github.com/justc2/worst-case-randomly-collected. We find that Algorithm 2 has comparable empirical performance to the original algorithm of [4], while being simpler and more computationally efficient. For more details of the experiments please see the supplementary material.
Importance Sampling Importance sampling, where elements are sampled independently but with different probabilities, is one of the simplest examples where data values may be correlated with sample membership. In this experiment the population size is n = 50 and the element i is included in the sample with probability 0.1 for i ≤ 25 and probability 0.5 for i > 25. We compare both Algorithm 2 and the sdp∞-based algorithm from [4] with two standard estimators for this setting. The first estimator is reweighting which computes the weighted mean of the xi, where the weight of each xi is the reciprocal of the probability that xi is included in the sample. The second is subgroup estimation which computes the sample means of xi for i ≤ 25 and xi for i > 25 separately, and then averages these two sample means.
We evaluate the results on three synthetic datasets: (1) Constant where xi = 1 for all i. (2) Intergroup variance where xi = 1 for i ≤ 25 and xi = −1 for i > 25. (3) Intragroup variance where xi = 1 for odd indices i and xi = −1 for even i. We also report worst-case `∞ and `2 error by solving sdp∞(a) and sdp2(a) for each estimator. The expected squared errors for each estimator and dataset in this experiment appear in Table 1. Note that on the three synthetic datasets the sdp∞ algorithm from [4] and Algorithm 2 have essentially identical expected error, while the former has lower sdp∞ error and the latter has lower sdp2 error as one might expect. Furthermore, the error on the synthetic datasets is consistently low for Algorithm 2, while the standard estimators each have large error for at least one setting where data values are correlated to sample membership.
Snowball Sampling In our snowball sampling experiment we randomly draw n = 50 points in the two dimensional unit square to be the population and let the target set B be the entire population. We construct a sample by first picking a random starting point, and then iteratively adding to the sample two of the five nearest neighbors of each point added so far, until k = 25 points are included in the sample. We compare Algorithm 2 and the sdp∞-based algorithm from [4] with the estimator that simply computes the sample mean. We evaluate these estimators on a dataset of spatially correlated values by setting xi equal to the sum of the two coordinates of the point in the unit square corresponding to element i. We also report worst-case `∞ and `2 error by solving sdp∞(a) and sdp2(a) for each estimator. The expected squared errors for each estimator and dataset in this experiment appear in Table 2. As in the previous example, we see that the sdp∞ algorithm and Algorithm 2 have equal error on the synthetic dataset (both clearly outperforming the sample mean estimator), and each algorithm does better than the other on the worst-case data values it was designed to optimize.
Selective Prediction In the selective prediction experiment the population consists of n = 32 timesteps, and the goal is to predict some future target time window {t, . . . , t+ w} given the sample consisting of the past {1, . . . t}. Here t < n is chosen uniformly at random, andw is chosen uniformly from {1, 2, 4, 8, 16}. In this setting we compare the sdp∞ algorithm from [4] and Algorithm 2 with the selective prediction estimator of [10, 18]. The selective prediction estimator simply computes the mean of the final w elements of the sample {1, . . . t}, and is known to achieve expected worst-case error O( 1logn ). We evaluate the worst-case `∞ and `2 error of these by solving sdp∞(a) and sdp2(a) and report the results in Table 3. As in the previous cases both the sdp∞ algorithm from [4] and Algorithm 2 outperform the selective prediction estimator, and each of the two algorithms outperforms the other on the worst-case data values of the appropriate type.
Funding Disclosure
This work was partially supported by the Wallenberg AI, Autonomous Systems and Software Program (WASP) funded by the Knut and Alice Wallenberg Foundation.

Review summaries:
The authors build on the work of Chen, Valiant and Valiant from NeurIPS 2020 on estimators minimizing the worst-case expected error by providing two new algorithms for the same problem formulation, but with explicit runtime guarantees. The first algorithm uses a standard semidefinite program within a gradient descent loop to arrive at a runtime that is (with some abusive simplification) $\tilde O(mn^5 \log n/\epsilon^6)$. The second relaxes the setting so that data values are bounded in $\ell_2$ norm (rather than $\ell_\infty$), and uses this new setting to reduce the runtime to $O(mn^3 \log n/\epsilon^3)$. As far as I can tell, these algorithms don't offer any guarantees on sample complexity, but the newly explicit runtime guarantee is (subject to questions below) an improvement. As an impossibility result, the authors provide a combinatorial condition on the sample-target distribution under which the worst-case expected error for all estimators is at least a constant. Finally, the authors verify their new algorithm experimentally, showing comparable performance to [4] but faster runtime.

This paper studies the _worst-case expected error framework_ of statistical estimation proposed by Chen, Valiant, and Valiant, where data values are worst-case, but the estimation algorithm exploits randomness in the data collection process. In this worst-case expected error framework, there are $n$ data values $x_1$ to $x_n$ and a probability distribution $P$ over pairs of subsets $A, B \subseteq {1, \ldots, n}$, where the variables $x_A$ over the sample set $A$ are revealed, and we want to estimate the value of some fixed function $g(X_B)$ of the variables $x_B$ over the target set $B$. Restricting attention to _semilinear_ estimators, which are linear combinations of the sample values whose weights could depend on the sample-target distribution $P$ or the indices $x_A$ and $x_B$ of the sample and target variables, Chen, Valiant, and Valiant gave a complicated semilinear estimator for the mean of the target set $B$ which is at most $\pi/2$ factor larger than the optimal semilinear estimator, when the data values are $\ell_\infty$-bounded.  This paper presents the following three results: 1. a simpler algorithm to compute a semilinear estimator within a multiplicative $\pi/2$ factor and an additive $\epsilon$ factor for $\ell_\infty$-bounded values, with explicit efficient runtime guarantees. 2. an efficient algorithm to compute a semilinear estimator within an additive $\epsilon$ error for $\ell_2$-bounded values. 3. a family of lower bound instances where all estimators (not necessarily semilinear) must have at least constant error, which generalizes the trivial lower bounds of disjoint $A$ and $B$.

This paper studies the recently introduced *worst-case expected error* framework of Chen-Valiant-Valiant [4]. In this setting the data values are worst case, but the statistician has knowledge of the data-generating process. The authors provide new efficient algorithms based on alternating minimization-maximization for estimating the mean in worst-case under the previously studied $\ell_\infty$ assumption [4] on the dataset as well as under a new $\ell_2$ assumption. The authors also provide a simple combinatorial condition denoting when the worst-case error is $\Omega(1)$, and thus consistent estimation is not possible in this setting. Experiments compare their $\ell_2$ algorithm with the algorithm proposed by [4]  on synthetic data in the settings of importance sampling, snowball sampling, and selective prediction. 

This work builds on a recent statistical estimation framework of Chen, Valiant, and Valiant [NeurIPS`20] that leverages knowledge about how samples are collected but makes no distributional assumptions on the data values. The authors explore three important questions in this space:  (1) Is is possible to compute a semilinear estimator for the $\ell_\infty$ case that achieves the same approximation quality as that of Chen-Valiant-Valiant, while also emitting explicit running time bounds?  (2) Is there an algorithm to compute a high-quality semilinear estimator in the $\ell_2$-bounded case?  (3) Are there simple properties of a data collection process that guarantees no estimator can achieve sub-constant worst-case expected error?  The results answer all three questions affirmatively. The authors show that online gradient descent can be used to design simpler SDP-based algorithms, answering (1) and (2). They also quantify a combinatorial property of joint sample-target distributions that, for any estimator, necessarily incur at least a constant error. The paper then presents experiments based on those in [Chen-Valiant-Valiant, NeurIPS`20] to demonstrate the effectiveness of their gradient descent-based approach.