Adaptive Sampling for Convex Regression

2018·Arxiv

Abstract

Abstract

In this paper, we introduce the first principled adaptive-sampling procedure for learning a convex function in the norm, a problem that arises often in the behavioral and social sciences. We present a function-specific measure of complexity and use it to prove that, for each convex function , our algorithm nearly attains the information-theoretically optimal, function-specific error rate. We also corroborate our theoretical contributions with numerical experiments, finding that our method substantially outperforms passive, uniform sampling for favorable synthetic and data-derived functions in low-noise settings with large sampling budgets. Our results also suggest an idealized “oracle strategy”, which we use to gauge the potential advance of any adaptive-sampling strategy over passive sampling, for any given convex function.

1 Introduction

Many functions that model individual economic utility, the output of manufacturing processes, and natural phenomena in the social sciences are either convex or concave. For example, convex functions are used to model utility functions that exhibit temporal discounting, a classic effect in behavioral economics where people value immediate rewards over delayed rewards [Frederick et al., 2002, Green and Myerson, 2004]. To measure such curves, it is common practice to manipulate a variable (e.g., delay) over a fixed, uniformly spaced grid of 5 points [Fisher, 1937], collect many repeated trials of data, and fit a function of assumed parametric form (e.g., exponential or hyperbolic) using maximum likelihood estimation. This approach can be brittle to model mismatch when the true function f lies outside the assumed class of functions. Moreover, non-linear parametric families can introduce challenges for constructing faithful and accurate confidence intervals when interpolating the estimator between measured design points.

Non-parametric convex regression (c.f. Dümbgen et al. [2004]) corrects for the shortcomings of parametric methods by making no assumptions other than that f is convex. In additional to faithfully modeling a large class of functions, non-parametric methods can also be employed to construct error bars at any (see Cai et al. [2013]). Unfortunately, even with shape restrictions, non-parametric methods may require prohibitively many samples for practical use.

In this paper, we propose a more parsimonious approach to non-parametric curve estimation by allowing the design points to be chosen sequentially and adaptively. Formally, we consider the problem of estimating an unknown convex function with an estimator which is close in the metric . The estimator is constructed from sequential, noisy evaluations from an oracle F at design points

and where represents zero-mean noise. We let denote the filtration generated by the design points and measurements up to time t, and assume that the number of samples is a stopping time with respect to , where is zero mean, -subgaussian. We refer to measurement allocation strategies for which does not depend on as passive, and adaptive sampling strategies for which may depend on

Our main contributions are the following:

• Inspired by Cai et al. [2013], we introduce the local approximation modulus, a new measure of local curvature for convex functions, , and a function-specific complexity measure , called the average approximation modulus. coin- cides with the average curvature of , up to logarithmic factors and endpoint considerations. We prove a function-specific lower bound on the sample complexity of actively estimating any convex function that scales at least as , up to logarithmic factors.

• The packing argument for constructing our lower bound explicitly describes a near-optimal, clairvoyant sampling allocation tailored to ; we call this the “oracle allocation” (Proposition 3.4). This allocation is instructive as an experimental benchmark when

• We introduce an active sampling procedure and an estimator whose sample complexity forany up to logarithmic factors, nearly matching our lower bounds.

• We show that for passive designs (e.g., sampled evenly on a grid), the sample complexity necessarily scales as , where coincides with the maximum curvature. We compare and for many natural classes of functions, including quadratic functions, exponential curves, and k-piecewise linear functions. For k-piecewise linear functions, the error of our active algorithm scales no slower than , whereas passive designs scale no faster than evaluations (see Remark 3.2).

Finally, we validate our theoretical claims with an empirical study using both synthetic functions and those derived from real data. We observe that in low-noise settings or when the sampling budget is large, active sampling can substantially outperform passive uniform sampling. Moreover, our algorithm constitutes the first theoretically justified algorithm (passive or active) that guarantees uniform accuracy, even at the boundaries of the interval [Cai et al., 2013, Dümbgen et al., 2004]. Even so, comparing the performance of our active algorithm to the oracle sampling strategy suggests room for modest but non-negligible improvements.

1.1 Related Work

Castro et al. [2005] and Korostelev [1999] studied the minimax rates of active non-parametric regression, showing that active and passive learning attain the same minimax rates of convergence for Holder smooth classes, but that active learning achieves faster rates when the function is known to be well approximated by a piecewise-constant function.

