Multi-Observation Regression
Rafael Frongillo Nishant A. Mehta Tom Morgan Bo Waggoner


On Euler Emulation of Observer-Based Stabilizers for Nonlinear Time-Delay Systemsfootnoteinfo
    [
    
=============================================================================================




  Recent work introduced loss functions which measure the error of a prediction based on multiple simultaneous observations or outcomes.
  In this paper, we explore the theoretical and practical questions that arise when using such multi-observation losses for regression on data sets of (x,y) pairs. 
When a loss depends on only one observation, the average empirical loss decomposes by applying the loss to each pair, but for the multi-observation case, empirical loss is not even well-defined, and the possibility of statistical guarantees is unclear without several (x,y) pairs with exactly the same x value. 
  We propose four algorithms formalizing the concept of empirical risk minimization for this problem, two of which have statistical guarantees in settings allowing both slow and fast convergence rates, but which are out-performed empirically by the other two.
  Empirical results demonstrate practicality of these algorithms in low-dimensional settings, while lower bounds demonstrate intrinsic difficulty in higher dimensions.
  Finally, we demonstrate the potential benefit of the algorithms over natural baselines that use traditional single-observation losses via both lower bounds and simulations.



  Information elicitation, empirical risk minimization, regression, high-order statistics





§ INTRODUCTION



The inference method of empirical risk minimization (ERM) selects a hypothesis by minimizing some loss function over a data set.
The question of which loss to use in which situation, however, is often a matter of debate.
Traditional loss functions have a familiar form: (f(x),y), where f is a hypothesis, x is a feature vector of some kind, and y is a label or observation.
Notably, such losses score the accuracy of a hypothesis's prediction f(x) based on a single observation y.

Recent work <cit.> proposed a natural generalization of ERM which depends on multiple observations: (f(x),y_1,y_2,…,y_m).
These multi-observation losses[Outside of the regression setting considered here, multi-observation losses have also appeared in the context of learning embeddings <cit.>, and some techniques from time series analysis can be viewed as a special case of multi-observation regression <cit.>.]
  allow more natural and direct capture of “higher-order” properties of the conditional distribution on y given x.
Their advantages can be formalized via elicitation complexity <cit.>.
For example, to regress on the variance of y given x, any statistically consistent loss function can be shown to require either a multi-dimensional range of f(x) or multiple observations.
Further, some statistics, such as the 2-norm of the conditional distribution, can be modeled much more simply with multiple observations than with multiple dimensions <cit.>.
In this paper, we give two more examples: the upper confidence bound (UCB), and MINVAR <cit.>.
See  <ref> for more on elicitation complexity.

Despite the groundwork laid by <cit.>, a very important methodological question remains: given a data set of (x,y) pairs, how should one use such a multi-observation loss function to actually perform ERM?
For the single-observation case, up to regularization, ERM is simply _f∈∑_(x,y)∈(f(x),y).
But the analogue in the multi-observation case is _f∈∑(f(x), y_1,…,y_m), requiring data of the form (x,y_1,y_2,…,y_m), i.e., multiple observations or labels y for every feature vector x.
While data of this form is sometimes obtainable, we would like to perform ERM in the typical case where each x has only one associated y value.



  
Contributions.
In this paper, we present several algorithms to perform multi-observation regression given data sets of (x,y) pairs. These algorithms reduce to ERM after constructing metasamples of the form (x,y_1,…,y_m).
The first two are unbiased algorithms, which take care to ensure that the x values are drawn independently from the underlying distribution.
We also present two biased algorithms which perform better in practice but whose statistical guarantees are much harder to analyze due to the lack of i.i.d. x values.

For our unbiased algorithms, we prove statistical convergence guarantees, with both slow and fast rates depending on the setting. 
Thus, we give the first excess risk bounds for ERM with multi-observation losses, aside from very preliminary (and much weaker) bounds of <cit.>.
The key technique is an analysis of ERM with corrupted samples, i.e. samples from distributions near but not equal to the same underlying distribution.
We consider both labeled and unlabeled sample complexity and discuss the importance of the distinction under the multi-observation paradigm. 
Our algorithmic and sample complexity contributions are most valuable and practically relevant in the low-dimensional regime; as we show in  <ref> and  <ref>, poor dependence on the dimension d of x is information-theoretically unavoidable when predicting higher-order statistics.

Finally, we demonstrate the advantage of multi-observation regression over the traditional single-observation approach in cases such as predicting conditional variance.
The intuition is that the single-observation approach requires (in some sense) fitting several surrogate statistics (e.g. first and second moment), which may follow some very complex relationships that are learning-theoretically difficult to fit.
We give both theoretical and empirical (simulation) examples. 
We often use the variance as an example due to its simplicity, but we stress that our results are fully general, encompassing the UCB, MINVAR, and countless other properties elicitable via multi-observation losses.



  
Applications.
Aside from augmenting the literature on fundamental properties of ERM, our results have applications to robust engineering design, uncertainty quantification, and finance.
In the design of an airfoil, building, truss, etc., engineers seek designs which minimize some objective (drag, cost, weight, displacement), but which are robust to changes in the environment or to manufacturing defects <cit.>.
To understand the quality of a design, expensive computer simulations are performed, which often have a stochastic component, and the goal is to minimize some risk measure of the objective, such as the 95% quantile of the drag for an airfoil <cit.>.
Our algorithms could increase the statistical efficiency of the popular practice of surrogate optimization, where a hypothesis f(x) is fit to the risk measure, for design parameters x, and then f is minimized directly <cit.>.
We show that two popular risk measures for robust design, the upper confidence bound <cit.> and MINVAR <cit.> are easily fit with multi-observation losses.

In finance, risk measures are used both in decision-making and to regulate banks and other financial institutions <cit.>.
Typically, banks must report estimates of their risk to authorities alongside their 10-day monetary losses <cit.>, and a statistical test is performed to judge the accuracy of these risk estimates <cit.>.
To perform this test, however, one essentially needs a loss function which is statistically consistent with the risk measure, i.e., which “elicits” it <cit.>; this has played a major role in academic responses to the recent Basel 3.5 risk regulation proposal <cit.>.
Previous work <cit.> has shown that multi-observation loss functions broaden the class of elicitable risk measures, but of course a bank cannot report “i.i.d.” monetary losses for the same day.
Our results could imply that using multiple monetary losses from nearby days would suffice.





§ BACKGROUND: ERM AND ELICITATION



Classically, in supervised learning an algorithm is presented with an i.i.d. sample of n labeled points (X_1, Y_1), …, (X_n, Y_n) with the objective of selecting some action f ∈ that obtains low expected loss _(X, Y) ∼×_X[ _f(X, Y) ], or risk, with respect to some loss function . 
Here, we use the notation _f(X, Y) = (f(X), Y), so that each f is a function mapping from  to , and we denote by f^* ∈ the risk minimizer over . 
We also always assume that the loss  is L-Lipschitz in its first argument. 
ERM, which returns any f̂∈ that minimizes the empirical risk 1/n∑_i=1^n _f(X, Y), is a natural choice for solving this problem. The performance of ERM is now known to be tightly characterized by the notion of Rademacher complexity. 


Let  be a class of functions mapping from a space  to , and let Z_1, …, Z_n be an i.i.d. sample from distribution P over . 
Let ϵ_1, …, ϵ_n be independent Rademacher random variables (distributed uniformly on {-1, 1}). The Rademacher complexity of  (with respect to P) is

    _n() = [ _ϵ_1, …, ϵ_n[ sup_g ∈1/n∑_i=1^n ϵ_i g(Z_i) ] ] .


In the above, we may for instance take the space = × and class = {_f : f ∈}.

The following uniform convergence result is well-known; a proof appears in  <ref>
for completeness.
 
Let (X_1, Y_1), …, (X_n, Y_n) be independent samples from distribution P and assume, for all f ∈, that (f(X), Y) ∈ [0, B] almost surely. Further assume that the loss  is L-Lipschitz in its first argument. 
With probability at least 1 - δ,

    sup_f ∈{ [ (f(X),Y) ] - 1/n∑_i=1^n (f(X_i),Y_i) }≤ 2 L _n() + B √(log(1/δ)/2 n) .


In particular, the upper bound holds for the empirical risk minimizer f̂. A straightforward argument then leads to a high probability bound on the excess risk of f̂ (the risk of f̂ minus the risk of f^*).



  
ERM for multi-observation losses.
To apply the above analysis to multi-observation loss functions, at first glance we would need samples of the form (X,Y_1,…,Y_m), with each Y_i drawn i.i.d. from the conditional distribution _X.
In reality, however, we will likely have samples of the form (X,Y) as above, and thus we cannot apply the traditional theory directly.
Nevertheless, there is hope to recover the kinds of guarantees of traditional ERM if nearby x values have similar conditional distributions on  (the domain of Y).
We will formalize this in  <ref> with a Lipschitz assumption on the total variation distance of the conditional distributions (see (<ref>)).



  
Elicitation and ERM.
How does the choice of loss function affect the resulting hypothesis f̂ chosen by ERM?
One way to understand this relationship is to ask what statistic a loss function elicits, meaning what minimizes the loss if given only a single x value and a distribution over the y values; specifically, what is _r (r,Y) as a function of the distribution of Y?
For example, it is well-known that squared loss elicits the mean [Y] in the sense that [Y] = _r (r - Y)^2.
It can easily be shown that the behavior of ERM can be completely characterized by what statistic the loss elicits, provided that conditional statistic as a function of x is in the function class  <cit.>.

Unfortunately, given a statistic, there do not always exist losses which elicit it; one such example is the variance <cit.>.
As a result, fitting a model to the conditional variances f(x) ≈[Y|X x] directly using ERM is impossible, and instead one would typically fit models to the first two moments ([Y|X x],[Y^2|X x]), from which the variance can be obtained.
This leads to the concept of elicitation complexity <cit.>: how many elicitable quantities does one need to obtain the desired statistic?
From previous work <cit.>, we now know that allowing multiple observations can vastly reduce this complexity: with 2 observations, the variance has complexity 1 (see  <ref>), and the complexity of the 2-norm of the distribution of Y drops from ||-1 to 1 again with only 2 observations.
This improvement in complexity can have real advantages in practice, as we show in  <ref>, and is the main motivation for studying multi-observation ERM.



§ ALGORITHMS



The setting we consider is slightly unusual in terms of the number and kinds of samples used.
Our algorithms will often take a number of different (x,y) pairs and “merge” them to create a data point of the form (x,y_1,…,y_m), i.e. multiple observations associated to a single x.
We refer to such a tuple as a metasample.
In this paper, n will always denote the number of metasamples constructed and used by a particular algorithm.
This differs from , the total number of data points drawn by the algorithm.
We focus on empirical risk minimization over the metasamples,

    _f∈∑_i=1^n _f(x_i, y_i,1, …, y_i,m) .

Therefore, the algorithmic questions are (1) how to draw or choose samples, and (2) how to construct metasamples.
(In  <ref>, we will prove theoretical risk guarantees for some of these algorithms.)



  
Learning paradigms.
Our algorithms will apply in two different paradigms.
In supervised learning, the algorithms draw  data points of the form (x,y) i.i.d., construct n metasamples, and run ERM.
The sample complexity is .
In pool-based active learning, the algorithms draw  unlabeled x points.
They may query up to one label y for each x, drawn independently from _x.
This results in a smaller number of labeled pairs (x,y), from which the algorithms construct the metasamples.
For those of our algorithms with theoretical guarantees, the label complexity will always be nm, because every label we draw is used in exactly one metasample.[This is not always true of other algorithms described below and tested in simulations, where the number of metasamples is larger compared to the number of labels because each y may appear in multiple metasamples.]

One paradigm we do not consider in this paper is (fully) active learning, where algorithms may repeatedly choose any x, query it to obtain an independent draw y ∼_x, and repeat.
Here, traditional algorithms and guarantees will generally carry over to the multi-observation setting: One can generally query as many i.i.d. observations from _x as desired, so multi-observation losses present no additional difficulty.
However, we note that multi-observation regression is well-motivated for active learning settings in practice, such as crowdsourcing settings where multiple crowd workers can label the same data point.



  
Algorithmic approach.
Our algorithms attempt to mimic an idealized ERM in the following way: (1) draw  data points; (2) choose n “clumps” of data points; (3) pick a “representative” x for each clump along with m labels from each clump.
For our “unbiased” algorithms with theoretical guarantees, it is important that each representative x be an i.i.d. draw from , so we first draw n data points, then draw some number of additional data points so as to produce a “clump” around each x.
This naturally leads to a label complexity of nm and a sample complexity of some  depending on the algorithm.
For our other “biased” algorithms, we may have “overlapping” clumps.
This can have the advantage of n becoming quite large, which may help ERM perform better in practice.
However, the metasamples become correlated because they share labels, making theoretical results challenging.



 §.§ Unbiased Algorithms


We first consider algorithms for which we will later be able to prove risk bounds, by ensuring that the x values in the metasamples are i.i.d. samples from .
The Naïve algorithm, Algorithm <ref>, starts by drawing n i.i.d. data points X_1^*,…,X_n^* and using them as the basis for a metasample.
For each X_i^*, it then draws many new data points, so that with high probability, enough points come close enough to form a good metasample.
It then moves on to the next X_i+1^*.



[t]0.46
  
  
    

[t]0.53
  

    






In  <ref> we prove a lemma stating that if we set =  m n^(d+3)/2 d^d/2log2 m n/δ, then with probability at least 1-δ, most of the points (X^*_j)_j ∈ [n] have their m nearest neighbors all within a proximity of 1/√(n).  The proof partitions  into “heavy” and “light” cells according to the probability mass of each and uses Hoeffding's inequality to bound the number of X^*_j points in light cells.

To improve on Naïve Sampling, Algorithm <ref> iteratively draws batches of new samples and finds a globally good assignment to the original base sample.[Note that metasamples constructed conditional on X still ensure that Y is drawn independently from _X.]  The following lemma guarantees this algorithm's performance for = O(mn + m(√(d)/ϵ)^d logmd^d/δϵ^d).

impsamplemcntr


[Improved Uniform Sampling Lemma]lemmaimpsamplem

Let ϵ∈ (0,1), d ∈  and let x^*_1, …, x^*_n, x_1, …, x_N be drawn independently from =Unif([0,1]^d). If N ≥ C m (n + (√(d)/ϵ)^d (logm/δ+d logd/ϵ)) for a universal constant C, then with probability at least 1 - δ, for each j ∈ [n], there are at least m points x_i_j,1, …, x_i_j,m satisfying x^*_j - x_i_j,m_2 ≤ϵ, and all the i_1,1, …, i_1,m, …, i_n,1, …, i_n,m∈ [N] are distinct.


The proof appears in  <ref>, the main idea being as follows.  We partition [0,1]^d into regions of diameter ϵ and, through a Poissonization argument, show that with sufficient probability, each region has at least m x_i's for every x_j^*.   <ref> also provides an analogous lemma for nonuniform distributions, where we set  = O(m n^(d+1)/2 d^d/2logm(nd)^d/δ) and ϵ=1/√(n).






 §.§ Biased Algorithms

We now briefly consider algorithms that we feel are likely to perform well in practice and indeed do so in our simulations. 
These algorithms construct a larger number of metasamples by reusing labels, thereby giving them access to more information but making theoretical guarantees very difficult. For simplicity, both algorithms are specified for single dimensional , but they can be generalized to higher dimensions at the expense of computational complexity.
The first algorithm, Sliding Window (Algorithm <ref>), simply iterates from left to right over the x-values on the real line and creates a metasample from each group of m adjacent points.[A similar algorithm, not formalized here, partitions the points into disjoint sets of adjacent m-tuples, treating each as a metasample. We found this algorithm equally difficult to analyze and less performant in practice.]
The second, dubbed “ϵ-Nearby” (Algorithm <ref>), sets a fixed upper distance limit ϵ and constructs a metasample from all m-tuples of data points whose x-values lie in an interval of diameter ϵ.
We demonstrate the performance of these algorithms in  <ref>.

[t]0.47
  

  

[t]0.52
  

  






§ RISK BOUNDS



Let  be an arbitrary probability distribution over ⊂^d, and, for each x ∈, let _x be the conditional distribution over  given x. 
For simplicity, we take = [0, 1]^d. 
In this section, we assume that loss values lie in [0, B].

As discussed above, to make headway we will need to relate the conditional distributions D_x for nearby x; formally, we make a Lipschitz assumption on their total variation distance,



    D_TV(_x, _x') ≤ K x - x'_2 . A1

Intuitively, this means that samples from nearby conditional distributions are almost interchangeable.



 §.§ Excess risk bounds for general situations


Consider an ideal setting where we are given an i.i.d. sample of n points X^*_1, …, X^*_n and, for each i ∈ [n], we have an i.i.d. sample of m points Y_i,1, …, Y_i,m drawn from _X^*_i.
Here, existing analyses of ERM still go through, modulo the label for each X^*_i now being ψ(Y_i,1, …, Y_i,m) for some known function ψ.
Our analysis relates the performance of our algorithms, which must construct their own “noisy” metasamples from imperfect data, to this ideal.
The analysis readily provides results under both the standard notion of sample complexity from supervised learning and the (much smaller) label complexity from pool-based active learning.


A key idea in the analysis is to view each metasample as being drawn in this idealized fashion (i.e. each Y_i,j∼_X^*_i), but with some chance of corruption.
We show this is possible by viewing any nearby distribution _X' as a mixture of _X^*_i with an arbitrary corruption.
We can then analyze ERM on a set of metasamples, most of which are ideal, but some of which are corrupted.

  
Assume that (<ref>) holds. Let 
= C m n^(d+1)/2 d^d/2logm(nd)^d/δ for some universal constant C, and let f̃ be the hypothesis returned by Algorithm <ref> on (n, m, , 1/√(n)).
Then, for n ≥ 2 log8/δ, with probability at least 1 - δ,

    [ _f̃(X, 𝐘) ] 
    ≤ [ _f^*(X, 𝐘) ] 
           + 2 L _n() 
           + 2 B ( 2 √(log4δ) + m K )1/√(n) ,

where X is drawn from , and, conditionally on X, 𝐘 = (Y_1, …, Y_m) is drawn from (_X)^m.

[Proof Sketch]
First, we appeal to Algorithm <ref> to obtain n metasamples {(X_i, Ỹ_i,1, …, Ỹ_i,m)}_i=1^n with each X_i an i.i.d. draw from  and each Ỹ_i,j an independent draw from some _X_i' with X_i' - X_i_2 ≤1/√(n).
This holds except for O(√(n)) metasamples.

Now, from Assumption (<ref>) we have, for each X_i, m independent samples Ỹ_i,j from distributions that are close to _X_i.
The key idea is to show that each metasample i's labels can be viewed as coming from _X_i^m with high probability and from an arbitrary distribution otherwise.
This argument is first made for each Ỹ_i,j: We can view it as a sample from a mixture that puts high probability on _X_i and small probability on some other distribution.
Under this view, with high probability, every Ỹ_i,1, …, Ỹ_i,m comes from the _X_i component of its mixture.
Of course, this fails to be true for some of the metasamples, which we show again number only O(√(n)) with high probability.

The final component is an analysis of ERM with corrupted samples.
Consider (even in the classical setting) running ERM on a set of n samples, of which O(√(n)) have been corrupted arbitrarily but the rest are drawn i.i.d. from the underlying distribution.
In this case, we show that standard generalization bounds continue to hold with an error loss of only O(1/√(n)).


Recall that nm is the label complexity (pool-based active learning paradigm) and  is the sample complexity (supervised learning paradigm).[The algorithm is the same in both paradigms, except for the timing of label draws: in supervised, labels are drawn together with each X, while in pool-based active learning, labels are only queried when X is added to a metasample.]
To illustrate exactly how our results translate, let us adopt the parametric setting where the Rademacher complexity term _n() decays at the rate of O(1/√(n)) with n (meta)samples. Because the sample complexity is =  C m n^(d+1)/2 d^d/2logm(nd)^d/δ, the excess risk decays as the rate Õ(1/N^1/(d+1)), where the Õ notation omits log factors.
Similarly, in pool-based active learning, we can write n' = nm for the label complexity and get excess risk decaying at a rate O(√(1/n')).
(We fix m here as it is inherent to the loss function.)





 §.§ Faster rates under strong convexity and the uniform distribution


