Distribution-Specific Hardness of Learning Neural Networks

arXiv:1609.01037v2 [cs.LG] 9 Mar 2017

Ohad Shamir
Weizmann Institute of Science
ohad.shamir@weizmann.ac.il

Abstract
Although neural networks are routinely and successfully trained in practice using simple gradientbased methods, most existing theoretical results are negative, showing that learning such networks is
difficult, in a worst-case sense over all data distributions. In this paper, we take a more nuanced view, and
consider whether specific assumptions on the “niceness” of the input distribution, or “niceness” of the
target function (e.g. in terms of smoothness, non-degeneracy, incoherence, random choice of parameters
etc.), are sufficient to guarantee learnability using gradient-based methods. We provide evidence that
neither class of assumptions alone is sufficient: On the one hand, for any member of a class of “nice”
target functions, there are difficult input distributions. On the other hand, we identify a family of simple
target functions, which are difficult to learn even if the input distribution is “nice”. To prove our results,
we develop some tools which may be of independent interest, such as extending Fourier-based hardness
techniques developed in the context of statistical queries [3], from the Boolean cube to Euclidean space
and to more general classes of functions.

1

Introduction

Artificial neural networks have seen a dramatic resurgence in recent years, and have proven to be a highly effective machine learning method in computer vision, natural language processing, and other challenging AI
problems. Moreover, successfully training such networks is routinely performed using simple and scalable
gradient-based methods, in particular stochastic gradient descent.
Despite this practical success, our theoretical understanding of the computational tractability of such
methods is quite limited, with most results being negative. For example, as discussed in [17], learning even
depth-2 networks in a formal PAC learning framework is computationally hard in the worst case, and even
if the algorithm is allowed to return arbitrary predictors. As common in such worst-case results, these are
proven using rather artificial constructions, quite different than the real-world problems on which neural
networks are highly successful. In particular, since the PAC framework focuses on distribution-free learning
(where the distribution generating the examples is unknown and rather arbitrary), the hardness results rely on
carefully crafted distributions, which allows one to relate the learning problem to (say) an NP-hard problem
or breaking a cryptographic system. However, what if we insist on “natural” distributions? Is it possible
to show that neural networks learning becomes computationally tractable? Can we show that they can be
learned using the standard heuristics employed in practice, such as stochastic gradient descent?
To understand what a “natural” distribution refers to, we need to separate the distribution over examples
(given as input-output pairs (x, y)) into two components:
• The input distribution p(x): “Natural” input distributions on Euclidean space tend to have properties
such as smoothness, non-degeneracy, incoherence etc.
1

• The target function h(x): In PAC learning, it is assumed that the output y equals h(x), where h
is some unkown target function from the hypothesis class we are considering. In studying neural
networks, it is common to consider the class of all networks which share some fixed architecture (e.g.
feedforward networks of a given depth and width). However, one may argue that the parameters of
real-world networks (e.g. the weights of each neuron) are not arbitrary, but exhibit various features
such as non-degeneracy or some “random like” appearance. Indeed, networks with a random structure
have been shown to be more amenable to analysis in various situations (see for instance [6, 2, 5] and
references therein).
Empirical evidence clearly suggest that many pairs of input distributions and target functions are computationally tractable, using standard methods. However, how do we characterize these pairs? Would appropriate
assumptions on one of them be sufficient to show learnability?
In this paper, we investigate these two components, and provide evidence that neither one of them alone
is enough to guarantee computationally tractable learning, at least with methods resembling those used in
practice. Specifically, we focus on simple, shallow ReLU networks, assume that the data can be perfectly
predicted by some such network, and even allow over-specification (a.k.a. improper learning), in the sense
that we allow the learning algorithm to output a predictor which is possibly larger and more complex than
the target function (this technique increases the power of the learner, and was shown to make the learning
problem easier in theory and in practice [17, 20, 21]). Even under such favorable conditions, we show the
following:
• Hardness for “natural” target functions. For each individual target function coming from a simple
class of small, shallow ReLU networks (even if its parameters are chosen randomly or in some other
oblivious way), we show that no algorithm invariant to linear transformations can successfully learn
it w.r.t. all input distributions in polynomial time (this corresponds, for instance, to standard gradientbased methods together with data whitening or preconditioning). This result is based on a reduction
from learning intersections of halfspaces. Although that problem is known to be hard in the worstcase over both input distributions and target functions, we essentially show that invariant algorithms
as above do not “distinguish” between worst-case and average-case: If one can learn a particular target
function with such an algorithm, then the algorithm can learn nearly all target functions in that class.
• Hardness for “natural” input distributions. We show that target functions of the form x 7→
ψ(hw, xi) for any periodic ψ are generally difficult to learn using gradient-based methods, even if
the input distribution is fixed and belongs to a very broad class of smooth input distributions (including, for instance, Gaussians and mixtures of Gaussians). Note that such functions can be constructed
by simple shallow networks, and can be seen as an extension of generalized linear models [18]. Unlike the previous result, which relies on a computational hardness assumption, the results here are
geometric in nature, and imply that the gradient of the objective function, nearly everywhere, contains virtually no signal on the underlying target function. Therefore, any algorithm which relies on
gradient information cannot learn such functions. Interestingly, the difficulty here is not in having a
plethora of spurious local minima or saddle points – the associated stochastic optimization problem
may actually have no such critical points. Instead, the objective function may exhibit properties such
as flatness nearly everywhere, unless one is already very close to the global optimum. This highlights
a potential pitfall in non-convex learning, which occurs already for a slight extension of generalized
linear models, and even for “nice” input distributions.
Together, these results indicate that in order to explain the practical success of neural network learning with
gradient-based methods, one would need to employ a careful combination of assumptions on both the input
2

distribution and the target function, and that results with even a “partially” distribution-free flavor (which
are common, for instance, in convex learning problems) may be difficult to attain here.
To prove our results, we develop some tools which may be of independent interest. In particular, the
techniques used to prove hardness of learning functions of the form x 7→ ψ(hw, xi) are based on Fourier
analysis, and have some close connections to hardness results on learning parities in the well-known framework of learning from statistical queries [14]: In both cases, one essentially shows that the Fourier transform
of the target function has very small support, and hence does not “correlate” with most functions, making
it difficult to learn using certain methods. However, we consider a more general and arguably more natural
class of input distributions over Euclidean space, rather than distributions on the Boolean cube. In a sense,
we show that learning general periodic functions over Euclidean space is difficult (at least with gradientbased methods), for the same reasons that learning parities over the Boolean cube is difficult in the statistical
queries framework.

Related Work
Recent years have seen quite a few papers on the theory of neural network learning. Below, we only briefly
mention those most relevant to our paper.
In a very elegant work, Janzamin et al. [13] have shown that a certain method based on tensor decompositions allows one to provably learn simple neural networks by a combination of assumptions on the input
distribution and the target function. However, a drawback of their method is that it requires rather precise
knowledge of the input distribution and its derivatives, which is rarely available in practice. In contrast, our
focus is on algorithms which do not utilize such knowledge. Other works which show computationallyefficient learnability of certain neural networks under sufficiently strong distributional assumptions include
[2, 17, 1, 22].
In the context of learning functions over the Boolean cube, it is known that even if we restrict ourself
to a particular input distribution (as long as it satisfies some mild conditions), it is difficult to learn parity
functions using statistical query algorithms [14, 3]. Moreover, it was recently shown that stochastic gradient
descent methods can be approximately posed as such algorithms [9]. Since parities can be implemented with
small real-valued networks, this implies that for “most” input distributions on the Boolean cube, there are
neural networks which are unlikely to be learnable with gradient-based methods. However, data provided
to neural networks in practice are not in the form of Boolean vectors, but rather vectors of floating-point
numbers. Moreover, some assumptions on the input distribution, such as smoothness and Gaussianity, only
make sense once we consider the support to be Euclidean space rather than Boolean cube. Perhaps these
are enough to guarantee computational tractability? A contribution of this paper is to show that this is not
the case, and to formally demonstrate how phenomena similar to the Boolean case also occurs in Euclidean
space, using appropriate target functions and distributions.
Finally, we note that [15] provides improper-learning hardness results, which hold even for a standard
Gaussian distribution on Euclidean space, and for any algorithm. However, unlike our paper, their focus is
on hardness of agnostic learning (where the target function is arbitrary and does not have to correspond to
a given class), the results are specific to the standard Gaussian distribution, and the proofs are based on a
reduction from the Boolean case.
The paper is structured as follows: In Sec. 2, we formally present some notation and concepts used
throughout the paper. In Sec. 4, we provide our hardness results for natural input distributions, and in Sec. 3,
we provide our hardness results for natural target functions. All proofs are presented in Sec. 5.

3

2

Preliminaries