Prior literature on convex and concave regression consider the passive design case, where the design points do not depend on measurements. Typically, the design points are chosen to be uniformly spaced on the unit interval, that is, for [Dümbgen et al., 2004]. If F is the set of Lipschitz, convex functions, then the of the least squares estimator decreases like , whereas if the convex function has Lipschitz gradients, the rate improves to [Dümbgen et al., 2004].

Recent work by Guntuboyina and Sen [2015] and Chatterjee [2016] has aimed at developing sharp errors bounds on the squared -norm of the least squares estimator, when samples are uniformly spaced on a grid. They show that even with this uniform allocation, the error to the true regression functions . For example, Chatterjee [2016] and Bellec et al. [2018] show that if is a k-piecewise linear function, then obtains the parametric error rate of . In a similar vein, Cai et al. [2013] proves sharp confidence intervals for for a fixed point

Our work draws heavily upon Cai et al. [2013] (who in turn build on Dümbgen et al. [2003]), whose aim was to characterize the function-specific sample complexity of estimating a convex f at a given point in the interior of [0, 1], from uniform measurements. We extend these tools to characterize the complexity of estimating f with uniform accuracy over the interval [0, 1], from measurements which may be chosen in an adaptive, function-dependent manner. We are thus able to obtain exceptionally granular, instance-specific results similar to those in the multi-arm bandit literature [Kaufmann et al., 2016], and in recent work studying the local minimax sample complexity of convex optimization [Zhu et al., 2016].

2 Eﬃciently Learning a Convex Function

We begin by establishing preliminary notation. The class of convex functions is denoted as . For an interval , define the left-, middle- and right-endpoints as . We define the secant approximation of f on an interval

and note that for a convex function, this approximation never underestimates f; that is, one has . We denote the error of the second approximation to f on I at the midpoint as

In addition, we overload notation so that for any x, t such that , we have We now state a remarkable fact about convex functions that is at the core of our analysis.

Lemma 2.1 is a special case of a more general lemma stated in Section 6 that upper bounds the supremum of the secant approximation error by a constant using only a single point within the interval. Convexity is critical to the proof of this lemma and such a property does not hold, for instance, on merely monotonic functions. We remark that the first inequality is trivial, and the second inequality is tight in the sense that it is achieved by on interval I = [0, 1] as

The above observations motivate our strategy of approximating f with secant approximations on disjoint intervals whose union is [0, 1]. The next definition relates the secant approximation error to the required sampling density.

Definition 1 (Local Approximation Modulus) We define the -approximation modulus of f at a point as the least t such that the midpoint secant approximation to bias

Note that , because convex functions are continuous on their domain. Intuitively, is the scale at which f “looks” linear around some x, up to a tolerance from the endpoints {0, 1}, smaller values of correspond to larger complexities, because they imply that f can only be approximated by a linear function on a small interval. But if x is the near will take small values, potentially overestimating the local complexity. We remedy this issue by defining the following left- and right-approximation points:

Within , we will show that the midpoint errors concisely describe how densely one would need to sample f in the neighborhood of in order to estimate it to the desired accuracy in the -norm. Moreover, we show that it suffices to sample at constant number of design points on the end-intervals and . At a high level, the main finding of this paper is as follows:

The sample complexity of learning a particular convex function accuracy in passive sampling is parametrized by the worst-case approximation modulus

In contrast, the sample complexity of active sampling algorithms is parametrized by the average approximation modulus

We emphasize that the algorithm presented in in this work guarantees accuracy on the whole interval [0, 1], whereas many passive algorithms pointwise and risk bounds [Cai et al., 2013, Dümbgen et al., 2004] only guarantee accuracy on a strictly smaller sub-interval.

2.1 Examples

Explicit parameterizations of provide intuition for when active sampling is advantageous. In this section, we describe different scalings of for various ; later, in Remark 3.2, we explain how these scalings can imply substantial differences in sample complexity.

2.1.1 Piecewise Linear Functions

Let be a Lipschitz, piecewise linear convex function with a constant number of pieces. It follows that where is the distance to the closest knot adjoining any two linear pieces of . It follows that

2.1.2 Bounded third-derivative:

Suppose . We may apply a Taylor series to find as , which makes an explicit connection between the curvature of the function and the differences between . This suggests that if the function has areas of high but localized curvature such as or then the difference between can be as vast as

2.1.3 Quadratic Functions:

Let for some real coefficients a, b, c. Ignoring the effect of endpoints, for all x due to the function having constant curvature, so

3 Main Results

In this section, we state a formal upper bound obtained by Algorithm 2, described in Section 4. Algorithm 2 takes in a confidence parameter , as well as a second parameter the degree to which the active sampling algorithm is ‘aggressive’; from simulations, we recommend setting . Lastly, at each round, Algorithm 2 maintains an estimator , whose performance is characterized by the following theorem:

Theorem 3.1 Let C > 0 be a universal constant, and for

Then, if Algorithm 2 is run with parameters and , with access to an oracle (1), the estimators and confidence estimates defined in Section 4.2 satisfy the following any-time guarantee:

In the case of the default parameter setting , we find that, for a possibly larger universal constant

Up to constants and logarithmic factors, the sample complexity is dominated by the term corresponds to the standard rate for estimating a scalar. The dependence on captures the number of points required to estimate f with a discretized proxy.

To better understand why is the appropriate quantity to consider, we now introduce a construction of local packings of , centered at a given . Recall that denotes the error of the secant approximation to f on the midpoint of I, constructed using the endpoints . We note that if any algorithm, even an active one, does not measure f on an interval , then one cannot distinguish between f and the alternative function . Thus, a key step to showing that approximately lower bounds the number of evaluations is to show that it approximately lower bounds the number of intervals . This is achieved in the following theorem proved in Section 5.3:

be a convex function, , and define

where is defined analogously. Then, there is an such that the points such that the intervals

have disjoint interiors, are contained in and satisfy . Moreover, the interval endpoints overlap so that

Note that corresponds to the actual size of the explict packing, and is a computable lower bound on . We now consider the class of alternative functions

We observe that , and by definition, if are distinct, then In particular, given any set of points , then there exist two convex functions in , such that for all i and . In Section 5.2, we formalize this argument to yield the following theorem:

be as in Lemma 3.2, and let and be as given by Equation (7). Let Alg be any active algorithm that returns an estimator a stopping time , and satisfies the correctness guarantee

Then the stopping time , under observations from f, is lower bounded by

and the average sample complexity over is at least

The above bounds hold when is replaced by

Remark 3.1 The additional logarithmic factor that arises in (9) is due to the fact that estimating a function corresponds to correctly performing simultaneous hypothesis tests, regarding the value of g on each of the intervals . However, for any fixed (and, in particular, f = g), one can devise an algorithm that does not suffer this logarithmic factor by ‘biasing’ the algorithm towards that function.

In addition to providing a lower bound, the packing of Theorem 3.2 defines a near-optimal covering as well, in the sense that it defines a sampling allocation that can be used to test the hypothesis that, for a given . Formally, we have the following:

Proposition 3.4 For every function , there exists a

a deterministic sampling allocation , and a test function constructed from the allocation

The design is explicitly constructed in Section 5.3 by augmenting the -intervals in Theorem 3.2 with at most three additional intervals to ensure coverage of all of [0, 1]. Crucially, we made use of the fact that intervals share endpoints, and have secant error exactly equal to . In light of Theorem 3.3, we see that the design is optimal for verifying that to scaling by constant factors. For this reason, we refer to this construction as the oracle allocation since it precisely characterizes the optimal sampling allocation taken if one . In general, this allocation may be too optimistic, since an algorithm which does not know the true cannot choose this allocation a fortiori.

3.1 Comparison between Upper and Lower Bounds

For the purpose of comparing upper and lower bounds, we will consider running Algorithm 2 with the setting ; any constant bounded away from zero will yield qualitatively similar results. We find that the upper bound of Theorem 3.1 and lower bound of Theorem 3.3 nearly match, with the following exceptions:

1. The upper bound involves a doubly logarithmic factor that depends on . This is a consequence of the law of the iterated logarithm, which Algorithm 2 uses to maintain uniform correctness of its confidence intervals over time.

2. Theorem 3.1 is given in terms of , whereas our lower bound is stated in terms of . The two quantities can be related by the following proposition, proved in Section 6.4.

cω

Hence, ignoring the contributions of the endpoints , rescaling by a multiplicative constant by at most c.

3. Lastly, the upper and lower bounds differs in that requires dividing through by . We conjecture that the lower bound more accurately reflects the true sample complexity; see Remark A.1.

3.2 Sample Complexity for Passive Designs

In this section, we show that the sample complexity for estimating a convex function f with an approximately uniform passive design up to error is governed by the parameter

Theorem 3.6 Consider a (possibly randomized) passive design , which is uniform in the sense that, for some , and any interval I = [a, b] with , one has that . Then, for a universal constant c, any and all such that , there exists an alternative

The proof for the above theorem is as follows. Let

which intuitively corresponds to the point with the highest local curvature. Further, let , so that . If we consider the alternative function