Let  be a class of linear predictors, so that it can be identified with a set ⊂^d. For w ∈, and fixed (x, 𝐲) ∈×^m, we assume that the loss has the generalized linear form
w ↦ c(⟨ w, ϕ(x) ⟩, 𝐲),
for some functions c and ϕ. 
The risk functional R is then

    R  w ↦[ c ( ⟨ w, ϕ(X) ⟩, 𝐘) ] ,

where X ∼ and, conditionally on X, 𝐘 = (Y_1, …, Y_m) ∼( _X )^m.

  
Assume that (<ref>) holds. 
Let  be a class of linear functionals as above, with the loss taking the generalized linear form. Suppose that ϕ(x)≤ B (the same B as for the upper bound on the loss). 
Let ε = (m K n)^-1 and =C m (n + (√(d)/ϵ)^d log3md^d/δϵ^d) for a universal constant C, and let f̃ be the hypothesis returned by Algorithm <ref> on (n,m,,ϵ). 
If  is the uniform distribution over [0, 1] and the risk functional R is σ-strongly convex, then, for any δ≤ 3 e^-4 and n ≥ 2 log8/δ, 
with probability at least 1 - δ

    [ _f̃(X, 𝐘) ] 
    ≤ [ _f^*(X, 𝐘) ] 
           + 3 B log3/δ/n 
           + 8 L^2 B^2 ( 32 + log3/δ)/σ n

