Optimizing Star-Convex Functions

arXiv:1511.04466v2 [] 11 May 2016

Jasper C.H. Lee

Paul Valiant

Department of Computer Science
Brown University
{jasperchlee,paul_valiant}@brown.edu
May 13, 2016

Abstract
We introduce a polynomial time algorithm for optimizing the class of star-convex functions,
under no Lipschitz or other smoothness assumptions whatsoever, and no restrictions except exponential boundedness on a region about the origin, and Lebesgue measurability. The algorithm’s
performance is polynomial in the requested number of digits of accuracy and the dimension of the
search domain. This contrasts with the previous best known algorithm of Nesterov and Polyak
which has exponential dependence on the number of digits of accuracy, but only nω dependence
on the dimension n (where ω is the matrix multiplication exponent), and which further requires
Lipschitz second differentiability of the function [13]. Star-convex functions constitute a rich class
of functions generalizing convex functions, including, for example: for any convex (or star-convex)

1/p
p
p
is
functions f , g, with global minima f (0) = g(0) = 0, their power mean h(x) = f (x) +g(x)
2
star-convex, for any real p, defining powers via limits as appropriate. Star-convex functions arise
as loss functions in non-convex machine learning contexts, including, for data points Xi , parameter
P
1/p
p
vector θ, and any real exponent p, the loss function hθ,X (θ̂) =
|(
θ̂
−
θ)
·
X
|
, significantly
i
i
generalizing the well-studied convex case where p ≥ 1. Further, for any function g > 0 on the
surface of the unit sphere (including discontinuous, or pathological functions that havedifferent

x
behavior at rational vs. irrational angles), the star-convex function h(x) = ||x||2 · g ||x||
extends
2
the arbitrary behavior of g to the whole space.
Despite a long history of successful gradient-based optimization algorithms, star-convex optimization is a uniquely challenging regime because 1) gradients and/or subgradients often do not
exist; and 2) even in cases when gradients exist, there are star-convex functions for which gradients
provably provide no information about the location of the global optimum. We thus bypass the
usual approach of relying on gradient oracles and introduce a new randomized cutting plane algorithm that relies only on function evaluations. Our algorithm essentially looks for structure at all
scales, since, unlike with convex functions, star-convex functions do not necessarily display simpler
behavior on smaller length scales. Thus, while our cutting plane algorithm refines a feasible region
of exponentially decreasing volume by iteratively removing “cuts”, unlike for the standard convex
case, the structure to efficiently discover such cuts may not be found within the feasible region: our
novel star-convex cutting plane approach discovers cuts by sampling the function exponentially far
outside the feasible region.
We emphasize that the class of star-convex functions we consider is as unrestricted as possible:
the class of Lebesgue measurable star-convex functions has theoretical appeal, introducing to the
domain of polynomial-time algorithms a huge class with many interesting pathologies. We view our
results as a step forward in understanding the scope of optimization techniques beyond the garden
of convex optimization and local gradient-based methods.

1

Introduction

Optimization is one of the most influential ideas in computer science, central to many rapidly developing areas within computer science, and also one of the primary exports to other fields, including
operations research, economics and finance, bioinformatics, and many design problems in engineering.
Convex optimization, in particular, has produced many general and robust algorithmic frameworks
that have each become fundamental tools in many different areas: linear programming has become
a general modeling tool, its simple structure powering many algorithms and reductions; semidefinite
programming is an area whose scope is rapidly expanding, yielding many of the best known approximation algorithms for optimizing constraint satisfaction problems [2]; convex optimization generalizes
both of these and has introduced powerful optimization techniques including interior point and cutting
plane methods. Our developing understanding of optimization has also led to new algorithmic design
principles, which in turn leads to new insights into optimization. Recent progress in algorithmic graph
theory has benefited enormously from the optimization perspective, as many recent results on max
flow/min cut [10, 6, 8, 3, 15], bipartite matching [10], and Laplacian solvers [20, 3, 17] have made
breakthroughs that stem from developing a deeper understanding of the characteristics of convex optimization techniques in the context of graph theory. These successes motivate the quest for a deeper
and broader understanding of optimization techniques: 1) to what degree can convex optimization
techniques be extended to non-convex functions; 2) can we develop general new tools for tackling
non-convex optimization problems; and 3) can new techniques from non-convex optimization yield
new insights into convex optimization?
As a partial answer to the first question, gradient descent—perhaps the most natural optimization
approach—has had enormous success recently in a variety of practically-motivated non-convex settings,
sometimes with provable guarantees. The method and its variants are the de facto standard for training
(deep) neural networks, a hot topic in high dimensional non-convex optimization, with many recent
practical results (e.g. [7, 16]). Gradient descent can be thought of as a “greedy” algorithm, which
repeatedly chooses the most attractive direction from the local landscape. The efficacy of gradient
descent algorithms relies on local assumptions about the function: for example, if the first derivative
is Lipschitz (slowly varying), then one can take large downhill steps in the gradient direction without
worrying that the function will change to going uphill along this direction. Thus the convergence of
gradient descent algorithms typically depends on a Lipschitz parameter (or other smoothness measure),
and conveniently does not depend on the dimension of the search space. Many of these algorithms
converge to within  of a local optimum in time poly(1/), with additional polynomial dependence
on the Lipschitz constant or other appropriate smoothness guarantee [12, 5]. In cases where all local
optima are global optima, then one has global convergence guarantees.
While one intuitively expects gradient descent algorithms to always converge to a local minimum,
in this paper we study the optimization of a natural generalization of convex functions where, despite the global optimum being the only stationary point/local minimum for every function in this
class, provably no variant of gradient descent converges in polynomial time. The class of star-convex
functions, which we define below, includes many functions of both practical and theoretical interest,
both generalizing common families of convex functions to wider parameter regimes, and introducing
new “pathologies” not found in the convex case. We show, essentially, how to make a gradient-based
cutting plane algorithm “robust” to many new pathologies, including lack of gradients or subgradients,
long narrow ridges and rapid oscillation in directions transverse to the global minimum, and, in fact,
almost arbitrary discontinuities. This challenging setting gives new insights into what fundamentally
enables cutting plane algorithms to work, which we view as progress towards answering questions 2
and 3 above.

1

(a) g defined on the unit circle

(b) Linear extension of g to f

Figure 1: An example star-convex function f defined by linearly extrapolating an arbitrary positive
function g defined on the unit circle (in white).

1.1

Star-convex functions

This paper focuses on the optimization of star-convex functions, a particular class of (typically) nonconvex functions that includes convex functions as a special case. We define these functions as follows,
based on the definition in Nesterov and Polyak [13].
Definition 1.1 (Star-convex functions). A function f : Rn → R is star-convex if there is a global
minimum x∗ ∈ Rn such that for all α ∈ [0, 1] and x ∈ Rn ,
f (αx∗ + (1 − α)x) ≤ αf (x∗ ) + (1 − α)f (x)
We call x∗ the star center of f . For convenience, in the rest of the paper we shall refer to the minimum
function value as f ∗ .
See Appendix A for a variety of basic constructions and examples of star-convex functions.
Intuitively, if we visualize the objective function as a landscape, star-convexity means that the
global optimum is “visible” from every point—there are no “ridges” on the way to the global optimum, but there could be many ridges in transverse directions. Since the global optimum is always
visible in a downhill direction from every point, gradient descent methods would seem to be effective.
Counterintuitively, these methods all fail in general. In Section 1.3 we demonstrate that standard
gradient descent and cutting plane methods fail in the absence of Lipschitz guarantees.
One broad class of star-convex functions that gives a good sense of the scope of this definition is
constructed by the following process: 1) pick an arbitrary positive function g(θ) on the unit circle (see
Figure 1a) which may be discontinuous, rapidly oscillating, or otherwise badly behaved; and 2) linearly
extend this function to a function f on the entire plane, about the value f (0) = 0 (see Figure 1b).
The name “star-convex” comes from the fact that each sublevel set (the set of x for which f (x) < c,
for some c) is “star-shaped”.
At any point at which a gradient or subgradient exists, the star center (global optimum) lies in
the halfspace opposite the (sub)gradient. However, even for the simple example of Figure 1, gradients
do not exist for angles θ at which g is discontinuous, and subgradients do not exist at most angles,
including those angles where g is a local maxima with respect to θ. Further, even for differentiable
star-convex functions, gradients can be misleading, since a rapidly oscillating g implies that gradients
typically point in the transverse direction, nearly orthogonal to the direction of the star center.
While it is often standard to design algorithms assuming one has access to an oracle that returns
both function values and gradients, we instead only assume access to the function value: even in the
case when gradients exist, it is unclear whether they are algorithmically helpful in our setting. (See
Section 1.3 for details.) Indeed, we pose this as an open problem: under what assumptions (short of
2

the Lipschitz guarantees on the second derivative of Nesterov and Polyak [13]), can one meaningfully
use a gradient oracle to optimize star-convex functions?
In this paper, we show that, assuming only Lebesgue measurability and an exponential bound on
the function value within a large ball, our algorithm optimizes a star-convex function in polylog(1/)
time where  is the desired accuracy in function value.
Theorem 1.2. (Informal) Given evaluation oracle access to a Lebesgue measurable star-convex function f : Rn → R and an error parameter , with the guarantee that x∗ is within radius R of the origin,
Algorithm 1 returns an estimate f0 of the minimum function value such that |f0 − f ∗ | ≤  with high
probability in time poly(n, log 1 , log R).
In previous work, Nesterov and Polyak [13] introduced a new adaptation of Newton’s method to
find local minima in functions that are twice-differentiable, with Lipschitz continuity of their second
derivatives. For the class of Lipschitz-twice-differentiable star-convex functions, they show that their
algorithm converges to within  of the global optimum in time O( √1 ), using star convexity to lower
bound the amount of progress in each of their optimization steps.
Our results are stronger in two significant senses: 1) as opposed to assuming Lipschitz twicedifferentiability, we make no continuity assumptions whatsoever, optimizing over an essentially arbitrary measurable function within the class of star-convex functions; 2) while the algorithm of [13]
requires exponential time to estimate the optimum to k digits of accuracy, our algorithm is polynomial
in the requested number of digits of accuracy.
The reader may note that the complexity of our algorithm does depend polynomially on the number
of dimensions n of the search space, whereas the Nesterov-Polyak algorithm uses a number of calls to
a Hessian oracle that is independent of the dimension n. However, our dependency on the dimension
n is necessary in the absence of Lipschitz guarantees, even for convex optimization [11].
We further point out an important distinction between star-convex optimization and star-shaped
optimization. A star-shaped function is defined similarly to a star-convex function, except the star
center is allowed to be located away from the global minimum. Certain NP-hard optimization problems,
including max-clique can be rephrased in terms of star-shaped optimization [1], and the problem is
in general impossible without continuity guarantees, for the global minimum may be hidden along a
single ray from the star center that is discontinuous from the rest of the function. In this sense, our
star-convex optimization algorithm narrowly avoids solving an NP-hard optimization problem.

1.2

Applications of Star-Convex Functions

We begin with the following concrete example: the function (|x|p + |y|p )1/p is convex for p ≥ 1, yet is
star-convex for all real p, positive or negative, taking limits as appropriate for p = 0. We generalize
this example significantly below.
Empirical risk minimization is a central technique in machine learning [19], where, given training
data xi with labels yi , we seek a hypothesis h so as to minimize
m
1 X
L(h(xi ), yi ),
m
i=1

where L is a loss function that describes the penalty for misprediction, on input our prediction h(xi )
and the true answer yi . We take the hypothesis h to be a linear function (this setting includes kernel
methods that preprocess xi first and then apply a linear function). If the loss function L is convex, then
finding the optimal hypothesis h is a convex optimization problem, which is widely used in practice:
for exponent p ≥ 1, we let L(ŷi , yi ) = |ŷi − yi |p and thus aim to optimize the convex function
arg min
h

m
X

|h · xi − yi |p .

i=1

3

(1)

In this paper we draw attention to the interesting regime of Equation 1 where p < 1. This regime
is not accessible with standard techniques. It is straightforward to verify, however, that the (1/p)th
power of the objective function of Equation 1 is star-convex, and thus the main result of our paper
yields an efficient optimization algorithm. Slightly more generally:
Corollary 1.3. For a loss function L(ŷi , yi ) such that there exists a real number p for which L(·, ·)p is
a star-convex function of its first argument, the following empirical risk minimization problem can be
optimized to accuracy  in time poly(n, log 1/) by inputting its (1/p)th power to the main algorithm
of this paper:
m
1 X
arg min
L(h(xi ), yi ),
m
h
i=1