then by construction, and f differ only on and . So if Alg can estimate f up to -norm error , then Alg can distinguish between f and . Consequently, standard information-theoretic arguments (Section 5.2) imply that any sampling algorithm must collect samples within . Theorem 3.6 then follows by the uniformity of the sampling procedure. In the case where the design is passive but not uniform, it is possible that the design performs well on particular functions . In Remark A.2, we show that nevertheless, if the design is not uniform, it will underperform on a ‘translation’ of f.

Remark 3.2 (Piecewise linear) If f is Lipschitz and piecewise linear with a constant number of pieces, then from Section 2.1 we have whereas . Theorem 3.6 implies that any passive sampling procedure requires measurements whereas Theorem 3.1 says that our active sampling procedure takes just . Thus, after n total samples the of passive sampling decays no faster than whereas active sampling decays like

4 Recursive Secant Approximation

We now introduce the recursive secant approximation algorithm for learning a convex function with noise. We begin by sampling each endpoint {0, 1} once. Subsequently, let t = 3, 4, . . . denote the number of samples taken, and let denote a binary tree of intervals contained in [0, 1], where

the children of an interval I are given by and . We let denote the set of leaves of . By construction, immediately satisfies the following properties stipulated in Lemma 4.1:

Lemma 4.1 For any , we have for any ; ; and

At each round t, we maintain three estimates of f. First, an estimator defined only at the points . Second, a secant-approximation estimator which extends the domain of

Note that is well defined, since by Lemma 4.1, for all , (a) there exists an such that and (b) if for , then x is a common endpoint of and , and thus the secant approximations coincide at Lastly, since is not guaranteed to be convex when measurements are noisy, we define an estimator projection onto

By definition so that by the triangle inequality. When not clear from context, we employ the use of a subscript t on to denote these functions once t samples have been taken.

4.1 Recursive Secant Approximation without Noise

To build intuition for Algorithm 2, we consider the following noiseless variant of our main algorithm, Algorithm 1, where the oracle returns noiseless queries F(x) = f(x). In this case, is set to be equal to f(x) at each point x that is queried, and a placeholder value of elsewhere. The algorithm maintains the invariant that, for all and have been measured and recorded in . Moreover, since the queries are noiseless, the secant approximation and no projection is required.

At each round, Algorithm 1 queries the interval for which the secant bias largest; note that if there is an interval has not been sampled, then will be queried, with ties broken arbitrarily. In preparation for the analysis of the noise-tolerant algorithm, we shall analyze the stopping time:

We shall prove the following proposition:

Proposition 4.2 For all . Moreover, for any we have

Proof Since for all queried points x, we have that, for . Moreover, since for any is constructed using secant approximations on a refinement of the intervals (see Lemma 6.1.)

It remains to bound . Let denote the set of points sampled at the start of round t; in the noiseless setting, , but bounding will be of broader interest for the noise-tolerant algorithm. Since are the endpoints of the intervals , which are adjacent, we have . Moreover, if denotes the parent-intervals of , we have . Thus, to bound , it suffices to bound . We adopt the shorthand

We now make a key observation about , which will allow us to relate to : for every , we have ; if not, then at the round at which I is bisected, we have , which implies that , a contradiction. The following lemma, proved in Section 4.4, shows that the inequality implies that the average modulus on each I cannot be too small.

Lemma 4.3 Let , and suppose that . Then for any ,

As a consequence, for any . To relate to the integral , we observe that the intervals are disjoint except at their endpoints, which yields

We remark that our bound on only used the fact that at time we had for each . This observation will be essential in generalizing to the setting with a noise oracle.

4.2 Recursive Secant Approximation with Noise

We now describe how to generalize Algorithm 1 to allow for noisy observations. Fix some time t and let be the collection of noisy function evaluation pairs. Recall that is independent, mean-zero -sub-Gaussian distributed noise, i.e. . In the algorithm, will denote the number of times the point been sampled so that if , and otherwise. Lastly, we let denote an anytime confidence interval such that

For example, suffices but we recommend using Kaufmann et al. [2016, Theorem 8]. In general can be chosen to be monotically decreasing in the t-argument, and increasing in the -argument. In addition to , we maintain a function such that

We shall let denote the event inside the probability operator in the above display. Finally, define confidence bounds

Crucially, our confidence bounds ensure the following sandwich relation, proved in Section 4.5:

Lemma 4.4 On , the following holds for all and

As a consequence, we find that

using the while loop of Line 10. This is to ensure that the stochastic variance always dominates the bias of the approximation. The parameter appears to have little effect on performance as long as it is smaller than 1; we recommend setting . The definition of in the algorithm is motivated by the sandwich relationship (Lemma 4.4) noted above. And in each case, is chosen in order to minimize the maximum confidence bound relevant to the interval . The values of