where X is drawn from , and, conditionally on X, 𝐘 = (Y_1, …, Y_m) is drawn from (_X)^m.


This result implies that the excess risk decays at the rate O(1/n), where nm is the label complexity.  Since the sample complexity is = C m (n + (mKn√(d))^d log3m(mKnd)^d/δ), this implies that the excess risk decay rate in terms of  is Õ(1/N^1/d).

[Strongly convex risk functional] 
Take the example of the variance with ⊂^d. 
For fixed (x, y_1, y_2), the loss and risk functional are, respectively,

    w ↦( ⟨ w, x ⟩ - 1/2 (y_1 - y_2)^2 )^2 ,      R  w ↦_X ∼𝒟
     Y_1 ∼𝒟_X 
     Y_2 ∼𝒟_X[ ( ⟨ w, X ⟩ - 1/2 (Y_1 - Y_2)^2 )^2 ] .

Then ∇^2_w R(w) = 2  [ X X^T ], and so if [ X X^T ] ≽σ I, then (2 σ)-strong convexity holds.
In the special case of d = 1, provided that X is non-trivial we clearly have strong convexity of the risk.




 §.§ Sample complexity and dimension



While our sample complexity results above apply for any dimension d of , they scale exponentially in d.
In this work, we are motivated by d = O(1) and show in the next section that the multi-observation approach can yield significant practical improvements in these settings.
Here, we briefly note that an exponential dependence on d is information-theoretically necessary for any algorithm to regress on such higher-order statistics, even, in particular, in the simple case of estimating the average (y  x) over the distribution.
Therefore, we leave investigation of strategies for the higher-dimensional case (for example, active learning approaches) to future work.

For example, in  <ref> we show:
booldimlowercntr

theorembooldimlower

 If  is in the d-dimensional hypercube and the Lipschitz constant is K = 1, no algorithm for regression on variance of y can have nontrivial accuracy with o(2^0.5d) samples.

Intuitively, the obstacle is that a subexponential number of samples can (and, under e.g. a uniform distribution, does) have all x values separated from each other by a constant distance.
So any given region will only have one (x,y) pair sampled, and information-theoretically, no knowledge can be gleaned about the variance in that region.
 <ref> also provides a similar lower bound for the case of the uniform distribution over [0,1]^d, where K = d.




§ COMPARISON TO SINGLE-OBSERVATION LOSSES


For completeness, we now compare multi-observation regression to the typical case of single-observation regression, where the loss depends on a single prediction and a single observation.
We first explore this comparison theoretically, and then give experimental results, both in the simple setting of fitting a parametric model to the variance [Y|X].
In our experiments, we also take the opportunity to compare the various algorithms presented in  <ref> and explore other statistics of interest besides variance.



 §.§ Lower bounds for single-observation losses





Consider the example of the variance. 
One method for regressing the variance is the following “two-estimator” approach: use a single-observation estimator to regress the first conditional moment [ Y | X ] and a separate single-observation estimator to regress the second conditional moment [ Y^2 | X ]; the variance can then be predicted in terms of the two learned hypotheses.

Intuitively, this approach can go very wrong in many settings.
Although the conditional variance might have a simple parametric form, we may not know a priori if these conditional moments have simple forms. Thus, we may either try to estimate them from classes that are too small, picking up a large approximation error that leads to poor estimators of the variance, or we may conservatively estimate them from a sufficiently rich class to match our Lipschitz assumption (<ref>), which may lead to 
overfitting. 
Indeed, 
an existing minimax lower bound of <cit.> for estimating the conditional mean [ Y | X] when the latter is 1 / (2 √(d))-Lipschitz with respect to x
(as per Assumption (<ref>)) with range in [0, 1] implies that the risk of estimating the variance (under squared loss) of any two-estimator approach f̂,

    [ ( f̂(X) - [Y | X] )^2 ] ,

cannot decay at a rate faster than n^-2/(2 + d). 
Note that the estimation setting is only made easier by the fact that the Lipschitz constant decreases with d; this lower bound applies a fortiori for K-Lipschitz functions when K does not decrease with the dimension. Lastly, since the bound is a minimax lower bound, it holds for all estimators, not just those that predict according to K-Lipschitz functions.