(We note that negative p is a somewhat peculiar regime in the above corollary, where running
a minimization algorithm on the (1/p)th power of the loss would actually yield a global maximum
instead of minimum; we include this case here for completeness’ sake.)
The regime where p ∈ (0, 1) has the structure that, when one is far from the right answer, small
changes to the hypothesis do not significantly affect performance, but the closer one gets to the true
parameters, the richer the landscape becomes. One of many interesting settings with this behavior
is biological evolution, where “fit” creatures have a rich landscape to traverse, while drastically unfit
creatures fail to survive. A paper by one of the authors proposed star-convex optimization as a
regime where evolutionary algorithms might be unexpectedly successful. This current work finds an
affirmative answer to an open question raised there by showing the first polynomial time algorithm
for optimizing this general class of functions. (The previous results by Nesterov and Polyak [13] fail
to apply to this
p setting because the objective function has regions that behave—for some parameter
regimes—like |x|, which is not Lipschitz.)
Corollary 1.4 (Extending Theorem 4.5 of [18]). There is a single mutation algorithm under which
for any p > 0, defining the loss function L(ŷ, y) = |ŷ − y|p , the class of constant-degree polynomials
from Rn → Rm with bounded coefficients is evolvable with respect to all distributions over the radius r
ball.

1.3

Need for Novel Algorithmic Techniques

Before we outline our approach to optimizing star-convex functions, we demonstrate the need for novel
algorithmic approaches by explaining how a wide variety of standard techniques fail on this
p challenging
p
class of functions. We start by explaining how simple star-convex functions such as ( |x| + |y|)2
confound gradient descent, cutting plane and related approaches, and end with a pathological example
with information-theoretic guarantees about its hardness to optimize, even given arbitrarily accurate
first-order (gradient) oracle access.
p
p 2
p
p Consider the star-convex function ( |x| + |y|) , which can be thought of essentially as |x| +
|y|. This function has unique global minimum (and star center) at the origin. Gradients of this
function go to infinity as either x or y goes to 0, thus the function essentially has deep canyons along
both the x and y axes. Intuitively, this function is challenging to optimize because, despite the origin
lying in a “downhill” direction from every point, gradient descent will typically converge to one of
the axes first, at which point the gradient has near-infinite magnitude in the direction of the nearest
axis, and this direction is thus near-orthogonal to the direction of the origin. In short, many variants
of gradient descent will have the following behavior: rapidly jump towards one of the two axes, and
then fail to make significant further progress as thep
search p
point oscillates around that axis. This is
in part a reflection of the fact that the gradient of ( |x| + |y|)2 is not Lipschitz, and in fact varies
arbitrarily rapidly as (x, y) converges towards either axis; since these axes are “canyons” in the search
4

space, typical algorithms will spend a disproportionate amount of their time stumbling around this
badly-behaved region.
There are many variants of gradient descent constructed to tackle optimization settings that are
challenging for “vanilla” gradient descent. One recent promising line of work concerns normalized
gradient descent [5], which is designed to work well despite regions where the magnitude of the gradient
is unusually small or large. However, the analysis techniques rely on the “strict local quasi-convexity”
property, which says that there is a global constant  such that the downhill halfspace induced by any
gradient contains a ball
p of radius
p 2 about the global optimum; this property does not apply even to
the simple function ( |x| + |y|) , where halfspaces cut arbitrarily close to the origin, the closer one
gets to either the x or y axis.
As an alternative to gradient descent methods, one might consider a gradient-based cutting plane
algorithm. Cutting plane methods work by iteratively cutting a large search space into ever-smaller
regions until they converge to a tiny region guaranteed to contain the global optimum. Cutting plane
algorithms, even for convex optimization, are known to typically lead to ill-conditioned behavior,
where one dimension of the search region becomes much smaller than the overall diameter of the
region, leading to bad numerical properties, difficulty finding and manipulating gradients,
etc.
p
p These
problems become even worse in our star-convex setting, even for the simple function ( |x|+ |y|)2 : as
the search space converges around one of the axes, the gradients become even more strongly oriented
towards that axis, yielding cuts that make the already thin region even thinner, instead of cutting in a
more productive transverse direction. Indeed, up to numerical precision, the computed gradients may
be found to point exactly towards the nearest axis, making cuts in the transverse direction impossible,
and thus halting the overall progress of the algorithm indefinitely.
Having seen simple examples showing how standard techniques cannot optimize star-convex functions, we now present a more sophisticated and pathological example which proves that, even with
access to a first-order oracle and the assumptions of infinite differentiability and boundedness on a
region around the origin, no deterministic algorithm can efficiently optimize star-convex functions
without Lipschitz guarantees, motivating the somewhat unusual randomized flavor of the algorithm
of this paper. In particular, we show that, for any deterministic polynomial time algorithm, there is
an infinitely differentiable star-convex function on R2 with a unique global minimum near the origin
such that for all queried points 1) the returned function value is always 1 and 2) the returned gradient depends only on the y-coordinate. The existence of such a function implies the inability of the
algorithm to find the x-coordinate of the star center.
Example 1.5. Given a deterministic polynomial time optimization algorithm A, we construct a starconvex function fA that algorithm A fails to optimize.
We assume that A will only query the function inside some disk of radius R about the origin, for
some (possibly exponentially large) R.
Now, we first choose a y-coordinate y0 (that A cannot specify exactly, for example an irrational
number) for the star center; we shall choose x0 later. Then, we simulate running the algorithm A,
assuming that in response to any query (x, y) that A makes to the function/gradient oracle, the return
value is fA (x, y) = 1 and the gradient ∇fA (x, y) is 0 in the x direction and 1/(y−y0 ) in the y direction.
Since A is a polynomial time algorithm, there must only be polynomially many queries before A
terminates. Therefore, we may choose x0 for the star center such that 1) (x0 , y0 ) is not exponentially
close to collinear with any pair of query points, and 2) the x coordinate of the star center, x0 , is chosen
to be significantly far from the x coordinates of all query points.
We thus construct a star-convex function fA subject to the constraints that 1) it passes through all
the points simulated above, with corresponding gradients, 2) it has star center at (x0 , y0 ), and 3) it
is infinitely differentiable, with each derivative being at most exponentially large in the degree of the
derivative and the running time of A. (Since no pair of query points are close to collinear with the
5

star center, standard interpolation lets us construct fA by first defining a smooth function gA on the
unit circle and linearly extending gA to the entire plane, except for an exponentially small region about
(x0 , y0 ) that we make smooth instead of “conically” shaped.)
Since the function and gradient oracles return zero information about the x-coordinate of the star
center, algorithm A clearly has a hopeless task optimizing this function.
In summary, gradient information is very hard to use effectively in star-convex optimization, even
for smooth functions. Our algorithm, outlined below, instead only queries the function value—not
because we object to gradients, but rather because they do not seem to help.

1.4

Our Approach

Our overall approach is via the ellipsoid method, which repeatedly refines an ellipsoidal region containing the star center (global optimum) by iteratively computing “cuts” that contain the star center,
while significantly reducing the overall volume of the ellipsoid.
As mentioned in Section 1.3, even for smooth functions with access to a gradient oracle, the cutting
planes induced by the gradients may yield no significant progress in some directions. Finding “useful”
cutting planes in the star-convex setting requires three novelties—see the introduction to Section 4 for
complete details. First, the star-convex function may be discontinuous, without even subgradients
defined at most points, so we instead rely on a sampling process involving the “blurred logarithm” of
the objective function. The blurred logarithm mitigates both the (potentially) exponential range of
the function, and the (potentially) arbitrary discontinuities in the domain. The blurred logarithm,
furthermore, is differentiable, and sampling results let us estimate and use its derivatives in our
algorithm. This technique is similar to that of randomized smoothing [4]. Second, the negative
gradient of the blurred logarithm might point away from the star center (global optimum)—despite
all gradients of the unblurred function (when they exist) pointing towards the star center—because
of the way blurring interacts with sharp wedge-shaped valleys. Addressing this requires an averaging
technique that repeatedly samples Gaussians-in-Gaussians, until it detects sufficient conditions to
conclude that the gradient may be used as a cut direction. Third and finally, the usual cutting plane
criterion—that the volume of the feasible region decreases exponentially—is no longer sufficient in the
star-convex setting, as an ellipsoid in two coordinates (x, y) might get repeatedly cut in the y direction
without any progress restricting the range of x. This is provably not a concern in convex optimization.
Our algorithm tackles this issue by “locking” axes of the ellipsoid smaller than some threshold τ , and
seeking cuts orthogonal to these axes; this orthogonal signal may be hidden by much larger derivatives
in other directions, and hence requires new techniques to expose. The counterintuitive approach is
that, in order to expose structure within an exponentially thin ellipsoid dimension we must search
exponentially far outside the ellipsoid in this dimension.

1.5

Stochastic Optimization

Our approach extends easily to the stochastic optimization setting. Here, the objective function is the
expectation over a distribution of star-convex functions which share the same star center and optimum
value. For example, in the context of empirical risk minimization with a star-convex loss function L
where we assume there is a true hypothesis h such that each labeled example is correctly predicted by
h, then for each example i, the loss L(ĥ(xi ), yi )—considered as a function of ĥ—is star-convex, with
global optimum 0 at ĥ = h. And the overall objective function is the expectation of this quantity, over
random examples: Ei [L(ĥ(xi ), yi )], which is thus also star-convex with global optimum 0 at ĥ = h.
The optimization algorithm for this new stochastic setting is actually identical to the general starconvex optimization algorithm. Since star-convex functions can take arbitrarily unrelated values on
nearby rays from the star center, our algorithm is already fully equipped to deal with the situation, and
6

adding explicitly unrelated values to the model by stochastically sampling from a family of functions
Li makes the problem no harder. The algorithm, analysis, and convergence are unchanged.

1.6

Outline

In Section 2 we introduce the oracle model by which we model access to star-convex functions. Because
we address such a large class of potentially badly-behaved functions (including, for example, starconvex functions which, in polar coordinates, behave very differently depending on whether the angle
is rational or irrational), we take great care to specify the input and output requirements of our
optimization algorithm appropriately, in the spirit of Lovász [9].
In Section 3 we present the overall star-convex optimization algorithm, a cutting-plane method
that at each stage tracks an ellipsoid that contains the star center (global optimum). The heart of the
algorithm consists of finding appropriate cutting planes at each step, which we present in Section 4.

2

Problem Statement and Overview

Our aim here is to discuss algorithms for optimizing star-convex functions in the most general setting
possible, and thus we must be careful about how the functions are specified to the algorithms. In
particular, we do not assume continuity, so functions with one behavior on the rational points and a
separate behavior on irrational points are possible. Therefore, it is essential that our algorithms have
access to the function values of f at inputs beyond the usual rational points expressible via standard
computer number representations. As motivation, see Example 2.1 which describes a star-convex
function that has value 1 on essentially all rational points an algorithm is likely to query, despite
having a rather different landscape on the irrational points, leading to a global minimum value of 0.
Example 2.1. We present the following example of a class of star-convex√functions√f in 2-dimensions
parameterized by integers (i, j) such that fi,j has global optimum at (1/ 2 + i, 1/ 3 + j). The class
has the property that, if i and j are chosen randomly from an exponential range, then, except with
exponentially small probability, any (probabilistic) polynomial time algorithm that accesses f only at
rational points will learn nothing about the location of the global optimum. It demonstrates the need
for access beyond the rationals.
We define a√star-convex
√ function fi,j parameterized by integers (i, j) that has a unique global
optimum at (1/ 2 + i, 1/ 3 + j). We evaluate f (x, y) via three cases:
√
√
1. If the ray from
√ (1/ 2 +√i, 1/ 3 + j) to (x, y) passes through no rational points, then let f (x, y) =
||(x, y) − (1/ 2 + i, 1/ 3 + j)||2 .
√
√
2. Otherwise, if the ray from (1/ 2 + i, 1/ 3 + j) to (x, y) passes through a rational point (p, q),
then:
√
√
(a) If bpc = i and bqc = j then let f (x, y) = ||(x, y) − (1/ 2 + i, 1/ 3 + j)||2 as above.
(b) Else, let

