Deep Semi-Random Features for Nonlinear Function Approximation
Kenji Kawaguchi,1∗ Bo Xie,2∗ Vikas Verma,3 Le Song2
1 Massachusetts

Institute of Technology
Institute of Technology
3 Aalto University

arXiv:1702.08882v7 [cs.LG] 21 Nov 2017

2 Georgia

Abstract
We propose semi-random features for nonlinear function approximation. The ﬂexibility of semi-random feature lies between the fully adjustable units in deep learning and the random features used in kernel methods. For one hidden layer
models with semi-random features, we prove with no unrealistic assumptions that the model classes contain an arbitrarily
good function as the width increases (universality), and despite non-convexity, we can ﬁnd such a good function (optimization theory) that generalizes to unseen new data (generalization bound). For deep models, with no unrealistic assumptions, we prove universal approximation ability, a lower
bound on approximation error, a partial optimization guarantee, and a generalization bound. Depending on the problems,
the generalization bound of deep semi-random features can
be exponentially better than the known bounds of deep ReLU
nets; our generalization error bound can be independent of the
depth, the number of trainable weights as well as the input dimensionality. In experiments, we show that semi-random features can match the performance of neural networks by using
slightly more units, and it outperforms random features by using signiﬁcantly fewer units. Moreover, we introduce a new
implicit ensemble method by using semi-random features.

Introduction
Many recent advances, such as human-level image classiﬁcation (Deng et al. 2009) and game playing (Bellemare et al.
2013; Silver et al. 2016) in machine learning are attributed to
large-scale nonlinear function models. There are two dominating paradigms for nonlinear modeling in machine learning: kernel methods and neural networks:
• Kernel methods employ pre-defined basis functions,
k(x, x ), called kernels to represent nonlinear functions
(Scholkopf and Smola 2001; Shawe-Taylor and Cristianini 2004). Learning algorithms that use kernel methods
often come with nice theoretical properties–globally optimal parameters can be found via convex optimization, and
statistical guarantees can be provided rigorously. However, kernel methods typically work with matrices that
are quadratic in the number of samples, leading to unfavorable computation and storage complexities. A popular
∗

Equal contribution
c 2018, Association for the Advancement of Artiﬁcial
Copyright 
Intelligence (www.aaai.org). All rights reserved.

approach to tackle such issues is to approximate kernel
functions using random features (Rahimi and Recht 2008;
Sindhwani, Avron, and Mahoney 2014; Pennington, Yu,
and Kumar 2015). One drawback of random features is
that its approximation powers suffer from the curse of
dimensionality (Barron 1993) because its bases are not
adaptive to the data.
• Neural networks use adjustable basis functions and learn
their parameters to approximate the target nonlinear function (LeCun, Bengio, and Hinton 2015). Such adaptive
nature allows neural networks to be compact yet expressive. As a result, they can be efﬁciently trained on
some of the largest datasets today. By incorporating domain speciﬁc network architectures, neural networks have
also achieved state-of-the-art results in many applications. However, learning the basis functions involves difﬁcult non-convex optimization. Few theoretical insights
are available in the literature and more research is needed
to understand the working mechanisms and theoretical
guarantees for neural networks (Choromanska, LeCun,
and Arous 2015; Swirszcz, Czarnecki, and Pascanu 2016;
Shamir 2016).
Can we have the best of both worlds? Can we develop a
framework for big nonlinear problems which has the ability to adapt basis functions, has low computational and storage complexity, while at the same time retaining some of
the theoretical properties of random features? Towards this
goal, we propose semi-random features to explore the space
of trade-off between ﬂexibility, provability and efﬁciency
in nonlinear function approximation. We show that semirandom features have a set of nice theoretical properties, like
random features, while possessing a (deep) representation
learning ability, like deep learning. More speciﬁcally:
• Despite the nonconvex learning problem, semi-random
feature model with one hidden layer has no bad local minimum;
• Depending on the problems, the generalization bound of
deep semi-random features can be exponentially better
than known bounds of deep ReLU nets;
• Semi-random features can be composed into multi-layer
architectures, and going deep in the architecture leads to
more expressive model than going wide;
• Semi-random features also lead to statistical stable func-

tion classes, where generalization bounds can be readily
provided.

Background
We brieﬂy review different ways of representing nonlinear
functions in this section.
Hand-designed basis. In a classical machine learning approach for nonlinear function approximation, users or domain experts typically handcraft a set of features φexpert :
X → H, a map from an input data space X to a (complete) inner product space H. Many empirical risk minimization algorithms then require us to compute the inner product of the features as φexpert (x), φexpert (x )H
for each pair (x, x ) ∈ X × X . Computing this inner
product can be expensive when the dimensionality of H
is large, and indeed it can be inﬁnite. For example, if
H is the space of square integrable functions, we need
to evaluate the integral as φexpert (x), φexpert (x )H =

φ
(x; ω)φexpert (x ; ω).
ω expert
Kernel methods. When our algorithms solely depend on
the inner product, the kernel trick avoids this computational
burden by introducing an easily computable kernel function
as kexpert (x , x) = φexpert (x), φexpert (x )H , resulting in
an implicit deﬁnition of the features φexpert (Scholkopf and
Smola 2001; Shawe-Taylor and Cristianini 2004). However,
the kernel approach typically scales poorly on large datasets.
Given a training set of m input points {xi }m
i=1 , evaluating a learned function at a new point x requires computing
m
fˆ(x) = i=1 αi kexpert (xi , x), the cost of which increases
linearly with m. Moreover, it usually requires computing (or
approximating) inverses of matrices of size m × m.

Random features. In order to scale to large datasets, one
can approximate the kernel by a set of random basis functions sampled according to some distributions. That is,
C
kexpert (x , x) ≈ C1 j=1 φrandom (x; rj )φrandom (x ; rj ),
where both the type of basis functions φrandom , and the
sampling distribution for the random parameter rj are determined by the kernel function. Due to its computational
advantage and theoretical foundation, the random feature approach has many applications and is an active research topic
(Rahimi and Recht 2008; Sindhwani, Avron, and Mahoney
2014; Pennington, Yu, and Kumar 2015).
Neural networks. Neural networks approximate functions
using weighted combination of adaptable basis functions
n
(2)
(1)
f (x) = k=1 wk φ(x; wk ), where both the combination
(2)
(1)
weights wk and the parameters wk in each basis function
φ are learned from data. Neural networks can be composed
into multilayers to express highly ﬂexible nonlinear functions.

Semi-Random Features
When comparing different nonlinear representaitons, we can
see that random features are designed to approximate a
known kernel, but not for learning these features from the
given dataset (i.e., it is not a representation learning). As

a result, when compared to neural network, it utilizes less
amount of information encoded in the dataset, which could
be disadvantageous. Neural networks, on the other hand,
pose a difﬁculty for theoretical developments due to nonconvexity in optimization.
This suggests a hybrid approach of random feature
and neural network, called semi-random feature (or semirandom unit), to learn representation (or feature) from
datasets. The goal is to obtain a new type of basis functions which can retain some theoretical guarantees via injected randomness (or diversity) in hidden weights, while at
the same time have the ability to adapt to the data at hand.
More concretely, semi-random features are deﬁned as


(1)
φs (x; r, w) = σs (x r) x w ,

where x = (1, x ) is assumed to be in R1+d , r =
(r0 , r  ) is sampled randomly, and w = (w0 , w  ) is
adjustable weights to be learned from data (hence, it is
“semi-random”). Furthermore, the family of functions σs for
s ∈ {0, 1, 2, . . . } is deﬁned as σs (z) = (z)s H(z), where H
is Heaviside step function (H(z) = 1 for z > 0 and 0 otherwise). For instance, σ0 is simply Heaviside step function,
σ1 is ramp function, and so on. We call the corresponding
semi-random features with s = 0 “linear semi-random features (LSR)” and with s = 1 “squared semi-random features (SSR)”. An illustration of example semi-random features can be found in Figure 4 in Appendix.
Unlike dropout, which uses a data independent random
switching mechanism (during training), the random switching in semi-random feature depends on the input data x
(during both training and testing), inducing highly nonlinear models with practical advantages. By further enhancing
this property, we additionally introduce linear semi-random
implicit-ensemble (LSR-IE) features in Section “Image classiﬁcation benchmarks”.
Intuitively, models with semi-random features have more
expressive power than those with random features because
of the learnable unit parameter w. Yet, these models are
less ﬂexible compared to neural networks, since the parameters in σs (x r) is sampled randomly. Depending on
the problems. this property of semi-random feature can result in exponential advantage over fully random feature in
expressiveness (as discussed in Section “Better than Random Feature?”), and exponential advantage over deep ReLU
models in generalization error bound (as discussed in Section “Generalization Guarantee”).

One Hidden Layer Model
With semi-random features φs in equation (1), we deﬁne one
hidden layer model for nonlinear function as
fˆns (x; w) =

n


(1)

(2)

φs (x; rk , wk )wk ,

(2)

k=1

where rk is sampled randomly for k ∈ {2, 3, · · · , n} as described in Section “Semi-Random Features”, and r1 is ﬁxed
to be the ﬁrst element of the standard basis as r1 = e1 =
(1, 0, 0, · · · 0) (to compactly represent a constant term in

x). We can think of this model as one hidden layer model
(1)
by considering φs (x; r, wk ) as the output of k-th unit of
(1)
the hidden layer, and wk as adjustable parameters associated with this hidden layer unit. This way of understanding
the model will become helpful when we generalize it to a
multilayer model in Section “Multilayer Model”. Note that
fˆns (x; w) is a nonlinear function of x. When it is clear, by the
notation w, we denote all adjustable parameters in the entire
model.
In matrix notation, the model in (2) can be rewritten as