We now describe how the lower bound of <cit.> implies the above result. 
As above, 
the goal is to minimize the L_2()-norm of the conditional variance of Y, i.e., to find an estimator f̂ for which the squared version thereof, (<ref>), 
is as small as possible.
In the two-estimator approach, f̂ takes the form f̂ = ĝ - ĥ for estimators ĝ and ĥ of [ Y^2 | X ] and (  [ Y | X] )^2 respectively. Then (<ref>) is equal to 

    [ ( ( ĝ(X) - [ Y^2 | X ] ) + ( [ Y | X ] - ĥ(x) ) )^2 ] .

The two-estimator approach needs to ensure that L_2() norm of ĝ(X) - [ Y^2 | X ] is small (and that the L_2() norm of the other term is small). Suppose that Y | X = x is Bernoulli for any x ∈. Then this latter necessary goal reduces to the familiar regression problem of minimizing

    [ ( ĝ(X) -  [ Y | X] )^2 ] .


The aforementioned minimax lower bound will apply to the above estimation problem in the following setting (we leave the full details to  <ref>):  take = [0, 1]^d and let  be the uniform distribution over . Suppose for all x ∈ that the distribution _x is a certain subclass of Bernoulli distributions
with x ↦ [ Y | X = x ] a K-Lipschitz function. Then for any estimator ĝ, there exists a law Y | X satisfying the aforementioned assumptions such that

    [ ( ĝ(X) -  [ Y | X ] )^2 ] = Ω( n^-2 / (2 + d)) ,

and there is a matching upper bound, so that this is the minimax optimal rate of convergence. 

To see how this rate compares to our results, suppose that the variance is captured by a generalized linear form which lies within our model , so that f^*  x ↦ =  [ Y | X = x]. Then a slight generalization of a standard exercise shows that the excess risk [ _f̃(X, 𝐘) ] -  [ _f^*(X, 𝐘) ] in Theorem <ref> takes the form (<ref>) with f̂ replaced by our direct, metasample-based estimator f̃. In the case of dimension d=1, the rate of Theorem <ref> is n^-1 and hence better than the rate of n^-2/3 of the two-estimator approach. Note that there is no contradiction with the minimax lower bound of <cit.>, because we assume that the variance takes a parametric form; were the variance to be an arbitrary K-Lipschitz function, our rates would degrade.[We believe in this case that our rate for d=1 would also be n^-2/3, based on existing fast rates results for classes whose metric entropy grows as ε^-1 (the class of K-Lipschitz functions exhibits such growth).]

One might wonder if it is possible to beat the above lower bound by using active learning. Unfortunately, the same lower bound holds in a similar setting even if one is allowed to use active learning strategies <cit.>.

In our experiments, we explore the performance of the two-estimator approach using ERM, which can be significantly worse than a direct multi-observation regression of the variance.







 §.§ Experiments



In our experiments we opted for synthetic data over real data.  Multi-observation loss functions help in learning higher order statistics about Y | X x; unfortunately with real data one generally does not know what the true value of those statistics are and thus has no objective way of comparing different algorithms.  By using synthetic data we can choose the underlying value of these statistics and evaluate algorithms by how closely they approximate it.

We consider three different statistics: variance, the upper confidence bound (UCB), and MINVAR.
The variance can be elicited, as observed in prior work <cit.>, by the two-observation loss function (r,y_1,y_2) = (r-1/2(y_1-y_2)^2)^2.
The UCB, often used in surrogate optimization for robust design in engineering, is defined by _λ(Y) = [Y] + λ√([Y]), for some fixed λ.
We show in  <ref> that it can be elicited by a two-observation loss function under some restrictions on the distribution of Y.
MINVAR (ρ_k), the expected minimum of k+1 i.i.d. draws of the distribution, is straightforwardly elicitable via squared loss (r,y_1,…,y_k+1) = (r-min{y_1,…,y_k+1})^2.


In all of our experiments we are trying to learn a statistic of Y | X x, where X ∼Unif(0,1).  For the variance, we tried different distributions Y | X x of the form f(x) + N(0,1).  We present here our results when f(x) is either a sine wave or a line. (Y|X x) = 1 in all cases.  For the _λ(Y) experiments, we chose λ = 8 and Y | X x ∼Γ(k(x), θ(x)) where k(x) and θ(x) were chosen such that [Y|X x] = 2 + sin(4π x) and _λ(Y|X x) = x+10.  For the ρ_k experiments we let k = 4 and Y | X x be an exponential distribution with mean 5(x+2).  Consequently, ρ_4(Y | X x) = x+2.


As a baseline for the variance and _λ, we compared our algorithms to the strategy of learning functions for E[Y|X x] and E[Y^2|X x] via ERM and then combining those functions to estimate [Y|X x] or _λ(Y|X x).  For ρ_k there is no such obvious single observation ERM strategy.  For _λ and ρ_k we also compared our algorithms against the “empirical approach” that is given the power to draw multiple i.i.d. labels for each random x it draws, compute the empirical statistic for that x, then fit a line to the results.  See  <ref> for more details about our experiments.




Our results are depicted in Figure <ref>.  Observe that in all of our experiments, Algorithms <ref> and <ref> performed the best.  Fitting lines to the moments never performed well, which is not surprising as in all cases at least one moment was non-linear.[We found that reducing the magnitude of the “non-linear part” led the two-moment approach to outperform our algorithms.  This occurred for variance and _λ when |f(x)| and λ were sufficiently small, respectively.]  However, even when fitting quadratics to moments which are quadratics (in the case of the variance when f(x) = 2x-1) our two observation algorithms still outperformed the two moment approach.  This demonstrates that for the two observation approach to be beneficial it is only necessary that the statistic is simpler than the underlying distribution, not that the underlying distribution is from an entirely unknown class.  Our algorithms only show very slight improvement over the empirical approach for _λ and ρ_k, but even being competitive is surprising given that the empirical approach has more power than our algorithms are afforded.




§ CONCLUSION AND FUTURE WORK



ERM with multi-observation loss functions presents challenges and complications as compared to traditional ERM, but also interesting opportunities.
With initial theoretical guarantees in hand, one next step is to explore  “risk” or “variance” regression problems encountered in practice for which the multi-observation approach may be useful.
Active learning settings may be among the most fruitful, as they are well-suited to collecting multiple labels at or near the same feature.

This investigation will interact with elicitation and the design of loss functions that are consistent for a desired property of the conditional distribution, such as the variance.
In particular, an open problem is to discover a loss function for the UCB property that does not require restrictions on the distribution of Y.

Another direction is to investigate the higher-dimensional setting further.
Our lower bounds show that, without additional assumptions, the number of samples required for regression on e.g. the variance require sample complexity exponential in the dimension of .
One approach could be investigation into settings with low intrinsic dimension, embedded in some higher-dimensional space.
Another could be investigating active learning applications, which should allow sidestepping these lower bounds.













§ PROOFS FROM BACKGROUND



[Proof of Lemma <ref>]

Let  be the probability measure operator with respect to P, and let  be the empirical measure operator with respect to the sample (X_1, Y_1), …, (X_n, Y_n); that is, for each f ∈, we have _f =  [ _f(X, Y) ] and _f = 1/n∑_i=1^n _f(X_i, Y_i).

From the bounded differences inequality (Lemma 1.2 of <cit.>) with bounded differences c_i = B/n,

    ( sup_f ∈ ( - ) _f > [ sup_f ∈ ( - ) _f ] + t  ) ≤ e^-2 n t^2 / B^2 .


Next, let (X'_i, Y'_i)_i ∈ [n] be an independent copy of (X_i, Y_i)_i ∈ [n]. 
Jensen's inequality implies that

    [ sup_f ∈ ( - ) _f ] 
       = [ sup_f ∈1/n∑_i=1^n (  [ _f(X'_i,Y'_i) ] - _f(X_i,Y_i) ) ] 
       ≤[ sup_f ∈1/n∑_i=1^n ( _f(X'_i,Y'_i) - _f(X_i,Y_i) ) ] ,

which, for independent Rademacher random variables ϵ_1, …, ϵ_n is in turn equal to

    [ _ϵ_1, …, ϵ_n[ sup_f ∈1/n∑_i=1^n ϵ_i ( _f(X'_i,Y'_i) - _f(X_i,Y_i) ) ] ] 
       ≤ 2 [ _ϵ_1, …, ϵ_n[ sup_f ∈1/n∑_i=1^n ϵ_i _f(X_i,Y_i) ] ] 
       = 2 _n({_f  f ∈} ) .

The result follows from Theorem 7 of <cit.>, which shows that under the L-Lipschitzness of the loss as a function of the prediction ŷ = f(x), we have the comparison inequality _n({_f  f ∈} )  ≤ L _n().





§ PROOFS OF SAMPLING METHODS