√
√
||(x, y) − (1/ 2 + i, 1/ 3 + j)||2
√
√
f (x, y) =
||(p, q) − (1/ 2 + i, 1/ 3 + j)||2

√
√
√
We √
note that since 1/ 2 and 1/ 3 are linearly independent over the rationals, no ray from (1/ 2 +
i, 1/ 3 + j) passes through more than one rational point. Hence the above three cases give a complete
definition of fi,j .
The function f is star-convex since each of the three cases above apply on the entirety of any ray
from the global optimum, and these cases each define a linear function on the ray. The derivative of
7

√
√
f along any of these rays is at most 1/(1 − 1/ 2), since (1 − 1/
√ 2) is the
√ closest that a point outside
the integer square at (i, j) can get to the global optimum (1/ 2 + i, 1/ 3 + j); thus the function is
linearly bounded in any finite region. Further, f evaluated at any rational point (p, q) will fall into
Case 2, and thus, unless (bpc, bqc) = (i, j), the function f will return 1.
Thus, no algorithm that queries f only at rational points can efficiently optimize this class of
functions, because unless the algorithm exactly guesses the pair (i, j)—which was drawn from an
exponentially large range—no information is gained (the value of f will always be 1 otherwise).
Since directly querying function values of general star-convex functions at rational points is so limiting, we instead introduce the notion of a weak sampling evaluation oracle for a star-convex function,
adapting the definition of a weak evaluation oracle by Lovász [9].
Definition 2.2. A weak sampling evaluation oracle for a function f : Rn → R takes as inputs a point
x ∈ Qn , a positive definite covariance matrix Σ, and an error parameter . The oracle first chooses a
random point y ← N (x, Σ), and returns a value r such that |f (y) − r| < .
The Gaussian sampling in Definition 2.2 can be equivalently changed to choosing a random point
in a ball of desired (small) radius, since any Gaussian distribution can be approximated to arbitrary
precision (in the total variation sense) as the convolution of itself with a small enough ball, and vice
versa.
Because inputs and outputs to the oracle must be expressible in a polynomial number of digits, we
consider “well-guaranteed” star-convex functions (in analogy with Lovász [9]), where the radius bound
R and function bound B below should be interpreted as huge numbers (with polynomial numbers of
digits), such as 10100 and 101000 respectively. Numerical accuracy issues in our analysis are analogous
to those for the standard ellipsoid algorithm and we do not discuss them further.
Definition 2.3. A weak sampling evaluation oracle for a function f is called well-guaranteed if the
oracle comes with two bounds, R and B, such that 1) the global minimum of f is within distance R of
the origin, and 2) within distance 10nR of the origin, |f (x)| ≤ B.
Such sampling gets around the obstacles of Examples 1.5 and 2.1 because even if the evaluation oracle is queried at a predictable rational point, the value it returns will represent the function evaluated
at an unpredictable, (typically) irrational point nearby.
At this point, our notion of oracle access may seem unnatural, in part because it allows access
to the function at irrational points that are not computationally expressible. However, there are two
natural and widely used justifications for this approach, of different flavors. First, the oracle represents
a computational model of a mathematical abstraction (the underlying star-convex function), where
even for mathematically pathological functions, we can actually in many cases implement simple code
to emulate oracle access in the manner described above. For example, it is in fact easy to implement
a weak sampling evaluation oracle for the function of Example 2.1, where rays from the star center
through rational points behave differently from the other rays. Since the rational points have Lebesgue
measure 0, the complexities introduced by case 2 of Example 2.1 happen with probability 0, and thus
we can write code to simulate a weak sampling
evaluation
oracle by only implementing the simpler
√
√
case 1, evaluating f (x, y) = ||(x, y) − (1/ 2 + i, 1/ 3 + j)||2 to the requested accuracy ±.
Second, setting oracle implementation issues aside, the model cleanly separates accessing a potentially pathological function, from computing properties of it, in this case, optimizing it. The more
pathological a function is, the more surprising it is that any efficient automated technique can extract
structure from it. Thus, in some sense, the unrealistic regime of pathological functions, for which we
do not even know how to write code emulating oracle access, yields the most surprising regime for
the success of the algorithmic results of this paper. This regime is the most insightful from a theory perspective almost exactly because it is the most unnatural and counterintuitive from a practical
perspective.
8

Having established weak sampling evaluation oracles as input to our algorithms, we must now take
the same care to define the form of their outputs. The output of a traditional (convex) optimization
problem is an x value such that f (x) is close to the global minimum. For star-convex functions,
accessed via a weak sampling evaluation oracle, one can instead ask for a small spherical region in
which f (x) remains close to its global optimum—given an input point x within distance τ of the star
4B
center, star convexity and boundedness of f imply that f (x) must be within value τ · 10nR−R−τ
of the
global optimum (see the proof of Lemma 3.2).
However, such an output requirement is too much to ask of a star-convex optimization algorithm:
see Example B.1 for an example of a star-convex function for which there is a large region in which the
distributions of f (x) in all small (Gaussian) balls are indistinguishable from each other except with
negligible probability, say 10−100 . That is, no algorithm can hope to distinguish a small ball around
the global optimum from a small ball anywhere else in this region, except by using close to 10100
function evaluations. Instead, since on this entire region, the function is close to the global optimum
except with 10−100 probability, the most natural option is for an optimization algorithm to return
any portion of this region, namely a region where the function value is within  of optimal except
with probability 10−100 . This is how we define the output of a star-convex optimization problem: in
analogy to the input convention, we ask optimization algorithms to return a Gaussian region, specified
by a mean and covariance, on which the function is near-optimal on all but a δ fraction of its points.
(See Example B.1 for details.)
Definition 2.4. The weak star-convex optimization problem for a Lebesgue measurable star-convex
function f , parameterized by δ, , F , is as follows. Given a well-guaranteed weak sampling evaluation
oracle for f , return with probability at least 1 − F a Gaussian G such that P r[f (x) ≤ f ∗ +  : x ←
G] ≥ 1 − δ.
We now formally state our main result.
Theorem 2.5. Algorithm 1 optimizes (in the sense of Definition 2.4) any Lebesgue measurable starconvex function f in time poly(n, 1/δ, log 1 , log F1 , log R, log B).
Observe that we require only Lebesgue measurability of our objective function, and we make no
further continuity or differentiability assumptions. The measurability is necessary to ensure that
probabilities and expectations regarding the function are well-defined, and is essentially the weakest
assumption possible for any probabilistic algorithm. The minimality in our assumptions contrasts that
of the work by Nesterov and Polyak, which assumes Lipschitz continuity in the second derivative [13].
It is mathematically interesting that the cardinality of the set of star-convex functions which we
optimize, even in the 2-dimensional case, equals the huge quantity 2|R| , while the cardinality of the
entire set of continuous functions—which is NP-hard to optimize, and strictly contains most standard
optimization settings—is only |R|.

3

Our Optimization Approach

In this section, we describe our general strategy of optimizing star-convex functions by an adaptation
of the ellipsoid algorithm. Our algorithm, like the standard ellipsoid algorithm, looks for the global
optimum inside an ellipsoid whose volume decreases by a fixed ratio for each in a series of iterations,
via algorithmically discovered cuts that remove portions of the ellipsoid that are discovered to lie in
the “uphill” direction from the center of the ellipsoid. In Section 4 we will explore the properties of
star-convex functions that will enable us to produce such cuts. However, it is crucial that, unlike for
standard convex optimization, such cuts are not enough. Consider, for example, the possibility that
for a star-convex function f (x, y), a cutting plane oracle only ever returns cuts in the x-direction, and
9

never reduces the size of the ellipsoid in the y direction, a situation which provably cannot occur in
the standard convex setting.
Thus, when constructing a cutting plane, our algorithm defines an exponentially small threshold
τ (defined in Definition 4.14), such that whenever a semi-principal axis of the ellipsoid is smaller
than τ , we guarantee a cut orthogonal to this axis. In this section, we make use of the following
characterization of our cutting plane algorithm, Algorithm 2, introduced and analyzed in Section 4:
Proposition (See Proposition 4.4). With negligible probability of failure, given an ellipsoid containing
the star center, centered at the origin with all semi-principal axes longer than τ scaled to 1, Algorithm
2 either a) returns a Gaussian region G such that P[f (x) ≤ f ∗ +  : x ← G] ≥ 1 − δ, or b) returns a
direction d⊥ , restricted to the subspace spanned by those ellipsoid semi-principal axes longer than τ ,
1
such that when normalized to a unit vector d̂⊥ , the cut {x : x · d̂⊥ ≤ 3n
} contains the global minimum.
Note that the above precondition on the ellipsoid input is satisfied by applying the appropriate
affine transformation to the current ellipsoid before supplying it to Algorithm 2. Also note that in
case a) above, the returned Gaussian region is a solution to the overall optimization problem.
Given these guarantees about the behavior of the cutting plane algorithm, Algorithm 2, we now
introduce our overall optimization algorithm, based on the ellipsoid method. h

1+

1

3n
For the following algorithm, we define the number of iterations m = 6(n+1) n log R
τ − (n − 1) log 2
according to a standard analysis of the multiplicative volume decrease of the ellipsoid method stated
in Lemma C.2 in Appendix C.

Algorithm 1 (Ellipsoid method)
Input: A radius R ball centered at the origin which is guaranteed to contain the global minimum.
Output: Either a) A Gaussian G such that P[f (x) ≤ f ∗ +  : x ← G] ≥ 1 − δ or b) An ellipsoid
E such that all function values in E are at most f ∗ + .
1. Let ellipsoid E1 be the input radius R ball.
2. For i ∈ [1, m + 1]
2a. If all the axes of Ei are shorter than τ , then Return Ei and Halt.
2b. Otherwise, execute Algorithm 2 with ellipsoid Ei .
2c. If it returns a Gaussian, then Return this Gaussian and Halt.
2d. Otherwise, it returns a cut direction d⊥ . Apply an affine transformation such that Ei
becomes the unit ball centered at the origin and d⊥ points in the negative x1 direction.
Use the construction in Lemma C.1 (see appendix) to construct a new ellipsoid Ei+1
that includes the intersection of ellipsoid Ei with the cut.
2e. If any semi-principal axis of the ellipsoid Ei+1 is larger than 3nR then apply Lemma C.3,
and if the center of the ellipsoid has distance > R from the origin then apply Lemma C.4,
to yield an ellipsoid of smaller volume, containing the entire intersection of Ei+1 with
the ball of radius R, that now has all semi-principal axes smaller than 3nR, and has
center in the ball of radius R.

The analysis of Algorithm 1 is a straightforward adaptation of standard techniques for analyzing
the ellipsoid algorithm, bounding the decrease in volume of the feasible ellipsoid at each step, until
either the algorithm returns an explicit Gaussian as the solution (in Step 2c), or terminates because
the ellipsoid is contained in an exponentially small ball (in Step 2a). See Appendix C for full details.
10

i

Lemma 3.1. Algorithm 1 halts within m + 1 iterations either through Step 2a or Step 2c.
The following lemma shows that, if the star center (global optimum) is found to lie within a ball of
exponentially small radius τ , then this entire ball has function value sufficiently close to the optimum
that any point within the ball can be returned.
Lemma 3.2. Given an ellipsoid E that a) has its center within R of the origin, b) is contained within
a ball of radius τ , and c) contains the star center x∗ , then f (x) ∈ [f ∗ , f ∗ + ] for all x ∈ E.
−−→
Proof. For any x ∈ E, let y be the intersection of the ray x∗ x with the sphere of radius 10nR about
the origin. The boundedness of f implies that f (y) ≤ B, and that f (x∗ ) ≥ −B. Thus star-convexity
∗ ||
2τ
2
∗
∗
yields that f (x) ≤ f ∗ + 2B ||x−x
||x∗ −y||2 ≤ f + 2B 10nR−R−τ < f +  since τ  .
Since the ellipsoid returned in Step 2a will always satisfy the preconditions of Lemma 3.2, the
entire ellipsoid has function value within  accuracy of the global minimum. In practice, we can just
return this region. However, if we do wish to conform to our problem statement (Definition 2.4), we
may simply return a Gaussian ball of sufficiently small radius and centered at the ellipsoid center,
such that there is only negligible probability of sampling a point more than τ from its center.
In summary, Lemmas 3.1 and 3.2 imply that Algorithm 1 optimizes a star-convex function in time
poly(n, 1/δ, log 1 , log F1 , log R, log B), proving Theorem 2.5.