We generally let bold-faced letters denote vectors. Given
a complex-valued number z = a + ib, we let
√
2
z = a − ib denote its complex conjugate, and |z| = a + b2 denote its modulus. Given a function f , we
let ∇f denote its gradient and ∇2 f denote its Hessian (assuming they exist).
Neural Networks. The focus of our results will be on learning predictors which can be described by
simple and shallow (depth 2 or 3) neural networks. A standard feedforward neural network is composed of
neurons, each of which computes the mapping x 7→ σ(hw, xi + b), where w, b are parameters and σ is a
scalar activation function, for example the popular ReLU function [z]+ = max{0, z}. These neurons are
arranged in parallel in layers, so the output of each layer can be compactly represented as x 7→ σ(W > x+b),
where W is a matrix (each column corresponding to the parameter vector of one of the neurons), b is a
vector, and σ applies an activation function on the coordinates of W > x. In vanilla feedforward networks,
such layers are connected to each other, so given an input x, the output equals
>
σk (Wk> σk−1 (Wk−1
. . . σ2 (W2> σ1 (W1> x + b1 ) + b2 ) . . . + bk−1 ) + bk ),

where Wi , bi , σi are parameter of the i-th layer. The number of layers k is denoted as the depth of the
network, and the maximal number of columns in Wi is denoted as the width of the network. For simplicity,
in this paper we focus on networks which output a real-valued number, and measure our performance with
respect to the squared loss (that is, given an input-output example (x, y), where x is a vector and y ∈ R, the
loss of a predictor p on the example is (p(x) − y)2 ).
Gradient-Based Methods. Gradient-based methods are a class of optimization algorithms for solving
problems of the form minw∈W F (w) (for some given function F and assuming w is a vector in Euclidean
space), based on computing ∇F (w) of approximations of ∇F (w) at various points w. Perhaps the simplest
such algorithm is gradient descent, which initializes deterministically or randomly at some point w1 , and
iteratively performs updates of the form wt+1 = wt − ηt ∇F (wt ), where ηt > 0 is a step size parameter. In
the context of statistical supervised learning problems, we are usually interested in solving problems of the
form minw∈W Ex∼D [`(f (w, x), h(x))], where {x 7→ f (w, x) : w ∈ W} is some class of predictors, h is a
target function, and ` is some loss function. Since the distribution D is generally unknown, one cannot compute the gradient of this function w.r.t. w directly, but can still compute approximations, e.g. by sampling
one x at random and computing the gradient (or sub-gradient) of `(f (w, P
x), h(x)). The same approach can
m
1
be used to solve empirical approximations of the above, i.e. minw∈W m
i=1 `(f (w, xi ), h(xi )) for some
m
dataset {(xi , h(xi ))}i=1 . These are generally known as stochastic gradient methods, and are one of the most
popular and scalable machine learning methods in practice.
PAC Learning. For the results of Sec. 3, we will rely on the following standard definition of PAC
learning with respect to Boolean functions: Given a hypothesis class H of functions from {0, 1}d to {0, 1},
we say that a learning algorithm PAC-learns H if for any  ∈ (0, 1), any distribution D over {0, 1}d ,
and any h? ∈ H, if the algorithm is given oracle access to i.i.d. samples (x, h? (x)) where x is sampled
according to D, then in time poly(d, 1/), the algorithm returns a function f : {0, 1}d 7→ {0, 1} such that
Prx∼D (f (x) 6= h? (x)) ≤  with high probability (for our purposes, it will be enough to consider any
constant close to 1). Note that in the definition above, we allow f not to belong to the hypothesis class H.
This is often denoted as “improper” learning, and allows the learning algorithm more power than in “proper”
learning, where f must be a member of H.
Fourier Analysis on Rd . In the analysis of Sec. 4, we will consider functions from Rd to the reals R or
complex numbers C, and view them as elements in the Hilbert space L2 (Rd ) of square integrable functions,

4

equipped with the inner product
Z
f (x) · g(x)dx

hf, gi =
x

p
and the norm kf k = hf, f i. We use f g or f ·g as shorthand for the function x 7→ f (x)g(x). Any function
f ∈ L2 (Rd ) has a Fourier transform fˆ ∈ L2 (Rd ), which for absolutely integrable functions can be defined
as
Z
fˆ(w) = exp (−2πi hx, wi) f (x)dx,
(1)
where exp(iz) = cos(z) + i · sin(z), i being the imaginary unit. In the proofs, we will use the following
well-known properties of the Fourier transform:
\
• Linearity: For scalars a, b and functions f, g, af
+ bg = afˆ + bĝ.
• Isometry: hf, gi = hfˆ, ĝi and kf k = kfˆk.
• Convolution: fcg = fˆ ∗ ĝ, where ∗ denotes the convolution operation: (f ∗ g)(w) =
z) dz.

3

R