fˆns (x; w) = σs (x R)  (x W(1) ) W (2) ,
(3)
where

(1)

(1)

W(1) = (w1 , w2 , . . . , wn(1) ) ∈ R(d+1)×n ,
R = (r1 , r2 , . . . , rn ) ∈ R(d+1)×n , and
(2)

W(2) = (w1 , . . . , wn(2) ) ∈ Rn+1 .
Here, (M1  M2 ) represents a Hadamard product of two
matrices M1 and M2 . Furthermore σs (M )ij = σs (Mij ),
given a matrix M of any size (with overloads of the symbol
σs ).
In the following subsections, we present our theoretical
results for one hidden layer model. All proofs in this paper
are deferred to the appendix.

Universal Approximation Ability
We show that our model class has universal approximation
ability. Given a ﬁnite s, our model class is deﬁned as
Fns = {x → fˆns (x; w) | w ∈ Rdw },

where dw = (d+1)n+n is the number of adjustable parameters. Let L2 (Ω) be the space of square integrable functions
on a compact set Ω ⊆ Rd . Then Theorem 1 states that we
can approximate any f ∈ L2 (Ω) arbitrarily well as we increase the number of units n. We discuss the importance of
the bias term r0 to obtain the universal approximation power
in Appendix.
Theorem 1 (Universal approximation) Let s be any fixed finite integer and let Ω = {0} be any fixed nonempty compact
subset of Rd . Then, for any f ∈ L2 (Ω), with probability one,
lim inf

n→∞ fˆ∈F s

n

f − fˆ

L2 (Ω)

= 0.

Optimization Theory
As we have conﬁrmed universal approximation ability of
our model class Fns in the previous section, we now want to
ﬁnd a good fˆ ∈ Fns via empirical loss minimization. More
speciﬁcally, given a dataset {(xi , yi )}m
i=1 , we will consider
the following optimization problem:
m
2
1 
yi − fˆns (xi ; w) .
minimize L(w) =
2m i=1
w∈Rdw

Let Y = (y1 , y2 , . . . , ym )
∈ R
and Ŷ =
(fns (x1 ; w), fns (x2 ; w), . . . , fns (xm ; w)) ∈ Rm . Given a


m

matrix M , let Pcol(M ) and Pnull(M ) be the projection matrices onto the column space and null space of M .
Our optimization problem turns out to be characterized by
the following m by nd matrix:
⎤
⎡


· · · σs (x
σs (x
1 r1 )x1
1 rn )x1
⎥
⎢
..
..
..
(4)
D=⎣
⎦,
.
.
.

σs (x
m r1 )xm

···


σs (x
m rn )xm


where σs (x
i rj )xi is a 1 × d block at (i, j)-th block entry. That is, at any global minimum, we have L(w) =
1
2
and Ŷ = Pcol(D) Y . Moreover, we can
2m Pker(DT ) Y
achieve a global minimum in polynomial time in dw based
on the following theorem.
Theorem 2 (No bad local minima and few bad critical
points) For any s and any n, the optimization problem of
L(w) has the following properties:
(i) it is non-convex (if D = 0)1 ,
(ii) every local minimum is a global minimum,
(2)

(iii) if wk = 0 for all k ∈ {1, 2, . . . , n}, every critical
point is a global minimum, and
1
Pnull(DT ) Y 2
(iv) at any global minimum, L(w) = 2m
and Ŷ = Pcol(D) Y .
Theorem 2 (optimization) together with Theorem 1 (universality) suggests that not only does our model class contain an arbitrarily good function (as n increases), but also we
can ﬁnd the best function in the model class given a dataset.
In the context of understanding the loss surface of neural networks (Kawaguchi 2016; Freeman and Bruna 2016), Theorem 2 implies that the potential problems in the loss surface
are due to the inclusion of r as an optimization variable.

Generalization Guarantee
In the previous sections, we have shown that our model class
has universal approximation ability and that we can learn the
best model given a ﬁnite dataset. A major remaining question is about the generalization property; how well can a
learned model generalize to unseen new observations? Theorem 3 bounds the generalization error; the difference between expected risk, 12 Ex (f (x) − fˆ(x; w∗ ))2 , and empirical risk, L(w). In Theorem 3, we use the following notations: σ, x = [σs (x r1 )x · · · σs (x rn )x] ∈ Rnd and
(2) (1)
(2) (1)
(2) (1)
w = (w1 w1 , w2 w2 , . . . , wn wn ) ∈ Rnd .
Theorem 3 (Generalization bound for shallow model) Let
s ≥ 0 and n > 0 be fixed. Consider any model class Fns
with w 2 ≤ CW and σ, x 2 ≤ Cσx almost surely.
Then, with probability at least 1 − δ, for any fˆ ∈ Fns ,
1
Ex (f (x) − fˆ(x; w))2 − L(w)
2

≤ (CY2 + CŶ2 )

log 1δ
C
+ 2(CY + CŶ ) √ Ŷ ,
2m
m

where CŶ = CW Cσx .
1
In the case of D = 0, L(w) is convex in a trivial way; our
model class only contains a single constant function x → 0.

By combining Theorem 2 (optimization) and Theorem 3
(generalization), we obtain the following remark.
Remark 1. (Expected risk bound) Let CW be values such
1
Pker(DT ) Y 2
that the global minimal value L(w) = 2m
is attainable in the model class (e.g., setting CW =
(D D)† D Y sufﬁces). Then, at any critical points with
(2)
wk = 0 for all k ∈ {1, 2, . . . , n} and any local minimum
such that w 2 < CW , we have
Ex (f (x) − fˆ(x; w∗ ))2 ≤

Pnull(DT ) Y 2
m

⎛

+O⎝

log

m

1
δ

⎞

⎠ , (5)

with probability at least 1 − δ. Here, O(·) notation simply
hides the constants in Theorem 3.
In the right-hand side of equation (5), the second term
goes to zero as m increases, and the ﬁrst term goes to zero
as n or d increases (because null(DT ) becomes a smaller
and smaller space as n or d increases, eventually containing
only 0). Hence, we can minimize the expected risk to zero.

Multilayer Model
We generalize one hidden layer model to H hidden layer
model by composing semi-random features in a nested fashion. More speciﬁcally, let nl be the number of units, or
width, in the l-th hidden layer for all l = 1, 2, . . . , H.
Then we will denote a model of fully-connected feedforward
semi-random networks with H hidden layers by
fˆs
(x; w) = h(H) W (H+1) ,
(6)
n1 ,...,nH

w

where for all l ∈ {2, 3, · · · , H},