Let d ∈ and let x^*_1, …, x^*_n, x_1, …, x_N be drawn independently from a distribution  on = [0,1]^d. If N ≥ m n^(d+3)/2 d^d/2log2 m n/δ, then with probability at least 1 - δ, there is a set 𝒥 of cardinality at least n - √((n log2/δ) / 2) for which, for each x^*_j with j ∈𝒥, there are at least m points x_i_j,1, …, x_i_j,m satisfying x^*_j - x_i_j,m_2 ≤1/√(n), and all the i_1,1, …, i_1,m, …, i_|𝒥|,1, …, i_|𝒥|,m∈ [N] are distinct.


Take C_1, …, C_r to be a partition of  for which every cell C_j has diameter sup_x,x' ∈ C_jx - x'_2 at most ε. Observe that we may always take r ≤(, ε/2), where (, ε) is the minimum number of radius-ε balls in the Euclidean norm ·_2 whose union contains .  Note that r ≤ (√(d)/ϵ)^d.

We take ε = 1/√(n) 
and partition the set of cells of  into light cells and heavy cells, where any light cell C_j satisfies (C_j) ≤1/r√(n). Since there are at most r cells, the aggregate probability measure among all the light cells is at most 1/√(n). Hoeffding's inequality implies that only with probability at most δ/2 will more than 
√((n log2/δ) / 2) 
samples of x_1^*, …, x_n^* fall into light cells. 
The remainder of the points therefore fall into heavy cells, each of which has probability measure at least 1/r√(n). Now, for some fixed x_j^* in a heavy cell, if we sample r√(n)log1/δ points, then with probability at least 1 - δ at least one of these latter points would fall into the same cell as x_j^*. Thus, if we sample 

    N = m n^3/2r log2 m n/δ≤ m n^(d+3)/2 d^d/2log2 m n/δ
 
points x_1, …, x_N, then with probability at least 1 - δ/2 every x_j^* in a heavy cell will have at least m samples falling into its cell.


temptheorem


*


We will prove this for the case when m = 1, from which the general result follows by setting δ = δ / m and repeating m times.

We partition [0,1]^d into r=(√(d)/ϵ)^d hypercubes of width w=ϵ / √(d), C_1,…,C_r.  S_i = |{i | x_i^* ∈ C_i}| is the number of x_i^*'s that lie in the i hypercube.  Similarly, 
let T_i = |{i | x_i ∈ C_i}| be the number of x_i's that fall in the i hypercube.  We will now show that for all i ∈ [r], T_i ≥ S_i with probability at least 1-δ.  The lemma follows from this as we can then simply match up points within each hypercube and the maximum distance between two points in a hypercube is w√(d)=ϵ.

In order to bound [∀ i   T_i ≥ S_i], we will first consider a slight alteration of our setting.  Instead of having fixed numbers of samples n and N, they will be random variables ñ∼Pois(2n) and Ñ∼Pois(N/2).  In this setting, we let S̃_̃ĩ and T̃_̃ĩ be the number of samples of each category falling into the i interval.  The key property we will use here is that 

    S̃_̃ĩ∼Bin(ñ,w^d) = Pois(2nw^d),
 
and analogously T̃_̃ĩ∼Pois(Nw^d/2).

Using this we have

    [T̃_̃ĩ < S̃_̃ĩ]    = [Pois(N w^d/2) < Pois(2n w^d)]
       ≤ e^-(√(Nw^d/2)-√(2nw^d))^2

where the final inequality is a standard bound on Poisson races which follows from a Chernoff bound <cit.>.  By a union bound, the probability that T̃_̃ĩ≥S̃_̃ĩ for all i ∈ [r] is lower bounded by

    1 - r · e^-(√(Nw^d/2)-√(2nw^d))^2.


Now observe that

    [∀ i   T_i ≥ S_i]    = [∀ i  T̃_̃ĩ≥S̃_̃ĩ|Ñ = N ∧ñ = n] 
       ≥[∀ i  T̃_̃ĩ≥S̃_̃ĩ|Ñ≤ N ∧ñ≥ n]

because increasing N or decreasing n only increases our chances of finding a matching.  Let I be the event that Ñ≤ N ∧ñ≥ n, and let I̅ be its negation.

    [∀ i  T̃_̃ĩ≥S̃_̃ĩ| I]    = [∀ i  T̃_̃ĩ≥S̃_̃ĩ] - [∀ i  T̃_̃ĩ≥S̃_̃ĩ|I̅][I̅]/[I]
       ≥ 1 -  r · e^-(√(Nw^d/2)-√(2nw^d))^2 - [I̅]
       ≥ 1 -  r · e^-(√(Nw^d/2)-√(2nw^d))^2 - e^-Θ(n+N).

The final inequality follows from standard Poisson tail bounds <cit.>.  Plugging in N = C (n + (√(d)/ϵ)^d ( log1/δ+dlogd/ϵ)), for a sufficiently large constant C, gives us the lemma.




Let d∈ and let x^*_1, …, x^*_n, x_1, …, x_N be drawn independently from a distribution  on = [0,1]^d. If N ≥ C m n^(d+1)/2 d^d/2(logm/δ + d log(nd)) for a sufficiently large constant C, then with probability at least 1 - δ, there is a set 𝒥 of cardinality at least n - √((n log2/δ) / 2) for which, for each x^*_j with j ∈𝒥, there are at least m points x_i_j,1, …, x_i_j,m satisfying x^*_j - x_i_j,m_2 ≤1/√(n), and all the i_1,1, …, i_1,m, …, i_|𝒥|,1, …, i_|𝒥|,m∈ [N] are distinct.


This follows from a combination of the proofs of Lemma <ref> and Lemma <ref>.  As in Lemma <ref>, we partition [0,1]^d into heavy and light cells of diameter 1/√(n) and see that with probability δ / 2 at most √((n log2/δ) / 2) x_i^*'s fall in light cells.  We then apply the argument of Lemma <ref> to these heavy cells, using 1/r√(n)  (for r =  (nd)^d/2) as their probability mass instead of w^d.  Thus, it suffices that 

    N≥ C m n^(d+1)/2 d^d/2(logm/δ + d log(nd)) ≥ C' m (n + r√(n)logmr/δ)
 
for a sufficiently large constant C.







§ PROOFS OF EXCESS RISK BOUNDS





 §.§ Proof of Theorem  <ref>



[Proof of Theorem <ref>]

For each i ∈ [n] and j ∈ [m], draw Y_i,j independently according to distribution _X_i^*. This “clean” sample is simply a theoretical device for the analysis.

We first set up some convenient notation. 
For each i ∈ [n], define 𝐘_i := (Y_i,1, …, Y_i,m) and Ỹ_i := (Ỹ_i,1, …, Ỹ_i,m). 
Let  be a probability measure operator, defined according to _f =  [ _f(X,𝐘) ]; here, X is drawn from , and, conditionally on X, 𝐘 = (Y_1, …, Y_m) is drawn from (_X)^m. For a fixed f,  takes _f to its expected value on a new draw from the distribution X and an m-tuple 𝐘 from _X.
We also define the empirical probability measure operators  and  via 

    _f = 1/n∑_i=1^n _f(X^*_j, 𝐘_i) 
        and    _f = 1/n∑_i=1^n _f(X^*_i, Ỹ_i) .


Now according to Lemma <ref> below, we have for any positive t_1,t_2,t_3,

    ( _ > _f^* + t_1 + 2 t_2 + t_3 ) 
       ≤( sup_f ∈ ( - ) _f > t_1 ) 
              + ( sup_f ∈ |( - ) _f| > t_2 ) 
              + ( ( - ) _f^* > t_3 ) .

From Lemma <ref>, the first probability is at most δ / 4 when t_1 = 2 L _n() + B √(log(4/δ)/2 n). 
From Hoeffding's inequality, the third probability is at most δ / 4 when t_3 = B √(log(4/δ)/2 n) (note that f^* is fixed). 
The remainder of the proof controls the second probability. As we will see, we will be able to take t_2 = O ( B √(log(1/δ)/n)) when the probability is at most δ/2.

 
First, under Lemma <ref>, 
with probability at least 1 - δ/4, 
there is a subset ℐ_G ⊂ [n] of cardinality at least 
n_G := n - √((n log8/δ) / 2) 
for which, for each X^*_i with i ∈ℐ_G, 
there are at least m points X_k_i,1, …, X_k_i,m within distance ε of x^*_i, 
and all the k_1,1, …, k_1,m, …, k_n,1, …, k_n,m∈ [N] are distinct.

Next, we make the observation that the observed sample can be obtained by the following corruption modifications to (𝐘_i)_i ∈ [n].


  * For i ∈ [n] ∖ℐ_G, draw Ỹ_i,j from distribution _X_i,j.

  * For i ∈ℐ_G, observe that Assumption (<ref>) implies that, without loss of generality, we can view each Y_i,j as drawn in the following way. First, set Ỹ_i,j to Y_i,j. Next, draw a Bernoulli random variable Z_i,j with success probability τ := K ε, and if Z_i,j = 1, we corrupt Ỹ_i,j by setting it (again) to a new draw from some distribution Q_i,j that can depend on both X^*_i and X_i,j.