f (z) · g(w −

Natural Target Functions

P
In this section, we consider simple target functions of the form x 7→ [ ni=1 [hwi , xi]+ ][0,1] , where [z]+ =
max{0, z} is the ReLU function, and [z][0,1] = min{1, max{0, z}} is the clipping operation on the interval
[0, 1]. This corresponds to depth-2 networks with no bias in the first layer, and where the outputs of the
first layer are simply summed and moved through a clipping non-linearity (this operation can also be easily
implemented using a second layer composed of two ReLU neurons). Letting W = [w1 , . . . , wn ], we can
write such predictors as x 7→ h(W > x) for an appropriate fixed function h. Our goal would be to show that
for such a target function, with virtually any choice of W (essentially, as long as its columns are linearly
independent), and any polynomial-time learning algorithm satisfying some conditions, there exists an input
distribution on which it must fail.
As the careful reader may have noticed, it is impossible to provide such a target-function-specific result
which holds for any algorithm. Indeed, if we fix the target function in advance, we can always “learn”
by returning the target function, regardless of the training data. Thus, imposing some constraints on the
algorithm is necessary. Specifically, we will consider algorithms which exhibit certain natural invariances
to the coordinate system used. One very natural invariance is with respect to orthogonal transformations:
For example, if we rotate the input instances xi in a fixed manner, then an orthogonally-invariant algorithm
will return a predictor which still makes the same predictions on those instances. Formally, this invariance
is defined as follows:
d
Definition 1. Let A be an algorithm which inputs a dataset ({xi , yi })m
i=1 (where xi ∈ R ) and outputs a
>
predictor x 7→ f (W x) (for some function f and matrix W dependent on the dataset). We say that A is
orthogonally-invariant, if for any orthogonal matrix M ∈ Rd×d , if we feed the algorithm with {M xi , yi }m
i=1 ,
>
the algorithm returns a predictor x 7→ f (WM x), where f is the same as before and WM is such that
> M x = W > x for all x .
WM
i
i
i

Remark 1. The definition as stated refers to deterministic algorithms. For stochastic algorithms, we will
understand orthogonal invariance to mean orthogonal invariance conditioned on any realization of the
algorithm’s random coin flips.
5

For example, standard gradient and stochastic gradient descent methods (possibly with coordinateoblivious regularization, such as L2 regularization) can be easily shown to be orthogonally-invariant1 . However, for our results we will need to make a somewhat stronger invariance assumption, namely invariance
to general invertible linear transformations of the data (not necessarily just orthogonal). This is formally
defined as follows:
Definition 2. An algorithm A is linearly-invariant, if it satisfies Definition 1 for any invertible matrix M ∈
Rd×d (rather than just orthogonal ones).
One well-known example of such an algorithm (which is also invariant to affine transformations) is the
Newton method [4]. More relevant to our purposes, linear invariance occurs whenever an orthogonallyinvariant algorithm preconditions or “whitens” the data so that its covariance has a fixed structure (e.g. the
identity matrix, possibly after a dimensionality reduction if the data is rank-deficient). For example, even
though gradient descent methods are not linearly invariant, they become so if we precede them by such a
preconditioning step. This is formalized in the following theorem:
−1 >
Theorem 1. Let A be any algorithm which given {xi , yi }m
i=1 , computes the whitening matrix P = D U
>
2
(where X = [x1 x2 . . . xm ], X = U DV is a thin SVD decomposition of X), feeds {P xi , yi }m
i=1 to
>
an orthogonally-invariant algorithm, and given the output predictor x 7→ f (W x), returns the predictor
x 7→ f ((P > W )> x). Then A is linearly-invariant.

It is easily verified that the covariance matrix of the transformed instances P x1 , . . . , P xm is the r × r
identity matrix (where r = Rank(X)), so this is indeed a whitening transform. We note that whitening is a
very common preprocessing heuristic, and even when not done explicitly, scalable approximate whitening
and preconditioning methods (such as Adagrad [8] and batch normalization [12]) are very common and
widely recognized as useful for training neural networks.
To show our result, we rely on a reduction from a PAC-learning problem known to be computationally hard, namely learning intersections of halfspaces. These are Boolean predictors parameterized by
w1 , . . . , wn ∈ Rd and b1 , . . . , bn ∈ R, which compute a mapping of the form
x →

n
^

(hwi , xi ≥ bi )

i=1

(where we let 1 correspond to ‘true’ and 0 to ‘false’). The problem of PAC-learning intersections of halfspaces over the Boolean cube (x ∈ {0, 1}d ) has been well-studied. In particular, two known hardness results
are the following:
• Klivans and Sherstov [16] show that under a certain well-studied cryptographic assumption (hardness
of finding unique shortest vectors in a high-dimensional lattice), no algorithm can PAC-learn intersection of nd = dδ halfspaces (where δ is any positive constant), even if the coordinates of wi and bi are
all integers, and maxi k(wi , bi )k ≤ poly(d).
• Daniely and Shalev-Shwartz [7] show that under an assumption related to the hardness of refuting
random K-SAT formulas, no algorithm can PAC-learn intersections of nd = ω(1) halfspaces (as
d → ∞), even if the coordinates of wi and bi are all integers, and maxi k(wi , bi )k ≤ O(d).
Essentially, this is because the gradient of any function g(W > x) = g(hw1 , xi , . . . , hwk , xi) w.r.t. any wi is proportional
to x. Thus, if we multiply x by an orthogonal M , the gradient also gets multiplied by M . Since M > M = I, the inner products
of instances x and gradients remain the same. Therefore, by induction, it can be shown that any algorithm which operates by
incrementally updating some iterate by linear combinations of gradients will be rotationally invariant.
2
That is, if X is of size d×m, then U is of size d×Rank(X), D is of size Rank(X)×Rank(X), and V is of size m×Rank(X).
1

6

In the theorem below, we will use the result of [7], which applies to an intersection of a smaller number of
halfspaces, and with smaller norms. However, similar results can be shown using [16], at the cost of worse
polynomial dependencies on d.
The main result of this section is the following:
P d
Theorem 2. Consider any network h(W?> x) = [ ni=1
[hwi? , xi]+ ][0,1] (where W? = [w1? , . . . , wn? ]), which
satisfies the following:
• nd ≥ ω(1) as d → ∞
• maxi kwi? k ≤ O(d)
• w1? . . . wn? are linearly independent, so the smallest singular value smin (W? ) of W? is strictly positive.
Then under the assumption stated in [7], there is no linearly-invariant
algorithm which for any  > 0 and
√
O(d dnd )
any distribution D over vectors of norm at most min{1,smin (W? )} , given only access to samples (x, h(W?> x))
where x ∼ D, runs in time poly(d, 1/) and returns with high probability a predictor x 7→ f (W > x) such
that

 
Ex∼D

f (W > x) − h(W?> x)

2

≤ .

Note that the result holds even if the returned predictor f (W > x) has a different structure than hW? (·),
and W is of a larger size than W? . Thus, it applies even if the algorithm is allowed to train a larger network
or more complicated predictor than hW? (·).
The proof (which is provided in Sec. 5) can be sketched as follows: First, the hardness assumption
for learning intersection of halfspaces is shown to imply hardness of learning networks x 7→ h(W > x) as
described above (and even if W has linearly independent columns – a restriction which will be important
later). However, this only implies that no algorithm can learn x 7→ h(W > x) for all W and all input
distributions D. In contrast, we want to show that learning would be difficult even for some fixed W? . To
do so, we show that if an algorithm is linearly invariant, then the ability to learn with respect to some W
and all distributions D means that we can learn with respect to all W and all D. Roughly speaking, we
argue that for linearly-invariant algorithms, “average-case” and “worst-case” hardness are the same here.
Intuitively, this is because given some arbitrary W, D, we can create a different input distribution D̃, so that
W, D̃ “look like” W? , D under some linear transformation (see Figure 1 for an illustration). Therefore, a
linearly-invariant algorithm which succeeds on one will also succeed on the other.
A bit more formally, let us fix some W? (with linearly independent columns), and suppose we have a
linearly-invariant algorithm which can successfully learn x 7→ h(W?> x) with respect to any input distribution. Let W, D be some other matrix and distribution with respect to which we wish to learn (where W has
full column rank and is of the same size as W? ). Then it can be shown that there is an invertible matrix M
such that W = M > W? . Since the algorithm successfully learns x 7→ h(W?> x) with respect to any input
distribution, it would also successfully learn if we use the input distribution D̃ defined by sampling x ∼ D
and returning M x. This means that the algorithm would succesfully learn from data distributed as
(x, h(W?> x) , x ∼ D̃ ⇐⇒ (M x, h(W?> (M x))) , x ∼ D ⇐⇒ (M x, h(W > x)) , x ∼ D.
Since the algorithm is linearly-invariant, it can be shown that this implies successful learning from (x, h(W > x))
where x ∼ D, as required.
In the sketch above, we have ignored some technical issues. For example, we need to be careful that M
has a bounded spectral norm, so that it induces a linear transformation which does not distort norms by too
7

1

1

0.5

0.5

0

0

−0.5

−0.5

−1

−1
−1

−0.5

0

0.5

1

−1

−0.5

0

0.5

1

Figure 1: Correspondence between W? , D (left figure) and W, D̃ (right figure). Arrows correspond to
columns of W? and W , and dots correspond to the support of D and D̃. D̃ is constructed so that the
same linear transformation mapping W? to W also maps D to D̃.
much (as all our arguments apply for input distributions supported on a bounded domain). A second issue
is that if we apply a linearly-invariant algorithm on a dataset transformed by M , then the invariance is only
with respect to the data, not necessarily with respect to new instances x sampled from the same distribution
(and this restriction is necessary for results such as Thm. 1 to hold without further assumptions). However,
it can be shown that if the dataset is large enough, invariance will still occur with high probability over the
sampling of x, which is sufficient for our purposes.

4

Natural Input Distributions

In this section, we consider the difficulty of gradient-based methods to learn certain target functions, even
with respect to smooth, well-behaved distributions over Rd . Specifically, we will consider functions of the
form x 7→ ψ(hw? , xi), where w? is a vector of bounded norm and ψ is a periodic function. Note that if ψ is
continuous and piecewise linear, then ψ(hw? , xi) can be implemented by a depth-2 neural ReLU network on
any bounded subset of the domain. More generally, any continuous periodic function can be approximated
arbitrarily well by such networks.
Our formal results rely on Fourier analysis and are a bit technical. Hence, we precede them with an
informal description, outlining the main ideas and techniques, and presenting a specific case study which
may be of independent interest (Subsection 4.1). The formal results are presented in Subsection 4.2.

4.1

Informal Description of Results and Techniques

Consider a target function of the form x 7→ ψ(hw? , xi), and any input distribution whose density function
can be written as the square ϕ2 of some function ϕ (the reason for this will become apparent shortly).
Suppose we attempt to learn this target function (with respect to the squared loss) using some hypothesis
class, which can be parameterized by a bounded-norm vector w in some subset W of an Euclidean space (not
necessarily of the same dimensionality as w? ), so each predictor in the class can be written as x 7→ f (w, x)
for some fixed mapping f . Thus, our goal is essentially to solve the stochastic optimization problem
h
i
min Ex∼ϕ2 (f (w, x) − ψ(hw? , xi))2 .
(2)
w:w∈W

In this section, we study the geometry of this objective function, and show that under mild conditions on f ,
and assuming the norm of w? is reasonably large, the following holds:
8

Figure 2: Graphical depiction of the objective function in Eq. (3), in 2 dimensions and where w? = (2, 2).
• For any fixed w, the value of the objective function is almost independent of w? , in the sense that
if we pick the direction of w? uniformly at random, the value is extremely concentrated around a
fixed value independent of w? (e.g. exponentially small in kw? k2 for a Gaussian or a mixture of
Gaussians).
• Similarly, the gradient of the objective function with respect to w is almost independent of w? , and is
extremely concentrated around a fixed value (again, exponentially small in kw? k2 for, say, a mixture
of Gaussians).
Therefore, assuming kw? k is reasonably large, any standard gradient-based method will follow a trajectory
nearly independent of w? . In fact, in practice we do not even have access to exact gradients of Eq. (2), but
only to noisy and biased versions of it (e.g. if we perform stochastic gradient descent, and certainly if we
use finite-precision computations). In that case, the noise will completely obliterate the exponentially small
signal about w? in the gradients, and will make the trajectory essentially independent of w? . As a result,
assuming ψ and the distribution is such that the function ψ(hw? , xi) is sensitive to the direction of w? , it
follows that these methods will fail to optimize Eq. (2) successfully. Finally, we note that in practice, it is
common to solve not Eq. (2) directly, but rather its empirical approximation with respect to some fixed finite
training set. Still, by concentration of measure, this empirical objective would converge to the one in Eq. (2)
given enough data, so the same issues will occur.
An important feature of our results is that they make virtually no structural assumptions on the predictors
x 7→ f (w, x). In particular, they can represent arbitrary classes of neural networks (as well as other predictor
classes). Thus, our results imply that target functions of the form x 7→ ψ(hw? , xi), where ψ is periodic,
would be difficult to learn using gradient-based methods, even if we allow improper learning and consider
predictor classes of a different structure.
To explain how such results are attained, let us study a concrete special case (not necessarily in the
context of neural networks). Consider the target function x 7→ cos(2π hw? , xi), and the hypothesis class
(parameterized by w) of functions x 7→ cos(2π hw, xi). Thus, Eq. (2) takes the form
h
i
min Ex∼ϕ2 (cos(2π hw, xi) − cos(2π hw? , xi)2 .
(3)
w

9

Furthermore, suppose the input distribution ϕ2 is a standard Gaussian on Rd . In two dimensions and for
w? = (2, 2), the objective function in Eq. (2) turns out to have the form illustrated in Figure 2. This
objective function has only three critical points: A global maximum at 0, and two global minima at w? and
−w? . Nevertheless, it would be difficult to optimize using gradient-based methods, since it is extremely flat
everywhere except close to the critical points. As we will see shortly, the same phenomenon occurs in higher
dimensions. In high dimensions, if the direction of w? is chosen randomly, we will be overwhelmingly likely
to initialize far from the global minima, and hence will start in a flat plateau in which most gradient-based
methods will stall3 .
We now turn to explain why Eq. (3) has the form shown in Figure 2. This will also help to illustrate our
proof techniques, which apply much more generally. The main idea is to analyze the Fourier transform of
Eq. (3). Letting cosw denote the function x 7→ cos(2π hw, xi), we can write Eq. (3) as
Z
(cos(2π hw, xi) − cos(2π hw? , xi)2 ϕ2 (x)dx = kcosw ·ϕ − cosw? ·ϕk2 ,
where k·k is the standard norm over the space L2 (Rd ) of square integrable functions. By standard properties
of the Fourier transform (as described in Sec. 2), this squared norm of a function equals the squared norm
of the function’s Fourier transform, which equals in turn
cos
\
\
w ·ϕ − cos
w ·ϕ

2

c[
osw ∗ ϕ̂ − c\
osw? ∗ ϕ̂

=

2

.

c[
osw (ξ) can be shown to equal 21 (δ(ξ − w) + δ(ξ + w)),Rwhere δ(·) is Dirac’s delta function (a “generalized” function which satisfies δ(z) = 0 for all z 6= 0, and δ(z)dz = 1). Plugging this into the above and
simplifying, we get
1
kϕ̂(· − w) + ϕ̂(· + w) − ϕ̂(· − w? ) − ϕ̂(· + w? )k2
4
Z
1
=
|ϕ̂(ξ − w) + ϕ̂(ξ + w) − ϕ̂(ξ − w? ) − ϕ̂(ξ + w? )|2 dξ.
4 ξ

(4)
2

If ϕ2 is a standard Gaussian, ϕ̂(ξ) can be shown to equal the Gaussian-like function (4π)d/2 a−kξk where
a = exp(4π 2 ). Plugging back, the expression above is proportional to
Z 
 
2
2
2
? 2
? 2
a−kξ−wk + a−kξ+wk − a−kξ−w k + a−kξ+w k
dξ.
(5)
ξ

The expression in each inner parenthesis can be viewed as a mixture of two Gaussian-like functions, with
centers at w, −w (or w? , −w? ). Thus, if w is far from w? , these two mixtures will have nearly disjoint
support, and Eq. (5) will have nearly the same value regardless of w – in other words, it is very flat. Since
this equation is nothing more than a re-formulation of the original objective function in Eq. (3) (up to a
constant), we get a similar behavior for Eq. (3) as well.
This behavior extends, however, much more generally than the specific objective in Eq. (3). First of all,
we can replace the standard Gaussian distribution ϕ2 by any distribution such that ϕ̂ has a localized support.
This would still imply that Eq. (4) refers to the difference of two functions with nearly disjoint support,
and the same flatness phenomenon will occur. Second, we can replace the cos function by any periodic
3

Although there are techniques to overcome flatness (e.g. by normalizing the gradient [19, 10]), in our case the normalization
factor will be huge and require extremely precise gradient information, which as discussed earlier, is unrealistic here.

10

function ψ. By properties of the Fourier transform of periodic functions, we still get localized functions in
the Fourier domain (more precisely, the Fourier transform will be localized around integer multiples of w,
up to scaling). Finally, instead of considering hypothesis classes of predictors x 7→ ψ(hw, xi) similar to the
target function, we can consider quite arbitrary mappings x 7→ f (w, x). Even though this function may no
longer be localized in the Fourier domain, it is enough that only the target function x 7→ ψ(hw? , xi) will
be localized: That implies that regardless how f looks like, under a random choice of w? , only a minuscule
portion of the L2 mass of f overlaps with the target function, hence getting sufficient signal on w? will be
difficult.
As mentioned in the introduction, these techniques and observations have some close resemblance to
hardness results for learning parities over the Boolean cube in the statistical queries learning model [3].
There as well, one considers a Fourier transform (but on the Boolean cube rather than Euclidean space), and
essentially show that functions with a “localized” Fourier transform are difficult to “detect” using any fixed
function. However, our results are different and more general, in the sense that they apply to generic smooth
distributions over Euclidean space, and to a general class of periodic functions, rather than just parities. On
the flip side, our results are constrained to methods which are based on gradients of the objective, whereas
the statistical queries framework is more general and considers algorithms which are based on computing
(approximate) expectations of arbitrary functions of the data. Extending our results to this generality is an
interesting topic for future research.

4.2

Formal Results

We now turn to provide a more formal statement of our results. The distributions we will consider consist of
arbitrary mixtures of densities, whose square roots have rapidly decaying tails in the Fourier domain. More
precisely, we have the following definition:
Definition 3. Let (r) be some function from [0, ∞) to [0, 1]. A function ϕ2 : Rd → R is (r) Fourierconcentrated if its square root ϕ belongs to L2 (Rd ), and satisfies
kϕ̂ · 1≥r k ≤ kϕ̂k (r),
where 1≥r is the indicator function of {x : kxk ≥ r}.
A canonical example is Gaussian distributions: Given a (non-degenerate, zero-mean) Gaussian density
function ϕ2 with covariance matrix Σ, its square root ϕ is proportional to a Gaussian with covariance 2Σ, and
its Fourier transform ϕ̂ is well-known to be proportional to a Gaussian with covariance (2Σ)−1 . By standard
Gaussian concentration results, it follows that ϕ2 is Fourier-concentrated with (r) = exp(−Ω(λmin r2 ))
where λmin is the minimal eigenvalue of Σ. A similar bound can be shown when the Gaussian has some
arbitrary mean. More generally, it is well-known that smooth functions (differentiable to sufficiently high
order with integrable derivatives) have Fourier transforms with rapidly decaying tails. For example, if we
consider the broad class of Schwartz functions (characterized by having values and all derivatives decaying
faster than polynomially in r), then the Fourier transform of any such function is also a Schwartz function,
which implies super-polynomial decay of (r) (see for instance [11], Chapter 11 and Proposition 11.25).
We now formally state our main result for this section. We consider any predictor of the form x 7→
f (w, x), where f is some fixed function and w is a parameter vector coming from some domain W, which
we will assume w.l.o.g. to be a subset of some Euclidean space4 (for example, f can represent a network
4

More generally, our analysis is applicable to any separable Hilbert space.

11

of a given architecture, with weights specified by w). When learning f based on data coming from an
underlying distribution, we are essentially attempting to solve the optimization problem
h
i
min Fw? (w) := Ex∼ϕ2 (f (w, x) − ψ(hw? , xi))2 .
w∈W

Assume that F is differentiable w.r.t. w, any gradient-based method to solve this problem proceeds by
computing (or approximating) ∇Fw? (w) at various points w. However, the following theorem shows that
at any w, and regardless of the type of predictor or network one is attempting to train, the gradient at w is
virtually independent of the underlying target function, and hence provides very little signal:
Theorem 3. Suppose that
• ψ : R → [−1, +1] is a periodic function of period 1, which has bounded variation on every finite
interval.
P
• ϕ2 is a density function on Rd , which can be written as a (possibly infinite) mixture ϕ2 = i αi ϕ2i ,
where each ϕi is an (r) Fourier-concentrated density function.
• At some fixed w, Ex∼ϕ2

2
∂
∂w f (w, x)

≤ Gw for some Gw .

Then for some universal positive constants c1 , c2 , c3 , if d ≥ c1 , and w? ∈ Rd is a vector of norm 2r chosen
uniformly at random, then
!
∞
X
2
Varw? (∇Fw? (w)) := Ew? k∇Fw? (w) − Ew? [∇Fw? (w)]k ≤ c2 Gw exp(−c3 d) +
(nr) .
n=1

We note that bounded variation is weaker than, say, Lipschitz continuity. Assuming (r) decays rapidly
with r – say, exponentially in r2 as is the case for a Gaussian mixture – we get that the bound in the theorem
is on the order of exp(−Ω(min{d, r2 })).
Overall, the theorem implies that if r, d are moderately large, the gradient of Fw? at any point w is
extremely concentrated around a fixed value, independent of w? . This implies that gradient-based methods,
which attempt to optimize Fw? via gradient information, are unlikely to succeed. One way to formalize
this is to consider any iterative algorithm (possibly randomized), which relies on an ε-approximate gradient
oracle to optimize Fw? : At every iteration t, the algorithm chooses a point wt ∈ W, and receives a vector gt
such that |∇Fw? − gt | ≤ ε. In our case, we will be interested in ε such that ε3 is on the order of the bound
in Thm. 3. Since the bound is extremely small for moderate d, r (say, smaller than machine precision), this
is a realistic model of gradient-based methods on finite-precision machines, even if one attempts to compute
the gradients accurately. The following theorem implies that if the number of iterations is not extremely
large (on the order of 1/ε, e.g. exp(Ω(d, r2 )) iterations for Gaussian mixtures), then with high probability,
a gradient-based method will return the same predictor independent of w? . However, since the objective
function Fw? is highly sensitive to the choice of w? , this means that no such gradient-based method can
train a reasonable predictor.
p
P
Theorem 4. Assume the conditions of Thm. 3, and let ε = 3 c2 (supw∈W Gw ) (exp(−c3 d) + ∞
n=1 (nr))
be the cube root of the bound specified there (uniformly over all w ∈ W). Then for any algorithm as above
and any p ∈ (0, 1), conditioned on an event which holds with probability 1 − p over the choice of w? , its
output after at most p/ε iterations will be independent of w? .

12

5
5.1

Proofs
Proof of Thm. 1

Let PM denote the whitening matrix employed if we transform the instances X by some invertible d × d
matrix M (that is, X becomes M X), and P the whitening matrix employed for the original data.
>,
Using the same notation as in the theorem, it is easily verified that P X = V > , and PM M X = VM
>
>
>
where UM DM VM is an SVD decomposition of the matrix M X. Since both V and VM are Rank(X) × m
matrices with rows consisting of orthonormal vectors, they are related by an orthogonal transformation (i.e.
> ). Therefore, R P X = P M X. Since the data
there is an orthogonal matrix RM such that RM V > = VM
M
M
> P M X = W > P X. This in turn
is fed to an orthogonally-invariant algorithm, its output WM satisfies WM
M
> R P X = W > P X, and hence W > R V > = W > V > . Multiplying both sides on the right by
implies WM
M
M M
>W
V and taking a transpose, we get that RM
M = W , and hence WM = RM W . In words, W and WM are
the same up to an orthogonal transformation RM depending on M . Therefore,
>
>
>
(PM
WM )> M X = WM
PM M X = W > RM
RM P X = W > P X = (P > W )> X,

so we see that the returned predictor makes the same predictions over the dataset, independent of the transformation matrix M .

5.2

Proof of Thm. 2

We start with the following auxiliary theorem, which reduces the hardness result of [7] to one about neural
networks of the type we discuss here:
Theorem 5. Under the assumption stated in [7], the following holds for any nd = ω(1) (as d → ∞):
There is no P
algorithm running in time poly(d, 1/), which for any distribution D on {0, 1}d , and any
>
d
[hwi , xi]+ ) (where W = [w1 , w2 , . . . , wnd ] and maxi kwi k ≤ O(d)), given only
h(W x) = σ( ni=1
access to samples (x, h(W > x)) where x ∼ D, returns with high probability a function f such that

2 
Ex∼D f (x) − h(W > x)
≤ .
Proof. Suppose by contradiction that there exists an algorithm
and hW as
h A which for any distribution
i
2
>
described in the theorem, returns a function f such that Ex∼D f (x) − h(W x)
≤  with high probability.
In particular, let us focus on distributions D supported on {0, 1}d−1 × {1}. For these distributions,
we argue that any intersection of halfspaces on Rd−1 specified by w1 , . . . , wnd ∈ Rd−1 with integer
coordinates, and integer b1 , . . . , bn , can be specified as x 7→ 1 − h(W > x) for some function h as described in the theorem statement. To see this, note that for any wi , bi and x = (x0 , 1) in the support of D,
h(−wi , bi ), xi = − hwi , x0 i + bi is an integer, hence
!
!
nd
nd
nd
X
X
_

σ
[h(−wi , bi ), xi]+ = σ
[− wi , x0 + bi ]+
=
− wi , x0 + bi > 0
i=1

i=1

=

nd
_

wi , x

i=1
0

< bi

i=1



= ¬

nd
^
i=1

13

!
wi , x

0

≥ bi



.

Therefore, for any distribution over examples labelled by an intersection of halfspaces x 7→ 1 − h(W > x)
(with integer-valued coordinates and bounded norms), by feeding A with {xi , 1 − yi }m
i=1 , the algorithm
2
>
returns a function f , such that with high probability, Ex f (x) − h(W x) ≤ , and therefore

2
Ex (1 − f (x)) − (1 − h(W > x)) ≤ .
In particular, if we consider the Boolean function f˜(x) = 1 − rnd(f (x)), where rnd(z) = 0 if z ≤ 1/2 and
rnd(z) = 1 if z > 1/2, we argue that Prx (f˜(x) 6= 1−h(W > x)) ≤ 8. Since  is arbitrary, and 1−h(W > x)
specifies an intersection of halfspaces, this would contradict the hardness result of [7], and therefore prove
the theorem. This argument follows from the following chain of inequalities, where 1 denotes the indicator
function:


Pr f˜(x)) 6= g(x) = Pr(f (x) > 1/2 ∧ g(x) = 1) + Pr(f (x) ≤ 1/2 ∧ g(x) = 0)
= E [1 (f (x) > 1/2 ∧ g(x) = 1)] + E [1 (f (x) ≤ 1/2 ∧ g(x) = 0)]
h
i
h
i
≤ E 4 ((1 − f (x)) − g(x))2 + E 4 ((1 − f (x)) − g(x))2
h
i
≤ 8 · E ((1 − f (x)) − g(x))2 ≤ 8.

Proposition 1. Thm. 5 holds even if we restrict w1 , . . . , wnd to be linearly independent, with smin (W ) ≥ 1.
Proof. Suppose by contradiction that there exists an algorithm A which succeeds for any W as stated above.
We will describe how to use A to get an algorithm which succeeds for any W as described in Thm. 5, hence
reaching a contradiction.
Specifically, suppose we have access to samples (x, h(W > x)), where x is supported on {0, 1}d , and
where W is any matrix as described in Thm. 5. We do the following: We map every x to x̃ ∈ {0, 1}d+nd by
x̃ = (x, 0, . . . , 0), run A on the transformed samples (x̃, h(W > x)) to get some predictor f˜ : {0, 1}d+nd 7→
R, and return the predictor f (x) = f˜((x, 0, . . . , 0)).
To see why this reduction works, we note that the mapping x 7→ x̃ we have defined, where x is distributed according to x, induces a distribution D̃ on {0, 1}d+nd . Let W̃ be the (d + nd ) × nd matrix
[W ; Ind ] (that is, we add another nd × nd unit matrix below W ). We have W̃ > W̃ = W > W + Ind , so
the minimal eigenvalue of W̃ > W̃ is at least 1, hence smin (W̃ ) ≥ 1, so W̃ satisfies the conditions in the
proposition. Moreover, the norm of each column of W̃ is larger than the norm of the corresponding column
in W by at most 1, so the norm constraint in Thm. 5 still holds. Finally, W̃ x̃ = W x for all x, and therefore hW̃ (x̃) = hW (x). Thus, the distribution of (x̃, h(W > x)) = (x̃, h(W̃ > x)) (which is used to feed the
algorithm A) is a valid distribution corresponding to the conditions of the proposition and Thm. 5 (only in
dimension d + nd ≤ 2d instead of d), so A returns with high probability a predictor f˜ such that

2 
Ex̃∼D̃ f˜(x̃) − h(W̃ > x̃)
≤ .
However, f˜(x̃) = f (x), h(W̃ > x̃) = h(W > x), so the returned predictor f satisfies

2 
>
Ex∼D f (x) − h(W x)
≤ .
This contradicts Thm. 5, which states that no efficient algorithm can return such a predictor for any sufficiently large dimension d and norm bound O(d).
14

In the definitions of orthogonal invariance and linear invariance, we only required the invariance to hold
with respect to instances xi in the dataset. A stronger condition is that the invariance is satisfied for any
x ∈ Rd . However, the following lemma shows that invariance w.r.t. a dataset sampled i.i.d. from some
distribution is sufficient to imply invariance w.r.t. “nearly all” x (under the same distribution):
d
Lemma 1. Suppose the dataset {xi , yi }m
i=1 is sampled i.i.d. from some distribution (where xi ∈ R ),
then the following holds with probability at least 1 − δ for any δ ∈ (0, 1): For any invertible M and
linearly-invariant algorithm (or orthogonal M and orthogonally-invariant algorithm), conditioned on the
algorithm’s internal randomness, the returned matrices W and WM (with respect to the original data and
the data transformed by M respectively) satisfy
>
Prx (WM
M x 6= W > x) ≤

d
.
δ(m + 1)

Proof. It is enough to prove that with probability at least 1 − δ over the sampling of x1 , . . . , xm ,
Prx (x ∈
/ span(x1 , . . . , xm )) ≤

d
.
δ(m + 1)

(6)

> M x = W > x for all i means that W > M x = W > x for any x in the span of
This is because the event WM
i
i
M
x1 , . . . , xm .
Let x1 , . . . , xm+1 be sampled i.i.d. according to D. Considering probabilities over this sample, we have


m+1
m+1
X
X
Pr (xj ∈
/ span(x1 , . . . , xj−1 )) = E 
1 (xj ∈
/ span(x1 , . . . , xj−1 )) ≤ d,
(7)
j=1

j=1

where the latter inequality is because each xj is a d-dimensional vector, hence the number of times we can
get a vector not in the span of the previous ones is at most d. Moreover, since the vectors are sampled i.i.d,
we have
Pr(xj+1 ∈
/ span(x1 , . . . , xj )) ≤ Pr(xj+1 ∈
/ span(x1 , . . . , xj−1 )) = Pr(xj ∈
/ span(x1 , . . . , xj−1 )),
so the probabilities in Eq. (7) monotonically decrease with j. Thus, Eq. (7) implies
(m + 1) Pr (xm+1 ∈
/ span(x1 , . . . , xm )) ≤ d ⇒ Pr (xm+1 ∈
/ span(x1 , . . . , xm )) ≤
This is equivalent to


Ex1 ,...,xm ∼D Prxm+1 (xm+1 ∈
/ span(x1 , . . . , xm )|x1 , . . . , xm ) ≤

d
,
m+1

so by Markov’s inequality, with probability at least 1 − δ over the sampling of x1 , . . . , xm ,
Prxm+1 (xm+1 ∈
/ span(x1 , . . . , xm )) ≤

d
.
δ(m + 1)

Since xm+1 is sampled independently, Eq. (6) and hence the lemma follows.

15

d
.
m+1

With these results in hand, we can finally turn to prove Thm. 2. Suppose by contradiction that there
exists an efficient
algorithm A, which for any distribution D supported on vectors of norm

√ linearly-invariant
at most O d 2dnd / min{1, smin (W? )}, returns w.h.p. a predictor x 7→ f (W̃ > x) such that

2 
>
>
≤ .
Ex∼D? f (W̃ x) − h(W? x)
We will show that the very same algorithm, if given poly(d, 1/) samples, can successfully learn w.r.t. any
d × nd matrix W and any distribution D satisfying Proposition 1 and Thm. 5, contradicting those results.
Indeed, let W and D be an arbitrary matrix and distribution as above. We first argue that there exists a
d × d invertible matrix M such that
√
O(d 2nd )
>
W = M W? , kM k ≤
.
(8)
min{1, smin (W? )}
To see this, note that W and W? are of the same size and our conditions imply that both of them have full
column rank. Thus, we can simply augment them to invertible d×d matrices [W Ŵ ] and [W? Ŵ? ], where the
columns of Ŵ (respectively Ŵ? ) are an orthonormal basis for the subspace orthogonal to the column space
of W (respectively W? ), and choosing M > = [W Ŵ ][W? Ŵ? ]−1 . Thus, kM k ≤ [W Ŵ ] · [W? Ŵ? ]−1 .
The spectral norm of [W Ŵ ] can be bounded by the Frobenius norm, whichpby the assumption on W from
2
Thm.
p 5 and the fact√that Ŵ consist of an orthonormal basis, is at most O(d) · nd + 1 · (d − nd ) =
O( d2 nd ) = O(d nd ). The spectral norm of [W? Ŵ? ]−1 can be bounded by 1/smin ([W? Ŵ? ]), where
smin ([W? Ŵ? ]) equals the square root of the smallest eigenvalue of [W? Ŵ? ]> [W? Ŵ? ], which can be easily
verified to be min{1, smin (W? )}.
Now, consider the following thought experiment: Suppose we would run the algorithm A with samples
xi is sampled from D (which by the assumptions of Thm. 5,
(M xi , h(W?> (M xi )), i = 1, 2, . . . , m, where
√
is supported on vectors of√norm at most
d).
By
Eq. (8), M xi is supported on the set of vectors of norm at
√
most kM k kxi k ≤ kM k d ≤ O(d dnd )/ min{1, smin (W? )}, and the outputs correspond to the network
specified by W? . Therefore, by assumption, the algorithm A would return w.h.p. a matrix W̃M such that

2 
>
Ex∼D f (W̃M
(M x)) − h(W?> (M x))
≤ .
By Eq. (8), W?> M = (M > W? )> = W > , so this is equivalent to

2 
>
>
≤ .
Ex∼D f ((W̃M M x) − h(W x))

(9)

Let W̃I be the matrix returned by A if we had fed it with the samples {(xi , h(W > xi ))}m
i=1 (or equivalently,
>
m
5
{(xi , h(W? (M xi )))}i=1 ) . Let Ex be the event (conditioned on the samples used by the algorithm) that
> M x = W̃ > x. By Lemma 1, w.h.p. over the samples fed to the
a freshly sampled x ∼ D satisfies W̃M
I
5

Note that if the algorithm is stochastic, both W̃M and W̃I are not fixed given the data, but also depend on the algorithm’s
internal randomness. However, the proof will still follow by conditioning on any possible realization of this randomness.

16

algorithm, Prx (Ex ) ≥ 1 − O(d/m). Therefore, w.h.p. over the samples x1 , . . . , xm ,

2 
>
>
Ex∼D f (W̃I x) − h(W x)




2
2
>
>
>
>
= Pr(Ex ) · Ex f (W̃I x) − h(W x) Ex + Pr(¬Ex ) · Ex f (W̃I x) − h(W x) ¬Ex
x
x


2
>
≤ Pr(Ex ) · Ex f (W̃M
M x) − h(W > x) Ex + Pr(¬Ex ) · 1
x
x


 
2
d
>
= Ex f (W̃M
M x) − h(W > x) 1(Ex ) + O
m

 
 
2 
d
d
>
>
≤ Ex f (W̃M M x) − h(W x)
+O
≤ +O
,
m
m
where we used the fact that h maps x to [0, 1], a union bound and Eq. (9). Now, recall that W̃I refers to
the output of the algorithm, given samples {(xi , h(W > xi ))}m
i=1 where m = poly(d, 1/). Thus, we have
shown that w.h.p., as long as the algorithm is fed with m ≥ d/ samples6 , the algorithm returns W̃I which
satisfies

2 
>
>
Ex∼D f (W̃I x) − h(W x)
≤ O().
This means that the algorithm succesfully learns the hypothesis x 7→ h(W > x) with respect to the distribution D. Since  is arbitrarily small and W, D were chosen arbitrarily, the result follows.

5.3

Proof of Thm. 3

P
We note that for any function q, Ex∼ϕ2 [q(x)] =
i αi · Ex∼ϕ2i [q(x)]. Thus, it is enough to prove the
2
bound in the theorem
ϕ consists of a single element whose square root is Fourier-concentrated. For
P when
2 , the result follows by applying the bound for each ϕ2 individually, and using
a mixture ϕ2 =
α
ϕ
i
i
i
i
Jensen’s inequality.
To simplify notation a bit, we let ψw (·) stand for the function ψ(hw, ·i), and let h(· − v) (where v is
some vector and h is a function on Rd ) stand for the function x 7→ h(x − v). Also, we will use several times
the fact that for any two L2 (Rd ) functions h1 , h2 ,
Z
Z
hh1 (· − v), h2 (· − v)i =
h1 (x − v)h2 (x − v)dx =
h1 (x)h2 (x)dx = hh1 , h2 i .
In other words, inner products (and hence also norms) in L2 (Rd ) are invariant to a shift in coordinates.
The proof is a combination of a few lemmas, presented below.
Lemma 2. For any w, it holds that ψw ϕ ∈ L2 (Rd ), and satisfies
X
[
ψ
az · ϕ̂(x − zw)
w ϕ(x) =
z∈Z

for any x, where Z is the set of integers and az are complex-valuedP
coefficients (corresponding to the Fourier
series expansion of ψ, hence depending only on ψ) which satisfy z∈Z |az |2 ≤ 1.
6

Even if the algorithm does not require that many samples, we can still artificially add more samples – these are merely used to
ensure that its linear invariance is with respect to a sufficiently large dataset.

17

Proof. First, we note that ψw ϕ ∈ L2 (Rd ), since both ψw and ϕ are locally integrable by the theorem’s
conditions, and satisfy
Z
Z
Z
2
2
2
2
2
kψw ϕk = ψw (x)ϕ (x)dx = ψ (hw, xi)ϕ (x)dx ≤ ϕ2 (x)dx = 1 < ∞.
2
d
[
As a result, ψ
w ϕ exists as a function in L (R ). Since ψ is a bounded variation function, it is equal
everywhere to its Fourier series expansion:
X
ψ(x) =
az exp (2πizx) ,
z∈Z

where i is the imaginary unit (note that since ψ is real-valued, the imaginary components eventually cancel
out, but it will be more convenient for us to represent the Fourier series in this compact form). By Parseval’s
R 1/2
P
identity, z |az |2 = −1/2 ψ 2 (x)dx, which is at most 1 (since ψ(x) ∈ [−1, +1]).
Based on this equation, we have
X
ψw (x) = ψ(hw, xi) =
az exp (2πiz hw, xi) .
z∈Z

We now wish to compute the Fourier transform of the above7 . First, we note that the Fourier transform of
exp(2πi hv, ·i) is given by δ(· − v), where δR is the Dirac delta function (a so-called generalized function
which satisfies δ(x) = 0 for all x 6= 0, and δ(x)dx = 1). Based on this and the linearity of the Fourier
transform, we have that
X
ψ̂w (x) =
az · δ (x − zw) ,
z∈Z

and therefore, by the convolution property of the Fourier transform, we have
Z


[
ψ
ψ̂w (z) · ϕ̂(x − z)dz
w ϕ(x) = ψ̂w ∗ ϕ̂ (x) =
Z
X
=
az · δ(z − zw) · ϕ̂(x − z)dz
z∈Z

=

X

az · ϕ̂(x − zw)

z∈Z

as required.
Lemma 3. For any distinct integers z1 6= z2 and any w such that kwk = 2r, it holds that
h |ϕ̂(· − z1 w)| , |ϕ̂(· − z2 w)| i ≤ 2 · (|z1 − z2 |r).
Proof. Let ∆ = |z2 − z1 |r, and v = (z2 − z1 )w, so v is a vector of norm 2∆. Since the inner product is
invariant to shifting the coordinates, we can assume without loss of generality that z1 = 0, and our goal is
to bound h|ϕ̂|, |ϕ̂(· − v)|i.
7

Strictly speaking, this function does not have a Fourier transform in the sense of Eq. (1), since the associated integrals do not
converge. However, the function still has a well-defined Fourier transform in the more general sense of a generalized function or
distribution (see e.g. [11] for a survey). In the derivation below, we will simply rely on some standard formulas from the Fourier
analysis literature, and refer to [11] for their formal justifications.

18

Using the convention that 1≤∆ is the indicator of {x : kxk ≤ ∆}, and 1>∆ is the indicator for its
complement, we have
h|ϕ̂| , |ϕ̂(· − v)|i = h|ϕ̂| , |ϕ̂(· − v)| 1≤∆ i + h|ϕ̂| , |ϕ̂(· − v)| 1>∆ i
= h|ϕ̂| , |ϕ̂(· − v)| 1≤∆ i + h|ϕ̂| 1>∆ , |ϕ̂(· − v)|i
≤ kϕ̂k kϕ̂(· − v)1≤∆ k + kϕ̂1>∆ k kϕ̂(· − v)k ,
where in the last step we used Cauchy-Schwartz. Using the fact that norms and inner products are invariant
to coordinate shifting, the above is at most
sZ
sZ
!
kϕ̂k kϕ̂1≤∆ (· + v)k + kϕ̂1>∆ k kϕ̂k = kϕ̂k
|ϕ̂(x)|2 1kx+vk≤∆ dx +
|ϕ̂(x)|2 1>∆ (x)dx .
By the triangle inequality and the assumption kvk = 2∆, the event kx + vk ≤ ∆ implies kxk ≥ ∆.
Therefore, the above can be upper bounded by
sZ
2 kϕ̂k
|ϕ̂(x)|2 1≥∆ (x)dx = 2 kϕ̂k · kϕ̂ · 1≥∆ k .
2
Since ϕ is Fourier-concentrated, this is at most 2(∆) kϕ̂k2 = 2(∆)
= 2(∆), where we use the
R kϕk
2
2
isometry of the Fourier transform and the assumption that kϕk = ϕ (x)dx = 1. Plugging back the
definition of ∆, the result follows.

Lemma 4. It holds that
X

|az1 | · |az2 | · (r|z1 − z2 |) ≤ 2

z1 6=z2 ∈Z

∞
X

(nr)

n=1

Proof. For simplicity, define 0 (x) = (x) for all x > 0, and (0) = 0. Then the expression in the lemma
equals
X
|az1 | · |az2 | · 0 (|z1 − z2 |r)
z1 ,z2 ∈Z

X 

=

|az1 |



p
p
0 (|z1 − z2 |r) |az2 | 0 (|z1 − z2 |r)

z1 ,z2 ∈Z

≤

s X

|az1 |2 0 (|z1 − z2 |r)

z1 ,z2 ∈Z

X

=

s X

|az2 |2 0 (|z1 − z2 |r)

z1 ,z2 ∈Z

|az1 |2 0 (|z1 − z2 |r)

z1 ,z2 ∈Z

where in the last step we used the fact P
that the two inner square roots are the same up to a different indexing.
0
Recalling the definition of  and that z |az |2 ≤ 1, the above is at most
X
X
X
|az1 |2
0 (|z1 − z2 |r) ≤ sup
0 (|z1 − z2 |r)
z1 ∈Z

=

z2 ∈Z
∞
X

0 (0) + 2

z1 ∈Z z ∈Z
2

!
0 (nr)

= 2

n=1

∞
X
n=1

19

(nr) .

Lemma 5. For any g ∈ L2 (Rd ), if d ≥ c0 (for some universal constant c0 ), and we sample w uniformly at
random from {w : kwk = 2r}, it holds that
!
D
∞
E
2 
X
[
E
g, ψ
≤ 10 kgk2 exp(−cd) +
(nr)
w ϕ − a0 hg, ϕ̂i
n=1

where a0 is the coefficient from Lemma 2 and c is a universal positive constant.
Proof. By symmetry, given any function f of w, the expectation Ew [f (w)] (where w is uniform on a
sphere) can be equivalently written as Ew∈W EU [f (U w)] = EU Ew∈W [f (U w)], where U is a rotation
matrix chosen uniformly at random (so that for any w, U w is uniformly distributed on the sphere of radius
kwk), and Ew∈W refers to a uniform distribution of w over some finite set W of vectors of norm 2r. In
particular, we will choose W = {w1 , . . . , wdexp(cd)e } (where c is some positive universal constant) which
satisfies the following:
∀i kwi k = 2r , ∀i 6= j | hwi , wj i | < 2r2 .
(10)
The existence of such a set follows from standard concentration of measure arguments (i.e. if we pick that
n
od
2r
2r
√
many vectors uniformly at random from − √
,
+
, and c is small enough, the vectors will satisfy the
d
d
above with overwhelming probability, hencesuch a set must exist).
D
E
2 
[
Thus, our goal is to bound EU Ew∈W E
g, ψw ϕ − a0 hg, ϕ̂i
. In fact, we will prove the bound
stated in the lemma for any U , and will focus on U = I without loss of generality (the argument for other
U is exactly the same). First, by applying Lemma 2, we have
 *
+
!2 
D
E
2 
X
[
Ew∈W
g, ψ
= Ew∈W 
g,
az ϕ̂(· − zw) − a0 hg, ϕ̂i 
w ϕ − a0 hg, ϕ̂i
z∈Z

*
= Ew∈W  g,

+2 
X

az ϕ̂(· − zw)

.

(11)

z∈Z\{0}

For any w ∈ W, let
Aw = {x ∈ Rd : ∃z ∈ Z \ {0} s.t. kx − zwk < r}.
In words, each Aw corresponds to the union of open balls of radius r around ±w, ±2w, ±3w . . .. An
important property of these sets is that they are disjoint: Aw ∩ Aw0 = ∅ for any distinct w, w0 ∈ W. To see
why, note that if there was some x in both of them, it would imply kx − z1 wk < r and kx − z2 w0 k < r
for some non-zero z1 , z2 ∈ Z, hence kz1 w − z2 w0 k < 2r by the triangle inequality. Squaring both sides
and performing some simple manipulations (using the facts that kwk = kw0 k = 2r and |z1 |, |z2 | ≥ 1), we
would get


z2
z1
0
2 2
2
2 2
2
0
2
+
≥ 2r2 ,
2|z1 z2 | · w, w > 4r (z1 + z2 − 1) ≥ 2r (z1 + z2 ) ⇒
w, w ≥ r
z2
z1
where we used the fact that x + 1/x ≥ 2 for all x > 0. This contradicts the assumption on W (see Eq. (10)),
and establishes that {Aw }w∈W are indeed disjoint sets.

20

We now continue by analyzing Eq. (11). Letting 1Aw be the indicator function to the set Aw , and 1ACw
be the indicator of its complement, and recalling that (a + b)2 ≤ 2(a2 + b2 ), we can upper bound Eq. (11)
by
*
*
+2 
+2 
X
X
az ϕ̂(· − zw)  . (12)
az ϕ̂(· − zw)  + 2 · Ew∈W  g, 1ACw
2 · Ew∈W  g, 1Aw
z∈Z\{0}

z∈Z\{0}

We consider each expectation separately. Starting with the first one, we have
*
+2 
D
E2 

X
[


az ϕ̂(· − zw)
= Ew∈W g, 1Aw ψw ϕ − a0 ϕ̂
Ew∈W
g, 1Aw
z∈Z\{0}



E2 
2
2 [
[
≤ Ew∈W k1Aw gk ψw ϕ − a0 ϕ̂
= Ew∈W 1Aw g, ψw ϕ − a0 ϕ̂



h

i
2
2
2
2
2
2
2
[
=
2
·
E
k1
gk
kψ
ϕk
+
|a
|
·
k
ϕ̂k
.
ψ
ϕ
+
ka
ϕ̂k
≤ 2 · Ew∈W k1Aw gk
w
0
w
0
w∈W
Aw
D

R 2
R
P
Since we have kϕ̂k = kϕk = 1, |a0 |2 ≤ z |az |2 ≤ 1, and kψw ϕk2 = ψw
(x)ϕ2 (x)dx ≤ ϕ2 (x)dx =
1, the above is at most
Z
h
i
4 X
2
4 · Ew∈W k1Aw gk ≤
1Aw (x)|g(x)|2 dx
|W|
w∈W
!
Z
X
4
=
1Aw (x) |g(x)|2 dx.
|W|
w∈W
P
Since Aw are disjoint sets, w∈W 1Aw (x) ≤ 1 for any x, so the above is at most
Z
4
|g(x)|2 dx ≤ 4 exp(−cd) kgk2 .
(13)
|W|
D
E2 
P
.
We now turn to analyze the second expectation in Eq. (12), namely Ew∈W g, 1ACw z∈Z\{0} az ϕ̂(· − zw)
We will upper bound the expression deterministically for any w, so we may drop the expectation. Applying
Cauchy-Schwartz, it is at most


2
D
E
X
X
kgk2 · 1ACw
az ϕ̂(· − zw)
= kgk2 
az1 az2 1ACw ϕ̂(· − z1 w), ϕ̂(· − z2 w)  .
z1 ,z2 ∈Z\{0}

z∈Z\{0}

(14)
We now divide the terms in the sum above to two cases:
• If z1 = z2 , then
Z
Z
D
E
1ACw ϕ̂(· − z1 w), ϕ̂(· − z2 w) =
1ACw (x)|ϕ̂(x − z1 w)|2 dx. =
1ACw (x + z1 w)|ϕ̂(x)|2 dx,
and by definition of AC
(x + z1 w) = 1 only if kxk ≥ r.
w and the assumption z1 6= 0, we have 1AC
w
Therefore, as ϕ is Fourier-concentrated, the above is at most
Z
|ϕ̂(x)|2 dx ≤ 2 (r) · kϕ̂k2 = 2 (r) · kϕk2 = 2 (r).
x:kxk≥r

21

• If z1 6= z2 , then by Lemma 3,
E
D
1ACw ϕ̂(· − z1 w), ϕ̂(· − z2 w) ≤ h|ϕ̂(· − z1 w)|, |ϕ̂(· − z2 w)|i ≤ 2(|z1 − z2 |r).
Plugging these two cases back into Eq. (14), we get the upper bound


X
X
kgk2 
|az |2 2 (r) + 2
|az1 | · |az2 | · (|z1 − z2 |r) .
z1 6=z2 ∈Z

z∈Z\{0}

Noting that

P

z

|az |2 ≤ 1, and applying Lemma 4, the above is at most
!
∞
∞
X
X
(nr),
(nr) ≤ 5 kgk2
kgk2 2 (r) + 4
n=1

n=1

P∞

where we used the fact that 2 (r) ≤ (r) ≤ n=1 (nr). Recalling this is an upper bound on the second
expectation in Eq. (12), and that the first expectation is upper bounded as in Eq. (13), we get that Eq. (12)
(and hence the expression in the lemma statement) is at most
!
∞
X
10 kgk2 exp(−cd) +
(nr)
n=1

as required.
With these lemmas in hand, we can now turn to prove the theorem. We have that
Varw? [∇Fw? (w)] = Ew? k∇Fw? (w) − Ew? [∇Fw? (w)]k2 ≤ Ew? k∇Fw? (w) − pk2
for any vector p which is not dependent of w? (this p will be determined later). Recalling the definition of
∂
the objective function F , and letting g(x) = (g1 (x), g2 (x), . . .) = ∂w
f (w, x), the above equals
2

Ew? Ex∼ϕ2 [(f (w, x) − ψ(hw? , xi)) g(x)] − p
X
2
=
Ew? Ex∼ϕ2 [f (w, x)gi (x) − ψ(hw? , xi)gi (x)] − pi
i

=

X

Ew? (hϕgi , ϕf (w, ·)i − hϕgi , ϕψw? i − pi )2

i

Let us now choose p so that pi = hϕgi , ϕf (w, ·)i − hϕgi , a0 ϕi (note that this choice is indeed independent
of w? ). Plugging back and applying Lemma 5 (using the L2 function ϕg
ci for each i), we get
D
E
2
X
X
\
Ew? (hϕgi , ϕψw? i − hϕgi , a0 ϕi)2 =
Ew? ϕg
ci , ϕψ
ci , a0 ϕi
b
w? − hϕg
i

i

≤ 10

X
i

kϕg
ci k2

exp(−cd) +

∞
X

!
(nr) ,

n=1

and since
d
X

kϕgi k

2

i=1

=

d Z
X

gi2 (x)ϕ2 (x)dx

Z
=

kg(x)k2 ϕ2 (x) = Ex∼ϕ2 kg(x)k2 ≤ G2w ,

i=1

the theorem follows.
22

5.4

Proof of Thm. 4

We will assume w.l.o.g. that the algorithm is deterministic: If it is randomized, we can simply prove the
statement for any possible realization of its random coin flips.
We consider an oracle which given a point w, returns Ew? [∇Fw? (w)] if |∇Fw? (w)−Ew? [∇Fw? (w)]| ≤
ε, and ∇Fw? (w) otherwise. Thus, it is enough to show that with probability at least 1−p, the oracle will only
return responses of the form Ew? [∇Fw? (w)], which is clearly independent of w? . Since the algorithm’s
output can depend on w? only through the oracle responses, this will prove the theorem.
The algorithm’s first point w1 is fixed before receiving any information from the oracle, and is therefore
independent of w? . By Thm. 3, we have that Varw? (∇Fw? (w1 )) ≤ ε, which by Chebyshev’s inequality,
implies that
Pr (|∇Fw? (w1 ) − Ew? [∇Fw? (w1 )]| > ε) ≤ ε,
where the probability is over the choice of w? . Assuming the event above does not occur, the oracle returns
Ew? [∇Fw? (w)], which does not depend on the actual choice of w? . This means that the next point w2
chosen by the algorithm is fixed independent of w? . Again by Thm. 3 and Chebyshev’s inequality,
Pr (|∇Fw? (w2 ) − Ew? [∇Fw? (w2 )]| > ε) ≤ ε.
Repeating this argument and applying a union bound, it follows that as long as the number of iterations T
satisfies T ε ≤ p (or equivalently T ≤ p/ε), the oracle reveals no information whatsoever on the choice of
w? all point chosen by the algorithm (and hence also its output) are independent of w? as required.
Acknowledgements
This research is supported in part by an FP7 Marie Curie CIG grant, Israel Science Foundation grant 425/13,
and the Intel ICRI-CI Institute.

References
[1] Alexandr Andoni, Rina Panigrahy, Gregory Valiant, and Li Zhang. Learning polynomials with neural
networks. In Proceedings of the 31st International Conference on Machine Learning (ICML-14), pages
1908–1916, 2014.
[2] Sanjeev Arora, Aditya Bhaskara, Rong Ge, and Tengyu Ma. Provable bounds for learning some deep
representations. In ICML, pages 584–592, 2014.
[3] Avrim Blum, Merrick Furst, Jeffrey Jackson, Michael Kearns, Yishay Mansour, and Steven Rudich.
Weakly learning dnf and characterizing statistical query learning using fourier analysis. In Proceedings
of the twenty-sixth annual ACM symposium on Theory of computing, pages 253–262. ACM, 1994.
[4] Stephen Boyd and Lieven Vandenberghe. Convex optimization. Cambridge university press, 2004.
[5] Anna Choromanska, Mikael Henaff, Michael Mathieu, Gérard Ben Arous, and Yann LeCun. The loss
surfaces of multilayer networks. In AISTATS, 2015.
[6] Amit Daniely, Roy Frostig, and Yoram Singer. Toward deeper understanding of neural networks: The
power of initialization and a dual view on expressivity. arXiv preprint arXiv:1602.05897, 2016.

23

[7] Amit Daniely and Shai Shalev-Shwartz. Complexity theoretic limitations on learning dnf’s. In 29th
Annual Conference on Learning Theory, pages 815–830, 2016.
[8] John Duchi, Elad Hazan, and Yoram Singer. Adaptive subgradient methods for online learning and
stochastic optimization. Journal of Machine Learning Research, 12(Jul):2121–2159, 2011.
[9] Vitaly Feldman, Cristobal Guzman, and Santosh Vempala. Statistical query algorithms for stochastic
convex optimization. arXiv preprint arXiv:1512.09170, 2015.
[10] Elad Hazan, Kfir Levy, and Shai Shalev-Shwartz. Beyond convexity: Stochastic quasi-convex optimization. In Advances in Neural Information Processing Systems, pages 1594–1602, 2015.
[11] John K. Hunter and Bruno Nachtergaele. Applied analysis. World Scientific Publishing Co., Inc., River
Edge, NJ, 2001.
[12] Sergey Ioffe and Christian Szegedy. Batch normalization: Accelerating deep network training by
reducing internal covariate shift. In Proceedings of the 32nd International Conference on Machine
Learning (ICML-15), pages 448–456, 2015.
[13] Majid Janzamin, Hanie Sedghi, and Anima Anandkumar. Beating the perils of non-convexity: Guaranteed training of neural networks using tensor methods. CoRR abs/1506.08473, 2015.
[14] Michael Kearns. Efficient noise-tolerant learning from statistical queries. Journal of the ACM (JACM),
45(6):983–1006, 1998.
[15] Adam R. Klivans and Pravesh Kothari. Embedding hard learning problems into gaussian space. In
APPROX/RANDOM, 2014.
[16] Adam R Klivans and Alexander A Sherstov. Cryptographic hardness for learning intersections of
halfspaces. Journal of Computer and System Sciences, 75(1):2–12, 2009.
[17] Roi Livni, Shai Shalev-Shwartz, and Ohad Shamir. On the computational efficiency of training neural
networks. In Advances in Neural Information Processing Systems, pages 855–863, 2014.
[18] Peter McCullagh and John A Nelder. Generalized linear models, volume 37. CRC press, 1989.
[19] Yurii E Nesterov. Minimization methods for nonsmooth convex and quasiconvex functions. Matekon,
29:519–531, 1984.
[20] Itay Safran and Ohad Shamir. On the quality of the initial basin in overspecified neural networks. In
ICML, 2016.
[21] Daniel Soudry and Yair Carmon. No bad local minima: Data independent training error guarantees for
multilayer neural networks. arXiv preprint arXiv:1605.08361, 2016.
[22] Yuchen Zhang, Jason D Lee, Martin J Wainwright, and Michael I Jordan. Learning halfspaces and
neural networks with random initialization. arXiv preprint arXiv:1511.07948, 2015.

24