(l)
(l−1)
(x) W (l) ) and
h(l)
w (x) = hr (x)  (hw

(l−1)
h(l)
(x) R(l) )
r (x) = σs (hr

is the output of the l-th semi-random hidden layer, and the
output of the l-th random switching layer respectively. Here,
(l)
(l)
(l)
W (l) = (w1 , w2 , . . . , wnl ) ∈ Rnl−1 ×nl and R(l) =
(l) (l)
(l)
(r1 , r2 , . . . , rnl ) ∈ Rnl−1 ×nl . Similarly to one hidden
(l)
layer model, rk is sampled randomly for k ∈ {2, 3, . . . , nl }
(l)
but r1 is ﬁxed to be e1 = (1, 0, 0, . . . , 0) (to compactly
write the effect of constant terms in x). The output of the ﬁrst
hidden layer is the same as that of one hidden layer model:
(1)
(1)
(1)
h(1)
) and h(1)
),
w (x) = hr (x)  (xW
r (x) = σs (xR

where the boldface notation emphasizes that we require the
bias terms at least in the ﬁrst layer. In other words, we keep
(l)
randomly updating the random switching layer hr (x), and
(l)
couple it with a linearly adjustable hidden layer hw (x)W (l)
(l+1)
to obtain the next semi-random hidden layer hw (x).
Convolutional semi-random feedforward neural networks
can be deﬁned in the same way as in equation (6) with
vector-matrix multiplication being replaced by c dimensional convolution (for some number c). In our experiments,
we will test both convolutional semi-random networks as
well as fully-connected versions. We will discuss further
generalizations of our architecture in Appendix.

Beneﬁt of Depth
We ﬁrst conﬁrm that our multilayer model class
Fns1 ,...,nH = {x → fˆns1 ,...,nH (x; w)|w ∈ Rdw }

preserves universal approximation ability.

Corollary 4 (Universal approximation with deep model)
Let s be any fixed finite integer and let Ω = {0} be any fixed
nonempty compact subset of Rd . Then, for any f ∈ L2 (Ω),
with probability one,
lim

inf

n1 ,...,nH →∞ fˆ∈F s

n1 ,...,nH

f − fˆ

L2 (Ω)

= 0.

We now know that both of one hidden layer models
and deeper models have universal approximation ability.
Then, a natural question arises: how can depth beneﬁt us?
To answer the question, note that H hidden layer model
only needs O(nH) number of parameters (by setting n =
n1 , n2 , · · · , nH ) to create around nH d paths, whereas one
hidden layer model requires O(nH ) number of parameters
to do so. Intuitively, because of this, the expressive power
would grow exponentially in depth H, if those exponential paths are not redundant to each other. The redundancy
among the paths would be minimized via randomness in the
switching and exponentially many combinations of nonlinearities σs .
We formalize this intuition by considering concrete degrees of approximation powers for our models. To do so,
we adopt a type of a degree of “smoothness” on the target functions from a previous work (Barron 1993). Consider Fourier representation of a target function as f (x) =


f˜(ω)eiω x . Deﬁne a class of smooth functions ΓC :
ω∈Rd



˜
ΓC = x → f (x) :
ω 2 |f (ω)| ≤ C .
ω∈Rd

Any f ∈ ΓC with ﬁnite C is continuously differentiable
in Rd , and the gradient of f can be written as ∇x f (x) =


iω f˜(ω)eiω x . Thus, via Plancherel theorem, we can
ω∈Rd
view C as the bound on ∇x f (x) L(Rd ) . See the previous
work (Barron 1993) for a more detailed discussion on the
properties of this function class ΓC . Theorem 5 states that a
lower bound on the approximation power gets better exponentially in depth H.
Theorem 5 (Lower bound on universal approximation
power) Let Ω = [0, 1]d . For every fixed finite integer s,
for any depth H ≥ 0, and for any set of nonzero widths
{n1 , n2 , . . . , nH },
 H −1/d

κC
sup
inf
f − fˆ L2 (Ω) ≥ 2
nl
,
s
d
f ∈ΓC fˆ∈Fn
···n
1

where κ ≥ (8πe

l=1

H

(π−1) −1

)

is a constant.

By setting n = n1 = · · · = nH , the lower bound in Theorem 5 becomes: supf ∈ΓC inf fˆ∈F s
f − fˆ L2 (Ω) ≥
κC −H/d
,
d2 n

H.

n1 ···nH

where we can easily see the beneﬁt of the depth

However, this analysis is still preliminary for the purpose
of fully understanding the beneﬁt of the depth. Our hope
here is to provide a formal statement to aid intuitions. To
help our intuitions further, we discuss about an upper bound
on universal approximation power for multilayer model in
Appendix.

Fns1 ,...,nH ,

Optimization Theory

where CŶ = CW Cσx .

Similarly to one hidden layer case, we consider the following optimization problem:

The generalization bound in Corollary 7 can be exponentially better than the known bounds of ReLU; e.g., see
(Sun et al. 2016; Neyshabur, Tomioka, and Srebro 2015;
Xie, Deng, and Xing 2015). The known generalization upper bounds of ReLU explicitly contain 2H factor, which
comes from ReLU nonlinearity and grows exponentially in
depth H. In contrast, the generalization bound in Corollary
7 does not explicitly depend on any of the depth, the number of trainable weights, and the dimensionality of the domain. Note that fn1 ,...,nH (x) = vec(w) vec(σ, x) =
vec(w) 2 vec(σ, x) 2 cos θ. Hence, even though the
dimensionality of vec(w) and vec(σ, x) grows exponentially in the depth, the norm bound CW Cσx would
stay near the norm of the output. Indeed, CW Cσx ≈
fn1 ,...,nH (x)/ cos θ. As a new related work, the generalization bounds of standard nets trained by a novel two-phase
training procedure in a recent paper (Kawaguchi, Kaelbling,
and Bengio 2017) have this desired qualitative property similarly to our generalization bounds of semi-random nets.

minimize L(H) (w) =
w∈Rdw

m
2
1 
yi − fns1 ,...,nH (x; w) .
2m i=1

Compared to one hidden layer case, our theoretical understanding of multilayer model is rather preliminary. Here,
given a function g(w1 , w2 , . . . , wn ), we say that w̄ =
(w̄1 , w̄2 , . . . , w̄n ) is a global minimum of g with respect to w1 if w̄1 is a global minimum of g̃(w1 ) =
g(w1 , w̄2 , . . . , w̄n ).
Corollary 6 (No bad local minima and few bad critical
points w.r.t. last two layers) For any s, any depth H ≥ 1,
and any set of nonzero widths {n1 , n2 , . . . , nH }, the optimization problem of L(H) (w) has the following property:
(i) every local minimum is a global minimum with respect
to (W (H) , W (H+1) ), and

1
Ex (f (x) − fˆ(x; w∗ ))2 − L(H) (w)
2

≤ (CY2 + CŶ2 )

(H+1)

= 0 for all k ∈ {1, 2, . . . , nH }, every criti(ii) if wk
cal point is a global minimum with respect to
(W (H) , W (H+1) ).
Future work is still required to investigate the theoretical
nature of the optimization problem with respect to all parameters. Some hardness results of a standard neural network
optimization come from the difﬁculty of learning activation
pattern via optimization of the variable r (Livni, ShalevShwartz, and Shamir 2014). In this sense, our optimization
problem is somewhat easier, and it would be interesting to
see if we can establish meaningful optimization theory for
semi-random model as a ﬁrst step to establish the theory for
neural networks in general.

Generalization Guarantee
The following corollary bounds the generalization error. In
the statement of the corollary, we can easily see that the generalization error goes to zero as m increases (as long as the
relevant norms are bounded). Hence, we can achieve generalization. In Corollary 7, we use the following notations:
σ, xk0 ,k1 ,...,kH = σk1 ,...,kH (x)xk0 and wk0 ,k1 ,...,kH =

H
(l)
(H+1)
. Let vec(M ) be a vectorization of
l=1 wkl−1 kl wkH
a tensor M .
Corollary 7 (Generalization bound for deep model) Let
s ≥ 0 and H ≥ 1 be fixed. Let nl > 0 be fixed for
l = 1, 2, . . . , H. Consider the model class Fns1 ,...,nH with
vec(w) 2 ≤ CW and vec(σ, x) 2 ≤ Cσx almost
surely. Then, with probability at least 1 − δ, for any fˆ ∈

log 1δ
C
+ 2(CY + CŶ ) √ Ŷ ,
2m
m

Experiments
We compare semi-random features with random features
(RF) and neural networks with ReLU on both UCI datasets
and image classiﬁcation benchmarks. We will study two
variants of semi-random features for s = 0 (LSR: linear semi-random features) and s = 1 (SSR: squared
semi-random features) in σs (·) from equation (1). Additional experimental details are presented in Appendix. The
source code of the proposed method is publicly available at:
http://github.com/zixu1986/semi-random.

Simple Test Function
We ﬁrst tested the methods with a simple sine function,
f (x) = sin(x), where we can easily understand what is
happening. Figure 1 shows the test errors with one standard deviations. As we can see, semi-random network (LSR)
performed the best. The problem of ReLU became clear
once we visualized the function learned at each iteration:
ReLU network had a difﬁculty to diversify activation units
to mimic the frequent oscillations of the sine function (i.e.,
it took long time to diversely allocate the breaking points of
its piecewise linear function). The visualizations of learned
functions at each iteration for each method are presented in
Appendix. On average, they took 54.39 (ReLU), 43.04 (random), 45.44 (semi-random) seconds. Their training errors
are presented in Appendix.

UCI datasets
We have comparisons on six large UCI datasets.The network
architecture used on this dataset is multi-layer networks with

Figure 1: Test error for a simple test function.
Table 1: Performance comparison on UCI datasets. RF for
random features, LSR and SSR for linear (s = 0) and
squared (s = 1) semi-random features respectively. mtr :
number of training data points; mte : number of test data
points; d dimension of the data.
Dataset
covtype
mtr = 522, 910
mte = 58, 102
d = 54

webspam
mtr = 280, 000
mte = 70, 000
d = 123

letter
mtr = 15, 000
mtr = 5, 000
d = 16

method err (%) Dataset

method err (%)

ReLU
RF
LSR
SSR

5.6 adult
ReLU
20.2 mtr = 32, 561
RF
5.7 mte = 16, 281 LSR
14.4 d = 123
SSR

15.0
14.9
14.8
14.9

ReLU
RF
LSR
SSR

1.1 senseit
ReLU
6.0 mtr = 78, 823
RF
1.1 mte = 19, 705 LSR
2.2 d = 100
SSR

13.6
16.0
13.9
13.3

ReLU
RF
LSR
SSR

13.1 sensor
ReLU
14.9 mtr = 48, 509
RF
6.5 mtr = 10, 000 LSR
5.6 d = 48
SSR

2.0
13.4
1.4
5.7

l = [1, 2, 4] hidden layers and k = [1, 2, 4, 8, 16] × d hidden
units per layer where d is the input data dimension.
Comparison of best performance. In Table 1, we listed
the best performance among different architectures for semirandom and fully random methods. For ReLU, we ﬁx
the number of hidden units to be d. On most datasets,
semi-random features achieve smaller errors compared with
ReLU by using more units. In addition, semi-random units
have signiﬁcant lower errors than random features.
Matching the performance of ReLU. The top row of Figure 2 demonstrates how many more units are required for
random and semi-random features to reach the test errors
of networks with ReLU. First, all three methods enjoy lower
test errors by increasing the number of hidden units. Second,
semi-random units can achieve comparable performance to
ReLU with slightly more units, around 2 to 4 times in Webspam dataset. In comparison, random features require many
more units, more than 16 times. These experiments clearly
show the beneﬁt of adaptivity in semi-random features.
Depth vs width. The bottom row of Figure 2 explores the
beneﬁt of depth. Here, “l-layer” indicates l hidden layer
model. To grow the number of total units, we can either use
more layers or more units per layer. Experiment results suggest that we can gain more in performance by going deeper.
The ability to beneﬁt from deeper architecture is an impor-

(a) Covtype dataset

(b) Webspam dataset

Figure 2: Top row: Linear semi-random features match the
performance of ReLU for two hidden layer networks in
two datasets. Bottom row: Depth vs width of linear semirandom features. Both plots show performance of semirandom units. Even the total number of units is the same,
deeper models achieve lower test error.
Table 2: Test error (in %) of different methods on three
image classiﬁcation benchmark datasets. 2×, 4× and 16×
mean the number of units used is 2 times, 4 times and 16
times of that used in neural network with ReLU respectively.
neuron type

MNIST

CIFAR10

SVHN

ReLU
RF
RF 2×
RF 4×
RF 16×
LSR
LSR 2×
LSR 4×
LSR-IE
LSR-IE 2×
LSR-IE 4×

0.70
8.80
5.71
4.10
2.69
0.97
0.78
0.71
0.59
0.47
0.54

16.3
59.2
55.8
49.8
40.7
21.4
17.4
18.7
20.0
16.8
14.9

3.9
73.9
70.5
58.4
37.1
7.6
6.9
6.4
6.9
5.9
4.8

tant feature that is not possessed by random features. The
details on how the test error changes w.r.t. the number of
layers and number of units per layer are shown in Figure
3. As we can see, on most datasets, more layers and more
units lead to smaller test errors. However, the adult dataset
is more noisy and it is easier to overﬁt. All types of neurons
perform relatively the same on this dataset, and more parameters actually lead to worse results. Furthermore, the squared
semi-random features have very similar error pattern to neural network with ReLU.

Image classiﬁcation benchmarks
We have also compared different methods on three image classiﬁcation benchmark datasets. Here we use publicly

4

0.03

1
2

# of layers

4
1

2

4

8

16

0.0105

0.144
0.136
2

4

8

16

1

# of units per layer

2

4

8

1

0.250

0.45
0.30
2

4

8

16

(a) ReLU

1

# of layers

2
4

8

16

# of units per layer

0.10
0.08
0.06

0.165

0.04

0.150
2

4

8

16

# of units per layer

1

2

4

1
2

0.012

0.180

1

# of layers

0.016

4

0.160

0.195

0.16
4

0.020

sensor

0.24

2

0.024

0.168

senseit

0.32

1

0.176

0.152

0.40

0.08

0.028

8

16

8

16

0.02

# of units per layer

(c) linear semi-random features

adult

0.20
0.19
0.18
0.17
0.16
0.15

letter

0.5

0.195

0.4

0.180

0.3

0.165

0.2
1

2

4

8

16

# of units per layer

0.10
2

4

8

16

# of units per layer

0.040
0.032
0.024

sensor

0.148
0.144

0.15
1

0.048

0.152

0.25
0.20

webspam

0.154
0.153
0.152
0.151
0.150
0.149

senseit

1

letter

0.184

covtype

2

4

2

# of layers

1

0.175
0.150
0.125
0.100
0.075

4

0.6

(b) random features
webspam

# of layers

adult

2

0.08

sensor

0.210

# of units per layer

4

covtype

1

# of units per layer

# of units per layer

0.10

0.150

senseit
0.75

1

0.12

0.156

0.225

0.60

16

0.14

0.162

letter

0.152

0.16

0.168

0.275

sensor
0.10
0.08
0.06
0.04
0.02

webspam
0.174

0.300

2

0.0120

0.160

1

# of units per layer

0.160
0.152

senseit

0.175
0.150
0.125
0.100
0.075
0.050

0.0135

adult

0.325

# of layers

0.06

0.168

covtype

1

2

0.09

0.0150

2

# of layers

0.12

letter

0.0165

# of layers

1

webspam
0.176

4

adult

0.15

4

covtype

0.140
0.136
1

2

4

8

16

# of units per layer

1

2

4

8

16

0.11
0.10
0.09
0.08
0.07
0.06

# of units per layer

(d) squared semi-random features

Figure 3: Detailed experiment results for all types of neurons and on all datasets. The heat map for each dataset shows how the
test error changes w.r.t. the number of layers and number of units per layer.
available and well-tuned neural network architectures from
tensorﬂow for the experiments. We simply replace ReLU by
random and semi-random units respectively. The results are
summarized in Table 2.
LSR-IE Unit. To improve the practical performance of LSR
unit while preserving its theoretical properties, we additionally introduced a new unit, called linear semi-random implicit-ensemble (LSR-IE) unit. Unlike an explicit ensemble,
LSR-IE trains only a single network and comparable to a
single standard network with dropout.
Recall that in LSR unit, we have a single random net (with
random weights r) and a single trainable net (with trainable weights w). In contrast, in LSR-IE unit, we increase
the number of the random nets (with random weights rk for
k = 1, . . . , K) while keeping the number of the trainable
net to be one. That is, by modifying equation (1) of semirandom unit, we implement the LSR-IE unit as follows: for
each k ∈ {1, . . . , K}, we have


φs (x; rk , w) = σs (x rk ) x w ,
(7)

with 5 × 5 ﬁlters and the number of channels is 32 and
64, respectively. Each convolution is followed by a maxpooling layer, then ﬁnally a fully-connected layer of 512
units with 0.5 dropout. Increasing the number of units for
semi-random leads to better performance. At four times
the size of the original network, semi-random feature can
achieve very close errors of 0.71%. In contrast, even when
increasing the number of units to 16 times more, random
features still cannot reach below 1%.
CIFAR10 dataset. CIFAR 10 contains internet images and
consists of 50,000 32 × 32 color images for training and
10,000 images for test. We use a convolutional neural network architecture with two convolution layers, each with
64 5 × 5 ﬁlters and followed by max-pooling. The fullyconnected layers contain 384 and 192 units. By using two
times more units, semi-random features are able to achieve
similar performance with ReLU. However, the performance
of random features lags behind by a huge margin.

and we randomly sample the index k at each step of SGD
during training. Then, during test, we average the preactivation of output units over k ∈ {1, . . . , K}. We set the
number of random nets K to be 2 in MNIST and CIFAR-10
and 4 in SVHN.

SVHN dataset. The Street View House Numbers (SVHN)
dataset contains house digits collected by Google Street
View. We use the 32 × 32 color images version and only predict the digits in the middle of the image. For training, we
combined the training set and the extra set to get a dataset
with 604,388 images. We use the same architecture as in the
CIFAR10 experiments.

MNIST dataset. MNIST is a popular dataset for recognizing handwritten digits. It contains 28 × 28 grey images,
60,000 for training and 10,000 for test. We use a convolution neural network consisting of two convolution layers,

In this paper, we proposed the method of semi-random features. For one hidden layer model, we proved that our model

Conclusion

class contains an arbitrarily good function as the width increases (universality), and we can ﬁnd such a good function (optimization theory) that generalizes to unseen new
data (generalization bound). For deep model, we proved universal approximation ability, a lower bound on approximation error, a partial optimization guarantee, and a generalization bound. Furthermore, we demonstrated the advantage of
semi-random features over fully-random features via empirical results and theoretical insights.
The idea of semi-random feature itself is more general
than what is explored in this paper, and it opens up several intersecting directions for future work. Indeed, we can
generalize any deep architecture by having an option to include semi-random units per unit level. We can also deﬁne a
more general semi-random feature as: given some nonconstant functions σ and g,


φ(x; r, w) = σ(x r)g x w ,

where x = (1, x) is assumed to be in R1+d , r = (r0 , r) is
sampled randomly, and w = (w0 , w) is adjustable weights
to be learned from data. This general formulation would lead
to a ﬂexibility to balance expressivity, generalization and
theoretical tractability.

References
Barron, A. R. 1993. Universal approximation bounds for
superpositions of a sigmoidal function. IEEE Transactions
on Information theory 39(3):930–945.
Bellemare, M. G.; Naddaf, Y.; Veness, J.; and Bowling, M.
2013. The arcade learning environment: An evaluation platform for general agents. Journal of Artificial Intelligence
Research 47:253–279.
Choromanska, A.; LeCun, Y.; and Arous, G. B. 2015. Open
problem: The landscape of the loss surfaces of multilayer
networks. In Proceedings of The 28th Conference on Learning Theory, 1756–1760.
Deng, J.; Dong, W.; Socher, R.; Li, L.-J.; Li, K.; and FeiFei, L. 2009. Imagenet: A large-scale hierarchical image database. In Computer Vision and Pattern Recognition,
2009, 248–255. IEEE.
Freeman, C. D., and Bruna, J. 2016. Topology and geometry of half-rectiﬁed network optimization. arXiv preprint
arXiv:1611.01540.
Huang, G.-B.; Chen, L.; Siew, C. K.; et al. 2006. Universal
approximation using incremental constructive feedforward
networks with random hidden nodes. IEEE Trans. Neural
Networks 17(4):879–892.
Kawaguchi, K.; Kaelbling, L. P.; and Bengio, Y.
2017. Generalization in deep learning. arXiv preprint
arXiv:1710.05468.
Kawaguchi, K. 2016. Deep learning without poor local minima. In Advances in Neural Information Processing Systems.
LeCun, Y.; Bengio, Y.; and Hinton, G. 2015. Deep learning.
Nature 521(7553):436–444.
Leshno, M.; Lin, V. Y.; Pinkus, A.; and Schocken, S. 1993.
Multilayer feedforward networks with a nonpolynomial ac-

tivation function can approximate any function. Neural networks 6(6):861–867.
Livni, R.; Shalev-Shwartz, S.; and Shamir, O. 2014. On
the computational efﬁciency of training neural networks. In
Advances in Neural Information Processing Systems, 855–
863.
Mohri, M.; Rostamizadeh, A.; and Talwalkar, A. 2012.
Foundations of machine learning. MIT press.
Neyshabur, B.; Tomioka, R.; and Srebro, N. 2015. Normbased capacity control in neural networks. In Proceedings
of The 28th Conference on Learning Theory, 1376–1401.
Pennington, J.; Yu, F.; and Kumar, S. 2015. Spherical random features for polynomial kernels. In Advances in Neural
Information Processing Systems, 1846–1854.
Rahimi, A., and Recht, B. 2008. Random features for largescale kernel machines. In Advances in Neural Information
Processing Systems, 1177–1184.
Rahimi, A., and Recht, B. 2009. Weighted sums of random
kitchen sinks: Replacing minimization with randomization
in learning. In Advances in neural information processing
systems, 1313–1320.
Scholkopf, B., and Smola, A. J. 2001. Learning with kernels: support vector machines, regularization, optimization,
and beyond. MIT press.
Shamir, O. 2016. Distribution-speciﬁc hardness of learning
neural networks. arXiv preprint arXiv:1609.01037.
Shawe-Taylor, J., and Cristianini, N. 2004. Kernel methods
for pattern analysis. Cambridge university press.
Silver, D.; Huang, A.; Maddison, C. J.; Guez, A.; Sifre, L.;
Van Den Driessche, G.; Schrittwieser, J.; Antonoglou, I.;
Panneershelvam, V.; Lanctot, M.; et al. 2016. Mastering
the game of go with deep neural networks and tree search.
Nature 529(7587):484–489.
Sindhwani, V.; Avron, H.; and Mahoney, M. W. 2014.
Quasi-monte carlo feature maps for shift-invariant kernels.
In International Conference on Machine Learning.
Sun, S.; Chen, W.; Wang, L.; Liu, X.; and Liu, T.-Y. 2016.
On the depth of deep neural networks: a theoretical view. In
Proceedings of the Thirtieth AAAI Conference on Artificial
Intelligence, 2066–2072. AAAI Press.
Swirszcz, G.; Czarnecki, W. M.; and Pascanu, R. 2016. Local minima in training of deep networks. arXiv preprint
arXiv:1611.06310.
Telgarsky, M. 2016. Beneﬁts of depth in neural networks. In
29th Annual Conference on Learning Theory, 1517–1539.
Xie, P.; Deng, Y.; and Xing, E. 2015. On the generalization error bounds of neural networks under diversityinducing mutual angular regularization. arXiv preprint
arXiv:1511.07110.
Xie, B.; Liang, Y.; and Song, L.
2016.
Diversity
leads to generalization in neural networks. arXiv preprint
arXiv:1611.03131.

Appendix
In practice, we can typically assume that x ∈ Ω ⊆ Rd
with some sufﬁciently large compact subspace Ω. Thus,
given such a Ω, we use the following sampling procedure in order for our method to be efﬁciently executable
in practice: r is sampled uniformly from a d dimensional
unit sphere Sd−1 , and r0 is sampled such that the probability measure on any open ball in R with its center in
[−radius(Ω), radius(Ω)] is nonzero.2 For example, uniform
distribution that covers [−radius(Ω), radius(Ω)] or normal
distribution with a nonzero ﬁnite variance sufﬁces the above
requirement.

Visualization of semi-random features
We visualize the surfaces of two semi-random features (a)
linear (LSR) (s = 0) and (b) square (SSR) (s = 1) in Figure 4. The semi-random features are deﬁned by two parts.
The random weights and nonlinear thresholding deﬁne a hyperplane, where on one side, it is zero and on the other side,
it is either an adjustable linear or adjustable square function. During learning, gradient descent is used to tune the
adjustable parts to ﬁt target functions.

Importance of Bias Terms in First Layer
We note the importance of the bias term r0 in the ﬁrst layer,
to obtain universal approximation power. Without the bias
term, Theorem 1 does not hold. Indeed, it is easy to see that
without the bias term, Heaviside step functions do not form a
function class with universal approximation ability for x ∈
Rd . The importance of the bias term can also be seen by
the binomial theorem. For example, let g(z) (z ∈ R) be a
polynomial of degree at least c for some constant c, then g
does not form a function class that contains the polynomial
of degree less than c in z. However, if z = x + r0 , it is
possible to contain polynomials of degrees less than c in x
by binomial theorem.

Proofs for One Hidden Layer Model
In this section, we provide the proofs of the theorems presented in Section “One Hidden Layer Model”. We denote a
constant function x → 1 by 1 (i.e., 1(x) = 1).

Proof of Theorem 1 (Universal Approximation)
To prove the theorem, we recall the following known result.
Lemma 8 (Leshno et al. 1993, Proposition 1) For any σ :
R → R, span({σ(r x + r0 ) : r ∈ Rd , r0 ∈ R}) is dense in
Lp (Ω) for every compact subset Ω ⊆ Rd and for every p ∈
[1, ∞),if and only if σ is not a polynomial of finite degree
(almost everywhere).
With this lemma, we ﬁrst prove that for every ﬁxed s, the
span of a set of our functions σs is dense in L2 (Ω).
2

The radius of a compact subspace Ω of Rd is deﬁned as
radius(Ω) = supx∈Ω x2 and a open ball is deﬁned as B(z̄, δ) =
{z : z − z̄ < δ} where z̄ is the center of the ball.

(a) linear semi-random

(b) square semi-random

Figure 4: Visualization of linear semi-random (LSR) (s = 0)
and squared semi-random (SSR) (s = 1) units. The random
weights deﬁne a hyperplane: on one side, it is zero while
on the other side, it is either adjustable linear or adjustable
square function.
Lemma 9 For every fixed finite integer s and for every compact subset Ω ⊆ Rd , span{σs (r x + r0 ), 1(x) : r ∈
Sd−1 , r0 ∈ R} is dense in L2 (Ω), where g(x) = 1 represents a constant function.
Proof Since σ0 is Heaviside step function, σ0 is not a polynomial of ﬁnite degree and span{σ0 (r x + r0 ) : r ∈
Rd , r0 ∈ R} is dense in L2 (Ω) (we know each of both
statements independently without Lemma 8). For the case
of s ≥ 1, σs (z) = z s σ0 (z) is not a polynomial of ﬁnite degree (i.e., otherwise, σ0 is a polynomial of ﬁnite degree). Thus, by Lemma 8, for ﬁxed ﬁnite integer s ≥ 1,
span{σs (r x + r0 ) : r ∈ Rd , r0 ∈ R} is dense in L2 (Ω).

Since σs ( rrx + r0 ) = r −s σs (r x + r r0 ) for all
r = 0 (because H(z) is insensitive to a positive scaling of
z),
span{σs (r x + r0 ) : r ∈ Rd , r0 ∈ R}

= span{σs (r x + r0 ) : r ∈ Sd−1 ∪ {0}, r0 ∈ R}

= span{σs (r x + r0 ), 1(x) : r ∈ Sd−1 , r0 ∈ R}.

This completes the proof.

In practice, we want to avoid wasting samples; we want
to have a choice to sample r0 from the interval that matters.
In order to do that, we admit the dependence of our function
class on Ω via the following lemma.
Lemma 10 For every fixed finite integer s and for every
compact subset Ω ⊆ Rd , span{σs (r x + r0 ), 1(x) : r ∈
Sd−1 , |r0 | ≤ radius(Ω)} is dense in L2 (Ω).

Proof For any r0 ∈ R \ [−radius(Ω), +radius(Ω)], since
|r x| ≤ r 2 x 2 ≤ radius(Ω), either σs (r x + r0 ) = 0
or = (r x + r0 )s , for all x ∈ Ω. Meanwhile, with
r0 = radius(Ω), σs (r x + r0 ) = (r x + r0 )s , for all
s
s s−k  k
x ∈ Ω. Since (r x + r0 )s =
(r x) ,
k=0 k r0
it is a polynomial in x, and r0 only affects the coefﬁcients of the polynomial. Thus, for any x ∈ Ω,
span{σs (r x + r0 ), 1(x) : r ∈ Sd−1 , r0 ∈ R} =
span{σs (r x + r0 ), 1(x) : r ∈ Sd−1 , |r0 | ≤ radius(Ω)},
because the difference in the nonzero coefﬁcients of the
polynomials can be taken care of by adjusting the coefﬁcients of the linear combination, unless the right-hand side

only contains r0 ∈ {0}. Having only r0 ∈ {0} implies
Ω = {0}, which is assumed to be false.
We use the following lemma to translate the statement
with uncountable samples in Lemma 10 to countable samples in Theorem 1.
Lemma 11 Let s be a fixed finite integer, Ω be any compact
subset of Rd , and (r̄, r̄0 ) is in
Sd−1 × [−radius(Ω), radius(Ω)].

Then, for any > 0, there exists a δ such that for any
(r, r0 ) ∈ B((r̄, r̄0 ), δ),

(σs (r x + r0 ) − σs (r̄ x + r̄0 ))2 < .

(ri , r0i ) in any δ-ball becomes one. This shows the last
line. The above inequality leads a contradiction |c| < |c| by
choosing < |c|/ f0 . This completes the proof.

Proof of Theorem 2 (No Bad Local Minima and
Few Bad Critical Points)
Since multiplying a constant in w does not change the optimization problem (in terms of optimizer), in the proof of
Theorem 2, we consider
L(w) :=
to be succinct.

m 

i=1

2
yi − fˆns (xi ; w) ,

x∈Ω

Proof For the case of s = 0, it directly follows the proof of
lemma II.5 in (Huang et al. 2006). For the case of s ≥ 1, we
can think of σs as a uniformly continuous function with a
compact domain that contains all the possible inputs of σs
for our choice of Ω and Sd−1 × [−radius(Ω), radius(Ω)].
The lemma immediately follows from its uniform continuity.

Proof of Theorem 2 (iii) and (iv) Deﬁne a m × (nd) matrix
as
⎤
⎡ (2)
(2)


· · · wn σs (x
w1 σs (x
1 r1 )x1
1 rn )x1
⎥
⎢
..
..
..
D = ⎣
⎦.
.
.
.
(2)


w1 σs (x
m r1 )xm

(2)


wn σs (x
m rn )xm

···

(2)

Note that if s → ∞, x → ∞, r → ∞, or r0 → ∞,
the proof of Lemma 11 does not work. We are now ready to
prove Theorem 1.
Proof of Theorem 1 Fix s and let σi : x → σs (ri x + r0i )
with ri and r0i being sampled randomly for i = 2, 3, . . .
as speciﬁed in Section “Semi-Random Features”. Let σ1 :
x → σs (r1 x + r01 ) = 1(x) = 1 as speciﬁed in Section
“One Hidden Layer Model”. Since span{σ1 , σ2 , ..., σn } ⊆
Fns (as we can get span{σ1 , σ2 , ..., σn } by only using
(1)
(2)
w0i and wi ), we only need to prove the universality of
span{σ1 , σ2 , ..., σn } as n → ∞. We prove the statement
by contradiction. Suppose that span{σ1 , σ2 , ...} is not dense
in L2 (Ω) (with some nonzero probability). Then, there exists a nonzero f0 ∈ L2 (Ω) such that f0 , σi  = 0 for all
i = 1, 2, . . .. However, from Lemma 10, either f0 , σ1  = 0
or there exists (r̄, r̄0 ) ∈ Sd−1 × [−radius(Ω), radius(Ω)]
such that | f0 , σr̄,r̄0  | = |c| > 0 with some constant c,
where σr̄,r̄0 = σs (r̄ x + r̄0 ). If f0 , σ1  = 0, we get the desired contradiction, and hence we considerer the latter case.
In the latter case, for any > 0,
|c| = | f0 , σr̄,r̄0  |
= | f0 , σr̄,r̄0  − f0 , σi  | ∀i ∈ {2, 3, ...}
= | f0 , σr̄,r̄0 − σi  | ∀i ∈ {2, 3, ...}
≤ σr̄,r̄0 − σi f0 ∀i ∈ {2, 3, ...}
< f0 ∃i ∈ {2, 3, ...},

where the last line follows from lemma 11 and our sampling
procedure; that is, for any > 0, if a sampled (ri , r0i )
is in a δ-ball, σi − σr̄,r̄0 < (from Lemma 11). Since
our sampling procedure allocates nonzero probability on
any such interval, as i → ∞, the probability of sampling


where wj σs (x
i rj )xi is a 1 × d block at (i, j)-th block
entry. Then, we can rewrite the objective function as

where

L(w) = D vec(w(1) ) − Y
(1)

vec(w(1) ) = (w1

(1)

, w2

2
2,

, . . . , wn(1) ) .

Therefore, at any critical point w.r.t. w(1) (i.e., by taking
gradient of L w.r.t. w(1) and setting it to zero), Ŷ is the projection of Y onto the column space of D .
n
(2)
 (1)
Since fˆn (xi ; w) = k=1 σs (x
i rk )(xi wk )wk ,
L(w) = D[[w]] − Y

where
(2)

(1)

[[w]] = (w1 w1

(2)

(1)

, w 2 w2

2
2,

, . . . , wn(2) wn(1) ) .

(8)

Therefore, we achieve a global minimum if Ŷ is the projection of Y onto the column space of D.
(2)
If wk = 0 for all k ∈ {1, 2, . . . , n}, the column space

of D is equal to that of D. Hence, Ŷ being the projection
of Y onto the column space of D is achieved by our
parameterization, which completes the proof of Theorem
2 (iv). Any critical point w.r.t. (w(1) , w (2) ) needs to be a
critical point w.r.t. w(1) . Hence, every critical point is a
(2)
global minimum if wk = 0 for all k ∈ {1, 2, . . . , n}. This
completes the proof of Theorem 2 (iii).
(2)

Proof of Theorem 2 (ii) From Theorem 2 (iii), if wk = 0
for all k ∈ {1, 2, . . . , n}, every local minimum is a global
minimum (because a set of local minima is included in a set
of critical points).

(2)

Consider any point w̄ where w̄k = 0 for some k ∈
{1, 2, . . . , n}. Without loss of generality, with some integer
(2)
(2)
c, let w̄k = 0 for all k ∈ {1, . . . , c − 1} and w̄k = 0 for
all k ∈ {c, . . . , n}.
Then, at the point w̄,
2
2,

L(w̄) = D1:c−1 [[w̄]]1:c−1 − Y

where D1:c−1 and [[w̄]]1:c−1 are the block elements of D
(equation 4) and [[w̄]] (equation 8) that correspond to the
nonzero w̄k . Then, from the proof of Theorem 2 (iii), if w̄ is
a critical point, then D1:c−1 [[w̄]]1:c−1 is the projection of Y
onto D1:c−1 .
If w̄ is a local minimum, for sufﬁciently small > 0, for
(1)
(2)
any k̄ ∈ {c, . . . , n}, for any (wk̄ , wk̄ ) ∈ Rd × R ,
2
2

D1:c−1 [[w̄]]1:c−1 − Y

≤ D1:c−1 [[w̄]]1:c−1 + Dk̄ [[w̄ ]]k̄ − Y

which is simpliﬁed to

+ Dk̄ [[w̄ ]]k̄

2

,

(9)

where [[w̄ ]]k̄ is a -perturbation of [[w̄]]k̄ as
[[w̄ ]]k̄ =

(2)
(1)
wk̄ (w̄k̄

+

(1)
wk̄ ).

Here, where D,k̄ is the corresponding k̄-th block of size
m × d. In equation (9), the ﬁrst term contains and 2 terms
whereas the second term contains 2 , 3 and 4 terms. With
sufﬁciently small, the term must be zero (as we can change
its sign by the direction of perturbation). Then, with sufﬁciently small, the 2 terms become dominant and must satisfy
(1)

(2)

0 ≤2 2 (D1:c−1 [[w̄]]1:c−1 − Y ) (Dk̄ wk̄ wk̄ )
2

+
(1)

(1)

(2) 2

Dk̄ w̄k̄ wk̄

,

(2)

for any (wk̄ , wk̄ ) ∈ Rd × R, which implies
(1)

(2)

(1)

for any |wk̄ | << minj∈{1,2,...,d} |(wk̄ )j | (see footnote3 ).
It implies that for any k̄ ∈ {c, . . . , n},
0 = (D1:c−1 [[w̄]]1:c−1 − Y ) Dk̄ .

Since D1:c−1 [[w̄]]1:c−1 is the projection of Y onto
D1:c−1 (as discussed above), this implies that
D (D[[w̄]] − Y ) = (D1:c−1 [[w̄]]1:c−1 − Y ) D = 0. Thus,
D[[w̄]] is the projection of Y onto D, which is a global
minimum.
3
For readers who have experienced deriving some conditions
of a local minimum, it may sound too good to be true that we can
ignore the 2 term from Dk̄ [[w̄ ]]k̄ 2 . However, this is true and
(2)
not that good since we are considering a spacial case of w̄k̄ = 0.
(2)

(2)

(2)

(1)

With a standard use of Rademacher complexity (e.g., see
section 3 in Mohri, Rostamizadeh, and Talwalkar for a clear
introduction of Rademacher complexity, and the proof of
lemma 12 in Xie, Liang, and Song for its concrete use),
1
Ex (f (x) − fˆ(x; w∗ ))2 − L(w)
2


(1)

Indeed, if w̄k̄ = 0, [[w̄ ]]k̄ = (w̄k̄ + wk̄ )(w̄k̄ + wk̄ )
creates the 2 term from Dk̄ [[w̄ ]]k̄ 2 that we cannot ignore.

(10)

log 1δ
+ 2(CY + CŶ )Rm (F),
2m

≤ (CY2 + CŶ2 )

with probability at least 1 − δ, where (CY2 + CŶ2 ) is an upper
bound on the value of L(w), and (CY + CŶ ) comes from
an upper bound on the Lipschitz constant of the squared loss
function. Note that f (x; w) = σ, x w ≤ CW Cσx .
We ﬁrst compute an upper bound on the empirical
Rademacher complexity of our model class F as follows:
with Rademacher variables ξi ,


m


ξi σ, xi  w
mR̂m (F) = E sup
w∈Sw i=1
m







≤ E sup 
ξi σ, xi  w 2

w∈Sw  i=1
 2
 m




= CW E 
ξi σ, xi 




0 ≤2 2 (D1:c−1 [[w̄]]1:c−1 − Y ) (Dk̄ wk̄ wk̄ ),
(2)

if Dk = 0. Thus, it is non-convex if D = 0.

Proof of Theorem 3 (Generalization Bound for
Shallow Model)

2
2,

0 ≤2(D1:c−1 [[w̄]]1:c−1 − Y ) (Dk̄ [[w̄ ]]k̄ )

Proof of Theorem 2 (i) It is sufﬁcient to prove non(1)
(2)
convexity w.r.t. the parameters (wk , wk ) for each k ∈
(k)
(2)
{1, 2, . . . , n}; that is, we prove that L (w1 , wk ) =
(1) (2)
L(w)|wi=k =0 = Dk wk wk − Y 22 is non-convex, where
D,k is the corresponding k-th block of size m × d. It is
indeed easy to see that L is non-convex. For example, at
(1)
(2)
(wk , wk ) = 0, its Hessian is


0
Dk
 0,
Dk
0

i=1

2

where the second line follows the Cauchy-Schwarz inequality. Furthermore,
  
 m

 m m


 


E 
ξi σ, xi  ≤ E
ξi ξj σ, xi  σ, xj 


i=1
i=1 j=1
2

m

σ, xi 22
≤
i=1

√
≤ Cσx m.

where the ﬁrst line uses Jensen’s inequality for the concave function, and the second line follows that for any



m 2
m m
≤
z, E
ξ
ξ
z
z
i
j
i
j
j=1
i=1
i=1 zi , where we used
the deﬁnition of Rademacher variables for E[ξi ξj ]. Hence,
σx
. Since the Rademacher complexity
R̂m (F) ≤ CW√C
m

Rm (F) is the expectation of R̂m (F) over samples (x, y),
σx
Rm (F) ≤ CW√C
.
m

Tensorial Structure
Since the output of our network is the sum of the outputs
of all the paths in the network, we observe the following
interesting structure:
fˆns1 ,...,nH (x)
=

n1
d



k0 =0 k1 =1

where

Proof of Theorem 5 (Lower Bound on Universal
Approximation Power)
Since the output of our network is the sum of the outputs of
all the paths in the network, we observed in Section “Tensorial Structure” that
fˆn1 ,...,nH (x)
=

n1
d



k0 =0 k1 =1

···

nH


[[σ]]k1 ,...,kH (x)xk0 [[w]]k0 ,k1 ,...,kH ,

kH =1

where
(1)

···

nH


[[σ]]k1 ,...,kH (x) = σs (x rk1 )

kH =1

σk1 ,...,kH (x)xk0 wk0 ,k1 ,...,kH ,
(1)

σk1 ,...,kH (x) = σs (x rk1 )


H

l=2

(l)

H



(l)

σs (h(l−1)
(x)rkl ),
r

l=2

and
[[w]]k0 ,k1 ,...,kH =

H


(l)
wkl−1 kl

l=1

(l)

σs (h(l−1)
(x) rkl ),
r

H




(H+1)

wk H

.

By re-writing gk0,k1,...,kH (x) = [[σ]]k1 ,...,kH (x)xk0 , this
means that

(H+1)

and wk0 ,k1 ,...,kH =
. This means
l=1 wkl−1 kl wkH
that the function is the weighted combination of (d + 1) ×
n1 × n2 . . . nH nonlinear basis functions, which is exponential in the number of layers H. Alternatively, the function
can also be viewed as the inner product between two tensors
σk1 ,...,kH (x)xk0 and wk0 ,k1 ,...,kH , which is nonlinear in
input x, but linear in the parameter tensor wk0 ,k1 ,...,kH .
We note that the parameter tensor is not arbitrary but highly
structured: it is composed using a collection of matrices with
H
dw = (d + 1)n1 + nH + l=2 nl−1 nl number of adjustable
parameters. Such special structure allows the function to
generate an exponential number of basis functions yet keep
the parameterization compact.

Fns1 ,...,nH ⊆ span({gk0,k1,...,kH : ∀k0 , ∀k1 , . . . , ∀kH }).
Then, Theorem 5 is a direct application of the following result:
Lemma 12 (Barron 1993, Theorem 6) For any fixed set of
M functions g1 , g2 , . . . , gM ,
sup

inf

f ∈ΓC fˆ∈span(g1 ,g2 ,...,dM )

f − fˆ

L2 (Ω)

≥κ

C −1/d
M
.
d

We apply Lemma 12 with our gk0,k1,...,kH to obtain the
statement of Theorem 5.

Proofs for Multilayer Model
In this section, we provide the proofs of the theorems and
corollary presented in Section “Multilayer Model”.

Proof of Corollary 4 (Universal Approximation
with Deep Model)
(1)

(2)

(3)

(H)

From the deﬁnition of r1 , r1 , r1 , . . . , r1 , the ﬁrst unit
in each hidden layer of the main semi-random net is not
affected by the random net; the output of the ﬁrst unit of
semi-random net is multiplied by one. By setting most of
the weights w to zeros, we can use only the ﬁrst unit at each
layer from the second hidden layer, creating a single path
that is not turned-off by the random activation net from the
second hidden layer. Then, by adjusting weights w in the
path, we can create an identity map from the second hidden
layer unit, from which the path starts. In theorem 1, we have
already shown that the output of any second hidden layer
unit has universal approximation ability. This completes the
proof.

Proof of Corollary 6 (No Bad Local Minima and
Few Bad Critical Points w.r.t. Two Last Layers)
It directly follows the proof of Theorem 2, because
any critical point (or any local minimum) is a critical
point (or respectively, a local minimum) with respect to
(W (H) , W (H+1) ).

Proof of Corollary 7 (Generalization Bound for
Deep Model)
By noticing that fn1 ,...,nH (x) = vec(w) vec(σ, x) (i.e.,
the tensorial structure as discussed in the proof of Theorem
5), it directly follows the proof of Theorem 3 (Generalization
Bound – shallow); in the proof of Theorem 3, replace σ, xi 
and w with the corresponding deﬁnitions for multilayer
models, which are vec(σ, xi ) and vec(w) of multilayer
models.

Discussion on an Upper Bound on
Approximation Error for Multilayer Model
Here, we continue the discussion in Section “Beneﬁt of
Depth”. An upper bound on approximation error can be decomposed into three illustrative terms as
f − fˆ2L2 (Ω)

=
f − fˆ2L2 (Ωi )

(11)

i∈S

≤


i∈S

f −

fi 2L2 (Ωi )

Figure 5: Training error for a simple test function.

+ fˆ − fˆi L2 (Ωi ) + fi − fˆi 2L2 (Ωi ) ,

where {Ωi }i∈S is an arbitrary partition of Ω, and fi = f (xi )
is a constant function with some xi ∈ Ωi . Then, the ﬁrst
term and the second term represent how much each of f and
fˆ varies in Ωi . In other words, we can bound these two terms
by some “smoothness” of f and fˆ (some discontinuity at a
set of measure zero poses no problem as we are taking the
norm of L2 (Ωi ) space). The last term in equation (11) represents the expressive power of fˆ on the ﬁnite points {xi }i∈S
(assuming that S is ﬁnite). Essentially, to obtain a good upper bound, we want to have a fˆ that is “smooth” and yet
expressive on some ﬁnite points.
For the ﬁrst term in equation (11), we can obtain a bound
via a simple calculation for any f ∈ ΓC as follows: by the
mean value theorem, there exists zx for each x such that

f − fi 2L2 (Ωi ) =
(f (x) − f (xi ))2
x∈Ωi

=
(∇f (zx ) (x − xi ))2 dx
x∈Ωi

2
L2 (Ωi ) x − xi
2
L2 (Rd ) x − xi
xi 2L2 (Ωi ) .

≤ ∇f (zx )

≤ ∇f (x)

2
L2 (Ωi )
2
L2 (Ωi )

≤ C2 x −

where x − xi 2L2 (Ωi ) = x∈Ωi x − xi 22 . Hence, the bound
goes to zero as the size of Ωi decreases.
We can use the same reasoning to bound the second term
in equation (11) with some “smoothness” of fˆ. A possible concern is that fˆ becomes less “smooth” as the width
n and depth H increase. However, from Bessel’s inequality
with Gram-Schmidt process, it is clear that such an effect is
bounded for a best fˆ with a ﬁnite s ≥ 1 (notice the difference from statistical learning theory, where we cannot focus
on the best fˆ). The last term in equation (11) represents the
error at points xi ∈ Ωi , which would be bounded via optimization theory, as the expressive power of fˆ increases as
depth H and width increase.
While this reasoning illustrates what factors may matter,
a formal proof is left to future work.

Additional Experimental Details
In this section, we provide additional Experimental details.

A simple test function
Figure 5 shows the training errors for a sine function experiment discussed in Section “Simple Test Function”. It shows
roughly the same patterns as in their test errors, indicating
that there is no large degree of overﬁtting.
Figures 6, 7, and 8 visualize the function learned at each
iteration for each method. We can see that semi-random features (LSR) learns the function quickly.
At each trial, 5000 points of the inputs x were sampled
uniformly from [−12π, 12π] for each of training dataset and
test dataset. We repeated this trial 20 times to obtain standard deviations. All hyperparameters were ﬁxed to the same
values for all methods, with the network architecture being
ﬁxed to 1-50-50-1. For all the methods, we used the learning rate 5×10−4 , momentum parameter 0.9, and mini-batch
size of 500. For all methods, the weights are initialized as
0.1 times random samples drawn from the standard normal
distribution.

UCI datasets
For all methods on all datasets, we use SGD with 0.9 momentum, a batch size of 128 to run 100 epochs (passes over
the whole dataset). The initial learning rate is set to 0.1 (for
some combinations, this leads to NaN and we use 0.02 as initial learning rate). We also use an exponential staircase decaying schedule for the learning rate where after one epoch
we decrease it to 0.95 of the original rate. The√parameters
are initialized as normal distributions times 1/ d where d
is the input dimension for the layer.
For ReLU and semi-random features (LSR and SSR),
their behaviors change similarly on most datasets: more layers and more units lead to better performance. Also, on most
datasets, adding more layers leads to better performance
than adding more units per layer. The dataset adult is an
exception where all methods have similar performance and
more parameters actually lead to worse performance. This
is likely due to that adult is quite noisy and using more parameters leads to overﬁtting. Overﬁtting is also observed on
senseit dataset for LSR where the training errors keep become smaller but test errors increase.

Image datasets
We train these convolution neural networks using SGD with
0.9 momentum with a batch size of 64 for MNIST and 128

1

1

0.5

0.5

1

0.5

0

0

0
-0.5

-0.5
-1

-0.5

-1
-1.5

-1

-1.5

-2
-40

-2

-30

-20

-10

0

10

20

30

40

-2.5
-40

-30

-20

(a) iteration 0

-10

0

10

(b) iteration 10

20

30

40

-1.5
-40

1.5

1.5

1

1

1

0.5

0.5

0.5

0

0

0

-0.5

-0.5

-0.5

-1

-1

-1

-30

-20

-10

0

10

20

30

40

-1.5
-40

-30

(d) iteration 104

-20

-10

0

-20

-10

0

10

(c) iteration 10

1.5

-1.5
-40

-30

2

10

20

30

40

-1.5
-40

-30

-20

(e) iteration 5 × 104

-10

0

20

30

40

3

10

20

30

40

(f) iteration 105

Figure 6: Function learned at each iteration (Semirandom - LSR).
1

1.5

1
0.8

1

0.5

0.6
0.4

0

0.5

0.2
0

0
-0.5

-0.2
-0.5

-0.4
-1

-0.6

-1

-0.8
-1.5
-40

-30

-20

-10

0

10

20

30

40

-1
-40

-30

-20

(a) iteration 0

-10

0

10

(b) iteration 10

20

30

40

-1.5
-40

-20

-10

0

10

(c) iteration 10

2

2.5

1.5

2

2

1.5

1.5

1

-30

2

20

30

40

20

30

40

3

2.5

1

1

0.5

0.5

0.5
0
-0.5

0

0

-0.5

-0.5

-1

-1

-1

-1.5

-1.5

-1.5

-2
-40

-30

-20

-10

0

10

(d) iteration 10

4

20

30

40

-2
-40

-30

-20

-10

0

10

(e) iteration 5 × 10

20

30

40

-2
-40

-30

4

Figure 7: Function learned at each iteration (fully random).

-20

-10

0

10

(f) iteration 10

5

1

1

1

0.8

0.8

0.8

0.6

0.6

0.6

0.4

0.4

0.4

0.2

0.2

0.2

0

0

0

-0.2

-0.2

-0.2

-0.4

-0.4

-0.4

-0.6

-0.6

-0.6

-0.8

-0.8

-0.8

-1
-40

-30

-20

-10

0

10

20

30

40

-1
-40

-30

-20

(a) iteration 0

-10

0

10

(b) iteration 10

20

30

40

-1
-40

1.5

1.5

1

1

1

0.5

0.5

0.5

0

0

0

-0.5

-0.5

-0.5

-1

-1

-1

-30

-20

-10

0

10

20

30

40

-1.5
-40

-30

(d) iteration 104

-20

-10

0

-20

10

-10

0

10

(c) iteration 10

1.5

-1.5
-40

-30

2

20

30

40

-1.5
-40

-30

-20

(e) iteration 5 × 104

-10

0

20

30

40

20

30

40

3

10

(f) iteration 105

Figure 8: Function learned at each iteration (ReLU).
for CIFAR10 and SVHN. On each dataset, we have tried initial learning rates of 0.1, 0.01, and 0.001 and used an exponential decaying schedule. For MNIST, the decay schedule
is decreasing to 0.95 of the previous rate after each epoch.
For CIFAR10 and SVHN, the schedule is decreasing to 0.1
of the previous rate after 120 epochs. The best performing
results on a validation set are then picked as the ﬁnal solution, and usually the learning rate of 0.1 worked best. The
parameters are also
√ randomly initialized from normal distributions times 1/ d where d is the input dimension for the
layer.

Better than Random Feature?
The experimental results veriﬁed our intuition that semirandom feature can outperform random feature with fewer
number of unites due to its learnable weights. We can also
strengthen this intuition via the following theoretical insights. Let f ∈ Frandom,n1 ···nH be a function that is a composition of any fully-random features with depth H where
the adjustable weights are only in the last layer. The following corollary states that a model class of any fully-random
features has a approximation power exponentially bad in the
dimensionality of x in the worst case.
Corollary 13 (Lower bound on approximation power for
fully-random feature) Let Ω = [0, 1]d . For any depth H ≥ 0,

and for any set of nonzero widths {n1 , n2 , . . . , nH },
sup
f ∈ΓC

inf

fˆ∈Frandom,n1 ···nH

f − fˆ

L2 (Ω)

≥

κC
(nH )−1/d ,
d2

where κ ≥ (8πe(π−1) )−1 is a constant.
Proof The model class of any random features is the span
of the features represented at the last hidden layer. Hence,
the statement of the corollary directly follows from the
proof of Theorem 5.
Corollary 13 (lower bound for fully-random feature) together with Theorem 5 (lower bound for semi-random feature) reﬂects our intuition that semi-random feature model
can potentially get exponential advantage over random feature by learning hidden layer’s weights. Again, because the
lower bound may not be tight, this is intended only to aid
our intuition.
We can also compare upper bounds on their approximation errors with an additional assumption. Assume that we
can represent a target function f using some basis as

f (x) =
σ(r x)(w x) p(r, w).
r∈Sd−1 ,w≤CW

Then, we can obtain the following results.
• If we have access to the true distribution p(r, w), f (x)
can be approximated as a ﬁnite sample average, obtaining
approximation error of O( √1n ).

• Without knowing the true distribution p(r, w), a purely
random feature approximation of f (x) incurs a large approximation error of O( q0 ·q1c·√n ). Here, q0 is inverse of
the hyper-surface area of a unit hypersphere and q1 is the
inverse of the volume of a ball of radius CW .
• Without knowing the true distribution p(r, w), semirandom feature approach with one hidden layer model can
√ ).
obtain a smaller approximation error of O( 35c
n

layer model), we obtain the following approximation: with


i.i.d.
ri ∼ q(r) := Uniform Sd−1 ,

A detailed derivation of this result is presented in Appendix.

By ﬁnding the best w, we can get at least as good as the case
i.i.d.
of wi ∼ p(w|r) where p(w|r) = p(w, r)/p(r) is the true
conditional distribution given r. If p(r) ∞ = c, then the
important weight βi can be much smaller, with a bound of
c
q0 ≤ 35c independent of the dimension. Accordingly, this
√ ).
approximation incurs an error of O( 35c
n

Derivation of Upper Bounds on
Approximation Errors
With an additional assumption, we compare upper bounds
on approximation errors for random features and semirandom features. Recall the additional assumption that we
can represent a target function f using some basis as

σ(r x)(w x) p(r, w),
f (x) =
r∈Sd−1 ,w≤CW

where we can write p(r, w) = p(r) p(w|r). If we have access to the true distribution p(r, w), then f (x) can be approximated as a ﬁnite sample average
n

1
f1 (x) =
σ(ri x)(wi x)
n i=1
i.i.d.

(ri , wi ) ∼ p(r, w).

Such an approximation will incur an error of O( √1n ).
On the other hand, without knowing the true distribution p(r, w), a purely random feature approximation of f (x)
will be sampling both r and w from uniform distribution:


i.i.d.
i.i.d.
with ri ∼ q(r) := Uniform Sd−1 and wi ∼ q(w) :=
Uniform ( wi ≤ CW ),
n

1  p(ri , wi )
σ(ri x)(wi x)
f2 (x) =
n i=1 q(ri )q(wi )
n

=

1
αi σ(ri x)(wi x)
n i=1
d

where q(r) = q0 := Γ( n2 )/(2π 2 ) (inverse of the hypersurface area of a unit hypersphere) and q(w) = q1 :=
d

d
Γ( d2 + 1)/(π 2 CW
) (inverse of the volume of a ball of radius CW ). Interestingly enough, for the unit hypersphere,
the hyper-surface area reaches a maximum and then decreases towards 0 as d increases. One can show that the
seven-dimensional unit hypersphere has a maximum hypersurface area that is less than 35. But the volume of a ball
of radius CW depends exponentially on the radius CW . If
p(r, w) ∞ = c, then the importance weight αi can be as
large as q0c·q1 . Due to the fact that q1 can be very small which
make the bound very large, this incurs a big approximation
error of O( q0 ·q1c·√n ).
In contrast, if we sample r from unit sphere and then optimize over w (i.e., semi-random approach with one hidden

n

1  p(ri )
σ(ri x)(wi x)
f3 (x) =
n i=1 q(ri )
n

=

1
βi σ(ri x)(wi x).
n i=1

Additional Future Work and Open Problems
We have discussed several open questions in the relevant
sections of the main text. Here, we discusses additional future work and open problems.
We emphasize that our theoretical understanding of multilayer model is of preliminary nature and it is not intended
to be complete. Future work may adopt an oscillation argument in the previous work by (Telgarsky 2016), in order
to directly compare deep models and shallow models. To
obtain a tight upper bound, future work may start with a
stronger assumption such as one used in the previous work
by (Rahimi and Recht 2009); an assumption on a relationship of random sampling and the true distribution. Indeed,
when compared with neural network, semi-random feature
would be advantageous when the target function is known
to be decomposed into the unknown part w and the random
part r with known distribution. However, we did not assume
such a knowledge in this paper, leaving the study of such a
case to feature work.