4

Computing Cuts for Star-Convex Functions

The general goal of this section is to explain how to compute a single cut, in the context of our
adapted ellipsoid algorithm explained in Section 3. Namely, given (weak sampling, in the sense of
Definition 2.2) access to a star-convex function f , and a bounding ellipsoid in which we know the
global optimum lies, we present an algorithm (Algorithm 2) that will either a) return a cut passing
close to the center of the ellipsoid and containing the global optimum, or b) directly return the answer
to the overall optimization problem.
There are several obstacles to this kind of algorithm that we must overcome. Star-convex functions may be very discontinuous, without any gradients or even subgradients defined at most points.
Furthermore, arbitrarily small changes in the input value x may produce exponentially large changes
in the output of f , for any x that is not already exponentially close to the global optimum. This
means that standard techniques to even approximate the function’s shape do not remotely suffice in
the context of star-convex functions. Finally, considering again the example of linearly extending an
arbitrary function of the unit sphere surface (see Figure 1 or item 3 in Appendix A), unlike regular
convex functions, star-convex functions do not become “locally flat” along thin dimensions of increasingly thin ellipsoids. Thus, the ellipsoid algorithm in general might repeatedly cut certain dimensions
while neglecting others.
The algorithm in this section (Algorithm 2) employs three key strategies that we intuitively explain
now.
1. We do not work directly with the star-convex function f , but instead work with a blurred version
of its logarithm (defined in Definition 4.1), which has a) well-defined derivatives, that b) we can
estimate efficiently. Given a badly-behaved (measurable) function f , and the pdf of a Gaussian,
g, blurring f by the Gaussian produces the convolution f ∗ g, which has well-behaved derivatives
since derivatives commute with convolution: namely, letting D stand for a derivative or second
derivative, we have D(f ∗ g) = f ∗ (Dg). Since the derivatives of a Gaussian pdf g are each
bounded, we can estimate derivatives of blurred versions of discontinuous functions by sampling.
Sampling bounds imply that if we have a random variable bounded to a range of size r, then
11

we can estimate its mean to within accuracy ± by taking the average of O(( r )2 ) samples. The
logarithm function (after appropriate translation so that its inputs are positive), maps a starconvex function f to a much smaller range, which enables accurate sampling (in polynomial
time). Therefore, we can efficiently estimate derivatives of the blurred logarithm of f .
2. While intuitively, the negative gradient of (the blurred logarithm of) f points towards the global
minimum, this signal might be overpowered by a confounding term in a different direction
(see Lemma 4.5), making the negative gradient point away from the global optimum in some
cases. To combat this, our algorithm repeatedly estimates the gradient at points sampled from a
distribution around the ellipsoid center, and for each gradient, estimates the confounding terms,
returning the corresponding gradient only once the confounding terms are found to be small.
Lemma 4.8 shows that the confounding terms are small in expectation, so Markov’s inequality
shows that our strategy will rapidly succeed.
3. The ellipsoid algorithm in general might repeatedly cut certain dimensions while neglecting
others, and our algorithm must actively combat this. If some axes of the ellipsoid ever become
exponentially small, then we “lock” them, meaning that we demand a cut orthogonal to those
dimensions, thus maintaining a global lower bound on the length of any axis of the ellipsoid.
Combined with the standard ellipsoid algorithm guarantee that the volume of the ellipsoid
decreases exponentially, this implies that the diameter of the ellipsoid can be made exponentially
small in polynomial time, letting us conclude our optimization. Finding cuts under this new
“locking” guarantee, however, requires a new algorithmic technique.
We must take advantage of the fact that the ellipsoid is exponentially small along a certain
direction b to somehow gain the new ability to produce a cut orthogonal to b. Convex functions
become increasingly well-behaved on increasingly narrow regions, however, star-convex functions
crucially do not (see item 3 in Section A). Thus, as the thickness of the ellipsoid in direction b
gets smaller, the structural properties of the function inside our ellipsoid do not give us any new
powers. Strangely, we can take advantage of one new parameter regime that at first appears
useless: relative to the thinness of the ellipsoid, we have exponentially much space in direction
b outside of the ellipsoid.
To implement the intuition of Step 1 above, we define Lz (x), a truncated and translated logarithm of the star-convex function f , which maps the potentially exponentially large range of f to
the (polynomial-sized) range [log 0 , log 2B], where 0 is defined below (Definition 4.14), and is slightly
smaller than our function accuracy bound . In the below definition, z intuitively represents our estimate of the global optimum function value, and will record, essentially, the smallest function evaluation
seen so far (see Algorithm 2).
Definition 4.1. Given an objective function f with bound |f (x)| ≤ B when ||x|| ≤ R0 and an offset
value z ≥ f ∗ , we define the truncated logarithmic version of f to be

0

f (x) − z ≤ 0
log 
Lz (x) = log 2B
f (x) − z ≥ 2B


log(f (x) − z) otherwise
While mapping to a small range, Lz (x) nonetheless gives us a precise view of small changes in the
function as we converge to the optimum. The next result shows that, if we “blur” Lz (x) by drawing
x from a Gaussian distribution, then not only can we efficiently estimate the expected value of the
“blurred logarithm of f ”, we can also estimate the derivatives of this expectation with respect to
changing either the mean or the variance of the Gaussian.
12

For an arbitrary bounded (measurable) function h, the derivative of its expected value over a
Gaussian of width σ with respect to either a) moving the center of the Gaussian or b) changing its
width σ, is bounded by O( σ1 ). Thus we normalize the estimates below in terms of the product of the
Gaussian width and the derivative, instead of estimating the derivative alone.
Proposition 4.2. Let N (µ, Σ) be a Gaussian with diagonal covariance matrix Σ consisting of elements
σ12 , σ22 , . . . , σn2 . For an error bound κ > 0 and a probability of error δ > 0, we can estimate each of
the following functions to within error κ with probability at least 1 − δ using poly(n, κ1 , log 1δ , log 2B/0 )
samples: 1) the expectation E[Lz (x) : x ← N (µ, Σ)]; 2) the (scaled) derivative σ1 · dµd 1 E[Lz (x) : x ←
N (µ, Σ)]; and 3) the derivative with respect to scaling σ1 · dσd 1 E[Lz (x) : x ← N (µ, Σ)]. A fortiori, these
derivatives exist.
0

Proof. Chernoff bounds yield the first claim, since Lz ∈ [log 0 , log 2B], so O(( log 2B/
)2 log 1δ ) samples
κ
suffice.
For the second claim, we note that the expectation can be rewritten as the integral over Rn of
Lz times the probability density function of the normal distribution: E[Lz (x) : x ← N (µ, Σ)] =
R
T −1
(2π)−n/2 |Σ|−1/2 Rn Lz (x) · e−(x−µ) Σ (x−µ)/2 dx. Thus the derivative of this expression with respect
to the first coordinate µ1 can be expressed by taking the derivative of the probability density function
inside the integral, which ends up scaling it by the vector (x1 − µ1 )/σ12 . Thus σ1 · dµd 1 E[Lz (x) : x ←
1
1
· Lz (x) : x ← N (µ, Σ)]. The multiplier x1σ−µ
is effectively bounded, because
N (µ, Σ)] = E[ x1σ−µ
1
1
x1 −µ1
σ1
1
0)
,
log
2B/
κ

for real numbers c, the probability P[

> c : x ← N (µ, Σ)] vanishes faster than exp( √cn ).

such that replacing the expression in the expectation,
Thus we can pick c = poly(n, log
x1 −µ1
κ
σ1 · Lz (x), by this expression clamped between ±c will change the expectation by less than 2 . We
κ
can thus estimate the expectation of this clamped quantity to within error 2 with probability at least
1 − δ via some sample of size poly(n, κ1 , log 1δ , log 2B/0 ), as desired.
The analysis for the third claim (the derivative with respect to σ1 ) is analogous: the derivative of
the Gaussian probability density function with respect to σ1 (as in the case of a univariate Gaussian
of variance σ12 ) scales the probability density function by |x1 − µ1 |2 /σ13 , and after scaling by σ1 , this
expression measures the square of the number of standard deviations from the mean, and as above, the
2
1|
product |x1 −µ
Lz (x) can be clamped to some polynomially-bounded interval [−c, c] without changing
σ12
the expectation by more than κ2 . The conclusion is analogous to the previous case.
As mentioned in the previous section, our guiding aim for the ellipsoid method is to prevent any
axis of the ellipsoid from getting too small. Therefore, in order to treat axes differently depending on
their length, we shall identify our basis as the unit vectors along the axes of the current ellipsoid and
distinguish between axes that are smaller than τ versus at least τ , where τ is an exponentially small
threshold for “thinness”, defined below in Definition 4.14.
Definition 4.3. Given an ellipsoid, consider an orthonormal basis parallel to its axes. Each semiprincipal axis of the ellipsoid whose length is less than τ , we call a “thin dimension”, and the rest are
“non-thin dimensions”. Given a vector µ, we decompose it into µ = µ⊥ + µ> where µ⊥ is non-zero
only in the non-thin dimensions, and µ> is non-zero only in the thin dimensions. Similarly, given the
identity matrix I, we decompose it into I = I⊥ + I> .
We apply a scaling to the non-thin dimensions so as to scale the non-thin semi-principal axes of
the ellipsoid to unit vectors (making the ellipsoid a unit ball in the non-thin dimensions). We keep
the thin dimensions as they are.
We present again Proposition 4.4, describing the guarantees on our cutting plane algorithm required
by Algorithm 1, and then state the cutting plane algorithm, Algorithm 2. Below, we make use of
13

constants defined in Definition 4.14 that may be interpreted as follows: k is a polynomial number
of mesh points; η is the mesh spacing; τ 0 is the minimum size of σ> , a Gaussian width in the thin
0 is a
dimensions that is somewhat larger than τ , the size of the ellipsoid in the thin dimensions; σ⊥
Gaussian width in the ⊥ (non-thin) dimensions, of inverse polynomial size; σ⊥ is polynomially smaller
0 , and s is a polynomial quantity. S is a polynomial number of samples defined in the proof
than σ⊥
of Proposition 4.4. Moreover, the function g in the algorithm below is a sum of a probability and
two derivatives of the blurred expectation of Lz (x) (the truncated logarithm in Definition 4.1), and is
defined in the statement of Lemma 4.10.
Proposition 4.4 (Correctness of Algorithm 2). With negligible probability of failure, Algorithm 2
either a) returns a Gaussian region G such that P[f (x) ≤ f ∗ +  : x ← G] ≥ 1 − δ, or b) returns a
direction d⊥ , restricted to the ⊥ dimensions, such that when normalized to a unit vector d̂⊥ , the cut
1
{x : x · d̂⊥ ≤ 3n
} contains the global minimum.
Algorithm 2 (Single cut with locked dimensions)
Take an orthonormal basis for the ellipsoid, as in Definition 4.3. We apply an affine transformation
so that a) the ellipsoid is centered at the origin, and b) the ellipsoid, when restricted to the ⊥
dimensions, is the unit ball.
Input: An ellipsoid containing the star center, under an affine transformation as above.
Output: Either a) A cut direction d⊥ or b) A Gaussian G.
1. For each i ∈ [0, k]
1a. Evaluate f at S samples from the Gaussian Gi of width τ 0 η i in the > dimensions and
0 in the ⊥ dimensions (that is, G = N (0, σ 02 I + τ 02 η 2i I )).
width σ⊥
i
>
⊥ ⊥
0
1b. If at least 1 − 31δ
32 fraction of the evaluations are within  of the minimum evaluation
(at this iteration i), then Return Gi and Halt.