4.3 Proof of Upper Bound, Theorem 3.1

Recall the definition set , and we shall assume that holds. Fix an be as in Algorithm 2 Line 5, and define the stopping time

The correctness guarantee is a direct consequence of (14) since on . Because is only updated using a decreasing sequence of values of , the guarantee immediately holds for all . In order to upper bound , we have the identity

Thus, a crucial part of bounding is showing that we do not ; this is accomplished by relating the stopping condition to the sampling rule.

As a consequence,

where the second line uses the fact that is monotone in its second argument, and . We can upper bound the inversion of to yield (see e.g. Kaufmann et al. [2016])

To wrap up, it suffices to prove that for some

Recalling the argument from Section 4.1, it suffices only to verify that, if the secant approximation error of its parent is lower bounded by . We prove this as follows: fix some for . If is the parent of I then there exists some previous time s < t such that

that is, . On the other hand, to split on we must also have that

Together, these two displays imply . Rearranging, we find which proves what we set out to verify.

4.4 Proof of Lemma 4.3

Fix . This proof relies on the following upper-continuity property of , whose proof is deferred to Section 6.2:

, and suppose that

which proves one side of the integral. A similar argument holds for

4.5 Proof of Lemma 4.4

Define . First note that

using the fact that the secant approximations are affine on Lemma 2.1. Adding and subtracting

whence we conclude . For the lower bound, we see

so that

Remark 4.1 In the proof of Lemma 4.4 we lower bound quite loose since this quantity is also lower bounded by and can be at least a factor of two smaller. Nevertheless, using matching upper and lower bounds for substantially simplifies clutter in the algorithm. It is straightforward to modify the algorithm to use these non-matching upper and lower bounds for superior empirical performance, and, indeed, our experiments implement this modification.

4.6 Proof of Lemma 4.5

Fix an , and let be the last round at which was sampled; note then that . It suffices to bound . If is a new, just-bisected interval then we must have that by the sampling rule (was sampled only a single time.

Otherwise, has been sampled more than once and This means that for one has that

where the last line follows by the sampling rule. It suffices for the right-hand side to be less than to meet the stopping condition.

5 Proof of Packing and Lower Bounds

5.1 Proof of Theorem 3.2

We construct the packing by choosing a sequence of interval midpoints and interval lengths , such that the intervals overlap only at their endpoints, and such that . To do this, we define . By definition of , we have the equality , and for each

One can think of as as the equivalent of , but starting at rather than zero. Note that is non-decreasing and continuous in t (Lemma 6.3), and thus, if there exists a such that , then the supremum in the definition of will be attained for a . Thus, we will terminate the construction at i = n, where n is the first number satisfying , or . Note that . Collecting what we have established thus far,

1. . By Lemma 6.4, it follows that

2. By definition, . And, by the stopping condition,

3. Hence, since , we have that , and that have disjoint interiors.

To conclude, we adopt an argument similar to the proof of Proposition 4.2. For ease of notation, define . We start off by showing thatfor

one has

and thus the number of intervals satisfies

Finally, we remove the last interval . Since the right endpoint the intervals are contained within , and we have

Note then that we may take

5.1.1 Proof of Lemma 5.1

We first need a technical lemma, which we prove in Section 6.5.

We shall now establish; the bound on the integral over is analogous. We can write as , where and . Now, set . Using Lemma 5.2, we can integrate

where (a) is precisely Lemma 5.2. Lastly, we can bound , since

5.2 Proof of Noisy Lower Bound, Theorem 3.3

It suffices to prove the theorem with N replaced by , since ; the case where is addressed at the end of the section. Let Alg be any algorithm satisfying the correctness guarantee (8) for some denote the law under g and Alg. Consider the local alternative class and intervals , where are defined Theorem 3.2 and (7), respectively. Let denote the random variable corresponding to the number of times Alg samples in the interior of , and observe that since the intervals in i have disjoint interiors, the stopping time

We can reduce to a multiple hypothesis testing problem by recalling that, for , . Hence, for , the events are pairwise disjoint. Further, by (8), one has . We also recall Birge’s inequality:

denote a family of probability distributions on a space denote pairwise disjoint events. If

To apply Birge’s inequality, we first compute for any such that g(x) = h(x) for all denote the and , which is equal to

where (i) uses the fact that g and h differ only on uses the fact that, on , one of {g, h} is equal to , one is equal to , and thus by Lemma 2.1, we have that