For each i, if Z_i,j = 0, we say that (i,j) are , and if (i,1), …, (i,m) all are , we say that i is . If some i is not , then it is . 
Clearly, for each i separately, with probability at least 1 - m K / √(n) over (Z_i,j)_j ∈ [m] it holds that i is  (recall that m τ = m K ε = m K / √(n)). 
Thus, from Hoeffding's inequality the probability (over (Z_i,j)_i ∈ [n], j ∈ [m]) that at least (C  +1) m K n_G/√(n) of the i's are  is at most e^-2 (n_G / n) (m K C)^2≤ e^-C^2 (recall that n ≥ 2 log8/δ) , so if C = √(log4/δ) then this probability is at most δ/4 (and our total probability of failure thus far is δ/2). We denote the (further diminished) good set of indices by ℐ'_G := {i ∈ℐ_G  i  is }; this set has cardinality at least n_G' := n_G ( 1 - √(log4/δ) + m K/√(n)) with probability at least 1 - δ/2.



From the above argument, we see that with probability at least 1 - δ/2, 
at most n'_B := n - n'_G corruption modifications occurred, and hence

    sup_f ∈|
        1/n∑_i=1^n _f(X^*_i, Ỹ_i)
        - 1/n∑_i=1^n _f(X^*_i, 𝐘_i) 
    | 
    ≤B (n - n'_G)/n .


Observe that

    n - n'_G 
       = n 
          - ( n - √((n log8/δ) / 2)) 
             ( 1 - √(log4/δ) + m K/√(n)) 
       = √(n)( √(log4/δ) + m K 
                                + √(log8/δ/2)) 
           - √(log8/δ/2)( √(log4/δ) + m K ) 
       ≤√(n)( 2 √(log4/δ) + m K ) ,

and thus we may take t_2 = B ( 2 √(log4/δ) + m K )/√(n).




 
  Under the hypotheses of Theorem <ref> we have the following, with probability taken over the random sample (i.e. , , and  are functions of the random sample):
  
    ( _ > _f^* + t_1 + 2 t_2 + t_3 ) 
       ≤( sup_f ∈ ( - ) _f > t_1 ) 
                + ( sup_f ∈ |( - ) _f| > t_2 ) 
                + ( ( - ) _f^* > t_3 ) .



First, observe that (using E for the 0-1 indicator function of a random event)

    ( sup_f ∈ ( - ) _f ≤ t_1 ) 
          ( sup_f ∈ |( - ) _f| ≤ t_2 ) 
          ( ( - ) _f^*≤ t_3 ) 
    
       ≤_≤_f^* + t_1 + 2 t_2 + t_3  .

To see this,

    _≤_ + t_1 
       ≤_f̂ +  (_ - _f̂) + t_1 
       (a)≤_f̂ + 2 t_2 + t_1 
       (b)≤_f^* + 2 t_2 + t_1 
       ≤_f^* + t_3 + 2 t_2 + t_1 ,

where (a) is from Lemma <ref> and (b) is from the optimality of ERM under .

By subtracting each side from one and rearranging, we get an implication on the negation of these events

    ( sup_f ∈ ( - ) _f > t_1 ) 
          ( sup_f ∈ |( - ) _f| > t_2 ) 
          ( ( - ) _f^* > t_3 ) 
    
       ≥_ > _f^* + t_1 + 2 t_2 + t_3

and we can use the union bound.



The following statement is true:

    _ - _f̂≤ 2 sup_f ∈| ( - ) _f | .



Observe that

    _ - _f̂   = ( _ - _) + ( _ - _f̂) 
       ≤( _ - _) + ( _f̂ - _f̂) 
       ≤ 2 sup_f ∈| ( - ) _f | ,

where the first inequality is from the optimality of  under .






 §.§ Proof of Theorem <ref>



[Proof of Theorem <ref>]
For each i ∈ [n] and j ∈ [m], draw Y_i,j independently according to distribution _X_i^*. This “clean” sample is simply a theoretical device for the analysis.

We first set up some convenient notation. 
For each i ∈ [n], define 𝐘_i := (Y_i,1, …, Y_i,m) and Ỹ_i := (Ỹ_i,1, …, Ỹ_i,m). 
Let  be a probability measure operator, defined according to _f =  [ _f(X,𝐘) ]; here, X is drawn from , and, conditionally on X, 𝐘 = (Y_1, …, Y_m) is drawn from (_X)^m. For a fixed f,  takes _f to its expected value on a new draw from the distribution X and an m-tuple 𝐘 from _X.
We also define the empirical probability measure operators  and  via 

    _f = 1/n∑_i=1^n _f(X^*_j, 𝐘_i) 
        and    _f = 1/n∑_i=1^n _f(X^*_i, Ỹ_i) .


Now according to Lemma <ref> below, we have for any positive t_1, t_2,

    ( _ > _f^* + t_1 + 2 t_2 ) 
       ≤( sup_f ∈{ (_f - _f^*) -  (_f - _f̂) } > t_1 ) 
              + ( sup_f ∈ |( - ) _f| > t_2 ) .

Now, since the risk is σ-strongly convex, the first probability is at most δ / 3 from Theorem 1 of <cit.> with a = 1 and λ = 2 σ, yielding the choice t_1 = 8 L^2 B^2 ( 32 + log(3/δ) ) / (σ n).

The remainder of the proof controls the second probability. As we will see, we will be able to take t_2 = O ( B log(1/δ)/n) when the probability is at most 2 δ/3.

First recall that by using Algorithm <ref> with our choice of parameters, Lemma <ref> gives us that with probability at least 1 - δ/3, for all i ∈ [n], 
there are at least m points X_k_i,1, …, X_k_i,m within distance ε of X^*_i, 
and all the k_1,1, …, k_1,m, …, k_n,1, …, k_n,m∈ [N] are distinct.

Next, we make the observation that the observed sample can be obtained by the following “corruption” modifications to (𝐘_i)_i ∈ [n].


  * For i ∈ [n] ∖ℐ_G, draw Ỹ_i,j from distribution _X_i,j.

  * For i ∈ℐ_G, observe that Assumption (<ref>) implies that, without loss of generality, we can view each Y_i,j as drawn in the following way. First, set Ỹ_i,j to Y_i,j. Next, draw a Bernoulli random variable Z_i,j with success probability τ := m K ε, and if Z_i,j = 1, we “corrupt” Ỹ_i,j by setting it (again) to a new draw from some distribution Q_i,j that can depend on both X^*_i and X_i,j.


For each i, if Z_i,j = 1, we say that (i,j) are , and if (i,1), …, (i,m) all are , we say that i is . If some i is not , then it is . 
Clearly, for each i separately, with probability at least 1 - 1/n over (Z_i,j)_j ∈ [m] it holds that i is  (recall that m τ = m K ε = 1/n). 
Thus, from a multiplicative Chernoff bound[The bound being used is (S_n ≥ R) ≤ 2^-R for R ≥ 6  [ S_n ], from equation (4.3) of <cit.>, where S_n is the sum of i.i.d. Bernoulli random variables with success probability 1/n.] the probability (over (Z_i,j)_i ∈ [n], j ∈ [m]) that at least 3/2log3/δ of the i's are  is at most δ/3 (for δ < 3 e^-4) (and our total probability of failure thus far is 2 δ/3). We denote the good set of indices by ℐ'_G := {i ∈ [n]  i  is }; this set has cardinality at least n'_G = n - 3/2log3/δ with probability at least 1 - 2 δ/3.




From the above argument, we see that with probability 1 - 2 δ/3, 
at most n'_B := n - n'_G corruption modifications occur, and hence

    sup_f ∈|
        1/n∑_i=1^n _f(X^*_i, Ỹ_i)
        - 1/n∑_i=1^n _f(X^*_i, 𝐘_i) 
    | 
    ≤B (n - n'_G)/n .


Thus, we may take t_2 = 3/2 B log3/δ/n.



 
  Under the hypotheses of Theorem <ref> we have the following, with probability taken over the random sample (i.e. , , and  are functions of the random sample):
  
    ( _ > _f^* + t_1 + 2 t_2 ) 
       ≤( 
    sup_f ∈{ (_f - _f^*) -  (_f - _f̂) } > t_1 ) 
                + ( sup_f ∈ |( - ) _f| > t_2 ) .



First, observe that (using E for the 0-1 indicator function of a random event)

    ( sup_f ∈{ (_f - _f^*) -  (_f - _f̂) }≤ t_1 ) 
          ( sup_f ∈ |( - ) _f| ≤ t_2 ) 
    
       ≤_≤_f^* + t_1 + 2 t_2  .

To see this,

    _≤_f^* +  ( _ - _f̂ ) + t_1 
    ≤_f^* + 2 sup_f ∈| ( - ) _f | + t_1 ,

where the second inequality is from Lemma <ref>.

By subtracting each side from one and rearranging, we get an implication on the negation of these events

    ( sup_f ∈{ (_f - _f^*) -  (_f - _f̂) } > t_1 ) 
          ( sup_f ∈ |( - ) _f| > t_2 ) 
    
       ≥_ > _f^* + t_1 + 2 t_2 + t_3

and we can use the union bound.








§ LOWER BOUNDS FOR HIGH DIMENSIONS
 

temptheorem2


*


  We consider the simplest nontrivial hypothesis class, constant functions (i.e. the set {f_c: c ∈} where each f_c: x ↦ c   (∀ x).
  The instances we construct will be realizable, i.e. in each instance, there will exist a constant c such that (y | x) = c   (∀ x).

  Consider the discrete uniform distribution on the boolean hypercube = {0,1}^d.
  We have = {0,1}.
  In instance A, [y=1  x] = 0.5 independently for all x, and (y  x) = 0.25 for all x.
  In the family of instances , we construct an instance B by drawing, for each x, θ_x ∈{0,1} uniformly and independently at random; then conditioned on x, the distribution of y is given by [y = 1  x] = θ_x.
  Notice that for all instances in the family , (y  x) = 0 for every x.
  Also, the Lipschitz constant satisfied by instances in  is 1, as all distinct x lie at a distance at least 1 from each other.

  Informal sketch.
  Any algorithm that is accurate must, with probability close to 1, produce a different output when given access to A than when given access to a uniform instance from B.
  However, by the principle of deferred decisions, we can rewrite the algorithm's behavior in the latter case as follows: Each time a uniformly random x is drawn, if the algorithm has already seen a sample from x, then set y consistent with that previous sample; otherwise, draw θ_x uniformly at random and set y = θ_x.

  Thus, the input to the algorithm is distributed exactly identically in both cases unless the algorithm obtains multiple samples from the same x.
  However, with o(2^0.5d) samples, the probability of this occurring is o(1) (by the “birthday paradox”), so the algorithm has the same distribution of outputs with probability 1 - o(1).

  Formal proof.
  Let M be an algorithm and write M(A) for the random variable which is M's hypothesis when run on samples from A, while M(B) is M's hypothesis when run on samples from a uniformly randomly chosen instance B from the family .
  Suppose that M satisfies that, with probability at least 2/3, its hypothesis (which is some constant c) is within ϵ of the correct variance, i.e. M(A) ≥ 0.25-ϵ and M(B) ≤ϵ each with probability at least 2/3.[One can also state this condition as an ≈ϵ generalization error guarantee for M with appropriate loss function.]
  Suppose ϵ < 0.125 and the number of samples drawn by M is o(2^0.5d); we show a contradiction.

  Use s to denote a set of samples each of the form (x,y), and let NR denote those sample-sets which have “no repeated x's”, i.e. NR = {s : each x in s is unique}.
  Use _A[s] to denote the probability of drawing a set of samples s given access to A, with _B[s] the probability of drawing s given access to a uniformly random instance from B, and so on.
  We have
  
    2/3   ≤[ M(A) ≥ 0.25-ϵ ] 
       =    ∑_s _A[s] [ M(s) ≥ 0.25-ϵ] 
       =    ∑_s ∈ NR_A[s] [ M(s) ≥ 0.25-ϵ ]    +  ∑_s ∉NR_A[s] [ M(s) ≥ 0.15]  
       ≤∑_s ∈ NR_A[s] [ M(s) ≥ 0.25-ϵ ]    +  ∑_s ∉NR_A[s] .

  Now, for all s ∈ NR, we claim _A[s] = _B[s], as each is equal to
  
    ∏_(x,y) ∈ s[x] [y  x]    = ∏_(x,y) ∈ s 2^-d(1/2).

  (In the case of A, this is immediate; in the case of B, by the principle of deferred decisions, we can construct B piece-by-piece; as each sample (x,y) ∈ s ∩ NR is drawn, we draw θ_x ∈{0,1} uniformly at random and set [y = 1] = θ_x, which results in a uniform distribution on y.)

  Meanwhile, ∑_s ∉NR_A[s] is the probability of drawing a sample with some repeated x value, which we claim is o(1) with o(2^0.5d) samples.
  The distribution is uniform on the 2^d possible x values.
  The probability of a repeat or “collision”, by Markov's inequality, is at most the expected number of collisions; with m samples, there are m  2 pairs each with a 2^-d chance of collision, so the expected number of collisions is O(m^2/2^d), which is o(1) for m = o(2^0.5d).

  So with o(2^0.5d) samples, we have by (<ref>) that
  
    2/3   ≤∑_s ∈ NR_B[s] [ M(s) ≥ 0.25-ϵ ]   +   o(1) 
       ≤[ M(B) ≥ 0.25-ϵ ] + o(1)

  which, if ϵ < 0.125, implies that [ M(B) ≤ϵ ] ≤1/3 + o(1).
  This contradicts the accuracy assumption that [ M(B) ≤ϵ] ≥2/3, so with this small number of samples, no such accurate M exists.

  
 
  With a uniform distribution on the unit hypercube [0,1]^d in d dimensions, with Lipschitz constant K = d, there is no algorithm for regression on variance with nontrivial accuracy drawing o(2^0.5d) samples.

  More precisely, estimating average variance over the hypercube to accuracy ϵ < 1/32 with success rate at least 2/3 requires Ω(2^0.5d) samples.


  The construction is very similar to the Boolean hypercube above.
  We have = {0,1}.
  On instance A, for every x, [ y = 1  x] = 0.5.
  Hence, (y  x) = 0.25 for all x.

  

  
Constructing .
  We now construct the family of instances  and show that each has Lipschitz constant K ≤ d and has average variance _x (y  x) ≤3/16.

  In each instance, the hypercube is divided into “corners” and “interior regions”.
  Let β = 1/21/2d and let α = 1/2 - β.
  Each “corner” C_v is a hypercube of side length α, inscribed in the unit hypercube and sharing the vertex v ∈{0,1}^d.
In other words, C_v = {x : x - v_∞≤α}.
  The portions of [0,1]^d not contained in any corner are considered the “interior regions”.

  To construct an instance in the family , draw θ_v ∈{0,1} i.i.d. uniformly for each v.
  For points x in some corner C_v, we have [y = 1  x ] = θ_v.
  For points x in the interior regions, let the notation x - C_v_2 denote min_x' ∈ C_vx' - x_2.
  If there exists a v such that x-C_v_2 < β, then, letting r = x-C_v_2, set [y = 1  x] = (1 - r/β)θ_v + (r/β)(1/2).
  (Note that this can only be true for at most one v, as this implies x - v_∞ < 1/2, i.e. x is contained within the “corner” of side length 1/2 touching v.)
  For all other x (those not within β of any C_v), we have [ y = 1  x] = 1/2.

  Now we show that any instance in family B has Lipschitz constant K = d.
  For shorthand, write p_x := [y = 1  x].
  Note that total variation distance between the distributions on y at x and at x' is |p_x - p_x'|.
  The Lipschitz constant is bounded by the maximal directional derivative of p_x with respect to x in any direction.
  This is zero if x lies within some C_v or if x is not within distance β of some C_v.
  Otherwise, if x' = _x”∈ C_vx” - x_2, then the absolute value of the directional derivative is maximized in the direction x' - x, where it is 1/2β.
  This gives a Lipschitz constant of 1/2β.

  Now we bound the average variance of y given x.
  The volume of each corner C_v is α^d and there are 2^d corners, so the total volume of the corners is
  
    2^d α^d
             = 2^d (1/2 - 1/21/d)^d  
       = (1 - 1/d)^d  
       ≥(1 - 1/2)^2  
       = 1/4

  assuming d ≥ 2.
  For any x ∈ C_v for any v, (y  x) = 0.
  For all other x, (y  x) ≤ 0.25; and the volume computation shows that they make up at most 3/4 of the hypercube.
  Hence, average variance in this instance is at most 0.25 (3/4) = 3/16.

  (We note that, by letting β = 1/21/C d for C ≥ 1, the same derivation gives Lipschitz constant K = C/d and an average variance bounded by 4C-1/16C^2≤1/4C.)

  

  
Indistinguishability.
  From here, the proof is almost identical to the Boolean case.
  Consider any algorithm M that is with ϵ < 1/32 of the correct average variance with probability at least 2/3.
  For each vertex v, let R_v = {x : x - v_∞ < 1/2}.
  There are 2^d disjoint regions R_v, each with volume 1/2^d, and with probability 1, each sample (x,y) has x ∈ R_v for some v.

  Let s denote a set of samples drawn by M and let NR = {s: s has no repeated xs}.
  If M draws o(2^d) samples, then [s ∈ NR] = o(1).
  If s ∉NR, then we claim _A[s] = _B[s], where the notation is shorthand for the probability of drawing the samples s given oracle access to A, or to a uniformly chosen member B of , respectively.
  The reason is that, by the principle of deferred decisions, we can in the case of B choose θ_v at the moment that a sample (x,y) is drawn with x ∈ R_v, which occurs at most once for each v because s ∉NR.
  Because the distribution on θ_v is uniform {0,1}, for any x ∈ R_v, the distribution of p_x = [y = 1  x] is uniform around 1/2, or in other words, the unconditional probability that this y = 1 is exactly 1/2.

  So,
  
    2/3   ≤[ M(A) ≥ 0.25 - ϵ ]  
       =    ∑_s ∈ NR_A[s] [M(s) ≥ 0.25 - ϵ]   +   o(1)  
       =    ∑_s ∈ NR_B[s] [M(s) ≥ 0.25 - ϵ]   +   o(1)  
       ≤[ M(B) ≥ 0.25 - ϵ ]   +   o(1) .

  For ϵ < 1/32, this contradicts the accuracy requirement that [M(B) ≤3/16 + ϵ] with probability 2/3.






§ SUPPORTING ARGUMENTS FOR LIPSCHITZ REGRESSION LOWER BOUND



In this section, we provide the full details linking the minimax lower bound of <cit.> to our application in  <ref>. All notation in this section is from <cit.>. 

Take = [0, 1]^d and let  be the uniform distribution over . For all x ∈, the distribution _x is a certain subclass of Bernoulli distributions[See Condition 2 of <cit.> for details; in our setting, we have t = θ(x) =  [ Y | X = x ], so the Bernoulli distribution can vary with x, as further explained on p. 1350 of <cit.>.] with x ↦ [ Y | X = x ] a K-Lipschitz function.
 Then for any estimator ĝ, there exists a law Y | X satisfying the aforementioned assumptions such that

    [ ( ĝ(X) -  [ Y | X ] )^2 ] = Ω( n^-2 / (2 + d)) ,

and there is a matching upper bound, so that this is the minimax optimal rate of convergence for this problem. 

To see how the above result follows from <cit.>, we take T(θ) = θ and observe that the Lipschitz condition is reflected in Stone's equation (1.2) by setting K_2 to our K, k to zero, and β = 1 (so that p = k + β = 1 and hence r = 1 / (2 + d)). Since θ is the parameter of a Bernoulli distribution, we need to verify from the lower bound construction that for all choices of θ∈Θ_n used in the lower bound, we have θ∈ (0, 1). We now do this verification.

As mentioned in the proof of Lemma 1 of <cit.>, for any binary sequence τ_n ∈{0,1}^V_n, the corresponding[We use the notation g^(τ)_n rather than Stone's notation g_n to make explicit the dependence on τ.] g^(τ)_n vanishes at the boundary of [0, 1]^d. Moreoever, by assumption each g_n^(τ) is K-Lipschitz for K = 1/(2 √(d)). Let us verify that, for all θ∈Θ_n and for all x ∈ [0,1]^d, it holds that θ(x) ∈ (0,1). Let θ_0 ≡1/2, and note that θ(x) = θ_0(x) + g^(τ)(x). It suffices to show that |g^(τ)_n(x)| ≤1/2 √(2) for all τ∈{0,1}^d and all x ∈ [0,1]^d. 
This is easily verified. First, observe that

    |g^(τ)_n(x)| 
    = | g^(τ)_n(x) - g^(τ)_n(0) | 
    = | g^(τ)_n(x) - g^(τ)_n(1) | ,

since g^(τ) vanishing at the boundary of [0,1]^d implies that g^(τ)_n(0) = g^(τ)_n(1) = 0. 
Therefore, 

    |g^(τ)_n(x)| 
       = min{| g^(τ)_n(x) - g^(τ)_n(1) |,
                             | g^(τ)_n(x) - g^(τ)_n(0) | }
       ≤ K min{x - 1,
                                     x - 0}
       = K 1/2·1
       = K √(d/2)
       = 1/2 √(2) .

The result then follows by applying Theorem 1 of <cit.> with q = 2.






§ ELICITING THE UPPER CONFIDENCE BOUND



Given a random variable Y, define _λ(Y) = [Y] + λσ[Y], where σ[Y] = √([Y]) = √([Y^2]-[Y]^2).
We show here the motivation behind the loss function used in  <ref>.

The level sets of _λ, as a function of the law of Y, i.e., the distribution, are given by

    r = _λ(Y) = [Y] + λσ[Y] .


We can rewrite this as follows:

    r - [Y]    = λσ[Y]
      
    
      (r - [Y])^2    = λ^2[Y]
    
      0    = λ^2[Y] - [Y]^2 + 2r[Y] - r^2 ,

though note that we have introduced another solution: both r = [Y]+λσ[Y] and r=[Y]-λσ[Y] now satisfy eq. (<ref>) but only the former satisfies eq. (<ref>).
Apart from this spurious solution, the following would be an identification function for _λ, meaning a distribution has zero expectation if and only if it is in the level set for r; see <cit.>.

    V(r,y_1,y_2) = λ^2/2 (y_1 - y_2)^2 - y_1y_2 + (y_1+y_2)r - r^2 ,

whence

    [V(r,Y_1,Y_2)] = λ^2 [Y] - [Y]^2 + 2[Y]r - r^2 ,

where of course Y_1,Y_2∼ Y are independent.

Despite the fact that V does not completely identify _λ, as for a given r, distributions where [Y]-σ[Y] = r also satisfy [V(r,Y_1,Y_2)]=0, we can still try to integrate -V with respect to r to get a loss function.  The result is

    (r,y_1,y_2) = -(λ^2/2 (y_1 - y_2)^2 - y_1y_2)r - 1/2 (y_1+y_2)r^2 + 1/3 r^3 .


As V was not a true identification function, we know that [(r,Y_1,Y_2)] has multiple extrema, at r=[Y]±λσ[Y].
What's worse, since  is cubic, we see that one can actually achieve arbitrarily negative loss as r → -∞.
Still, we can impose conditions on r and Y so that r = _λ(Y) is the unique minimizer of [].
In particular, if we restrict r to the range [Y]-2λσ[Y] ≤ r ≤∞, then the loss will elicit _λ.
If we further restrict r ≥[Y]-λσ[Y], then [] is quasi-convex.
Finally, if we further restrict r ≥[Y], then [] is convex.

Another possible loss is obtained by integrating -rV, that is, r times -V.  This gives us,

    (r,y_1,y_2) = -1/2(λ^2/2 (y_1 - y_2)^2 - y_1y_2) r^2 - 1/3 (y_1+y_2) r^3 + 1/4 r^4 .

This removes the problem of unbounded negative loss for incorrect reports, but can still have a local optimum at r = [Y] - λσ[Y].




§ SIMULATION DETAILS




  
Algorithm implementation. For Algorithms <ref> and <ref> we used ϵ = 1/2√(n), which we found to generally perform the best.[For ρ_k larger values of ϵ performed better for Algorithm <ref>, but caused it to be very computationally expensive.  This is because the number of metasamples is roughly (n ϵ)^m and for ρ_4, m=5.] There is some question about how to apply Algorithm <ref> in the setting where one is given a fixed set of samples, rather than an oracle for drawing samples.  The strategy that we found to work best, which we used here, was given n̂ samples, to use the first n for X_1^*,…,X_n^* and the remaining N = n̂ - n for X^(1)_1,…,X^(1)_N.  We then binary searched to find the largest n which allowed a perfect matching.[Note that this modification may introduce some bias so our theoretical results about Algorithm <ref> no longer directly apply.  Nevertheless, we found this modification to be effective in practice.]  In all of our experiments, our two-observation methods used the linear functions as their hypothesis class.  In all cases, the true statistic is within this class. 



  
Empirical approach. For _λ we compared against the standard strategy used in practice.  This strategy is to sample n / k random points x_1,…,x_n/k from X, and then for each x_i sample k values y_i,1,…,y_i,k∼ Y | X x_i.  For each i, we compute the empirical _λ u_i of y_i,1,…,y_i,k and then fit a line to (x_1,u_1),…,(x_n/k,u_n/k) via least-squares regression.  We found that best results were achieved by letting k = √(n).  We used an analogous strategy for ρ_k where we empirically estimated ρ_k using k+1 y samples for each x.  



  
Evaluation framework. As we used simple underlying statistics (1, x+10 and x+2) and simple hypothesis classes, we were able to compute the mean squared error between an algorithm's reported hypothesis and the true statistic via closed form expressions.  The only exception to this was for the two moment method of learning _λ, where the instead we estimated the mean squared error via 1000 sample Monte Carlo integration.  Each data point in Figure <ref> is the median of 1000 independent trials.