2. Otherwise, let z be the minimum of all samples in Step 1.
δ
3. Repeatedly sample the following, estimating g to within ± 32
each time
2
02
− σ⊥
)I⊥ ) and (σ> = eX ; X ← Unif[log τ 0 , log R/s])
g(µ0⊥ , σ> ) : µ0⊥ ← N (0, (σ⊥

3a. Accept the first pair (µ0⊥ , σ> ) such that g(µ0⊥ , σ> ) >

7
32 δ.

2 I + σ 2 I )]]
4. Return the gradient d⊥ = ∇⊥ [E[Lz (x) : x ← N (µ0⊥ , σ⊥
⊥
> >
0
µ⊥ changes, computed via Proposition 4.2).

(the derivative as

We note that a special case of the above algorithm is when there are no > (very thin, thinner than
τ ) dimensions. This applies, for example, at the beginning of the optimization process.
The rest of this section develops the mathematical analysis leading to the proof of Proposition 4.4.
The following lemma (Lemma 4.5) analyzes how Lz (x) (the truncated logarithm of f ) decreases as
x moves towards the global optimum, or equivalently, how Lz (x) increases as x moves away. A crucial
complicating factor is that we always average Lz (x) over x drawn from a Gaussian, and moving the
mean of a Gaussian away from the global optimum is not exactly the same thing as moving every
point in the Gaussian away from the global optimum. Lemma 4.5 expresses the effect of moving
the mean away from the global optimum—restricted to the non-thin ⊥ dimensions—in terms of a
14

positive probability, minus three confounding derivatives. If we can show that the left hand side of the
expression in the following lemma is positive, this means that a cut in the direction of the gradient of
(a Gaussian-blurred) Lz (x) in the ⊥ dimensions is guaranteed to contain the global minimum.
Lemma 4.5. For a star-convex function f with star center (global minimum) at the origin satisfying
2 , σ2 ,
f (0) ≤ z, and any mean µ = µ⊥ + µ> , and variances σ⊥
>

d 
2
2
E Lz (x) : x ← N (αµ⊥ + µ> , σ⊥
I⊥ + σ>
I> )
dα

2
2
≥ P[f (x) − z ∈ (0 , 2B) : x ← N (µ, σ⊥
I⊥ + σ>
I> )]
α=1


d 
2
2
E Lz (x) : x ← N (µ⊥ + αµ> , σ⊥
I⊥ + σ >
I> )
dα
α=1

d 
2
2
−
E Lz (x) : x ← N (µ⊥ + µ> , α2 σ⊥
I⊥ + σ>
I> )
dα
α=1

d 
2
2 2
−
E Lz (x) : x ← N (µ⊥ + µ> , σ⊥ I⊥ + α σ> I> )
dα
α=1
−

Proof. Let f ∗ ≤ z be the function value at the global minimum x = 0. By the star-convexity of
f with the origin as the star center, (f (αx) − f ∗ ) ≥ α(f (x) − f ∗ ) for all α > 1, which implies
(f (αx) − z) ≥ α(f (x) − z), and thus log(f (αx) − z) ≥ log α + log(f (x) − z). Thus the corresponding
inequality holds for Lz (x), provided, f (x) − z ∈ (0 , 2B) and α is close enough to 1, so that Lz behaves
like log f (x):


2
2
2
2
I> )] · log α
I⊥ + σ>
I> ) ≥ P[f (x) − z ∈ (0 , 2B) : x ← N (µ, σ⊥
I⊥ + α2 σ>
E Lz (x) : x ← N (αµ, α2 σ⊥


2
2
+ E Lz (x) : x ← N (µ, σ⊥ I⊥ + σ> I> )
Consider the left hand side as a function g(α); rearranging to put the rightmost term on the left
hand side, the inequality says that
2
2
I> )] · log α
I⊥ + σ>
g(α) − g(1) ≥ P[f (x) − z ∈ (0 , 2B) : x ← N (µ, σ⊥

(2)

By Proposition 4.2, g has a derivative at α = 1, which equals limα→1 g(α)−g(1)
= limα→1 g(α)−g(1)
α−1
log α
(by L’Hôpital’s rule), which by Equation 2 is thus greater than or equal to the probability P[f (x) − z ∈
2 I + σ 2 I )].
(0 , 2B) : x ← N (µ, σ⊥
⊥
> >
Therefore

d 
2
2
I⊥ + α2 σ>
I> )
E Lz (x) : x ← N (αµ, α2 σ⊥
dα

2
2
≥ P[f (x) − z ∈ (0 , 2B) : x ← N (µ, σ⊥
I⊥ + σ>
I> )]
α=1

Rewriting the left hand side of the above inequality into the sum of the four derivatives in the lemma
statement gives the result required.
As explained above, if we are able to lower bound the left hand side of Lemma 4.5, then we will
be able to make a cut to the ellipsoid that contains the global optimum, while being perpendicular to
the thin dimensions. Therefore, we need to upper bound the three derivatives on the right hand side
of the inequality of Lemma 4.5. Lemmas 4.6, 4.7, and 4.9 bound√these three terms respectively.

2B 2√ 2
In the following, we define the radius τ 0 = τ 16
to be slightly larger than τ (see
δ · log 0
π
Definition 4.14 below for details). Lemma 4.6 considers widths σ> ≥ τ 0 , in line with the operation of
our Algorithm 2.
The next two lemmas directly bound the first two derivative terms on the right in the expression
of Lemma 4.5, via a direct calculation of how fast the average of an arbitrary function can change
with respect to the Gaussian parameters.
15

Lemma 4.6. For all z, µ⊥ , σ⊥ and all ||µ> ||2 ≤ τ and σ> ≥ τ 0 we have:

d 
2
2
E Lz (x) : x ← N (µ⊥ + αµ> , σ⊥
I⊥ + σ>
I> )
dα

≤
α=1

δ
16

Proof. Given µ⊥ , σ⊥ and σ> ≥ τ 0 , define Pµ> to be the probability measure of the distribution
2 I + σ 2 I ). We have
N (µ⊥ + µ> , σ⊥
⊥
> >

Z
1
L.H.S. = lim
Lz (x)(dP(1+∆α )µ> − dPµ> )
∆α →0 |∆α |
n

Z R
1
≤ lim
|Lz (x)| dP(1+∆α )µ> − dPµ>
∆α →0 |∆α |
Rn
Z



1
2B
≤ lim
log
dP(1+∆α )µ> − dPµ> since [log 0 , log 2B] contains both Lz (x) and 0
0
∆α →0 |∆α |

n
R
 √

1
2B 2 2 ||∆α µ> ||2
√
≤ lim
log
(total variation distance between Gaussians)
∆α →0 |∆α |
0
σ>
π
 √

δ
2B 2 2 ||µ> ||2
√
≤
by the definition of τ 0
= log
0
16
π σ>

Lemma 4.7. For any z, µ, σ⊥ and σ> we have the following inequality:

d 
2
2
I> )
I⊥ + σ>
E Lz (x) : x ← N (µ, α2 σ⊥
dα

α=1



p
2B
≤ 2 dim(⊥) log
0


√
2B
≤ 2n log
0

where dim(⊥) is the number of ⊥ dimensions.
Proof. This proof is analogous to that of Lemma 4.6, making use of Pinsker’s inequality to bound the
total variation distance between Gaussians of different variances, via their KL divergence.
For our overall strategy of bounding each of the three terms in the right hand side of the inequality
of Lemma 4.5: the bound of Lemma 4.6 we use as is, however the bound of Lemma 4.7 is somewhat
larger than 1, and needs to be improved (in order to ultimately compare it with δ). We accomplish
this with an averaging argument, where we sample the mean of our Gaussian from a somewhat larger
Gaussian, which will effectively decrease the expectation of the left hand side of the expression in
Lemma 4.7 by the ratio of the variances of the two Gaussians. We analyze the effect of this “doublesampling” process on the relevant quantities in Lemma 4.8.
0 we have the following two equalities:
Lemma 4.8. For all z, µ> , σ> , µ⊥ and σ⊥ < σ⊥
2
2
02
2
E[P[f (x) − z ∈ (0 , 2B) : x ← N (µ0⊥ + µ> , σ⊥
I⊥ + σ>
I> )] : µ0⊥ ← N (µ⊥ , (σ⊥
− σ⊥
)I⊥ )]
02
2
= P[f (x) − z ∈ (0 , 2B) : x ← N (µ⊥ + µ> , σ⊥
I⊥ + σ>
I> )]

and



d 
2
2
02
2
E Lz (x) : x ← N (µ0⊥ + µ> , α2 σ⊥
I⊥ + σ>
I> )
: µ0⊥ ← N (µ⊥ , (σ⊥
− σ⊥
)I⊥ )
dα
α=1
 2

σ⊥
d 
02
2
=
E Lz (x) : x ← N (µ⊥ + µ> , α2 σ⊥
I⊥ + σ>
I> )
0
σ⊥
dα
α=1
E

16



Proof. Because adding two Gaussian random variables produces a Gaussian random variable, the
distribution
2
2
02
2
(x ← N (µ0⊥ + µ> , α2 σ⊥
I⊥ + σ>
I> ) : µ0⊥ ← N (µ⊥ , (σ⊥
− σ⊥
)I⊥ ))
02 +(α2 −1)σ 2 )I +σ 2 I ), which gives the first equality.
is equivalent to the distribution N (µ⊥ +µ> , (σ⊥
⊥ ⊥
> >
To prove the second equality, we move the derivative outside the expectation, and combine the two
expectations into the equivalent expression


02
2
2
E Lz (x) : x ← N (µ⊥ + µ> , (σ⊥
+ (α2 − 1)σ⊥
)I⊥ + σ>
I> )

Further observe that the derivatives of the α dependencies in the two expressions we are comparing
d
d 2 02
0 2
02 + (α2 − 1)σ 2 )
are dα
(σ⊥
⊥ α=1 and dα α σ⊥ α=1 , which have ratio exactly (σ⊥ /σ⊥ ) , giving the second
equality.
Having bounded the first two “confounding derivatives” from the right hand side of Lemma 4.5,
we now bound the third. We cannot directly bound this derivative, so we again employ an averaging
argument. Intuitively, this term records how Lz increases as the width σ> increases; however, since
Lz is bounded between 0 and 2B, the derivative cannot stay large over a large range of σ> . Crucially,
we let σ> vary over a huge range outside the ellipsoid, between τ 0 > τ and R/s (with s a polynomial
factor, defined in Definition 4.14). In order to give σ> a large enough range for the following bound
to be meaningful, τ 0 (and hence τ , the threshold for “thinness”) has to be exponentially small, as
specified below in Definition 4.14.
Lemma 4.9. For all z, µ, σ⊥ ,


d 
2
2
E
I> )
I⊥ + α 2 σ >
E Lz (x) : x ← N (µ⊥ + µ> , σ⊥
dα
log 2B − log 0
≤
log(R/s) − log τ 0


: (σ> = eX ; X ← Unif[log τ 0 , log R/s])
α=1

Proof. The probability density function of σ> is fσ> (y) =

1
y

·

1
log(R/s)−log τ 0 .

Let



2
2
I> ) .
I⊥ + σ>
h(σ> ) = E Lz (x) : x ← N (µ⊥ + µ> , σ⊥
The following quantity thus equals the left hand side of the inequality in the lemma statement without
1
the scaling factor log(R/s)−log
τ0 :
Z

R/s

τ0

1 d
h(αy)
y dα

Z

R/s

dy =
α=1

τ0

d
h(y) dy = h(R/s) − h(τ 0 ) ≤ log 2B − log 0
dy

since Lz (x) and hence h(σ> ) is bounded between log 0 and log 2B. The lemma statement follows.
Recalling that in Lemma 4.6 we satisfactorily bounded the second term from the right hand side
of the inequality of Lemma 4.5, we now combine the results of the previous lemmas so that we may
use Markov’s inequality to bound the sum of the remaining three terms. The following lemma makes
an assumption that the probability that f (x) − z ∈
/ (0 , 2B) is bounded away from 1; in the context
of our Algorithm 2, the case where this assumption is false turns out to imply that our algorithm
has already successfully optimized the function—despite this algorithm being intended merely to seek
another cut. This result is shown in Lemma 4.12.

17

Lemma 4.10. Given z, µ> and σ⊥ , define an auxiliary function
2
2
g(µ0⊥ , σ> ) = P[f (x) − z ∈ (0 , 2B) : x ← N (µ0⊥ + µ> , σ⊥
I⊥ + σ>
I> )]


d
2
2
−
E Lz (x) : x ← N (µ0⊥ + µ> , α2 σ⊥
I⊥ + σ>
I> )
dα
α=1

d 
2
2
−
E Lz (x) : x ← N (µ0⊥ + µ> , σ⊥
I⊥ + α 2 σ >
I> )
dα
α=1
0 > σ such that for all σ ∈ [τ 0 , R/s] we have P[f (x) − z ∈ (0 , 2B) : x ←
Suppose we have µ⊥ and σ⊥
⊥
>
02
2
N (µ⊥ + µ> , σ⊥ I⊥ + σ> I> )] ≥ 22
δ.
Thus,
letting
D
be the joint distribution of independent random
32
0
02
2
X
variables µ⊥ ← N (µ⊥ , (σ⊥ − σ⊥ )I⊥ ) and (σ> = e ; X ← Unif[log τ 0 , log R/s]), we have that

1
E
√
P[g(µ0⊥ , σ> ) > E : (µ0⊥ , σ> ) ← D] ≥
2
2(1 + 2 2n log 2B
0 )
where
22
E= δ−
32



σ⊥
0
σ⊥

2

√


2n log

2B
0


−

log 2B − log 0
log(R/s) − log τ 0

Proof. We bound the expectations of the individual terms in g(µ⊥ , σ> ) when we draw (µ⊥ , σ> ) ← D.
For the first term, we make use of Lemma 4.8 to simplify the double expectations.
2
2
I> )] : (µ⊥ , σ> ) ← D]
I⊥ + σ>
E[P[f (x) − z ∈ (0 , 2B) : x ← N (µ, σ⊥
02
2
= E[P[f (x) − z ∈ (0 , 2B) : x ← N (µ0⊥ + µ> , σ⊥
I⊥ + σ>
I> )] : (σ> = eX ; X ← Unif[log τ 0 , log(R/s)])]
22
≥ δ by assumption.
32

The other two terms are bounded by Lemmas 4.7 (via Lemma 4.8) and 4.9 respectively, yielding
that E[g(µ⊥ , σ> ) : (µ⊥ , σ⊥ ) ← D] ≥ E. Having bounded the expectation, we now upper bound g
so that we can apply Markov’s inequality. The first term of g is a probability and
√ hence is upper
bounded by 1; the next two terms are each upper bounded via Lemma 4.7 by 2n log 2B
0 . Thus
√
2B
g(µ⊥ , σ> ) ≤ 1 + 2 2n log 0 .
E
Markov’s inequality yields that P[g(µ⊥ , σ> ) > E2 ] > 2(1+2√2n
.
log 2B )
0

Using Lemma 4.10 and the previous lemmas, we have successfully upper bounded the confounding
derivatives and hence lower bounded the left hand side of Lemma 4.5, but only when the probability
term is bounded away from 0 for all σ> . As we mentioned above, in fact, if the probability term is
too small for some σ> , then it turns out we already have a Gaussian region that can be returned
as the optimization output. Intuitively, if f is “flat” (except on a small fraction of the points) on a
large Gaussian region that is known to be near the global minimum (“near”, relative to the size of the
Gaussian), then the function value on this Gaussian is essentially the global minimum, and thus this
Gaussian region may be returned as the overall answer to the optimization problem. The next two
lemmas formalize this result.
Lemma 4.11. Given a probability distribution on the positive real line with pdf proportional to f (x) =
2
e−(x−c) /2 xn−1 , for c ≥ 0 then any set with probability at least 0.9 under this distribution contains two
√
points with ratio at least 1 + 1/(6 max{c, n}).
√
c+ c2 +4(n−1)
Proof. The pdf has maximum at x value m =
(as can be found by differentiating the
2
2 /2
−(x−m)
logarithm of f ). We upper bound f by f (x) ≤ f (m) · e
, since, taking the logarithm of both
18

sides, the function value and first derivative match at x = m (the derivative is 0), while the second
derivative of the (logarithm of the) right hand side equals −1 while the left hand side yields
d2
d2
n−1
log
f
(x)
=
−1
−
log xn−1 = −1 −
dx2
dx2
x2

(3)

which is smaller than −1 for all x. Thus the integral of f√from 0 to ∞ is at most f (m) times the
2
integral of the Gaussian e−(x−m) /2 , which is at most f (m) 2π.
We correspondingly lower bound f (x).
the value m (specifying the x coordinate of the
√ We note √
c+ c2 +4(n−1)
≥ n − 1. Thus for x ≥ m, Equation 3 yields the
maximum of f ) is bounded as m =
2
2
d2
bound dx2 log f (x) ≥ −2. Analogously to above this yields the lower bound that f (x) ≥ f (m)e−(x−m)
√
for x ≥ m—where the standard deviation of this Gaussian is now 1/ 2 instead of 1.
√
2
Given this lower bound f (x) ≥ f (m)e−(x−m) for x ≥ m and the upper bound f (m) 2π on the
overall integral of f , the probability of x lying in the interval [m, m + 13 ] is at least
1
√
f (m) 2π

Z

m+1/3

2

f (m)e−(x−m) dx > 0.1

m

Similarly, the probability of x lying in the interval [m + 32 , ∞) is greater than 0.1.
Thus, any set S of probability at least 0.9 on the distribution with pdf proportional to f must
contain points from both intervals [m, m + 13 ] and [m + 23 , ∞). Hence S must contain two points with
1
ratio at least m+2/3
m+1/3 = 1 + 3m+1 .
√
√
We conclude the lemma by bounding 3m + 1 in terms ofpmax{c, n}. Let y = max{c, pn}, which
is at least 1 for n ≥ 1. We have 3m + 1 = 1 + 32 (c + c2 + 4(n − 1)) ≤ 1 + 32 (y + 5y 2 − 4),
which is easily seen to be at most 6y for y ≥ 1. Thus S contains two points with ratio at least
√
1
1 + 3m+1
≥ 1 + 1/(6 max{c, n}), as desired.
Lemma 4.12. Given a star-convex function f with global optimum at the origin, if for some location
µ we have P[f (x) ∈ [z, z + 0 ] : x ← N (µ, I)] > 0.95 then the function value at the global optimum,
√
f ∗ = f (0), satisfies f ∗ ≥ z − 60 max{||µ||2 , n}.
Proof. Consider a ray through the global minimum (at the origin), and let the ray be defined as
all positive multiples of a unit vector x. Consider the set Cx ⊂ R+ of scaling factors r such that
f (rx) ∈ [z, z + 0 ]. We note that if Cx contains two values r < r0 with some ratio κ = r0 /r then by
the star-convexity of f (on this ray), the global minimum of f must have value at least z − 0 /(κ − 1).
Thus we want to show that there exists a direction x with a set Cx that contains two widely–spaced
elements; namely, it is impossible for each Cx to only contain values within small ratio of each other.
We note that we may express P r[f (x) ∈ [z, z + 0 ] : x ← N (µ, I)] in terms of the sets Cx as an
integral in polar coordinates. Letting Sn−1 denote the (n − 1)-dimensional sphere, this probability
equals the ratio
Z
Z
Z
Z
|rx−µ|2
|rx−µ|2
n−1 −
2
r
e
dr dx
rn−1 e− 2 dr dx
(4)
x∈Sn−1

x∈Sn−1

r∈Cx

r∈R+

Consider those directions x with a positive component in the direction of µ; this halfspace H comprises
at least half the probability mass of the Gaussian N (µ, I), and thus, within this halfspace, P[f (x) ∈
[z, z + 0 ] : x ← N (µ, I) | x ∈ H] > 0.9. Suppose for the sake of contradiction that for each x ∈ H,
defining c < ||µ||2 to be the component of µ in the direction of x (where c ≥ 0), the set Cx does
√
not contain any points of ratio at least 1 + 1/(6 max{c, n}). Thus by Lemma 4.11, we have the
19

.R
R
2
2
bound r∈Cx rn−1 e−(r−c) /2 dr r∈R+ rn−1 e−(r−c) /2 dr < 0.9. For each x, with c as defined, the ratio
between these integrals
the corresponding inner integrals of Equation 4,
 equals the ratio between
R
R
|rx−µ|2
|rx−µ|2
−
−
n−1
n−1
2
2
dr
dr, yielding that the average value of this ratio over x in
e
e
r∈Cx r
r∈R+ r
the halfspace H is at most 0.9. This contradicts the fact derived earlier that P[f (x) ∈ [z, z + 0 ] : x ←
N (µ, I) | x ∈ H] > 0.9.
Therefore there must exist a direction x such that Cx contains two points of ratio at least 1 +
√
√
1/(6 max{c, n}), where since c ≤ ||µ||2 , this ratio is thus at least 1 + 1/(6 max{||µ||2 , n}). We then
√
conclude f ∗ ≥ z − 60 max{||µ||2 , n}.
We now wish to combine Lemmas 4.10 and 4.12, in the sense that we reason either one or the other
applies. The naı̈ve interpretation of the above would involve algorithmically checking the condition of
Lemma 4.12 for each and every Gaussian induced by all the σ> in the continuous (exponential) range
as in Lemma 4.10, which is obviously impossible. Therefore, we need to choose a mesh in that range
such that the total variation distance between a Gaussian with a width σ> in the > dimensions and
the closest Gaussian in the mesh is upper bounded by some quantity proportional to δ. This allows
us to reason that after sampling all points in the mesh, either Lemma 4.12 applies, in which case we
can conclude our optimization, or we can apply Lemma 4.10 to produce a cut. The following lemma
shows that a geometric spacing in the mesh suffices, thus explaining Step 1 of Algorithm 2.
2I +
Lemma 4.13. Given µ and σ⊥ , let Pσ> be the probability measure of the distribution N (µ, σ⊥
⊥
2 I ). For σ < σ 0 , we have the inequalities
σ>
>
>
>
s
r
 0  s
 0 
σ>
σ>
1
1
dim(>)
n
0 | ≤
0 ) ≤
|Pσ> − Pσ>
DKL (Pσ> ||Pσ>
log
≤
log
2
2
2
σ>
2
σ>
0 | denotes the total variation distance and DKL denotes the KL divergence. That
where 12 |Pσ> − Pσ>

σ0

>
is, if σ>
≤e
is at most p.

2p2
n

, then the difference in probabilities of event E happening under the two distributions

Proof. The first inequality is Pinsker’s inequality. A direct calculation from the standard expression
for the KL-divergence of multivariate Gaussians shows the second inequality. The third inequality is
self-evident.
We now have all the theoretical tools to reason about Algorithm 2. In the following definition, we
choose all the quantities we have used in the paper.
Definition 4.14.
!
r r
4
1
1
1
Let s = n 1 +
n + + log + log B + log R + log
3
δ

F
s
r


12 −1
δ/8
1
0
0
 = 1+ 0
σ⊥ = σ⊥
σ⊥
log 2B − log 0 2n
√
δ2
1
π
0 δ
8n
√

τ =τ
·
η
=
e
16 log 2B
0 2 2
√

0
σ⊥
=

R
τ =
s
0

k=

1
3ns
 16
2B − δ
0
 16
2B δ



log 0
log η

Finally, we present the proof of Proposition 4.4 (restated below for convenience), which establishes
the correctness of Algorithm 2.

20

Proposition 4.4 (Correctness of Algorithm 2). With negligible probability of failure, Algorithm 2
either a) returns a Gaussian region G such that P[f (x) ≤ f ∗ +  : x ← G] ≥ 1 − δ, or b) returns a
direction d⊥ , restricted to the ⊥ dimensions, such that when normalized to a unit vector d̂⊥ , the cut
1
{x : x · d̂⊥ ≤ 3n
} contains the global minimum.
Proof. To keep the proof simpler, we now translate the coordinate system to place the global optimum
1
1
at the origin. If the probability error parameter δ ≥ 20
, then make δ = 21
, which can only improve
the results.
If the algorithm halts in Step 1b, this means that for the returned Gaussian Gi , at least 1 − 31
32 δ
0
fraction of the samples were within a range [zi , zi +  ], denoting by zi the smallest observed sample
from Gi . Chernoff bounds imply that (except with probability of failure F ), the true probability
P[f (x) ∈ [zi , zi + 0 ] : x ← Gi ] is at least 1 − δ, provided we take S = poly(1/δ, | log F |) samples.
√
1
Thus we apply Lemma 4.12, since δ < 20
, and conclude that f ∗ ≥ zi − 60 max{||µ||2 , n}, where
we must now bound “µ”, which from the notation of Lemma 4.12 means the “distance between the
global minimum and the center of the ellipsoid, measured by the number of standard deviations of
Gaussian Gi ”. Specifically, in the ⊥ dimensions, the ellipsoid has radius 1 and the Gaussian Gi =
02 I + τ 02 η 2i I ) has radius σ 0 ; in the > dimensions, the Gaussian has radius larger than the
N (0, σ⊥
⊥
>
⊥
ellipsoid, so the distance in these directions is less than 1, which is certainly less than σ10 . Thus in total
⊥
√
we bound ||µ||2 ≤ σ20 . Plugging this bound into the conclusion of Lemma 4.12, where ||µ||2 ≥ n, we
⊥

have f ∗ ≥ zi − 60 σ20 ≥ zi + 0 − , by the definition of 0 . Thus, the returned Gaussian satisfies the
⊥
desired properties of the output of our optimization algorithm: P[f (x) ∈ [f ∗ , f ∗ + ] : x ← Gi ] ≥ 1 − δ.
We now analyze the situation when the algorithm does not halt in Step 1b. In this case, for each
i, the proportion of samples larger than zi + 0 was observed to be at least 31
32 δ. Since z = mini zi ,
31
0
the proportion larger than z +  must also be at least 32 δ. As above, by Chernoff bounds, for each i
(except with probability of failure F ), the true probability P[f (x) > z + 0 : x ← Gi ] is at least 30
32 δ,
provided we take S = poly(1/δ, | log F |) samples.
Since as i ranges from 0 to k, the Gaussians Gi vary exponentially in their width in the > dimensions,
as σ>,i = τ 0 η i , these widths form a fine (exponentially spaced) mesh over this entire region from τ 0
0 ∈ [τ 0 , R ], there is a an i for which the Gaussian G has width σ
to Rs . Thus, for any σ>
i
>,i such that
s
0
σ>
σ>,i

02 I + σ 02 I ) we invoke Lemma 4.13 to conclude that
≤ η. Thus for the Gaussian G 0 = N (0, σ⊥
⊥
> >
p
0
the difference in the probability of f (x) > z +  between x ← Gi and x ← G 0 is at most n2 log η = 4δ .
0 ∈ [τ 0 , R ], letting G 0 = N (0, σ 02 I +σ 02 I ), the probability P[f (x) > z +0 : x ← G 0 ]
Thus for any σ>
> >
⊥ ⊥
s
22
is at least 32 δ. (Namely, given that the algorithm did not halt in Step 1b, we have a guarantee that
0 , despite only taking a polynomial (k + 1) number
holds over an exponentially wide range of widths σ⊥
of iterations to check, and a polynomial (S) number of samples from f per iteration.)
In order to apply Lemma 4.10, we need the tiny variant of the above claim, where instead of
bounding the probability that f (x) > z +0 , we instead need a bound on the probability that f (x)−z ∈
(0 , 2B). However, B was chosen to be a truly huge number, such that we have the global guarantee
of Definition 2.3 that for all x within distance 10nR of the origin, |f (x)| ≤ B. We consider any
case where our algorithm evaluates the function outside this ball to be a failure of the algorithm.
Because by the condition of Step 2e of Algorithm 1, the ellipsoid under consideration in Algorithm 2
has semi-principal axes of length at most 3nR, and its center lies within R of the (original) origin,
each point in the ellipsoid has distance at most (3n + 1)R from the origin. Further, by construction,
each Gaussian has standard deviation in the ⊥ dimensions bounded by the size of the ellipsoid over
s, and standard deviation in the > dimensions at most Rs , where s is chosen so that, over the entire
course of the algorithm, no Gaussian sample will ever be more than s standard deviations from its
mean, except with negligible probability. Thus, except with negligible probability, all samples from the
algorithm are in the region within 10nR of the origin (as a very loose bound), and thus have function

1≤

21

value |f (x)| ≤ B. Having analyzed this negligible probability of failure, we assume for the rest of this
proof that all function evaluations have magnitude less than B, and condition all probabilities on the
assumption that no failure has occurred.
0 ∈ [τ 0 , R ] we have P[f (x) − z ∈ (0 , 2B) : x ← G 0 ] ≥ 22 δ, and we may now
Thus for any σ>
s
32
invoke Lemma 4.10, in order to describe the properties of the function g(µ0⊥ , σ> ) at the center of Step
3. Using the parameters defined in Definition 4.14, the bound E of Lemma 4.10 is found to equal
1
2 δ. Lemma 4.10 thus yields (with distribution D as defined in the lemma, and as used in Step 3 of
Algorithm 2) that
1
δ
√
P[g(µ0⊥ , σ> ) > δ : (µ0⊥ , σ> ) ← D] ≥
4
4(1 + 2 2n log 2B
0 )
In Step 3, we estimate each of the three terms of g via Proposition 4.2, and take enough samples
1
to ensure that, except with negligible probability, our estimate of g is accurate to within 32
δ. Thus
1
7
0
0
0
for each sampled µ⊥ , σ> such that g(µ⊥ , σ> ) > 4 δ, our estimate of g(µ⊥ , σ> ) will be at least 32
δ, and
thus the condition in Step 3 will succeed.
δ
Each iteration of Step 3 of the algorithm thus succeeds with probability at least 4(1+2√2n
,
log 2B )
0

and Chernoff bounds imply that, except with failure probability F , Step 3 will successfully terminate
in poly(| log F |, n, 1/δ, log 2B
0 ) many iterations.
7
Given that our observed value of g(µ0⊥ , σ> ) when Step 3 terminates is at least 32
δ and our estimates
1
6
0
are accurate to within 32 δ, the true value of g(µ⊥ , σ> ) must be at least 32 δ.
2 I + σ 2 I )] estimated in Step 4
We now show that the gradient d⊥ = ∇⊥ [Lz (x) : x ← N (µ0⊥ , σ⊥
⊥
> >
has positive component in the direction away from the global optimum, which will enable us to make
a cut. Lemma 4.5 shows that the component of the gradient at location µ0⊥ in the direction away from
the global optimum is at least
g(µ0⊥ , σ> ) −


d 
2
2
I> )
I⊥ + σ>
E Lz (x) : x ← N (µ⊥ + αµ> , σ⊥
dα

α=1

6
From the above bound g(µ0⊥ , σ> ) ≥ 32
δ and the bound of Lemma 4.6 the term we subtract from g
δ
is at most 16 , we conclude that the component of the gradient in the direction away from the global
optimum is at least 81 δ.
We thus estimate the gradient in Step 4 by estimating the derivative in each of the dimensions in
1
⊥ to within error 16n
δ, which guarantees that the total gradient vector has component in the direction
1
away from the global optimum at least 16
δ. (As an implementation detail, we note that Proposition 4.2
computes a version of the derivative scaled by σ⊥ , to any desired accuracy κ, using samples scaling
1
polynomially with κ1 . Thus for the unscaled derivative to have accuracy 16n
δ, the scaled derivative
1
must have accuracy 16n δ · σ⊥ , which requires time polynomial to the inverse of this quantity. From the
definition of σ⊥ in Definition 4.14, its inverse is polynomial in the overall parameters of the algorithm.)
Thus our estimate of the gradient, normalized to a unit vector d̂⊥ , defines a halfspace {x : (x −
µ0⊥ ) · d̂⊥ ≤ 0} that contains the global optimum. Further, since µ0⊥ came from a sample µ0⊥ ←
02 − σ 2 )I ), where σ 0 = 1 (in Definition 4.14), and s was chosen so that no Gaussian
N (0, (σ⊥
⊥ ⊥
⊥
3ns
1
sample will ever be more than s standard deviations from its mean, we have that ||µ0⊥ ||2 ≤ 3n
.
1
0
Thus µ⊥ · d̂⊥ ≤ 3n , from which we conclude that the global optimum is contained in the halfspace
1
{x : x · d̂⊥ ≤ 3n
}, as desired.

22

References
[1] Karthekeyan Chandrasekaran, Daniel Dadush, and Santosh Vempala. Thin partitions: Isoperimetric inequalities and a sampling algorithm for star shaped bodies. In Proc. SODA’10, pages
1630–1645, 2010.
[2] Eden Chlamtac and Madhur Tulsiani. Handbook on Semidefinite, Conic and Polynomial Optimization, chapter Convex Relaxations and Integrality Gaps, pages 139–169. Springer US, Boston,
MA, 2012.
[3] Paul Christiano, Jonathan A. Kelner, Aleksander Madry, Daniel A. Spielman, and Shang-Hua
Teng. Electrical flows, Laplacian systems, and faster approximation of maximum flow in undirected graphs. In Proc. STOC’11, pages 273–282, 2011.
[4] John C Duchi, Peter L Bartlett, and Martin J Wainwright. Randomized smoothing for stochastic
optimization. SIAM J. Optimiz., 22(2):674–701, 2012.
[5] Elad Hazan, Kfir Levy, and Shai Shalev-Shwartz. Beyond convexity: Stochastic quasi-convex
optimization. In Proc. NIPS’15, pages 1594–1602, 2015.
[6] Jonathan A. Kelner, Yin Tat Lee, Lorenzo Orecchia, and Aaron Sidford. An almost-linear-time
algorithm for approximate max flow in undirected graphs, and its multicommodity generalizations.
In Proc. SODA’14, pages 217–226, 2014.
[7] Alex Krizhevsky, Ilya Sutskever, and Geoffrey E Hinton. Imagenet classification with deep convolutional neural networks. In Proc. NIPS’12, pages 1097–1105, 2012.
[8] Yin Tat Lee and Aaron Sidford. Path finding methods for linear programming: Solving linear
programs in Õ(vrank) iterations and faster algorithms for maximum flow. In Proc. FOCS’14,
pages 424–433, 2014.
[9] László Lovász. An algorithmic theory of numbers, graphs and convexity. SIAM, 1987.
[10] Aleksander Madry. Navigating central path with electrical flows: From flows to matchings, and
back. In Proc. FOCS’13, pages 253–262, 2013.
[11] A.S. Nemirovskii and D.B. Yudin. Problem Complexity and Method Efficiency in Optimization.
Wiley, 1983.
[12] Yurii Nesterov. A method of solving a convex programming problem with convergence rate
o(1/k 2 ). Soviet Mathematics Doklady, 27(2):372–376, 1983.
[13] Yurii Nesterov and B.T. Polyak. Cubic regularization of Newton method and its global performance. Math. Program., 108(1):177–205, 2006.
[14] Walter Rudin. Well-distributed measurable sets. Am. Math. Monthly, 90(1):41–42, 1983.
[15] Jonah Sherman. Nearly maximum flows in nearly linear time. In Proc. FOCS’13, pages 263–269,
2013.
[16] David Silver, Aja Huang, Chris J Maddison, Arthur Guez, Laurent Sifre, George van den Driessche, Julian Schrittwieser, Ioannis Antonoglou, Veda Panneershelvam, Marc Lanctot, et al. Mastering the game of go with deep neural networks and tree search. Nature, 529(7587):484–489,
2016.
23

[17] Daniel A. Spielman and Shang-Hua Teng. Nearly-linear time algorithms for graph partitioning,
graph sparsification, and solving linear systems. In Proc. STOC’04, pages 81–90, 2004.
[18] Paul Valiant. Evolvability of real functions. ACM Trans. Comput. Theory, 6(3):12:1–12:19, 2014.
[19] V. Vapnik. The Nature of Statistical Learning Theory. Information Science and Statistics. Springer
New York, 2013.
[20] Nisheeth K Vishnoi. Lx = b. Found. Trends Theoretical Computer Science, 8(1-2):1–141, 2012.

24

A

The Scope of Star-Convex Functions

In this section we introduce some basic building blocks to construct large families of star-convex
functions, with the aim of demonstrating the richness of the definition. More technical examples of
star-convex functions are introduced in the rest of the paper to illustrate particular points (including
Examples 1.5, 2.1, and Example B.1 below).
Basic Star-Convex Functions:
1. Any convex function is star-convex.
2. Even in 1 dimension, star-convex functions can be non-convex:
(
|x|
if |x| < 1
f (x) =
|2x| otherwise
Or, for a Lipschitz twice-differentiable example, f (x) = |x|(1 − e−|x| ) (from Nesterov and
Polyak [13]).
3. In n > 1 dimensions, take an arbitrary positive function g(x) on the unit sphere, and extend it
to the origin in a way that is star-convex on each ray through the origin, for example, extending
g linearly to define


x
f (x) = ||x||2 · g
||x||2
Ways to Combine Star-Convex Functions: Given any star-convex functions f, g that have star
centers at f (0) = g(0) = 0, we can combine them to generate a new star convex function in the
following ways:
4. A star-convex function can be shifted in the x and f senses so that its global optimum is at
an arbitrary x and f value; in general, for any affine transformation A and real number c, the
function h(x) = f (A(x)) + c is star-convex.
5. For any positive power p ≥ 1, the function h(x) = f (x)p is star-convex.
6. The sum h(x) = f (x) + g(x) is star-convex.
7. The product f (x) · g(x) is star-convex.


f (x)p +g(x)p
2

1/p

is star-convex, defining
p
powers via limits as appropriate. (The case p = 0 corresponds to h(x) = f (x)g(x).)

8. The power mean: for real number p, the function h(x) =

Practical Examples of Star-Convex Function Classes:
We combine the basic examples above to yield the following general classes of star-convex functions.
9. Sums of squares of monomials: the square of any monomial (xi y j z k · · · )2 , with i, j, k, . . . ∈ Z+ ,
will be star-convex about the origin, and thus sums of such terms will be star-convex despite
typically being non-convex, for example h(x, y) = x2 y 2 + x2 + y 2 . Sums-of-squares arise in many
different contexts.
10. More generally, any polynomial of |x|, |y|, |z|, . . . with nonnegative coefficients is star-convex.

25

11. In the standard machine learning setting where one is trying to learn a parameter vector θ: for a
given hypothesis θ̂, and each training example Xi , the hypothesis gives an error LXi (θ, θ̂) that is
(typically) a convex function of its second argument, θ̂. Averaging these convex losses over all i,
via the power mean of exponent p leads to a star-convex loss, as a function of θ̂, with minimum
value 0 at θ̂ = θ. The paper by one of the authors [18] discusses a class of such star-convex loss
functions which have the form
!1/p
X
hθ,X (θ̂) =
|(θ̂ − θ) · Xi |p
i

Note that, for p ≥ 1, the loss function is convex.

B

Example of a Pathological Star-Convex Function

See the discussion before Definition 2.4.
Example B.1. We present a class of minimization problems where we argue that no algorithm will be
able to return a (rational) point with objective function value near the global minimum, and instead,
in order to be effective, an algorithm must return a region on which “the function value is low with
high probability”.
For a given δ ∈ (0, 1) we define a class of minimization problems in the unit square, parameterized
by three uniformly random real numbers X, Y, θ ← [0, 1]. The star-convex function f (x, y) will have
unique global minimum at f (X, Y ) = 0, and is defined as
(
100
|x − X| if y = Y or f rac(angle(X − x, Y − Y ) · 1010 + θ) < δ
f (x, y) = |y − Y | +
,
0
otherwise
where f rac(x) denotes the fractional part of x, and angle(x, y) denotes the angle of location (x, y)
about the origin.
We first note that the “if y = Y ” condition is only possible if Y is rational (since the input y must
be rational), and hence this condition occurs with probability 0; this condition is part of the definition
of f only to avoid making the entire line y = Y a global minimum of f . We thus ignore the y = Y
possibility in what follows.
Fixing x, y, X, Y , the probability with respect to the random choice of θ that the evaluation of f (x, y)
falls into case 1 equals δ. Since θ only affects f (x, y) in case 1, if 1/δ is significantly larger than the
runtime of the algorithm, then with high probability the algorithm will never observe an evaluation via
case 1.
Given that the algorithm never observes an evaluation via case 1, and the only dependence of f on
X is via case 1, the algorithm will have no information about X, with high probability.
For any x that is -far from X (for some ), the function f will be at least , with probability at
100
least δ over choices of θ. Further, since regions of high function value appear every 10−10 radians
99
about the origin, for any ball of radius 10−10 centered at x and any θ, the the function f will take
value at least  with probability close to δ on a random point in the ball.
Thus, any optimization algorithm that runs in time o( 1δ ) and outputs a set that can be decomposed
99
into 10−10 -radius balls can expect the function value on its output set to be low only with 1 − 1/δ
probability. In short, we cannot expect any reasonable algorithm to return a set on which the function
value is always near-optimal, or even near-optimal with probability 1. The best we can hope for is a
polynomial relation between the error probability and the runtime, as we achieve in Theorem 2.5.

26

Letting the optimization algorithm specify output regions of double-exponentially small size does not
help, since the set of angles for which case 1 applies could be replaced by the construction by Rudin [14]
of a measurable set on [0, 1], whose intersection with any subinterval 0 < a < b < 1 has measure strictly
between 0 and (b − a) · δ Furthermore, letting the optimization algorithm return rational points, instead
of full-dimension sets, still does not help, since we could modify f in the style of Example 2.1 so that
its values on all the rationals with x coordinate more than 10−100 far from X are high.

C

Details for Adapting the Ellipsoid Algorithm

In this section we fill in many standard adaptations of the ellipsoid algorithm omitted from the body
of Section 3, supporting the analysis of Algorithm 1.
In analogy with the standard ellipsoid algorithm, the following lemma provides a means to construct
a new ellipsoid from the current ellipsoid and a cutting plane. Such scenarios can always be affinely
transformed into cutting an origin-centered unit ball along a basis direction.
Lemma C.1. Consider the unit ball Bn−1 in n-dimensional space centered at the origin, and the
1
half-space H = {x ∈ Rn : x1 ≥ − 3n
}, where xi denotes the coordinate of x in the ith dimension. Then
the ellipsoid E contains the intersection Bn−1 ∩ H, where


2



2−1X
1
n
2
2
x
≤
1
E = x ∈ Rn :
+
x
−
1
i
2


3(n + 1)
n2
(1 − 3(n+1)
)2
i6=1

P
Proof. Consider an arbitrary point x ∈ Bn−1 , then we have i6=1 x2i ≤ 1 − x21 . Thus, it suffices to
show that the function

2
1
2
n2 − 1
f (x1 ) =
x
−
+
(1 − x21 )
1
2
2
2
3(n
+
1)
n
(1 − 3(n+1) )
1
, 1].
satisfies f (x1 ) ≤ 1 for all x1 ∈ [− 3n
The proof is simple. The two roots of f (x1 ) = 1 are

−(5n+1)
12n2 +5n+1

and 1. Also, f 0 (1) > 0, and since
h
i
−(5n+1)
f is quadratic in x1 , it must be the case that f (x1 ) ≤ 1 for all x1 ∈ 12n
2 +5n+1 , 1 . Note that since,
1
3n

≤

5n+1
12n2 +5n+1

1
, 1].
for all n ≥ 1, we have f (x1 ) ≤ 1 for all x ∈ [− 3n

We also want to upper bound the volume of the new ellipsoid relative to the original one, in order
to show a fixed ratio of volume decrease at each round.
Lemma C.2. The volume of the ellipsoid constructed in Lemma C.1 is less than e
of the unit ball.

1
− 6(n+1)

times that

Proof. The ratio between the volumes is


2
1−
3(n + 1)



n2
n2 − 1

 n−1
2

2
− 3(n+1)

<e

n−1

1
− 6(n+1)

e 2(n2 −1) = e

The next two lemmas demonstrate the straightforward constructions required by Step 2e in Algorithm 1 to keep all semi-principal axes of the feasible ellipsoid bounded by 3nR, and keep the ellipsoid
centered within distance R of the origin.
27

P
i 2
Lemma C.3. Let E be an ellipsoid in Rn consisting of the points x = (x1 , . . . , xn ) such that ni=1 ( xia−c
) ≤
i
1, namely with center c and semi-principal axis lengths specified by the vector a. If the first j semiprincipal axes have lengths ai ≥ 3n (for i ≤ j), then replacing these axes with n, multiplying the
remaining axes (j + 1, . . . , n) by (n + 1)/n, and replacing the first j elements of the center c with 0
yields an ellipsoid E 0 defined by
j 
X
x i 2
i=1

n

+


n 
X
n
xi − ci 2
≤1
·
n+1
ai

(5)

i=j+1

that has smaller volume than E and contains the entire intersection of E with the unit ball.
Proof. Let E>j be the (n − j)-dimensional
defined as E restricted to dimensions > j, namely,
Pnellipsoid,
xi −ci 2
points x>j = (xj+1 , . . . , xn ) satisfying i=j+1 ( ai ) ≤ 1, and let Bj be the j-dimensional unit ball
centered at the origin. Then the intersection of E with the n-dimensional unit ball is contained in the
cartesian product Bj × E>j . We now show that E 0 contains this cartesian product, by showing that
points in this cartesian product satisfy Equation 5.
Consider a point x ∈ Bj × E>j . The first j coordinates of x, corresponding to Bj , satisfy
Pj
Pj xi 2
1
2
i=1 xi ≤ 1, and thus contribute to Equation
i=1 ( n ) ≤ n2 . The remaining coordiPn 5 xati −cmost
2
nates of x, corresponding to E>j , satisfy i=j+1 ( ai i ) ≤ 1 and thus contribute to Equation 5 at
P
n xi −ci 2
n 2
most ni=j+1 ( n+1
ai ) ≤ ( n+1 ) . Adding these two bounds yields that, for points x ∈ Bj × E>j , the
n 2
left hand side of Equation 5 is at most n12 + ( n+1
) ≤ 1, as desired, proving the lemma.
Lemma C.4. The intersection of a ball B and an ellipsoid E is contained within the intersection of
the ellipsoid and a version of the ball translated to be centered at the nearest point of E to the center
of B.
Lemma C.4, interpreted under an affine transformation, and combined with Lemma C.3, shows
that we can run the ellipsoid method with every ellipsoid guaranteed to have axes bounded by 3nR,
and always centered inside the ball of radius R.
The following two lemmas give important invariants that the algorithm maintains throughout the
iterations.
Lemma C.5. The global minimum is contained in the ellipsoid Ei for all i.
Proof. We prove by induction. The base case is trivial.
For the inductive case, assume the lemma is true for the ellipsoid Ei . Apply an affine transformation
such that the Ei is a unit ball. By assumption, the cut d⊥ satisfies the property that, taking into
1
account the affine transformation, the half-space H = {x ∈ Rn : x · d̂⊥ ≥ − 3n
} contains the global
minimum. Since the global minimum is also contained in Ei by assumption, it must be contained in
Ei+1 by Lemma C.1.
 1 
1+ 3n
Lemma C.6. No semi-principal axis of any ellipsoid Ei is ever less than
τ.
2
Proof. We prove by induction. The base case is again trivial.
For the inductive case, assume the lemma is true for the ellipsoid Ei . For simplicity, when we refer
to thin (less than τ ) and non-thin (at least τ ) directions, they are always with respect to the lengths
of the axes of Ei . Observe that, by construction, all the thin directions in the axes of Ei are also axes
directions of Ei+1 . Also, by Step 2(c) in Algorithm 1, Ei+1 always contain the center of Ei . Therefore,
the semi-principal axes of Ei+1 in the thin directions must be longer than that of Ei , which is at least
the quantity in the lemma statement by the induction hypothesis. Now consider an arbitrary axis v
of Ei+1 that is not in any of the thin directions. The direction of this axis is a linear combination of
28

the non-thin directions in the axes of Ei . Therefore, the diameter ofEi in 
the direction of v is at least
1
1+
1
3n
2τ . By construction of the cut, 2||v||2 ≥ (1 + 3n
)τ , and so ||v||2 ≥
τ , completing the proof of
2
the lemma.
We restate and prove Lemma 3.1, showing that Algorithm 1 halts after the designated number of
iterations.
Lemma (Lemma 3.1). Algorithm 1 halts within m + 1 iterations either through Step 2a or Step 2c.
Proof. By Lemma C.2 and induction, after m iterations, the volume of ellipsoid Em is at most
 1  n−1
1+ 3n
τ times the volume of the input radius R ball. Since, by Lemma C.6, all its semiτ
2
 1 
1+ 3n
principal axes have length at least
τ , all the axes must also be of length at most τ .
2
Therefore, the algorithm will definitely halt in Step 2a of iteration m + 1, if it has not already done
so in the previous iterations.

29

