# Efficient Algorithms for Learning Depth-2 Neural Networks with General ReLU
Activations

Pranjal Awasthi
Google Research
<EMAIL_ADDRESS>Alex Tang11footnotemark: 1
Northwestern University
<EMAIL_ADDRESS>Aravindan Vijayaraghavan
Northwestern University
<EMAIL_ADDRESS>The last two authors are supported by the National
Science Foundation (NSF) under Grant No. CCF-1652491, CCF-1637585 and CCF
1934931.

###### Abstract

We present polynomial time and sample efficient algorithms for learning an
unknown depth-2 feedforward neural network with general ReLU activations,
under mild non-degeneracy assumptions. In particular, we consider learning an
unknown network of the form
$f(x)={a}^{\mathsf{T}}\sigma({W}^{\mathsf{T}}x+b)$, where $x$ is drawn from
the Gaussian distribution, and $\sigma(t)\coloneqq\max(t,0)$ is the ReLU
activation. Prior works for learning networks with ReLU activations assume
that the bias $b$ is zero. In order to deal with the presence of the bias
terms, our proposed algorithm consists of robustly decomposing multiple higher
order tensors arising from the Hermite expansion of the function $f(x)$. Using
these ideas we also establish identifiability of the network parameters under
minimal assumptions.

##  1 Introduction

The empirical success of deep learning in recent years has led to a flurry of
recent works exploring various theoretical aspects of deep learning such as
learning, optimization and generalization. A fundamental question in the
theory of deep learning is to identify conditions under which one can design
provably time-efficient and sample-efficient learning algorithms for neural
networks. Perhaps surprisingly, even for the simplest case of a depth-2
feedforward neural network, the learning question remains unresolved. In this
work we make progress on this front by studying the problem of learning an
unknown neural network of the form

$y=f(x)=a^{\mathsf{T}}\sigma(W^{\mathsf{T}}x+b).$ (1)

We are given access to a finite amount of samples of the form $(x_{i},y_{i})$
drawn i.i.d. from the data distribution, where each $x_{i}$ is comes from the
standard Gaussian distribution $\mathcal{N}(0,I)$, and $y_{i}=f(x_{i})$. The
goal is to design an algorithm that outputs an approximation of the function
$f$ up to an arbitrary error measured in the expected mean squared metric
(squared $\ell_{2}$ loss). An efficient learning algorithm has running time
and sample complexity that are polynomial in the different problem parameters
such as the input dimensionality, number of hidden units, and the desired
error.

Without any further assumptions on the depth-2 network, efficient learning
algorithms are unlikely. The recent work of [DKKZ20] provides evidence by
proving exponential statistical query lower bounds (even when $x$ is Gaussian)
that rule out a broad class of algorithms.

Several recent works have designed efficient algorithms for depth-2 neural
networks in the special setting when the bias term $b=0$. One prominent line
of work [GLM18, BJW19] give polynomial time algorithms under the non-
degeneracy assumption that the matrix $W$ has full-rank. Another body of work
relaxes the full-rank assumption by designing algorithm that incur an
exponential dependence on the number of hidden units [DKKZ20, CKM20, GK19], or
a quasipolynomial dependence when the coefficients $\\{a_{i}:i\in[m]\\}$ are
all non-negative [DK20]. There is little existing literature on learning
neural networks in the presence of the bias term. A notable exception is an
approach based on computing “score functions” [JSA15] that applies to certain
activations with bias and requires various assumptions that are not satisfied
by the ReLU function. The diminished expressivity of neural networks without
the bias terms leads to the following compelling question:

Can we design polynomial time algorithms even in the presence of bias terms in
the ReLU units?

We answer the question affirmatively by designing efficient algorithms for
learning depth-2 neural networks with general ReLU activations, under the
assumption that $W_{d\times m}$ has linearly independent columns (hence $m\leq
d$). In fact, our algorithms can be extended to work under much weaker
assumptions on $W$, that allow for $m\gg d$ ($m\leq O(d^{\ell})$ for any
constant $\ell\geq 1$) in a natural smoothed analysis setting considered in
prior works [GKLW19] (see Theorem 3.2 and Corollary 3.3). An important
consequence of our techniques is the fact that the network parameters are
identifiable up to signs, as long as no two columns of $W$ are parallel, and
all the $\\{a_{i}:i\in[m]\\}$ are non-zero. Furthermore we show that this
ambiguity in recovering the signs is inherent unless stronger assumptions are
made.

Conceptual and technical challenges with bias terms. Similar to prior works
[GKLW19, JSA15], our techniques rely on the use of tensor decomposition
algorithms to recover the parameters of the network. In the absence of any
bias, it can be shown that the $4$th Hermite coefficient of the function
$f(x)$ takes the form $\hat{f}_{4}=\sum_{i=1}a_{i}w^{\otimes 4}_{i}$ where
$w_{i}$ are the columns of $W$ and $a=(a_{1},a_{2},\dots,)$. When the columns
of $W$ are linearly independent, existing algorithms for tensor decompositions
in the full-rank setting can be used to recover the parameters
[Har70].111Tensor decompositions will in fact recover each ReLU activation up
to an ambiguity in the sign. However, in the full-rank setting, the correct
sign can also be recovered (as we demonstrate later in Theorem 3.1).  However,
when bias terms are present, there are several challenges that we highlight
below.

In the presence of biases the $k$th Hermite expansion of $f(x)$ takes the form
$\hat{f}_{k}=\sum_{i=1}^{m}a_{i}g_{k}(b_{i})w^{\otimes k}_{i}$, where $g_{k}$
is a function that may vanish on some of the unknown $b_{i}$ parameters. This
creates a hurdle in recovering the corresponding $w_{i}$. A simple example
where the above approach fails is when some of the $b_{i}=\pm 1$, since the
corresponding rank-1 terms vanish from the decomposition of $\hat{f}_{4}$. To
overcome this obstacle, we first show give a precise expression for the
function $g_{k}(b_{i})$ involving the $(k-2)$th Hermite polynomial (see Lemma
3.5). We then design an algorithm that decomposes multiple tensors obtained
from Hermite coefficients to recover the parameters. We use various properties
of Hermite polynomials to analyze the algorithm e.g., the separation of roots
of consecutive Hermite polynomials is used to argue that each $w_{i}$ is
recovered from decomposing at least one of the tensors.

Secondly, in the presence of the bias terms, recovery of all the parameters
(even up to sign ambiguities) may not even be possible from polynomially many
samples. For instance, consider a particular hidden node with output
$\sigma(w^{\top}_{i}x+b_{i})$. If $b_{i}$ is a large positive number then it
behaves like a linear function (always active). Hence if multiple $b_{i}$s are
large positive numbers then one can only hope to recover a linear combination
of their corresponding weights and biases. On the other hand if $b_{i}$ is a
large negative constant then the activation is $0$ except with probability
exponentially small in $|b_{i}|$. We cannot afford a sample complexity that is
exponential in the magnitude of the parameters. Furthermore, when the columns
of $W$ are not linearly independent, the tensor decomposition based method
will only recover good approximations up to a sign ambiguity for the terms
whose bias does not have very large magnitude i.e., we recover
$\pm(w_{i}^{\top}x+b_{i})$ if $|b_{i}|$ is not very large.

To handle the above issue we proceed in two stages. In the first stage we
recover the network parameters (up to signs) of all the “good” terms, i.e.,
hidden units with biases of small magnitude. To handle the “bad” terms (large
magnitude bias) we show that a linear functions is a good approximation to the
residual function comprising of the bad terms. Based on the above, we show
that one can solve a truncated linear regression problem to learn a function
$g(x)$ that achieves low mean squared error with respect to the target $f(x)$.
The output function $g(x)$ is also a depth-2 ReLU network with at most two
additional hidden units than the target network.

There are several other technical challenges that arise in the analysis
sketched above, when there are sampling errors due to using only a polynomial
amount of data (for example, the tensors obtained from $\hat{f}_{k}$ may have
some rank-$1$ terms that are small but not negligible, that may affect the
robust recovery guarantees for tensor decompositions). We obtain our robust
guarantees by leveraging many useful properties of Hermite polynomials, and a
careful analysis of how the errors propagate.

The rest of the paper is organized as follows. We present preliminaries in
Section 1.1 followed by related work in Section 2. We then formally present
and discuss our main results in Section 3. In Section 4 we present our main
algorithm and analysis in the population setting, i.e., under the assumption
that one has access to infinite data from the distribution. We then present
the finite sample extension of our algorithm in Section 5 that achieves
polynomial runtime and sample complexity.

###  1.1 Model Setup and Preliminaries

We consider the supervised learning problem with input $x\in\mathbb{R}^{d}$
drawn from a standard $d$-dimensional Gaussian distribution
$\mathcal{N}(0,I_{d\times d})$ and labels $y$ generated by a neural network
$y=f(x)=a^{\mathsf{T}}\sigma(W^{\mathsf{T}}x+b)$, where
$a,b\in\mathbb{R}^{m}$, $W\in\mathbb{R}^{d\times m}$ and $\sigma$ is the
element-wise ReLU activation function, i.e., $\sigma(t)=\max(t,0)$. We denote
the column vectors of $W$ as $w_{i}\in\mathbb{R}^{d}$ with $i\in[m]$ and
$a_{i}$ as the $i$’th element of vector $a$, similarly for $b$ and $x$. We
pose a constraint on magnitudes of $a,b,W$ such that they are all $B$-bounded
for some $1\leq B\leq\mathrm{poly}(m,d)$, i.e. $\lVert a\rVert_{\infty},\lVert
b\rVert_{\infty},\lVert W\rVert_{\infty}\leq B$, and
$\min_{i\in[m]}|a_{i}|\geq 1/B$. Furthermore, we assume $\|w_{i}\|_{2}=1$
without loss of generality. If $w_{i}$ are not unit vectors, we can always
scale $a_{i}$ and $b_{i}$ to $\|w_{i}\|a_{i}$ and $\frac{b_{i}}{\|w_{i}\|}$
respectively so that $w_{i}$ are normalized. We will denote by $\Phi()$ the
cumulative density function (CDF) of the standard Gaussian distribution.
Finally, for a matrix $M$, we will use $s_{k}(M)$ to denote the $k$th largest
singular value of $M$.

For some parameters $n_{1},n_{2},\dots$ we will say that a quantity
$q\leq\mathrm{poly}(n_{1},n_{2},\dots)$ if and only if there exists constants
$c_{0}>0,c_{1}>0,c_{2}>0,\dots$ such that $q\leq
c_{0}n_{1}^{c_{1}}n_{2}^{c_{2}}\dots$. If these constant depend on another
parameter $\ell>0$ which is also a constant, then we will denote this by
$\mathrm{poly}_{\ell}(n_{1},n_{2},\dots)$. We will say that an probabilistic
event occurs with high probability if and if it occurs with probability
$1-(mdB)^{-\omega(1)}$ i.e., the failure probability is smaller that any
inverse polynomial in the parameters $m,d,B$. Finally, we will use sign
variables of the form $\xi_{i}\in\\{1,-1\\}$; they will typically capture an
ambiguity in the sign of the parameters of the $i$th unit.

##### Hermite Polynomials and Hermite Coefficients

The $d$-dimensional or $d$-variate Hermite polynomials of the $k$th order
evaluated at a point $x$ represented by the tensor
$He_{k}(x)\in(\mathbb{R}^{d})^{\otimes k}$ with $k\in\mathbb{N}$ (here
$\mathbb{N}$ includes $0$) form a complete orthogonal basis for the weighted
$L^{2}(\mathbb{R}^{d},\gamma)$ space with inner product defined as
$\mathbb{E}[f(x)g(x)]=\int_{\mathbb{R}^{d}}f(x)g(x)\gamma(x)dx$ for two
functions $f,g$, where $\gamma(x)=\exp(-\|x\|^{2}/2)/\sqrt{2\pi}$. We can thus
define the $k$’th Hermite coefficient of $f(x)$ by
$\hat{f}_{k}=\mathbb{E}[f(x)He_{k}(x)]$, similar to how the Fourier
coefficients are constructed. Throughout the context we wll also use
$He_{k}^{\alpha}(x)\in\mathbb{R}$ to denote a specific entry of $He_{k}(x)$
with $\alpha\in[d]^{k}$ as a multi-index if $d>1$. We now briefly introduce
the definition of Hermite polynomials along with some useful facts below.

######  Definition 1.1.

Let $D_{t}^{(k)}$ be the total differential operator taken $k$ times with
respect to $t$. For a function $g:\mathbb{R}^{d}\rightarrow\mathbb{R}$,
$D_{t}^{(k)}g(t)|_{t=t_{0}}\in(\mathbb{R}^{d})^{\otimes k}$, where the
$\alpha$’th element of $D_{t}^{(k)}g(t)$, for
$\alpha=(i_{1},i_{2},\dots,i_{k})$, is
$\frac{d}{dt_{i_{1}}}...\frac{d}{dt_{i_{k}}}g(t)$ with $\alpha\in[d]^{k}$
being a multi-index. Note that the above is invariant to permutations, i.e.,
for any permutation $\alpha^{\prime}$ of the indices in $\alpha$, the
$\alpha^{\prime}$th element of $D_{t}^{(k)}g(t)$ is the same as the $\alpha$th
element.

######  Definition 1.2.

Let $x\in\mathbb{R}^{d}$, $\gamma(x)=\exp(-\frac{\|x\|^{2}}{2})$, the
(probabilist’s) $k$’th $d$-dimensional Hermite polynomial
$He_{k}(x)\in(\mathbb{R}^{d})^{\otimes k}$ is given by

$He_{k}(x)=\frac{(-1)^{k}}{\gamma(x)}\cdot D_{x}^{(k)}\gamma(x)$ (2)

A particularly useful fact for 1-dimensional Hermite polynomials is their
relation with derivatives of a standard univariate Gaussian function.

######  Fact 1.3.

The $k$’th order derivative of $\gamma(x)$ can be written in terms of
1-dimensional Hermite polynomials as

$\frac{d^{k}}{dx^{k}}\gamma(x)=(-1)^{k}\cdot
He_{k}(x)\cdot\gamma(x)=He_{k}(-x)\gamma(x)$ (3)

We will utilize this fact to express the Hermite coefficients of
$f(x)=a^{\mathsf{T}}\sigma(W^{\mathsf{T}}x+b)$. Finally, for 1-dimensional
Hermite polynomials $He_{k}(x)$, and a sign variable $\xi\in\\{\pm 1\\}$, we
have $He_{k}(\xi x)=\xi^{k}He_{k}(x)$ from the odd/even function properties of
Hermite polynomials.

##### Tensor Decomposition

The tensor product $u\otimes v\otimes
w\in\mathbb{R}^{d_{1}}\otimes\mathbb{R}^{d_{2}}\otimes\mathbb{R}^{d_{3}}$ of
vectors $u\in\mathbb{R}^{d_{1}},v\in\mathbb{R}^{d_{2}},w\in\mathbb{R}^{d_{3}}$
is a rank-$1$ tensor. Similarly, we will use
$u^{\otimes\ell}\in(\mathbb{R}^{d})^{\otimes\ell}$ to denote the tensor
product of $u$ with itself $\ell$ times. An order-$t$ tensor
$T\in\mathbb{R}^{d_{1}}\otimes\mathbb{R}^{d_{2}}\otimes\dots\otimes\mathbb{R}^{d_{t}}$
is represented using a $t$-way array $\mathbb{R}^{d_{1}\times
d_{2}\times\dots\times d_{t}}$ that has $t$ modes corresponding to the $t$
different indices. Given two matrices $U,V$ with $k$ columns each given by
$U=(u_{i}:i\in[k])$ and $V=(v_{u}:i\in[k])$, the Khatri-Rao product $M=U\odot
V$ is a matrix formed by the $i$th column being $u_{i}\otimes v_{i}$. We will
also use $U^{\odot 2}=U\odot U$ (and similarly for higher orders). A claim
about preserving the full-column-rank property (and analogously minimum
singular value) under the Khatri-Rao product is included below.

Flattening or Reshaping:  Given an order-$t$ tensor $T$, for
$t_{1},t_{2},t_{3}\geq 0$ such that $t_{1}+t_{2}+t_{3}=t$ define
$T^{\prime}=\text{flatten}(T,t_{1},t_{2},t_{3})$ as the order-$3$ tensor
$T^{\prime}\in\mathbb{R}^{d^{t_{1}}\times d^{t_{2}}\times d^{t_{3}}}$,
obtained by flattening and combining the first $t_{1}$ modes, the next $t_{2}$
and last $t_{3}$ modes respectively. When $t_{3}=0$, the output is a matrix in
$\mathbb{R}^{d^{t_{1}}\times d^{t_{2}}}$.

Tensor decompositions of order $3$ and above, unlike matrix decompositions
(which are of order $2$) are known to be unique under mild conditions. While
tensor decompositions are NP-hard in the worst-case, polynomial time
algorithms for tensor decompositions are known under certain non-degeneracy
conditions (see e.g., [JGKA19, Vij20]). In particular, Jennrich’s algorithm
[Har70] provides polynomial time guarantees for recovering all the rank-$1$
terms of a decomposition of a tensor $T=\sum_{i=1}^{k}u_{i}\otimes
v_{i}\otimes w_{i}$, when the $\\{u_{i}:i\in[k]\\}$ are linearly independent,
the $\\{v_{i}:i\in[k]\\}$ are linearly independent, and no two of the vectors
$\\{w_{i}:i\in[k]\\}$ are parallel. This algorithm and its guarantee can also
be made robust to some noise (of an inverse polynomial magnitude), when
measured in Frobenius norm. In this paper, the following claims are especially
vital to formulate our main results.

######  Claim 1.4 (Implication of Lemma A.4 of [BCV14]222In fact one can show
that a stronger statement that a certain quantity called Kruskal-rank
increases, see [BCV14]. ).

Let $U\in\mathbb{R}^{d_{1}\times k},V\in\mathbb{R}^{d_{2}\times k}$, and
suppose the smallest column length $\min_{i\in[k]}\lVert
v_{i}\rVert_{2}\geq\kappa$. Then the Khatri-Rao product $U\odot V$ has rank
$k$ and satisfies $s_{k}(U\odot V)\geq\kappa\cdot s_{k}(U)/\sqrt{2k}$.

A robust guarantee we will use for Jennrich’s algorithm [Har70] is given below
(see also [GVX14, Moi18] for robust analysis).

######  Theorem 1.5 (Theorem 2.3 of [BCMV14]).

Suppose $\varepsilon_{\ref{thm:jennrich}}>0$ we are given tensor
$\widetilde{T}=T+E\in\mathbb{R}^{m\times n\times p}$, where $T$ has a
decomposition $T=\sum_{i=1}^{k}u_{i}\otimes v_{i}\otimes w_{i}$ satisfying the
following conditions:

  1. 1.

Matrices $U=(u_{i}:i\in[k]),V=(v_{i}:i\in[k])$ have condition number (ratio of
the maximum singular value $\sigma_{1}$ to the least singular value
$\sigma_{k}$) at most $\kappa$,

  2. 2.

For all $i\neq j$, the submatrix $W_{\\{i,j\\}}$ has
$s_{2}(W_{\\{i,j\\}})\geq\delta$.

  3. 3.

Each entry of $E$ is bounded by
$\eta_{\ref{thm:jennrich}}(\varepsilon,\kappa,\max\\{n,m,p\\},\delta)=\frac{\mathrm{poly}(\varepsilon)}{\mathsf{poly}(\kappa,\max\\{n,m,p\\},\tfrac{1}{\delta})}$.

Then there exists a polynomial time algorithm that on input $\widetilde{T}$
returns a decomposition
$\\{(\widetilde{u}_{i},\widetilde{v}_{i},\widetilde{w}_{i}):i\in[k]\\}$ s.t.
there is a permutation $p:[k]\to[k]$ with

$\forall
i\in[k],~{}~{}\|\widetilde{u}_{i}\otimes\widetilde{v}_{i}\otimes\widetilde{w}_{i}-u_{p(i)}\otimes
v_{p(i)}\otimes w_{p(i)}\|_{F}\leq\varepsilon_{\ref{thm:jennrich}}.$ (4)

##  2 Related Work

By now there is a vast literature exploring various aspects of deep learning
from a theoretical perspective. Here we discuss the works most relevant in the
context of our results. As discussed earlier, the recent works of [GLM18,
BJW19, GKLW19] provide polynomial time algorithms for learning depth-2
feedforward ReLU networks under the assumption that the input distribution is
Gaussian and that the matrix $W$ is full rank. Some of these works consider a
setting where the output is also a high dimensional vector [BJW19, GKLW19],
and also consider learning beyond the Gaussian distribution. However, these
works do not extend to the case of non-zero bias.

The work of [JSA15] proposed a general approach based on tensor decompositions
for learning an unknown depth-2 neural network that could also handle the
presence of the bias terms. The tensor used in the work of [JSA15] is formed
by taking the weighted average of a “score” function evaluated on each data
point. In this way their approach generalizes to a large class of
distributions provided one has access to the score function. However, for most
data distributions computing the score function is a hard task itself. When
the input distribution is Gaussian, then the score functions correspond to the
Hermite coefficients of the target function $f(x)$ and can be evaluated
efficiently. However, the analysis in [JSA15] does not extend to the case of
ReLU activations for several reasons. Their technique needs certain smoothness
and symmetry assumptions on the activations that do not hold for ReLU. These
assumptions also ensure that all the terms in the appropriate tensor are non-
zero. We do not make such assumptions, and tackle one of the main challenges
by showing that one can indeed recover a good approximation to the network by
analyzing multiple higher order tensors. Furthermore, the authors in [JSA15]
assume that the biases, and the spectral norm of $W$ are both bounded by a
constant. As a result they do not handle the case of biases of large magnitude
where some of the ReLU units mostly function as linear functions (with high
probability).

There have also been works on designing learning algorithms for neural
networks without assumptions on the linear independence of columns of $W$.
These results incur an exponential dependence on either the input
dimensionality or the number of parameters in the unknown network [DKKZ20,
CKM20], or quasipolynomial dependence when the coefficients
$\\{a_{i}:i\in[m]\\}$ are all non-negative [DK20]. In particular, the result
of [CKM20] provides a learning algorithm for arbitrary depth neural networks
under the Gaussian distribution with an “FPT” guarantee; its running time is
polynomial in the dimension, but exponential in the number of ReLU units.
Given the recent statistical query lower bounds on learning deep neural
networks [DKKZ20, GGJ+20], getting a fully polynomial time algorithm without
any assumptions is a challenging open problem, even under Gaussian marginals.

Polynomial time algorithms with fewer assumptions and beyond depth-2 can be
designed if the activation functions in the first hidden layer are sigmoid
functions [GK19]. Finally, there is also extensive literature on analyzing the
convergence properties of gradient descent and stochastic gradient descent for
neural networks. The results in this setting implicitly or explicitly assume
that the target function is well approximated in the Neural Tangent Kernel
(NTK) space of an unknown network. Under this assumption these results show
that gradient descent on massively overparameterized neural networks can learn
the target [APVZ14, Dan17, DFS16, AZLS19, DZPS19, ADH+19, LXS+19, CB18,
JGH21].

##  3 Main Results

There are two related but different goals that we consider in learning the
ReLU network:

  * •

Achieves low error: Output a ReLU network
$g(x)=a^{\prime\top}\sigma(W^{\prime\top}x+b^{\prime})$ such that the $L_{2}$
error is at most $\varepsilon$ for a given $\varepsilon>0$ i.e.,
$\operatorname*{\mathbb{E}}_{x\sim\mathcal{N}(0,I_{d\times
d})}[(f(x)-g(x))^{2}]\leq\varepsilon^{2}$.

  * •

Parameter recovery: Output $\widetilde{W},\widetilde{a},\widetilde{b}$, such
that each parameter is $\varepsilon$-close (up to permuting the $m$ co-
ordinates of $a,b\in\mathbb{R}^{m}$ and reordering the corresponding columns
of $W$).

We remark that the second goal is harder and implies the first; in particular,
when $\varepsilon=0$, the second goal corresponds to identifiability of the
model. However in some cases, parameter recovery may be impossible to achieve
(see later for some examples) even though we can achieve the goal of achieving
low error. As we have seen earlier, given $N$ samples if $b_{i}\gg\sqrt{\log
N}$, then $\sigma(w_{i}^{\top}x+b_{i})$ will be indistinguishable from the
linear function $w_{i}^{\top}x+b$ w.h.p.; hence if there are multiple such
$b_{i}\in[m]$ with large magnitude, the best we can hope to do is recover the
sum of all those linear terms. Our first result shows that this is the only
obstacle when we are in the full-rank or undercomplete setting i.e.,
$\\{w_{i}:i\in[m]\\}$ are linearly independent (in a robust sense).

######  Theorem 3.1 (Full-rank setting).

Suppose $\varepsilon\in(0,1)$ and
$N\geq\mathrm{poly}(m,d,1/\varepsilon,1/s_{m}(W),B)$ samples be generated by a
ReLU network $f(x)=a^{\top}\sigma(W^{\top}x+b)$ that is $B$-bounded, and
$|b_{i}|<c\sqrt{\log(1/\varepsilon mdB)}$ for all $i\in[m]$. Then there exists
an algorithm that runs in $\mathrm{poly}(N,m,d)$ time and with high
probability recovers $\widetilde{a}_{i}$, $\widetilde{b}_{i}$,
$\widetilde{w}_{i}$ such that $\lVert
w_{i}-\widetilde{w}_{i}\rVert_{2}+|a_{i}-\widetilde{a}_{i}|+|b_{i}-\widetilde{b}_{i}|<\varepsilon$
for all $i\in[m]$.

The above theorem recovers all the parameters when the biases
$\\{b_{i}:i\in[m]\\}$ of each ReLU unit does not have very large magnitude.
Moreover even when there are $b_{i}$ of large magnitude, we can learn a
depth-$2$ ReLU network $g$ that achieves low error, and simultaneously recover
parameters for the terms that have a small magnitude of $b_{i}$ (up to a
potential ambiguity in signs). In fact, our algorithm and guarantees are more
general, and can operate under the much milder condition that
$\\{w_{i}^{\otimes\ell}:i\in[m]\\}$ are linearly independent for any constant
$\ell\geq 1$; the setting when $\ell>1$ corresponds to what is often called
the overcomplete setting. In what follows, for any constant
$\ell\in\mathbb{N}$ we use $\mathrm{poly}_{\ell}(n_{1},n_{2},\dots)$ to denote
a polynomial dependency on $n_{1},n_{2},\dots$, and potentially exponential
dependence on $\ell$.

######  Theorem 3.2.

Suppose $\ell\in\mathbb{N}$ be a constant, and $\varepsilon>0$. If we are
given $N$ i.i.d. samples as described above from a ReLU network
$f(x)=a^{\top}\sigma(W^{\top}x+b)$ that is $B$-bounded then there is an
algorithm that given
$N\geq\mathrm{poly}_{\ell}(m,d,1/\varepsilon,1/s_{m}(W^{\odot\ell}),B)$ runs
in $\mathrm{poly}(N,m,d)$ time and with high probability finds a ReLU network
$g(x)={a^{\prime}}^{\top}\sigma({W^{\prime}}^{\top}x+{b^{\prime}})$ with at
most $m+2$ hidden units, such that the $L_{2}$ error
$\operatorname*{\mathbb{E}}_{x\sim\mathcal{N}(0,I_{d\times
d})}[(f(x)-{g}(x))^{2}]\leq\varepsilon^{2}$. Furthermore there are constants
$c=c(\ell)>0,c^{\prime}>0$ and signs $\xi_{i}\in\\{\pm 1\\}~{}\forall
i\in[m]$, such that in $\mathrm{poly}(N,m,d)$ time, for all $i\in[m]$ with
$|b_{i}|<c\sqrt{\log(1/(\varepsilon\cdot mdB))}$, we can recover
$(\widetilde{a}_{i},\widetilde{w}_{i},\widetilde{b}_{i})$, such that
$|a_{i}-\widetilde{a}_{i}|+\lVert
w_{i}-\xi_{i}\widetilde{w}_{i}\rVert_{2}+|b_{i}-\xi_{i}\widetilde{b}_{i}|<c^{\prime}\varepsilon/(mB)$.

In the special case of $\ell=1$ in Theorem 3.2, we need the least singular
value $s_{m}(W)>0$ (this necessitates that $m\leq d$). This corresponds to the
full-rank setting considered in Theorem 3.1. In contrast to the full-rank
setting, for $\ell>1$ we only require that the set of vectors
$w_{1}^{\otimes\ell},w_{2}^{\otimes\ell},\dots,w_{m}^{\otimes\ell}$ are
linearly independent (in a robust sense), which one can expect for much larger
values of $m$ typically. The following corollary formalizes this in the
smoothed analysis framework of Spielman and Teng [ST04], which is a popular
paradigm for reasoning about non-worst-case instances [Rou20]. Combining the
above theorem with existing results on smoothed analysis [BCPV19] implies
polynomial time learning guarantees for non-degenerate instances with
$m=O(d^{\ell})$ for any constant $\ell>0$. Below, $\widehat{W}$ denotes the
columns of $W$ are $\tau$-smoothed i.e., randomly perturbed with standard
Gaussian of average length $\tau$ that is at least inverse polynomial (See
Section 6 for the formal smoothed analysis model and result).

######  Corollary 3.3 (Smoothed Analysis).

Suppose $\ell\in\mathbb{N}$ and $\varepsilon>0$ are constants in the smoothed
analysis model with smoothing parameter $\tau>0$, and also assume the ReLU
network $f(x)=a^{\top}\sigma(\widehat{W}^{\top}x+b)$ is $B$-bounded with
$m\leq 0.99\binom{d+\ell-1}{\ell}$. Then there is an algorithm that given
$N\geq\mathrm{poly}_{\ell}(m,d,1/\varepsilon,B,1/\tau)$ samples runs in
$\mathrm{poly}(N,m,d)$ time and with high probability finds a ReLU network
$g(x)={a^{\prime}}^{\top}\sigma({W^{\prime}}^{\top}x+{b^{\prime}})$ with at
most $m+2$ hidden units, such that the $L_{2}$ error
$\operatorname*{\mathbb{E}}_{x\sim\mathcal{N}(0,I_{d\times
d})}[(f(x)-{g}(x))^{2}]\leq\varepsilon^{2}$. Furthermore there are constants
$c,c^{\prime}>0$ and signs $\xi_{i}\in\\{\pm 1\\}~{}\forall i\in[m]$, such
that in $\mathrm{poly}(N,m,d)$ time, for all $i\in[m]$ with
$|b_{i}|<c\sqrt{\log(1/(\varepsilon\cdot mdB))}$, we can recover
$(\widetilde{a}_{i},\widetilde{w}_{i},\widetilde{b}_{i})$, such that
$|a_{i}-\widetilde{a}_{i}|+\lVert
w_{i}-\xi_{i}\widetilde{w}_{i}\rVert_{2}+|b_{i}-\xi_{i}\widetilde{b}_{i}|<c^{\prime}\varepsilon/(mB)$.

While our algorithm and the analysis give guarantees that are robust to
sampling errors and inverse polynomial error, even the non-robust analysis has
implications and, implies identifiability of the model (up to ambiguity in the
signs) as long as no two rows of $W$ are parallel. Note that in general
identifiability may not imply any finite sample complexity bounds.

######  Theorem 3.4 (Partial Identifiability).

Suppose we are given samples from a ReLU network
$f(x)=a^{\top}\sigma(W^{\top}x+b)$ where $\min_{i\in[m]}|a_{i}|>0$ and no two
columns of $W$ are parallel to each other. Then given samples, the model
parameters $(a_{i},b_{i},w_{i}:i\in[m])$ are identified up to ambiguity in the
signs and reordering indices i.e., we can recover
$\\{(a_{i},\xi_{i}b_{i},\xi_{i}w_{i}):i\in[m]\\}$ for some
$\xi_{i}\in\\{+1,-1\\}~{}\forall i\in[m]$.
Moreover given any $(\xi_{i}\in\\{+1,-1\\}:i\in[m])$ such that

$\sum_{i=1}^{m}a_{i}\Phi(\xi_{i}b_{i})\xi_{i}w_{i}=\sum_{i=1}^{m}a_{i}\Phi(b_{i})w_{i},~{}~{}\text{
and
}~{}\sum_{i=1}^{m}a_{i}b_{i}\Phi(\xi_{i}b_{i})=\sum_{i=1}^{m}a_{i}b_{i}\Phi(b_{i}),$
(5)

we have that the set of parameters
$((a_{i},\xi_{i}b_{i},\xi_{i}w_{i}:i\in[m])$ also gives rise to the same
distribution.

The above theorem shows that under a very mild assumption on $W$, the
parameters can be identified up to signs. However, this ambiguity in the signs
may be unavoidable – the second part of the Theorem 3.4 shows that any
combination of signs that match the zeroth and first Hermite coefficient gives
rise to a valid solution (this corresponds to the $d+1$ equations in (5)).
Even in the case when all the $b_{i}=0$, we have non-identifiability due to
ambiguities in signs whenever the $\\{w_{i}:i\in[m]\\}$ are not linearly
independent for an appropriate setting of the $\\{a_{i}\\}$; see Claim 4.8 for
a formal statement. On the other hand, Theorem 4.1 gives unique
identifiability result in the full-rank setting (as there is only one setting
of the signs that match the first Hermite coefficient in the full-rank
setting).

Our results rely on the precise expressions for higher order Hermite
coefficients of $f(x)$ given below.

######  Lemma 3.5.

Let $\hat{f}_{k}=\mathbb{E}_{x\sim\mathcal{N}(0,I)}[f(x)He_{k}(x)]$ (with
$k\in\mathbb{N}$) be the $k$’th Hermite coefficient (this is an order-$k$
tensor) of $f(x)=a^{\top}\sigma(W^{\top}x+b)$. Then

$\displaystyle\hat{f}_{0}$
$\displaystyle=\sum_{i=1}^{m}a_{i}\Big{(}b_{i}\Phi(b_{i})+\frac{\exp(-\frac{b_{i}^{2}}{2})}{\sqrt{2\pi}}\Big{)}~{},~{}\hat{f}_{1}=\sum_{i=1}^{m}a_{i}\Phi(b_{i})w_{i}$
(6) $\displaystyle\forall k\geq 2,~{}~{}\hat{f}_{k}$
$\displaystyle=\sum_{i=1}^{m}(-1)^{k}\cdot a_{i}\cdot
He_{k-2}(b_{i})\cdot\frac{\exp(\frac{-b_{i}^{2}}{2})}{\sqrt{2\pi}}\cdot
w_{i}^{\otimes k}$ (7)

We prove this by considering higher order derivatives and using properties of
Hermite polynomials. A key property we use here is that the $k$’th derivative
of a standard Gaussian function is itself multiplied by the $k$’th Hermite
polynomial (with sign flipped for odd $k$). This significantly simplifies the
expression for the coefficient $g(b_{i})$ of $w_{i}^{\otimes k}$.

We remark that the above lemma may also be used to give an expression for the
training objective for depth-2 ReLU networks, analogous to the result of
[GLM18] for ReLU activations with no bias, that provides an expression as a
combination of tensor decomposition problems of increasing order. The authors
in [GLM18] crucially use the form of the decomposition to design a new
regularized objective on which the convergence of gradient descent can be
analyzed. The decomposition presented below for general ReLU activations opens
the door for analyzing gradient descent in the non-zero bias setting.

######  Proposition 3.6.

Let
$\widetilde{f}(x)=\widetilde{a}^{\top}\sigma(\widetilde{W}^{\top}x+\widetilde{b})$
be the model trained using samples generated by the ground-truth ReLU network
$f(x)=a^{\top}\sigma(W^{\top}x+b)$. Then the statistical risk with respect to
the $\ell_{2}$ loss function can be expressed as follows

$\displaystyle L(\widetilde{a},\widetilde{b},\widetilde{W})$
$\displaystyle=\sum_{k\in\mathbb{N}}\frac{1}{k!}\Big{\|}T_{k}-\hat{f}_{k}\Big{\|}_{F}^{2}$
$\displaystyle~{}~{}\text{ where }~{}~{}T_{0}$
$\displaystyle=\sum_{i=1}^{m}\widetilde{a}_{i}(\widetilde{b}_{i}\Phi(\widetilde{b}_{i})+\frac{\exp(-\widetilde{b}_{i}^{2}/2)}{\sqrt{2\pi}}),~{}\text{
and
}~{}T_{1}=\sum_{i=1}^{m}\widetilde{a}_{i}\Phi(\widetilde{b}_{i})\widetilde{w}_{i}$
$\displaystyle\forall k\geq 2,~{}~{}T_{k}$
$\displaystyle=\sum_{i=1}^{m}(-1)^{k}\cdot\widetilde{a}_{i}\cdot
He_{k-2}(\widetilde{b_{i}})\cdot\frac{\exp(-\widetilde{b_{i}}^{2}/2)}{\sqrt{2\pi}}\cdot\widetilde{w}_{i}^{\otimes
k}$

Please refer to Appendix A for the proof.

Observe that by setting $b_{i}=\widetilde{b}_{i}=0$ in our above expression,
when $k\geq 2$, we immediately recover the objective function given in Theorem
2.1 of [GLM18] as

$\sum_{\begin{subarray}{c}k\geq 2\end{subarray}}\frac{He_{k-2}(0)^{2}}{2\pi
k!}\Big{\|}\sum_{i=1}^{m}\widetilde{a}_{i}\widetilde{w}_{i}^{\otimes
k}-\sum_{i=1}^{m}a_{i}w_{i}^{\otimes
k}\Big{\|}_{F}^{2}=\sum_{\begin{subarray}{c}k\geq 2\\\ k~{}\text{is
even}\end{subarray}}\frac{((k-3)!!)^{2}}{2\pi
k!}\Big{\|}\sum_{i=1}^{m}\widetilde{a}_{i}\widetilde{w}_{i}^{\otimes
k}-\sum_{i=1}^{m}a_{i}w_{i}^{\otimes k}\Big{\|}_{F}^{2}$ (8)

for $k=0,1$, we also have

$\Big{|}T_{0}-\hat{f}_{0}\Big{|}^{2}+\Big{\|}T_{1}-\hat{f}_{1}\Big{\|}_{2}^{2}=\frac{1}{2\pi}\Big{|}\sum_{i=1}^{m}\widetilde{a}_{i}-\sum_{i=1}^{m}a_{i}\Big{|}^{2}+\frac{1}{4}\Big{\|}\sum_{i=1}^{m}\widetilde{a}_{i}\widetilde{w}_{i}-\sum_{i=1}^{m}a_{i}w_{i}\Big{\|}_{2}^{2}$
(9)

Furthermore, note that the Hermite coefficients of the ReLU function
$\sigma(x)$ are $\hat{\sigma}_{0}=1/\sqrt{2\pi}$, $\hat{\sigma}_{1}=1/2$ and
$\hat{\sigma}_{k}=(-1)^{\frac{k-2}{2}}(k-3)!!/\sqrt{2\pi k!}$ for $k\geq 2$
and $k$ being even; otherwise $\hat{\sigma}_{k}=0$.

##  4 Non-robust Algorithm and Analysis

Our algorithms for learning the parameters of
$f(x)=a^{\top}\sigma(W^{\top}x+b)$ decompose tensors obtained from the Hermite
coefficients $\\{\hat{f}_{t}\in(\mathbb{R}^{d})^{\otimes t}\\}$ of the
function $f$. In this section, we design an algorithm assuming that we have
access to all the necessary Hermite coefficients exactly (no noise or sampling
errors). This will illustrate the basic algorithmic ideas and the
identifiability result. However with polynomial samples, we can only hope to
estimate these quantities up to inverse polynomial accuracy. In Section 5 we
describe how we deal with the challenges that arise from errors.

Our first result is a polynomial time algorithm in the full-rank setting that
recovers all the parameters exactly.

######  Theorem 4.1 (Full-rank non-robust setting).

Suppose the parameters $\\{(a_{i},b_{i},w_{i}):i\in[m]\\}$ satisfies: (i)
$a_{i}\neq 0$ for all $i\in[m]$, (ii) $\\{w_{i}:i\in[m]\\}$ are linearly
independent. Then given $\\{\hat{f}_{t}:0\leq t\leq 4\\}$ exactly, Algorithm 1
recovers (with probability $1$) the unknown parameters $a$, $b$ and $W$ in
$\mathrm{poly}(m,d)$ time.

See Theorem 3.1 for the analogous theorem in the presence of errors in
estimating the Hermite coefficients $\\{\hat{f}_{t}\\}$. Our algorithm for
recovering the parameters estimates different Hermite coefficient tensors
$\\{\hat{f}_{t}:0\leq t\leq 4\\}$ and uses tensor decomposition algorithms on
these tensors to first find the $\\{w_{i}:i\in[m]\\}$ up to some ambiguity in
signs. We can also recover all the coefficients $\\{b_{i}:i\in[m]\\}$ up to
signs (corresponding to the signs of $w_{i}$ ), and all the
$\\{a_{i}:i\in[m]\\}$ (no sign ambiguities). This portion of the algorithm
extends to higher order $\ell$, under a weaker assumption on the matrix $W$.

######  Theorem 4.2.

Suppose the parameters $\\{(a_{i},b_{i},w_{i}):i\in[m]\\}$ satisfies: (i) no
two $\\{w_{i}:i\in[m]\\}$ are linearly dependent and, (ii) for a constant
$\ell\in\mathbb{N}$, $\\{w_{i}^{\otimes\ell}:i\in[m]\\}$ are linearly
independent, (iii) $a_{i}\neq 0$ for all $i\in[m]$. Then given
$\\{\hat{f}_{t}:0\leq t\leq 2\ell+2\\}$ exactly, Algorithm 1 in
$\mathrm{poly}_{\ell}(m,d)$ time outputs (with probability $1$)
$\\{\hat{a}_{i},\hat{w}_{i},\hat{b}_{i}:i\in[m]\\}$ such that we can recover
the parameters up to a reordering of the indices $[m]$ and up to signs i.e.,
for some $\\{\xi_{i}\in\\{-1,1\\}:i\in[m]\\}$ we have $\hat{a}_{i}=a_{i}$,
$\hat{w}_{i}=\xi_{i}w_{i}$ and $\hat{b}_{i}=\xi_{i}b_{i}$. Furthermore, given
exact statistical query access to the distribution $\mathcal{N}(0,I_{d\times
d})$,333This means that for any function $h(x,y)$ that can be computed in
polynomial time, one can obtain $\mathbb{E}_{x,y}[h(x,y)]$ exactly. there
exists an algorithm that runs in time $\mathrm{poly}(m,d,B)$ and outputs a
function $g(x)$ such that $\mathbb{E}_{x\sim\mathcal{N}(0,I_{d\times
d})}\big{(}f(x)-g(x)\big{)}^{2}=0$.

We now describe the algorithm for general $\ell\geq 1$ (this specializes to
the full-rank setting for $\ell=1$).

Input:
$\hat{f}_{\ell},\hat{f}_{\ell+1},\hat{f}_{\ell+2},\hat{f}_{\ell+3},\hat{f}_{2\ell+1},\hat{f}_{2\ell+2}$;

1\. Let
$T^{\prime}=\text{flatten}(\hat{f}_{2\ell+1},\ell,\ell,1)\in\mathbb{R}^{d^{\ell}\times
d^{\ell}\times d}$ and
$T^{\prime\prime}=\text{flatten}(\hat{f}_{2\ell+2},\ell,\ell,2)\in\mathbb{R}^{d^{\ell}\times
d^{\ell}\times d^{2}}$ be order-3 tensors obtained by flattening
$\hat{f}_{2\ell+1}$ and $\hat{f}_{2\ell+2}$.

2\. Set
$k^{\prime}=\text{rank}(\text{flatten}(\hat{f}_{2\ell+1},\ell,\ell+1,0))$. Run
Jennrich’s algorithm [Har70] on $T^{\prime}$ to recover rank-1 terms
$\\{\alpha^{\prime}_{i}u_{i}^{\otimes\ell}\otimes u_{i}^{\otimes\ell}\otimes
u_{i}~{}|~{}i\in[k^{\prime}]\\}$, where $\forall
i\in[k^{\prime}],~{}u_{i}\in\mathbb{S}^{d-1}$ and
$\alpha^{\prime}_{i}\in\mathbb{R}$.

3\. Set
$k^{\prime\prime}=\text{rank}(\text{flatten}(\hat{f}_{2\ell+2},\ell,\ell+1,0))$.
Run Jennrich’s algorithm [Har70] on $T^{\prime\prime}$ to recover rank-1 terms
$\\{\alpha^{\prime\prime}_{i}v_{i}^{\otimes\ell}\otimes
v_{i}^{\otimes\ell}\otimes v_{i}^{\otimes 2}~{}|~{}i\in[k^{\prime\prime}]\\}$,
where $\forall i\in[k^{\prime\prime}],~{}v_{i}\in\mathbb{S}^{d-1}$ and
$\alpha^{\prime\prime}_{i}\in\mathbb{R}$.

4\. Remove duplicates and negations (i.e., antipodal pairs of the form $v$ and
$-v$) from
$\\{u_{1},u_{2},\dots,u_{k^{\prime}}\\}\cup\\{v_{1},v_{2},\dots,v_{k^{\prime\prime}}\\}$
to get $\widetilde{w}_{1},\widetilde{w}_{2},\dots,\widetilde{w}_{m}$.

5\. Run subroutine
RecoverScalars$(m,\ell,\\{\widetilde{w}_{i}:i\in[m]\\},\hat{f}_{\ell},\hat{f}_{\ell+1},\hat{f}_{\ell+2},\hat{f}_{\ell+3})$
to get $\\{\widetilde{a}_{i},\widetilde{b}_{i}:i\in[m]\\}$.

6\. If $\ell=1$ (full-rank setting), run Algorithm 3 (FixSigns) on parameters
$m,\hat{f}_{1}$ and
$(\widetilde{a}_{i},\widetilde{b}_{i},\widetilde{w}_{i}:i\in[m])$ to get
$(a^{\prime}_{i}=\widetilde{a}_{i},b^{\prime}_{i},w^{\prime}_{i}:i\in[m])$.

Result: Output
$\\{\widetilde{w}_{i},\widetilde{a}_{i},\widetilde{b}_{i}:i\in[m]\\}$

Algorithm 1 Algorithm for order $\ell$: recover $a$, $b$, $W$ given
$\\{\hat{f}_{t}:0\leq t\leq 2\ell+2\\}$

Subroutine Algorithm 2 finds the unknown parameters
$a_{1},\dots,a_{m}\in\mathbb{R}$ and $b_{1},\dots,b_{m}\in\mathbb{R}$ given
$w_{1},w_{2},\dots,w_{m}$. While Algorithm 1 changes a little when we have
errors in the estimates, the subroutine Algorithm 2 remains the same even for
the robust version of the algorithm.

Input: $m,\ell,(\widetilde{w}_{i}:i\in[m])$, and tensors
$T_{\ell},T_{\ell+1},T_{\ell+2},T_{\ell+3}$ which are tensors of orders
$\ell,\ell+1,\ell+2,\ell+3$ respectively;

if  _$\ell=1$_ then

1\. For $j\in\\{2,3\\}$, solve the system of linear equation
$\sum_{i=1}^{m}\zeta_{j}(i)\mathsf{vec}(\widetilde{w}_{i}^{\otimes
j})=\mathsf{vec}(T_{j})$ to recover unknowns
$\\{\zeta_{j}(i)~{}|~{}i\in[m]\\}$;

2\. For each $i\in[m]$, set $b_{i}=-\frac{\zeta_{3}(i)}{\zeta_{2}(i)}$, and
$a_{i}=\zeta_{2}(i)\cdot\sqrt{2\pi}e^{b_{i}^{2}/2}$.

else

1\. For $j\in\\{\ell,\ell+1,\ell+2,\ell+3\\}$, solve the system of linear
equation $\sum_{i=1}^{m}\zeta_{j}(i)\mathsf{vec}(\widetilde{w}_{i}^{\otimes
j})=\mathsf{vec}(T_{j})$ to recover unknowns
$\\{\zeta_{j}(i)~{}|~{}i\in[m]\\}$;

2\. For each $i\in[m]$,
$q_{i}:=\operatorname*{argmax}_{\begin{subarray}{c}j\in\\{\ell+1,\ell+2\\}\end{subarray}}|\zeta_{j}(i)|$,
set
$\widetilde{b}_{i}=-\frac{\zeta_{q+1}(i)+q\cdot\zeta_{q-1}(i)}{\gamma_{q}(i)}$,
and
$\widetilde{a}_{i}=\sqrt{2\pi}(-1)^{q}\zeta_{q}(i)e^{b_{i}^{2}/2}/He_{q}(b_{i})$,
as described in (12) (Lemma 4.5).

end if

Result: Output $(\widetilde{a}_{i},\widetilde{b}_{i}:i\in[m])$

Algorithm 2 Subroutine RecoverScalars to recover $b_{i}$ (up to signs) and
$a_{i}$ given $w_{i}$ (up to signs) for all $i\in[m]$.

The above two algorithms together recover for all $i\in[m]$, the $a_{i}$ and
up to a sign the $w_{i}$ and $b_{i}$. In the special case of $\ell=1$ which we
refer to as the full-rank setting, we can also recover the correct signs, and
hence recover all the parameters.

Input: $m$, $\hat{f}_{1}\in\mathbb{R}^{d}$ and estimates
$\tilde{a}_{i},\tilde{b}_{i},\tilde{w}_{i}$ for each $i\in[m]$;

1\. Solve the system of linear equation
$\sum_{i=1}^{m}z_{i}\tilde{a}_{i}\tilde{w}_{i}=\hat{f}_{1}$ to recover
unknowns $\\{z_{i}~{}|~{}i\in[m]\\}$;

2\. Set $\tilde{\xi}_{i}=\text{sign}(z_{i})$ for each $i\in[m]$.

Result: Output
$(\tilde{a}_{i},\tilde{\xi}_{i}\tilde{b}_{i},\tilde{\xi}_{i}\tilde{w}_{i}:i\in[m])$

Algorithm 3 Algorithm Fix signs in full-rank setting.

Algorithm 1 decomposes two different tensors obtained from consecutive Hermite
coefficients $\hat{f}_{2\ell+1},\hat{f}_{2\ell+2}$ to obtain the $\\{w_{i}\\}$
up to signs. We use two different tensors because the bias $b_{i}$ could make
the coefficient of the $i$th term in the decomposition $0$ (e.g.,
$He_{2\ell-1}(b_{i})=0$ for $\hat{f}_{2\ell+1}$); hence $w_{i}$ cannot be
recovered by decomposing $\hat{f}_{2\ell+1}$. Hence $\hat{f}_{2\ell+1}$ can
degenerate to a rank $m^{\prime}<m$ tensor, and Jennrich’s algorithm will
return only $m^{\prime}<m$ eigenvectors that correspond to non-zero
eigenvalues.

The following lemma addresses this issue by showing that two consecutive
Hermite polynomials can not both take small values at any point
$x\in\mathbb{R}$. This implies a separation between roots of consecutive
Hermite polynomials or , and establishes a “robust” version that will be
useful in Section 5. Moreover, this lemma also shows that when $|x|$ is not
close to $0$, at least one out of every two consecutive odd Hermite
polynomials takes a value of large magnitude at $x$.

######  Lemma 4.3 (Separation of Roots).

For all $k\in\mathbb{N},x\in\mathbb{R}$,
$\max\\{|He_{k}(x)|,|He_{k+1}(x)|\\}\geq\sqrt{k!/2}$.

###### Proof.

First, $\forall k\in\mathbb{N}$, by Turán’s inequality [Tur50] we have

$He_{k+1}^{2}(x)-He_{k}(x)He_{k+2}(x)=k!\cdot\sum_{i=0}^{k}\frac{He_{i}(x)^{2}}{i!}>0$
(10)

Set $\varepsilon=\sqrt{k!/2}$ and assume for contradiction that
$|He_{k+1}(x)|,|He_{k+2}(x)|<\varepsilon$. The LHS of (10) is at most
$\varepsilon^{2}+\varepsilon|He_{k}(x)|$, and the RHS of (10) is at least

$k!\cdot(1+x^{2}+...+\frac{He_{k}(x)^{2}}{k!})>k!\cdot(1+\frac{He_{k}(x)^{2}}{k!})$
(11)

Therefore, if $|He_{k}(x)|=t$, combining both sides we get
$\varepsilon^{2}+\varepsilon t\geq k!+t^{2}$. This implies on the one hand
that $\varepsilon t\geq k!/2$, and on the other hand that $\varepsilon t\geq
t^{2}$. However for our choices of $\varepsilon=\sqrt{k!/2}$, no value of $t$
is feasible. This yields the required contradiction for the first claim.

∎

The following claim shows that Jennrich’s algorithm for decomposing a tensor
successfully recovers all the rank-$1$ terms whose appropriate
$He_{k}(b_{i})\neq 0$. This claim along with Lemma 4.3 shows that Steps 2-3 of
Algorithm 1 successfully recovers all the $\\{w_{i}:i\in[m]\\}$ up to signs.

######  Claim 4.4.

Let $\ell_{1},\ell_{2}\geq\ell$ and $T\in\mathbb{R}^{d^{\ell_{1}}\times
d^{\ell_{2}}\times d^{\ell_{3}}}$ have a decomposition
$T=\sum_{i=1}^{m}\alpha_{i}w_{i}^{\otimes\ell_{1}}\otimes
w_{i}^{\otimes\ell_{2}}\otimes w_{i}^{\otimes\ell_{3}}$, with
$\\{w_{i}^{\otimes\ell}:i\in[m]\\}$ being linearly independent. Consider
matrix
$M=\text{flatten}(T,\ell_{1},\ell_{2}+\ell_{3},0)\in\mathbb{R}^{d^{\ell_{1}}\times
d^{\ell_{2}+\ell_{3}}}$, and let $r:=rank(M)$. Then Jennrich’s algorithm
applied with rank $r$ runs in
$\mathrm{poly}_{\ell_{1}+\ell_{2}+\ell_{3}}(m,d)$ time recovers (w.p. $1$) the
rank-$1$ terms corresponding to $\\{i\in[m]:|\alpha_{i}|>0\\}$. Moreover for
each $i$ with $|\alpha_{i}|>0$, we have $\tilde{w}_{i}=\xi_{i}w_{i}$ for some
$\xi_{i}\in\\{+1,-1\\}$.

###### Proof.

Let $Q=\\{i\in[m]:|\alpha_{i}|>0\\}$. Firstly $r=rank(M)=|Q|$, since $M$ has a
decomposition

$M=\sum_{i\in
Q}\alpha_{i}\big{(}w_{i}^{\otimes\ell_{1}}\big{)}\big{(}w_{i}^{\otimes\ell_{2}+\ell_{3}}\big{)}^{\top}=M_{1}\text{diag}(\alpha_{Q})M_{2}^{\top}$

where $M_{1},M_{2},\text{diag}(\alpha_{Q})$ all have full column rank $|Q|$.
Secondly, from assumption (ii) of Theorem 4.2 and Claim 1.4 applied with
$\ell_{1},\ell_{2}\geq\ell$, we have that
$\\{w_{i}^{\otimes\ell^{\prime}}:i\in[m]\\}$ are linearly independent for
every $\ell^{\prime}\geq\ell$. Hence the the factor matrices
$U=(w_{i}^{\otimes\ell_{1}}:i\in Q)$ and $V=(w_{i}^{\otimes\ell_{2}}:i\in Q)$
also have full column rank. Similarly from (iii) no two vectors in
$\\{\alpha_{i}w_{i}^{\otimes\ell_{3}}:i\in Q\\}$ are parallel. Hence, they
satisfy the conditions of Jennrich’s algorithm. Since there is no error in the
tensor, Jennrich’s algorithm (Theorem 1.5) succeeds with probability $1$ (see
[Vij20]). Finally since each rank-1 term is recovered exactly when
$\alpha_{i}\neq 0$, the vector in $\mathbb{R}^{d}$ obtained from the term will
correspond to either $w_{i}$ or $-w_{i}$ as required. ∎

The above claim was useful in recovering $w_{i}$ up to a sign ambiguity. The
following lemma is useful for recovering $a_{i}$ parameters (no sign
ambiguities) and the $b_{i}$ parameters up to sign ambiguity, once we have
recovered the $w_{i}$ up to sign ambiguity. It uses various properties of
Hermite polynomials along with Lemma 3.5 and Lemma 4.3.

######  Lemma 4.5.

Suppose $k\in\mathbb{N},k\geq 2$. Suppose for some unknowns
$\beta,z\in\mathbb{R}$ with $\beta\neq 0$, we are given values of
$\gamma_{j}=(-1)^{j}\xi^{j}\beta He_{j}(z)~{}\forall j\in\\{k,k+1,k+2,k+3\\}$
for some $\xi\in\\{+1,-1\\}$. Then $z,\beta$ are uniquely determined by

$\text{For
}q:=\operatorname*{argmax}_{\begin{subarray}{c}j\in\\{k+1,k+2\\}\end{subarray}}|\gamma_{j}|,~{}~{}\xi
z=-\frac{\gamma_{q+1}+q\cdot\gamma_{q-1}}{\gamma_{q}},~{}~{}~{}~{}\beta=(-1)^{q}\frac{\gamma_{q}}{He_{q}(\xi
z)}$ (12)

###### Proof.

We use the following fact about Hermite polynomials:

$He_{r+1}(z)=zHe_{r}(z)-r\cdot He_{r-1}(z).$ (13)

From Lemma 4.3, we know that $\max\\{|He_{k+1}(z)|,|He_{k+2}(z)|\\}>0$ and
hence $\gamma_{q}\neq 0$. Substituting in the recurrence (13) with $r=q$,

$\displaystyle z$ $\displaystyle=\frac{He_{q+1}(z)+q\cdot
He_{q-1}(z)}{He_{q}(z)}=\frac{\beta He_{q+1}(z)+q\cdot\beta He_{q-1}(z)}{\beta
He_{q}(z)}=-\xi\Big{(}\frac{\gamma_{q+1}+q\cdot\gamma_{q-1}}{\gamma_{q}}\Big{)},$

where we used the fact that the Hermite polynomials are odd functions for odd
$q$ and even polynomials for even $q$. The $\beta$ value is also recovered
since $\gamma_{q}=(-1)^{q}\beta(\xi^{q}He_{q}(z))=(-1)^{q}\beta\cdot
He_{q}(\xi z)$.

∎

A robust version of this lemma (see Section B.2.2) will be important in the
robust analysis of Section 5. The following claim applies the above lemma for
each $i\in[m]$ with $\xi_{i}z=b_{i}$ and
$\beta=a_{i}e^{-b_{i}^{2}/2}/\sqrt{2\pi}$, to show that Step 5 of the
algorithm recovers the correct $\\{(a_{i},\xi_{i}b_{i})\\}_{i\in[m]}$ given
the $\\{\xi_{i}w_{i}\\}_{i\in[m]}$.

######  Claim 4.6.

Given $\\{\widetilde{w}_{i}=\xi_{i}w_{i}:i\in[m]\\}$ where
$\xi_{i}\in\\{+1,1\\}~{}\forall i\in[m]$, Step 5 of Alg. 1 recovers
$\\{(a_{i},\xi_{i}b_{i},\xi_{i}w_{i}):i\in[m]\\}$.

###### Proof.

We first prove for $\ell\geq 2$. For each of the
$j\in\\{\ell,\ell+1,\ell+2,\ell+2\\}$, we have from Lemma 3.5 and the Hermite
polynomials $He_{j}$ being odd functions for odd $j$ and even functions for
even $j$,

$\hat{f}_{j}=\sum_{i=1}^{m}(-1)^{j}\cdot
a_{i}\cdot\cdot\frac{He_{j-2}(b_{i})\exp(-b_{i}^{2}/2)}{\sqrt{2\pi}}\cdot
w_{i}^{\otimes j}=\sum_{i=1}^{m}(-1)^{j}\cdot
a_{i}\cdot\frac{He_{j-2}(\xi_{i}b_{i})\exp(-b_{i}^{2}/2)}{\sqrt{2\pi}}\cdot\widetilde{w}_{i}^{\otimes
j}.$

Moreover the vectors $\\{\xi_{i}^{j}w_{i}^{\otimes j}:i\in[m]\\}$ are linearly
independent by assumption (and from Claim 1.4 for $j>\ell$). Hence the linear
system for each $j$ has a unique solution

$\forall\ell\leq j\leq\ell+3,~{}\forall i\in[m],~{}\text{ we have
}\zeta_{j}(i)=a_{i}\cdot\frac{(-1)^{j}}{\sqrt{2\pi}}\exp(-b_{i}^{2}/2)~{}He_{j-2}(\xi_{i}b_{i}).$

Lemma 4.5 applied with $\beta=\frac{1}{\sqrt{2\pi}}a_{i}e^{-b_{i}^{2}/2}$ and
$\gamma_{j-2}=(-1)^{j}He_{j-2}(\xi_{i}b_{i})$ (note $\beta\neq 0$) proves that
Alg. 2 recovers $a_{i},\xi_{i}b_{i}$.

For $\ell=1$, we note $\forall z\in\mathbb{R},~{}He_{0}(z)=1$, and
$He_{1}(z)=z$. From Lemma 3.5, we see that one set of solutions to the linear
system is

$\forall
i\in[m],~{}~{}\zeta_{2}(i)=a_{i}\frac{\exp(-b_{i}^{2}/2)}{\sqrt{2\pi}},~{}~{}\text{
and
}~{}~{}\zeta_{3}(i)=-a_{i}\xi_{i}b_{i}\cdot\frac{\exp(-b_{i}^{2}/2)}{\sqrt{2\pi}}.$

Moreover the vectors $\\{\xi_{i}w_{i}:i\in[m]\\}$ are linearly independent.
Hence, $\zeta_{2},\zeta_{3}$ are the unique solutions to the system. Hence
Algorithm 2 recovers $a_{i},\xi_{i}b_{i}$ as claimed. ∎

We now complete the proof of the non-robust analysis for any constant
$\ell\geq 0$.

###### Proof of Theorem 4.2.

The proof follows by combining Claim 4.4 and Claim 4.6, along with Lemma 4.3.
Let $Q_{1}=\\{i:|He_{2\ell-1}(b_{i})|>0\\}$ and
$Q_{2}=\\{i:|He_{2\ell+1}(b_{i})|>0\\}$. From Claim 4.4, Step 2 recovers all
the rank-$1$ terms in $Q_{1}$ with probability $1$; hence we obtain in
particular $\\{\xi_{i}w_{i}\in\mathbb{S}^{d-1}~{}|~{}i\in Q_{1}\\}$ for some
signs $\xi_{i}\in\\{+1,-1\\}$. Similarly, in Step 3 we recover w.p. $1$, the
$\\{\xi_{i}w_{i}\in\mathbb{S}^{d-1}~{}|~{}i\in Q_{2}\\}$ for some
$\xi_{i}\in\\{1,-1\\}$.

From Lemma 4.3, we know that no $x\in\mathbb{R}$ is a simultaneous root of
$He_{2\ell+1}(x),He_{2\ell+2}(x)$. Hence $Q_{1}\cup Q_{2}=\\{1,2,\dots,m\\}$.
Thus we obtain $\\{\xi_{i}w_{i}\in\mathbb{S}^{d-1}:i\in[m]\\}$ in Step 4 for
some signs $\xi_{i}\in\\{1,-1\\}$ for all $i\in[m]$. Finally using Claim 4.6,
we recover for each $i\in[m]$, the $a_{i},\xi_{i}b_{i}\in\mathbb{R}$
corresponding to $\xi_{i}w_{i}$.

Next, in order to recover a function $g(x)$ of zero $L_{2}$ error we set up a
linear regression problem. Given $x\in\mathbb{R}^{d}$ consider a $2m$
dimensional feature space $\phi(x)$ where
$\phi(x)_{2i}=a_{i}\sigma(\xi_{i}w^{\top}_{i}x+\xi_{i}b_{i})$ and
$\phi(x)_{2i+1}=a_{i}\sigma(-\xi_{i}w^{\top}_{i}x-\xi_{i}b_{i})$. Then is is
easy to see that the target network $f(x)$ can be equivalently written as
$f(x)={\beta^{*}}^{\top}\phi(x)$ for some vector $\beta^{*}$. Hence we can
recover another vector $\beta$ of zero $L_{2}$ error by solving ordinary least
squares, i.e,
$\beta=\mathbb{E}[\phi(x)^{\top}\phi(x)]^{\dagger}\mathbb{E}[\phi(x)y]$.444We
remark that using Claim B.5, we can further consolidate the terms to get a
ReLU network with at most $m+2$ hidden units.  Notice that both the
expectations can be calculated exactly given exact statistical query access to
the data distribution. In Section 5.3 we provide a more general analysis of
the above argument with finite sample analysis that will let us approximate
$f(x)$ up to arbitrary accuracy in the presence of sampling errors. ∎

We now complete the proof of recovery in the full-rank setting.

###### Proof of Theorem 4.1.

We first apply Theorem 4.2 (and its above proof) with $\ell=1$. We note that
the conditions are satisfied since $s_{\min}(W)>0$ and all the $a_{i}\neq 0$.
Theorem 4.2 guarantees that the first 5 steps of Algorithm 1 recovers (with
probability $1$) for each $i\in[m]$, $a_{i},\widetilde{b}_{i}=\xi_{i}b_{i}$
and $\widetilde{w}_{i}=\xi_{i}w_{i}$ for some $\xi_{i}\in\\{+1,1\\}$. From
Lemma 3.5, we have that

$\hat{f}_{1}=\sum_{i=1}^{m}a_{i}\Phi(b_{i})w_{i}=\sum_{i=1}^{m}z^{*}_{i}a_{i}\widetilde{w}_{i},~{}~{}\text{for}~{}~{}z^{*}_{i}=\xi_{i}\Phi(b_{i})~{}\forall
i\in[m],$

and $\Phi(b_{i})$ is the Gaussian CDF and restricted to $(0,1)$. Moreover the
$\\{w_{i}:i\in[m]\\}$ are linearly independent. Hence there is a unique
solution $(z^{*}_{i}:i\in[m])$ to the system of linear equations in the
unknowns $(z_{i}:i\in[m])$ in step 6 of Algorithm 1, and
$\widetilde{\xi}_{i}=\xi_{i}$ as required. Hence
$(a_{i},\tilde{\xi}_{i}\widetilde{b}_{i},\tilde{\xi}_{i}\widetilde{w}_{i}:i\in[m])$
are the true parameters of the network (up to reordering indices). ∎

Proof of Identifiability (Theorem 3.4): Theorem 3.4 follows by verifying that
the conditions of Theorem 4.2 hold for $\ell=m$. Conditions (i) and (iii)
follow from the conditions of Theorem 3.4.

We now verify condition (ii). For a matrix $U$ with columns
$\\{u_{i}:i\in[m]\\}$, the $\mathsf{krank}(U)$ (denoting the Kruskal-rank) is
at least $k$ iff every $k$ of the $m$ columns of $U$ are linearly independent.
Note that $\mathsf{krank}(U)\leq\mathsf{rank}(U)\leq m$. The $\mathsf{krank}$
increases under the Khatri-Rao product.

######  Fact 4.7 (Lemma A.4 of [BCV14]).

For two matrices $U,V$ with $m$ columns, $\mathsf{krank}(U\odot
V)=\min(\mathsf{krank}(U)+\mathsf{krank}(V)-1,m)$.

Let $U=W^{\top}$ (with $i$th column $w_{i}$). Note that since no two columns
are parallel, $\mathsf{krank}(U)\geq 2$. By applying the above fact on matrix
$M=U^{\odot m}$ with $i$th column $w_{i}^{\otimes m}$, we get that
$\mathsf{krank}(M)=m$, as required. Hence, Theorem 4.2 can be applied to
recover for all $i\in[m]$, all the unknown $a_{i}$, and up to ambiguities in
signs given by (unknown) $\xi_{i}\in\\{1,-1\\}$ the $b_{i}$ and $w_{i}$ as
well (we recover $\xi_{i}b_{i},\xi_{i}w_{i}$).

For the second half of the claim, let $\xi_{i}\in\\{1,-1\\}~{}\forall i\in[m]$
be any combination of signs. Consider the solution
$a^{\prime}_{i}=a_{i},w^{\prime}_{i}=\xi_{i}w_{i},b^{\prime}_{i}=\xi_{i}b_{i}$,
and let $g(x)$ represent the corresponding ReLU function given by these
parameters. Note that the Hermite polynomial $He_{t}(\xi z)=\xi^{t}He_{t}(z)$
for all $\xi\in\\{\pm 1\\}$ and $z\in\mathbb{R}$ . Hence, the Hermite
coefficients of order at least $2$ are equal for $f$ and $g$ i.e., for all
$t\geq 2$

$\hat{g}_{t}=\sum_{i=1}^{m}(-1)^{t}a_{i}~{}He_{t-2}(\xi_{i}b_{i})\cdot\frac{e^{-b_{i}^{2}/2}}{\sqrt{2\pi}}\cdot(\xi_{i}w_{i})^{\otimes
t}=\sum_{i=1}^{m}(-1)^{t}a_{i}~{}He_{t-2}(b_{i})\cdot\frac{e^{-b_{i}^{2}/2}}{\sqrt{2\pi}}\cdot
w_{i}^{\otimes t}=\hat{f}_{t}.$

Condition 5 also implies that zeroth and first Hermite coefficients of $f,g$
are also equal. All the Hermite coefficients are hence equal (and the
functions are squared-integrable w.r.t. the Gaussian measure for bounded
$a_{i}$). Thus the two functions $f$ and $g$ being identical follows since
Hermite polynomials form a complete orthogonal system.

∎

##### Non-identifiability of signs when $\\{w_{i}:i\in[m]\\}$ are not linearly
independent.

Theorem 4.1 shows that when the $\\{w_{i}:i\in[m]\\}$ are linearly
independent, the model is identifiable. The following claim shows that even in
the special setting when the biases $b_{i}=0~{}\forall i\in[m]$, whenever the
$(w_{i}:i\in[m])$ are linearly dependent, the model is non-identifiable (for
appropriate $(a_{i}:i\in[m])$) because of ambiguities in the signs (Theorem
4.2 also shows it is identifiable up to this sign ambiguity).

######  Claim 4.8.

Suppose $w_{1},\dots,w_{m}$ are linearly dependent. Then there exists
$a_{1},\dots,a_{m}$ (not all $0$) and signs
$\xi_{1},\xi_{2},\dots,\xi_{m}\in\\{\pm 1\\}$ with not all $+1$ such that the
ReLU networks $f$ and $g$ defined as

$f(x)\coloneqq\sum_{i=1}^{m}a_{i}\sigma(w_{i}^{\top}x),~{}~{}g(x)\coloneqq\sum_{i=1}^{m}a_{i}\sigma(\xi_{i}w_{i}^{\top}x)\text{
satisfy }f(x)=g(x)~{}\forall x\in\mathbb{R}^{d}.$

###### Proof.

Since $\\{w_{i}:i\in[m]\\}$ are linearly dependent, there exists
$(\beta_{i}:i\in[m])$ which are not all $0$ such that
$\sum_{i=1}^{n}\beta_{i}w_{i}=0$. Define $a_{i}=\beta_{i}$ for each $i\in[m]$,
and let $\xi_{i}=-1$ if $\beta_{i}\neq 0$ and $\xi_{i}=1$ otherwise. Let
$f(x)=\sum_{i=1}^{m}a_{i}\sigma(w_{i}^{\top}x)$ and
$g(x)=\sum_{i=1}^{m}a_{i}\sigma(\xi_{i}w_{i}^{\top}x)$.

From Lemma 3.5, it is easy to verify that all the even Hermite coefficients
are equal, and the odd Hermite coefficients for $\ell\geq 3$ are all $0$ since
$b_{i}=0$. Moreover the $\ell=1$ order Hermite coefficients are equal since

$\sum_{i=1}^{m}a_{i}w_{i}-\sum_{i=1}^{m}a_{i}\xi_{i}w_{i}=\sum_{i=1}^{m}a_{i}(1-\xi_{i})w_{i}=\sum_{i=1}^{m}2\beta_{i}w_{i}=0.$

All the Hermite coefficients of $f$ and $g$ are equal (and the functions are
also squared-integrable w.r.t. the Gaussian measure when the $a_{i}$ are
bounded). As the Hermite polynomials form a complete orthogonal basis, the two
ReLU network functions $f(x)$ and $g(x)$ are also equal. This concludes the
proof. ∎

##  5 Robustness Analysis

In this section, we prove Theorem 3.1 and Theorem 3.2 which give polynomial
time and sample complexity bounds for our algorithms. In the previous section
we showed that given oracle access to $\\{\hat{f}_{k}\\}$, we can recover the
exact network parameters $a,b,W$ (or at least up to signs). In reality, we can
only access polynomially many samples in polynomial time, and we will have
sampling errors when estimating $\\{\hat{f}_{k}\\}$. Therefore, given data
generated from the target network $(x_{1},y_{1}),...,(x_{N},y_{N})$, we will
approximate $\hat{f}_{k}$ through the empirical estimator

$T_{k}=\frac{1}{N}\sum_{i=1}^{N}y_{i}He_{k}(x_{i})$ (14)

Observe that $T_{k}$ is an unbiased estimator for $\hat{f}_{k}$. We first show
using standard concentration bounds that for any $\eta>0$, with
$N\geq\mathrm{poly}(d^{k},m,B,1/\eta)$ samples, the empirical estimates with
high probability satisfies
$\forall\ell\in[k],~{}\lVert\xi_{\ell}\rVert_{F}\coloneqq\|T_{\ell}-\hat{f}_{\ell}\|_{F}\leq\eta$
(see Appendix B.1). Hence for any constant $k$, with polynomial samples, we
can obtain with high probability, estimates for the tensors
$\\{T_{0},T_{1},\dots,T_{k}\\}$ that are accurate up to any desired inverse-
polynomial error.

The main algorithm in the robust setting is Algorithm 4 described below, which
approximately recovers the parameters for the activation units (up to signs)
that do not have large positive bias. The guarantees are given in the
following Theorem 5.1.

######  Theorem 5.1.

Suppose $\ell\in\mathbb{N}$ be a constant, and $\varepsilon>0$. If we are
given $N$ i.i.d. samples as described above from a ReLU network
$f(x)=a^{\top}\sigma(W^{\top}x+b)$ that is $B$-bounded Then there are
constants $c=c(\ell)>0,c^{\prime}>0$, signs $\xi_{i}\in\\{\pm 1\\}~{}\forall
i\in[m]$ and a permutation $\pi:[m]\to[m]$ such that Algorithm 4 given
$N\geq\mathrm{poly}_{\ell}(m,d,1/\varepsilon,1/s_{m}(W^{\odot\ell}),B)$ runs
in $\mathrm{poly}_{\ell}(N,m,d)$ time and with high probability outputs
$\\{\widetilde{a}_{i},\widetilde{b}_{i},\widetilde{w}_{i}:i\in[m^{\prime}]\\}$
such that for all $i\in[m]$ with $|b_{i}|<c\sqrt{\log(1/(\varepsilon\cdot
mdB))}$ we have that $\lVert
w_{i}-\xi_{\pi(i)}\widetilde{w}_{\pi(i)}\rVert_{2}+|a_{i}-\widetilde{a}_{\pi(i)}|+|b_{i}-\xi_{\pi(i)}\widetilde{b}_{\pi(i)}|\leq\varepsilon$.

In fact, the analysis just assumes that $\lVert T_{k}-\hat{f}_{k}\rVert_{F}$
are upper bounded up to an amount that is inverse polynomial in the different
parameters (this could also include other sources of error) to approximate the
2-layer ReLU network that approximates $f$ up to desired inverse polynomial
error $\varepsilon$.

Input: Estimates $T_{0},\dots,T_{2\ell+2}$ for
$\hat{f}_{0},\hat{f}_{1},\dots,\hat{f}_{2\ell+2}$;

Parameters: $\eta_{0},\eta_{1},\eta_{2},\eta_{3}>0$.;

1\. Let order-3 tensors
$T^{\prime}=\text{flatten}(T_{2\ell+1},\ell,\ell,1)\in\mathbb{R}^{d^{\ell}\times
d^{\ell}\times d}$ and let
$T^{\prime\prime}=\text{flatten}(T_{2\ell+2},\ell,\ell,2)\in\mathbb{R}^{d^{\ell}\times
d^{\ell}\times d^{2}}$.

2\. Set $k^{\prime}=\max_{r\leq
m}s_{r}(\text{flatten}(T_{2\ell+1},\ell,\ell+1,0))>\eta_{1}$. Run Jennrich’s
algorithm on $T^{\prime}$ to recover rank-1 terms
$\\{\alpha^{\prime}_{i}u_{i}^{\otimes\ell}\otimes u_{i}^{\otimes\ell}\otimes
u_{i}~{}|~{}i\in[k^{\prime}]\\}$, where $\forall
i\in[k^{\prime}],~{}u_{i}\in\mathbb{S}^{d-1}$ and
$\alpha^{\prime}_{i}\in\mathbb{R}$.

3\. Set $k^{\prime\prime}=\max_{r\leq
m}s_{r}(\text{flatten}(T_{2\ell+2},\ell,\ell+2,0))>\eta_{1}$. Run Jennrich’s
algorithm on $T^{\prime\prime}$ to recover rank-1 terms
$\\{\alpha^{\prime\prime}_{i}v_{i}^{\otimes\ell}\otimes
v_{i}^{\otimes\ell}\otimes v_{i}^{\otimes 2}~{}|~{}i\in[k^{\prime\prime}]\\}$,
where $\forall i\in[k^{\prime\prime}],~{}v_{i}\in\mathbb{S}^{d-1}$ and
$\alpha^{\prime\prime}_{i}\in\mathbb{R}$.

4\. Remove all the rank-$1$ terms in steps 2 and 3 with Frobenius norm
$<\eta_{2}$ i.e., $\alpha^{\prime\prime}_{i}$ or
$\alpha^{\prime}_{i}<\eta_{2}$. Also remove all duplicates from
$\\{u_{1},u_{2},\dots,u_{k^{\prime}}\\}\cup\\{v_{1},v_{2},\dots,v_{k^{\prime\prime}}\\}$
even up to signs i.e., remove iteratively from the above set vectors $v$ if
either of $+v,-v$ are within $\eta_{3}$ in $\ell_{2}$ distance of the other
vectors in the set, to get
$\widetilde{w}_{1},\widetilde{w}_{2},\dots,\widetilde{w}_{m^{\prime}}$.

5\. Run the subroutine
RecoverScalars$(\ell,\\{\widetilde{w}_{i}:i\in[m^{\prime}]\\},T_{\ell},T_{\ell+1},T_{\ell+2},T_{\ell+3})$
(i.e., Alg. 2) to get
$\\{\widetilde{a}_{i},\widetilde{b}_{i}:i\in[m^{\prime}]\\}$.

Result: Output $\\{\widetilde{w}_{i},\widetilde{a}_{i},\widetilde{b}_{i}:1\leq
i\leq m^{\prime}\\}$.

Algorithm 4 for order $\ell$: recover $\\{a_{i}\\}$, and (up to signs)
$\\{b_{i},w_{i}\\}$ given estimates $\\{T_{0},\dots,T_{2\ell+2}\\}$.

The following algorithm (Algorithm 5) shows how to find a depth-2 ReLU network
that fits the data i.e., achieves arbitrarily small $L_{2}$ error. The
algorithm uses Algorithm 4 as a black-box to first approximately recover the
unknown parameters of the activation units (with not very large bias) up to
signs, and then setup an appropriate linear regression problem to find a
network that fits the data.

Input: $N$ i.i.d. samples of the form $(x_{i},y_{i})$;

Parameters: $\varepsilon,\eta_{0},\eta_{1},\eta_{2},\eta_{3}>0$.;

1\. Construct estimates $T_{0},\dots,T_{2\ell+2}$ for
$\hat{f}_{0},\hat{f}_{1},\dots,\hat{f}_{2\ell+2}$ using the first
$\frac{N}{2}$ samples.

2\. Let $S=\\{(\tilde{w}_{i},\tilde{a}_{i},\tilde{b}_{i})\\}$ be the output of
Algorithm 4 on inputs $T_{0},\dots,T_{2\ell+2}$ when run with parameters
$\eta_{0},\eta_{1},\eta_{2},\eta_{3}>0$.

3\. For each $(x_{i},y_{i})$ and $i\in[N/2+1,N]$ construct the feature mapping
$\phi(x_{i})=(Z(x_{i}),Z^{\prime}(x_{i}))$ as described in the proof of Lemma
5.8.

4\. Set $\tau=20m(8|S|+d)B\sqrt{\log(\frac{mdB|S|}{\varepsilon})}$ and find,
via projected gradient descent, a vector $\hat{\beta}$ such that
$\hat{L}_{\tau}(\hat{\beta})\leq\min_{\beta:\|\beta\|\leq\sqrt{8|S|}+m(1+B)}\hat{L}_{\tau}(\beta)+\frac{\varepsilon^{2}}{100}.$
Here $\hat{L}_{\tau}(\beta)$ is defined as
$\hat{L}_{\tau}(\beta)=\frac{2}{N}\sum_{i=\frac{N}{2}+1}^{N}(y_{i}-\beta^{\top}\phi(x_{i}))^{2}\mathds{1}\big{(}\|\phi(x_{i})\|<\tau\big{)}.$

Result: The function $g(x)=\hat{\beta}^{\top}\phi(x)$.

Algorithm 5 Outputs a function $g(x)$ that approximates the target network
$f(x)$ in mean squared error.

The error parameters $\eta_{0},\eta_{1},\eta_{2},\eta_{3}$ can be set with
appropriate polynomial dependencies on
$\varepsilon,d^{\ell},m,B,s_{m}(W^{\odot\ell})$ to obtain the recovery
guarantees in Theorem 5.1 and Theorem 3.2. See Section 5.4 for details.

#####  Overview of Analysis.

The error in the tensors $T_{k}$ introduces additional challenges that we
described in Section 1. The analysis is technical and long, but we now briefly
describe the main components.

_(i)_ Recall from Section 1, that when there are errors, it may not even be
possible to recover the parameters of some ReLU units! In particular when the
bias $b_{i}$ is large in magnitude, the ReLU unit will be indistinguishable
from a simple linear function. It will contribute negligibly to any of the
higher order Hermite coefficients, and hence will be impossible to recover
them individually (especially if there are multiple such units). For a desired
recovery error $\varepsilon>0$, the $m$ hidden ReLU units are split into
groups (for analysis)

$G=\\{i\in[m]~{}|~{}|b_{i}|<O\big{(}\sqrt{\log(1/(\varepsilon\cdot
mdB))}\big{)}\\},\text{ and }P=\\{1,2,\dots,m\\}\setminus G.$

We aim to recover all of the parameters of the units corresponding to $G$ up
to signs. For the terms in $P$, we will show the existence of a linear
function that approximates the total contribution from all the terms in $P$.

_(ii)_ The tensor decomposition steps (steps 2-3) are simpler in the no-noise
setting: the parameter $w_{i}$ of the $i$th ReLU unit can be recovered (up to
sign ambiguity) as long its bias $b_{i}$ is not a root of $He_{2\ell-1}$. When
there is noise, there could be terms $i\in[m]$ for which $b_{i}$ are not roots
of $He_{k-2}(x)$, and yet their signal can get swamped by the sampling error
in the tensor. We can only hope to recover those $k\leq m$ components whose
corresponding coefficient is above some chosen threshold $\eta_{1}$ (the other
terms are considered as part of the error tensor). However a technical issue
that arises is that the robust recovery guarantees for tensor decomposition
algorithms lose polynomial factors in different parameters including the least
singular value ($s_{k}(\cdot)$) of the factor matrices. Hence, for each of
step 2 and 3, we argue that recovery is possible only if the coefficient of
the corresponding term is significantly large, and this may give reasonable
estimates for only a subset of these $m$ terms with coefficients
$>\varepsilon$.

_(iii)_ When decomposing two consecutive tensors $T_{2\ell+1}$ and
$T_{2\ell+2}$, we use Lemma 4.3 to argue that each $i\in G$ will have a large
coefficient in at least one of these two tensors. Hence we can stitch together
estimates $\\{\widetilde{w}_{i}:i\in G\\}$ which are accurate up to a sign and
small error. This will in turn be used to recover $a_{i},b_{i}$ for $i\in G$,
with properties of Hermite polynomials used to ensure that the errors do not
propagate badly.

_(iv)_ We argue that the other ReLU units in $P=[m]\setminus G$ can be
approximated altogether using a linear function. This is obtained by
subtracting from estimates $T_{0},T_{1}$ with the corresponding terms from
$G$.

_(v)_ The above arguments let us compute good approximations to the parameters
for the units in $G$, but only up to signs. In order to use this to learn a
good predictor for $f(x)$ we consider solving a truncated linear regression
problem in an expanded feature space. At a high level, for each $i\in G$,
given estimates $(\widetilde{a}_{i},\widetilde{w}_{i},\widetilde{b}_{i})$ we
consider an expanded feature representation for this unit into an
$8$-dimensional vector where each coordinate is of the form
$\xi_{i_{3}}\widetilde{a}_{i}\sigma(\xi_{i_{1}}\widetilde{w}_{i}\cdot
x+\xi_{i_{2}}\widetilde{b}_{i})$ for
$\xi_{i_{1}},\xi_{i_{2}},\xi_{i_{3}}\in\\{-1,+1\\}$.555While this portion of
the algorithm works more generally with ambiguities in the sign of $a_{i}$,
$b_{i}$, and $w_{i}$, in our case, the sign ambiguity of $w_{i}$ and $b_{i}$
are coordinated, and the sign of the $a_{i}$ are also recovered correctly;
hence a $2$-dimensional vector suffices in this case. Moreover the terms can
be consolidated to get an equivalent ReLU network with at most $|G|+2$ hidden
units (see Claim B.5 and Lemma 5.8). Repeating this for every $i\in G$ it is
easy to see that there is a linear function in the expanded space that
approximates the part of the function $f(x)$ that depends on units in $G$.
Combining with the previous argument that the units in $P=[m]\setminus G$ can
be approximated by a linear function in the original feature space, we deduce
that there is an $O(d+m)$ dimensional feature space where $f(x)$ admits a good
linear approximation. We then solve a truncated least squares problem in this
space to obtain our final function $g(x)$ that approximates $f(x)$ in $L_{2}$
error.

In the following sections, we will state the main claims and intermediate
steps that prove the robustness of the algorithm and establish Theorems 5.1,
3.1, 3.2.

###  5.1 Estimating the Hermite Coefficients

First, we derive concentration bounds on $\xi_{k}$, which will be followed by
error bounds of the recovered parameters
$\widetilde{a},\widetilde{b},\widetilde{W}$ in terms of $\xi_{k}$.

######  Lemma 5.2.

For any $\eta>0$, if $T_{k}$ is estimated from $N\geq
c_{k}d^{k}m^{2}B^{4}\cdot\mathrm{polylog}(mdB/\eta)/(\eta^{2})$ samples, then
for some constant $c_{k}>0$ that depends only on $k$, we have with probability
at least $1-(mdB)^{-\log(md)}$,

$\|T_{k}-\hat{f}_{k}\|_{F}\leq\eta.$ (15)

The above statement follows by applying Rosenthal inequality [Ros70, P+15]
along with Markov’s inequality. Please refer to Section B for the proof of the
theorem.

###  5.2 Recovering the Parameters under Errors

Suppose $\varepsilon>0$ is the desired recovery error. The $m$ hidden units
are split into groups

$G=\Big{\\{}i\in[m]~{}|~{}|b_{i}|<c_{\ell}\sqrt{\log\Big{(}\frac{1}{\varepsilon
mdBs_{m}(W^{\odot\ell})}\Big{)}}\Big{\\}},\text{ and
}P=\\{1,2,\dots,m\\}\setminus G.$ (16)

where $c_{\ell}$ is an appropriate constant that depends only on the constant
$\ell>0$. Note that under the assumption that $s_{m}(W^{\odot\ell})\geq
1/\mathrm{poly}(m,d,B)$ in Theorem 5.1, this reduces to

$G=\Big{\\{}i\in[m]~{}|~{}|b_{i}|<c^{\prime}_{\ell}\sqrt{\log(1/\varepsilon
mdB)}\Big{\\}}.$

We aim to recover all of the parameters of the units corresponding to $G$. For
the terms in $P$, we will learn a linear function that approximates the total
contribution from all the terms in $P$.

##### Recovery of Weight Vectors $w_{i}$ for the Terms in $G$.

We first state the following important lemma showing that Jennrich’s algorithm
run with an appropriate choice of rank $k$ will recover each large term up to
a sign ambiguity.

######  Lemma 5.3.

Suppose $\varepsilon_{2}\in(0,\tfrac{1}{4})$, and
$\ell_{1},\ell_{2}\geq\ell,\ell_{3}>0$ be constants for some fixed $\ell$, and
$T=\mathsf{flatten}(\hat{f}_{\ell_{1}+\ell_{2}+\ell_{3}},\ell_{1},\ell_{2},\ell_{3})$
have decomposition $T=\sum_{i=1}^{m}\lambda_{i}(u_{i}\otimes v_{i}\otimes
z_{i})$ with $\lambda_{i}\in\mathbb{R}$ and unit vectors
$u_{i}=w_{i}^{\otimes\ell_{1}}\in\mathbb{R}^{d^{\ell_{1}}},v_{i}=w_{i}^{\otimes\ell_{2}}\in\mathbb{R}^{d^{\ell_{2}}},z_{i}=w_{i}^{\otimes\ell_{3}}\in\mathbb{R}^{d^{\ell_{3}}}$.
There exists
$\eta_{1}=\mathrm{poly}(\varepsilon_{2},s_{m}(W^{\odot\ell}))/\mathrm{poly}(m,d,B)>0$
and
$\varepsilon_{1}\coloneqq\max\\{2\varepsilon_{2},1/\mathrm{poly}(1/\varepsilon_{2},1/s_{m}(W^{\odot\ell}),d^{\ell_{1}+\ell_{2}+\ell_{3}},B)\\}$
such that if

$\displaystyle\|T-\widetilde{T}\|_{F}$
$\displaystyle\leq\eta^{\prime}_{1}\coloneqq\min\Big{\\{}\mathrm{poly}(\varepsilon_{2},s_{m}(W^{\odot\ell}))/\mathrm{poly}(m,d,B,1/\eta_{1}),\frac{\eta_{1}}{2}\Big{\\}},$
(17)

then Jennrich’s algorithm runs with rank
$k^{\prime}\coloneqq\operatorname*{argmax}_{r\leq
m}s_{r}(\text{flatten}(\widetilde{T},\ell_{1},\ell_{2}+\ell_{3},0))>\eta_{1}$
and w.h.p. outputs666 Note that one can also choose to pad the output with
zeros to output $m$ sets of parameters instead of $k^{\prime}$ if required.
$\\{\widetilde{\lambda}_{i},\widetilde{w}_{i}\\}_{i\in[k^{\prime}]}$ such that
there exists a permutation $\pi:[m]\to[m]$ and signs
$\xi_{i}\in\\{1,-1\\}~{}\forall i\in[m]$ satisfying:

$\displaystyle(i)~{}~{}\forall i\in[m],$
$\displaystyle~{}~{}|\lambda_{i}-\widetilde{\lambda}_{\pi(i)}|\leq\varepsilon_{2}^{2},\text{
and }$ (18) $\displaystyle(ii)~{}~{}\forall i\in[m],$ $\displaystyle\text{
s.t. }|\lambda_{i}|>\varepsilon_{1},~{}\text{we have}~{}\lVert w_{i}^{\otimes
t}-\xi_{\pi(i)}^{t}\widetilde{w}_{\pi(i)}^{\otimes
t}\rVert_{2}\leq\varepsilon_{2},~{}~{}~{}\forall t\in[2\ell].$ (19)

A direct application of Lemma 5.3 establishes the following claim, showing
that we can recover all the weight vectors $w_{i}$ for each term $i\in G$ up
to a sign ambiguity.

######  Lemma 5.4.

For any $\varepsilon_{2}>0$, there exists an
$\eta_{2}^{\prime}=\frac{\mathrm{poly}(\varepsilon_{2},s_{m}(W^{\odot\ell}))}{\mathrm{poly}_{\ell}(m,d,B)}>0$
such that if the estimates $\lVert
T_{k}-\hat{f}_{k}\rVert_{F}\leq\eta_{2}^{\prime}$ for all
$k\in\\{0,1,\dots,2\ell+2\\}$, then steps 1-4 of Algorithm 4 finds a set
$\\{\widetilde{w}_{i}:i\in[m^{\prime}]\\}$ such that there exists a one-to-one
map $\pi:[m^{\prime}]\to[m]$ satisfying (i) every $i\in G$ has a pre-image in
$\pi$ (i.e., every term in $G$ is recovered), and for appropriate signs
$\\{\xi_{i}\in\\{1,-1\\}:i\in[m^{\prime}]\\}$,

$\forall i\in[m^{\prime}],\forall
t\in[2\ell],~{}~{}\lVert\xi_{i}^{t}\widetilde{w}_{i}^{\otimes
t}-w_{\pi(i)}^{\otimes t}\rVert_{F}\leq\varepsilon_{2}.$ (20)

In particular $\forall i\in[m^{\prime}]$, we have
$\lVert\xi_{i}\widetilde{w}_{i}-w_{\pi(i)}\rVert_{2}\leq\varepsilon_{2}$.

##### Recovering Error for the Parameters $a_{i},b_{i}$ with $i\in G$.

The following lemmas now proves the recovery for each $i\in G$, the $a_{i}$
(no sign ambiguity) and the $b_{i}$ up to the same sign ambiguity as in
$w_{i}$ (and in fact, this holds for all the terms output in steps 1-5 of
Algorithm 4).

######  Lemma 5.5.

For $\varepsilon>0$ in the definition of $G$ in (16), there exists
$\eta_{3}^{\prime}=\frac{\mathrm{poly}(\varepsilon,s_{m}(W^{\otimes\ell}))}{\mathrm{poly}_{\ell}(m,d,B)}>0$,
and
$\varepsilon_{3}^{\prime}=\frac{\mathrm{poly}(\varepsilon,s_{m}(W^{\otimes\ell}))}{\mathrm{poly}_{\ell}(m,d,B)}>0$
such that for some $\xi_{i}\in\\{\pm 1\\}~{}\forall i\in[m^{\prime}]$

$\text{if}~{}~{}\lVert
T_{k}-\hat{f}_{k}\rVert_{F}\leq\eta_{3}^{\prime}~{}~{}\forall k\leq
2\ell+2,~{}~{}~{}~{}\text{ and }~{}~{}~{}\lVert\widetilde{w}_{i}^{\otimes
t}-\xi^{t}w_{i}^{\otimes t}\rVert_{F}\leq\varepsilon_{3}^{\prime},\forall
i\in[m^{\prime}],\forall\ell\leq t\leq\ell+3.$

then steps 5-6 of the algorithm finds
$(\widetilde{a}_{i},\widetilde{b}_{i}:i\in[m^{\prime}])$ such that

$|\widetilde{a}_{i}-a_{i}|\leq\varepsilon,\text{ and
}|\widetilde{b}_{i}-\xi_{i}b_{i}|\leq\varepsilon.$ (21)

Note that in Lemma 5.4 we showed that $G$ is contained in the $m^{\prime}$
terms output in steps 1-4 (and hence step 5 as well). The above uses the
following two lemmas which gives a robust version of Lemma 4.5 when there are
errors in the estimates. We remark that $\beta=\alpha e^{-z^{2}/2}$ in the
notation of Lemma 4.5.

######  Lemma 5.6 (Robust version of Lemma 4.5 for $k\geq 2$).

Suppose $k\in\mathbb{N},k\geq 2,B\geq 1$, and $\alpha,z\in\mathbb{R}$ be
unknown parameters. There exists a constant $c_{k}=c(k)\geq 1$ such that for
any $\varepsilon\in(0,\tfrac{1}{4})$ satisfying (i)
$|\alpha|\in[\tfrac{1}{B},B]$ and $|z|\leq B$, and (ii)
$|z|<2\sqrt{\log(c_{k}/\varepsilon^{1/4}(1+B)^{3}))}$ , if we are given values
$\gamma_{k},\gamma_{k+1},\gamma_{k+2},\gamma_{k+3}$ s.t. for some
$\xi\in\\{\pm 1\\}$,

$\Big{|}\gamma_{j}-\alpha\cdot\frac{(-1)^{j}}{\sqrt{2\pi}}e^{-z^{2}/2}He_{j}(\xi
z)\Big{|}\leq\varepsilon^{\prime}=\frac{\varepsilon^{4}}{2(1+B)^{2}}~{}~{}~{}\forall
j\in\\{k,k+1,k+2,k+3\\},$

then the estimates $\widetilde{z},\widetilde{\alpha}$ obtained as:

$\displaystyle\widetilde{z}=-\frac{\gamma_{q+1}+q\cdot\gamma_{q-1}}{\gamma_{q}}~{}~{}\text{where
}q:=\operatorname*{argmax}_{\begin{subarray}{c}j\in\\{k+1,k+2\\}\end{subarray}}|\gamma_{j}|,~{}~{}\text{and}~{}~{}\widetilde{\alpha}$
$\displaystyle=(-1)^{q}\frac{\sqrt{2\pi}\gamma_{q}}{e^{-\frac{\widetilde{z}^{2}}{2}}He_{q}(\widetilde{z})}$
$\displaystyle\text{satisfy}~{}~{}~{}\big{|}\widetilde{z}-\xi
z\big{|}\leq\frac{\varepsilon|z|}{B+1}\leq\varepsilon,~{}~{}\text{ and
}~{}~{}\big{|}\widetilde{\alpha}-\alpha\big{|}$
$\displaystyle\leq\varepsilon.$ (22)

The simpler variant of the above lemma (Lemma 5.6) for $k=1$ which is used in
the full-rank setting, follows a very similar analysis and is stated below.

######  Lemma 5.7 (Robust version of Lemma 4.5 for $k=1$).

Suppose $B\geq 1$, and $\alpha,z\in\mathbb{R}$ be unknowns. There exists a
constant $c\geq 1$ such that for any $\varepsilon\in(0,\tfrac{1}{4})$
satisfying (i) $|\alpha|\in[\tfrac{1}{B},B]$ and $|z|\leq B$, and (ii)
$|z|<2\sqrt{\log(c/\varepsilon^{1/4}(1+B)^{3}))}$, if we are given values
$\gamma_{0},\gamma_{1}$ s.t. for some $\xi\in\\{\pm 1\\}$,

$\Big{|}\gamma_{j}-\alpha\cdot\frac{(-1)^{j}}{\sqrt{2\pi}}e^{-z^{2}/2}He_{j}(\xi
z)\Big{|}\leq\varepsilon^{\prime}=\frac{\varepsilon^{4}}{2(1+B)^{2}}~{}~{}~{}\forall
j\in\\{0,1\\},$

then the estimates $\widetilde{z},\widetilde{\alpha}$ obtained as:

$\displaystyle\widetilde{z}$
$\displaystyle=-\frac{\gamma_{1}}{\gamma_{0}},~{}~{}\text{ and
}~{}~{}\widetilde{\alpha}=\frac{\sqrt{2\pi}\gamma_{0}}{e^{-\frac{\widetilde{z}^{2}}{2}}}$
$\displaystyle\text{satisfy}~{}~{}~{}\big{|}\widetilde{z}-\xi
z\big{|}\leq\frac{\varepsilon|z|}{B+1}$
$\displaystyle\leq\varepsilon,~{}~{}\text{ and
}~{}~{}\big{|}\widetilde{\alpha}-\alpha\big{|}\leq\varepsilon.$ (23)

###  5.3 Learning Guarantees via Linear Regression

In the previous sections we designed algorithms based on tensor decompositions
that, given i.i.d. samples from a network
$f(x)=\sum_{i=1}^{m}a_{i}\sigma(w^{\top}_{i}x+b_{i})$, can recover
approximations (up to signs) for “good units”, i.e,
$G=\\{i\in[m]:|b_{i}|<O\big{(}\sqrt{\log(\frac{1}{\varepsilon
mdB})}\big{)}\\}$. In this section we will show how to use these
approximations to perform improper learning of the target network $f(x)$ via a
simple linear regression subroutine. Our algorithm will output a functions of
the form
$g(x)=\sum_{i=1}^{m^{\prime}}a^{\prime}_{i}\sigma(w^{\prime\top}_{i}x+b^{\prime}_{i})+w^{\prime\prime\top}x+C$,
where $m^{\prime}\leq 8m$. In particular we will prove the following.

######  Lemma 5.8.

Let $\varepsilon>0$ and $f(x)=\sum_{i=1}^{m}a_{i}\sigma(w^{\top}_{i}x+b_{i})$
be an unknown target network. Let $S$ be a given set of tuples of the form
$(\widetilde{w}_{i},\widetilde{b}_{i},\widetilde{a}_{i})$ with
$\|\tilde{w}_{i}\|=1$, such that for each $i\in G$, there exists $j\in S$, and
$\xi_{j_{1}},\xi_{j_{2}},\xi_{j_{3}}\in\\{-1,+1\\}$, such that
$\|w_{i}-\xi_{j_{1}}\widetilde{w}_{j}\|\leq O(\frac{\varepsilon}{mdB})$,
$|b_{i}-\xi_{j_{2}}\widetilde{b}_{j}|\leq O(\frac{\varepsilon}{mdB})$, and
$|a_{i}-\xi_{j_{3}}\widetilde{a}_{j}|\leq O(\frac{\varepsilon}{mdB})$. Then
for any $\delta\in(0,1)$, given
$N=\mathrm{poly}(m,d,B,\frac{1}{\varepsilon},\log(\frac{1}{\delta}))$ i.i.d.
samples of the form $(x,y=f(x))$ where $x\sim N(0,I)$, there exists an
algorithm (Algorithm 5) that runs time polynomial in $N$ and with probability
at least $1-\delta$ outputs a network $g(x)$ of the form
$g(x)=\sum_{i=1}^{m^{\prime}}a^{\prime}_{i}\sigma(w^{\prime\top}_{i}x+b^{\prime}_{i})+w^{\prime\prime\top}x+C$,
where $m^{\prime}\leq 8|S|$, such that

$\mathbb{E}_{x\sim\mathcal{N}(0,I_{d\times
d})}\big{(}f(x)-g(x)\big{)}^{2}\leq\varepsilon^{2}.$

Furthermore, when $\xi_{j_{1}}=\xi_{j_{2}}$ and $\xi_{j_{3}}=+1$ for all $j\in
S$ (i.e., the sign ambiguity of $w_{i}$ and $b_{i}$ are the same, and there is
no ambiguity in the sign of $a_{i}$ for all $i\in G$), then the number of
hidden units in $g(x)$ is at most $|S|+2$.

While the above lemma is more general, when it is applied in the context of
Theorem 3.2 it satisfies the conditions of the “furthermore” portion of the
lemma. Our algorithm for recovering $g(x)$ will set up a linear regression
instance in an appropriate feature space. In order to do this we will need the
lemma stated below that shows that there is a good linear approximation for
the units not in $G$, i.e., $P=[m]\setminus G$.

######  Lemma 5.9 (Approximating $f_{P}$).

Let $c>2$ be a fixed constant. Consider
$f_{P}(x)=\sum_{i=1}^{m}a_{i}\sigma(w^{\top}_{i}x+b_{i})\mathds{1}\big{(}|b_{i}|\geq
c\sqrt{\log(\frac{1}{\varepsilon mdB})}\big{)}$. Then there exists a function
$g_{P}(x)=\beta^{\top}_{P}x+C_{P}$ where $\|\beta_{P}\|\leq mB$ and
$|C_{P}|\leq mB^{2}$ such that for a constant $c^{\prime}>0$ that depends on
$c$,

$\mathbb{E}_{x\sim\mathcal{N}(0,I)}[f_{P}(x)-g_{P}(x)]^{2}=c^{\prime}\varepsilon^{2}.$

We first establish the main result assuming the lemma above and provide a
proof of the lemma at the end of the subsection.

###### Proof of Lemma 5.8.

In order to find the approximate network $g(x)$ we will set up a linear
regression problem in an appropriate feature space. We begin by describing the
construction of the feature space and showing that there does indeed exist a
linear function in the space that approximates $f(x)$. We first focus on the
terms in the set $G$, i.e.,

$f_{G}(x)=\sum_{i\in G}a_{i}\sigma(w^{\top}_{i}x+b_{i}).$

In order to approximate $f_{G}(x)$ we create for each
$(\widetilde{w}_{j},\widetilde{b}_{j},\widetilde{a}_{j})\in S$, eight features
$Z_{j,1},\dots,Z_{j,8}$ where each feature is of the form
$\xi_{j_{3}}\widetilde{a}_{j}\sigma(\xi_{j_{1}}\widetilde{w}^{\top}_{j}x+\xi_{j_{2}}\widetilde{b}_{j})$
for $\xi_{j_{1}},\xi_{j_{2}},\xi_{j_{3}}\in\\{-1,+1\\}$. Consider a particular
$i\in G$. Since the set $S$ consists of a good approximation
$(\widetilde{w}_{j},\widetilde{b}_{j},\widetilde{a}_{j})$ for the unit $i$, it
is easy to see that one of the eight features corresponding to $Z_{j,:}$
approximates the $i$th unit well (by matching the signs appropriately). In
other words we have that there exists $r\in[8]$ such that

$\displaystyle|Z_{j,r}-a_{i}\sigma(w^{\top}_{i}x+b_{i})|$
$\displaystyle=|{a^{\prime}_{j}}\sigma({w^{\prime}_{j}}^{\top}x+{b^{\prime}}_{j})-a_{i}\sigma(w^{\top}_{i}\cdot
x+b_{i})|$ (24)
$\displaystyle\leq|(a^{\prime}_{j}-a_{i})\sigma(w^{\top}_{i}x+b_{i})|+|a^{\prime}_{j}\big{(}\sigma({w^{\prime}}^{\top}_{j}x+b^{\prime}_{j})-\sigma(w^{\top}_{i}x+b_{i})\big{)}|$
(25) $\displaystyle\leq
O\Big{(}\frac{\varepsilon}{mdB}\Big{)}(B+\|x\|)+O\Big{(}\frac{\varepsilon}{mdB}\Big{)}+O\Big{(}\frac{\varepsilon}{mdB}\Big{)}\|x\|.$
(26)

Noting that $\mathbb{E}[\|x\|^{2}]=d$ we get that there exists a vector
$\beta^{*}_{1}$ with $\|\beta^{*}_{1}\|_{2}\leq\sqrt{8|S|}$ in the feature
space $Z(x)$ defined as above such that

$\displaystyle\mathbb{E}\big{[}f_{1}(x)-{\beta^{*}_{1}}^{\top}Z(x)\big{]}^{2}\leq\frac{\varepsilon^{2}}{100}.$
(27)

To approximate terms not in $G$, i.e., $f_{P}(x)=\sum_{i\notin
G}a_{i}\sigma(w^{\top}_{i}x+b_{i})$, we use Lemma 5.9 to get that there exists
a vector $\beta^{*}_{2}$ with $\|\beta^{*}_{2}\|\leq m(1+B)$ in the
$Z^{\prime}=(x,1)$ feature space such that

$\displaystyle\mathbb{E}\big{[}f_{P}(x)-{\beta^{*}_{2}}^{\top}Z^{\prime}(x)\big{]}^{2}\leq\frac{\varepsilon^{2}}{100}.$
(28)

Combining the above and noting that $y=f(x)=f_{G}(x)+f_{P}(x)$, we get that
there exists a vector $\beta^{*}$ in the $\phi(x)=(Z(x),Z^{\prime}(x))$
feature space with $\|\beta^{*}\|\leq\sqrt{8|S|}+m(1+B)$such that

$\displaystyle\mathbb{E}\big{[}y-{\beta^{*}}^{\top}\phi(x)\big{]}^{2}\leq\frac{\varepsilon^{2}}{20}.$
(29)

In order to approximate $\beta^{*}$ we solve a truncated least squares
problem. In particular, define the truncated squared loss
$L_{\tau}(\beta)=\mathbb{E}[(y-\beta^{\top}\phi(x))^{2}\mathds{1}({\|\phi(x)\|<\tau})]$.
Furthermore we define the empirical counter part $\hat{L}_{\tau}(\beta)$ based
on $N$ i.i.d. samples drawn from the distribution of $\phi(x)$. For an
appropriate value of $\tau$ we will output $\hat{\beta}$ such that

$\displaystyle\hat{L}_{\tau}(\hat{\beta})$
$\displaystyle\leq\min_{\beta:\|\beta\|\leq\sqrt{8|S|}+m(1+B)}\hat{L}_{\tau}(\beta)+\frac{\varepsilon^{2}}{100}.$
(30)

In particular we will set
$\tau=20m(8|S|+d)B\sqrt{\log(\frac{mdB|S|}{\varepsilon})}$. Notice that the
empirical truncated loss above is convex and for the chosen value of $\tau$,
has gradients bounded in norm by
$\mathrm{poly}(m,d,B,|S|,\frac{1}{\varepsilon})$. Hence we can use the
projected gradient descent algorithm [BBV04] to obtain a $\hat{\beta}$ that
achieves the above guarantee in
$N\cdot\mathrm{poly}(m,d,B,|S|,\frac{1}{\varepsilon})$ time. Furthermore using
standard uniform convergence bounds for bounded loss functions [MRT18] we get
that if
$N=\mathrm{poly}(m,d,B,|S|,\frac{1}{\varepsilon},\log(\frac{1}{\delta}))$ then
with probability at least $1-\delta$ we have

$\displaystyle L_{\tau}(\hat{\beta})$
$\displaystyle\leq\min_{\beta:\|\beta\|\leq\sqrt{8|S|}+m(1+B)}{L}_{\tau}(\beta)+\frac{\varepsilon^{2}}{50}$
(31) $\displaystyle\leq{L}_{\tau}(\beta)+\frac{\varepsilon^{2}}{50}.$ (32)

Finally, it remains to relate the truncated loss $L_{\tau}(\beta)$ to the true
loss $L(\beta)=\mathbb{E}[y-\beta^{\top}\phi(x)]^{2}$. We have that for any
$\beta$ such that $\|\beta\|_{2}\leq\sqrt{8|S|}+m(1+B)$,

$\displaystyle|L_{\tau}(\beta)-L(\beta)|$
$\displaystyle=\mathbb{E}[(y-\beta^{\top}\phi(x))^{2}\mathds{1}(\|\phi(x)\|\geq\tau)].$
(33)

Next notice that if $\|\phi(x)\|\geq 2^{j}\tau$ then we must have that either
$|\tilde{a}_{j}\sigma(\tilde{w}^{\top}_{j}\cdot
x+\tilde{b}_{j})|\geq\frac{2^{j}\tau}{8|S|+d}$ or that for some $i\in[d]$,
$|x_{i}|\geq\frac{2^{j}\tau}{8|S|+d}$. For our choice of $\tau$, this
probability is bounded by
$(8|S|+d)e^{-2^{2j}\Omega(\log(\frac{mdB|S|}{\varepsilon}))}$. Hence we get
that

$\displaystyle|L_{\tau}(\beta)-L(\beta)|$
$\displaystyle=\mathbb{E}[(y-\beta^{\top}\phi(x))^{2}\mathds{1}(\|\phi(x)\geq\tau\|)]$
(34)
$\displaystyle=\sum_{j=0}^{\infty}\mathbb{E}[(y-\beta^{\top}\phi(x))^{2}\mathds{1}(\|\phi(x)\|\in[2^{j}\tau,2^{j+1}\tau))]$
(35)
$\displaystyle\leq\sum_{j=0}^{\infty}O(2^{2j}m^{2}(8|S|+d)^{2}\tau^{2})(8|S|+d)e^{-2^{2j}\Omega(\log(\frac{mdB|S|}{\varepsilon}))}$
(36) $\displaystyle\leq\frac{\varepsilon^{2}}{50}.$ (37)

Hence, the output network
$g(x)=\hat{\beta}\cdot\phi(x)=\sum_{i=1}^{m^{\prime}}a^{\prime}_{i}\sigma(w^{\prime\top}_{i}x+b^{\prime}_{i})+w^{\prime\prime\top}x+C$
satisfies with probability at least $1-\delta$ that

$\mathbb{E}_{x\sim\mathcal{N}(0,I_{d\times
d})}\big{(}f(x)-g(x)\big{)}^{2}\leq\varepsilon^{2}.$

Notice that since a linear function can be simulated via two ReLU units (see
Claim B.5), our output function $g(x)$ is indeed a depth-2 neural network with
$m^{\prime}+2\leq 8m$ hidden units.

Furthermore, while the statement of Lemma 5.8 assumes that the signs of units
in $G$ are completely unknown, the output of the tensor decomposition
procedure from Theorem 3.2 in fact recovers, for each $i\in G$, the signs of
$a_{i}$ exactly and the signs of the corresponding $(w_{i},b_{i})$ are either
both correct or both incorrect. Hence when applying Lemma 5.8 to our
application we only need to create two features for each unit in $G$. In other
words we can output a network of the form
$g(x)={w^{{}^{\prime\prime}}}^{\top}x+C+\sum_{i=1}^{m^{\prime}}a^{\prime}_{i}\sigma(w^{\prime\top}_{i}x+b^{\prime}_{i})+a^{\prime\prime}_{i}\sigma(-w^{\prime\top}_{i}x-b^{\prime}_{i})$,
where $m^{\prime}\leq m$. Finally, from Claim B.5, the above network can be
written as a depth-2 network with ReLU activations and at most $m+2$ hidden
units. ∎

We end the subsection with the proof of Lemma 5.9.

###### Proof of Lemma 5.9.

Consider a particular unit $i$ such that
$b_{i}>c\sqrt{\log(\frac{1}{\varepsilon mdB})}$. Then notice that
$z_{i}=w^{\top}_{i}x+b_{i}\sim\mathcal{N}(b_{i},1)$. By using standard
properties of the Gaussian cdf, we get that by approximating $\sigma(z_{i})$
by the linear term $z_{i}$ we incur the error

$\displaystyle\mathbb{E}_{z_{i}\sim\mathcal{N}(b_{i},1)}(z_{i}-\sigma(z_{i}))^{2}$
$\displaystyle=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{0}z^{2}_{i}e^{-(z_{i}-b_{i})^{2}/2}dz_{i}$
(38) $\displaystyle\leq O(b^{2}_{i})e^{-b^{2}_{i}/2}\leq
O\Big{(}\frac{\varepsilon^{2}}{m^{2}}\Big{)}.$ (39)

Similarly for a unit with $b_{i}<-c\sqrt{\log(\frac{1}{\varepsilon mdB})}$, by
approximating $\sigma(z_{i})$ with the constant zero function we incur the
error

$\displaystyle\mathbb{E}_{z_{i}\sim\mathcal{N}(b_{i},1)}(\sigma(z_{i}))^{2}$
$\displaystyle=\frac{1}{\sqrt{2\pi}}\int_{0}^{\infty}z^{2}_{i}e^{-(z_{i}-b_{i})^{2}/2}dz_{i}$
(40) $\displaystyle\leq O(b^{2}_{i})e^{-b^{2}_{i}/2}\leq
O\Big{(}\frac{\varepsilon^{2}}{m^{2}}\Big{)}.$ (41)

Hence, each unit $i\in P$ with $z_{i}=w^{\top}_{i}x+b_{i}$ has a good linear
approximation $\widetilde{z}_{i}=\widetilde{w}^{\top}_{i}x+\widetilde{b}_{i}$
of low error. Combining the above we get that

$\displaystyle\mathbb{E}_{x\sim\mathcal{N}(0,I_{d\times d})}\Big{(}\sum_{i\in
P}a_{i}(z_{i}-\widetilde{z}_{i})\Big{)}^{2}$ $\displaystyle\leq
c^{\prime}\varepsilon^{2},$ (42)

for a constant $c^{\prime}$ that depends on $c$. Furthermore it is easy to see
that the linear approximation $\sum_{i\in P}a_{i}\widetilde{z}_{i}$ is of the
form $\beta^{\top}_{P}x+C_{P}$ where $\|\beta_{P}\|\leq\sum_{i\in
P}|a_{i}|\|w_{i}\|\leq mB$ and $|C_{P}|\leq\sum_{i\in P}|a_{i}b_{i}|\leq
mB^{2}$. ∎

###  5.4 Wrapping up the proofs

With the lemmas above, we can now complete the proof of Theorem 3.1, Theorem
5.1 and Theorem 3.2.

##### Proof of Theorem 5.1

We first set the parameters according to the polynomial bounds from the
different lemmas in this section.

For the final error $\varepsilon$ in approximating $f$, we will set
$\varepsilon^{\prime}\coloneqq\varepsilon/(4mB)$. Also set
$\varepsilon^{\prime}_{3},\eta^{\prime}_{3}$ according to Lemma 5.5 with the
$\varepsilon$ in Lemma 5.5 set to $\varepsilon^{\prime}$. Then set
$\varepsilon_{2}=\varepsilon^{\prime}_{3}$, and also set
$\varepsilon_{1}=\sqrt{\varepsilon_{2}}$. Now we can set the algorithm
parameters $\eta_{3}\coloneqq 2\varepsilon_{2}$, and $\eta_{2}\coloneqq
4\varepsilon_{1}$, and
$\eta_{0}=\min\\{\eta_{3}^{\prime},\eta^{\prime}_{2}\\}$, where
$\eta^{\prime}_{2}$ is given by Lemma 5.4. Moreover $\eta_{1}$ (and
$\eta_{1}^{\prime}$) are set according to Lemma 5.3.

First by using Lemma 5.2 we see that with
$\mathrm{poly}_{\ell}(d,m,B,1/s_{m}(W^{\odot\ell}),1/\varepsilon)$ we can
estimate all the Hermite coefficients up to $2\ell+2$ up to $\eta_{0}$ error
in Frobenius norm. Then, for our setting of parameters we have from Lemma 5.4
that for every $\widetilde{w}_{i}$ for $i\in[m^{\prime}]$ output by steps 1-4
of Algorithm 4, we have that there exists a $w_{i}$ (up to relabeling $i$)
such that $\lVert w_{i}^{\otimes t}-\widetilde{w}_{i}^{\otimes
t}\rVert_{F}\leq\varepsilon_{2}=\varepsilon^{\prime}_{3}<\varepsilon^{\prime}$.
Then we can apply Lemma 5.5 to conclude that for all such terms
$i\in[m^{\prime}]$ that are output we get estimates
$\widetilde{a}_{i},\widetilde{b}_{i}$ with
$|\widetilde{a}_{i}-a_{i}|+|\widetilde{b}_{i}-b_{i}|\leq\varepsilon^{\prime}$.
Moreover using Lemma 5.4 and Lemma 5.5 also show that every $i\in G$ is also
one of the $m^{\prime}$ terms that are output. Hence for each
$i\in\widetilde{G}$, we have recovered each parameter up to error
$\varepsilon^{\prime}$. This completes the proof.

##### Proof of the full-rank setting: Theorem 3.1

The guarantees for Theorem 3.1 hold for the following Algorithm 6, which is a
robust variant of Algorithm 1 in the special case of $\ell=1$. It first uses
Algorithm 4 to approximately recover for each $i\in[m]$, the $a_{i}$, and up
to an ambiguity in a sign (captured by unknown $\xi_{i}\in\\{1,-1\\}$) close
estimates of $w_{i}$ and $b_{i}$. Then it runs Algorithm 3 to disambiguate the
sign by recovering $\xi_{i}$.

Input: Estimates $T_{0},T_{1},T_{2},T_{3},T_{4}$;

Parameters: $\eta_{0},\eta_{1},\eta_{2},\eta_{3}>0$.;

1\. Run Algorithm 4 on parameters $(\eta_{0},\eta_{1},\eta_{2},\eta_{3})$ with
inputs $T_{0},T_{1},T_{2},T_{3},T_{4}$ to receive results
$(\widetilde{w}_{i},\widetilde{a}_{i},\widetilde{b}_{i})_{i\in[m]}$. Note that
$\widetilde{b}_{i},\widetilde{w}_{i}$ are only recovered up to signs.

2\. Run Algorithm 3 (FixSigns) on parameters $m,T_{1}$ and
$(\widetilde{a}_{i},\widetilde{b}_{i},\widetilde{w}_{i}:i\in[m])$ to recover
$(\widetilde{a}_{i},\widetilde{b}_{i},\widetilde{w}_{i}:i\in[m])$.

Result: Output $\\{\widetilde{w}_{i},\widetilde{a}_{i},\widetilde{b}_{i}:1\leq
i\leq m\\}$.

Algorithm 6 Robust full-rank algorithm: recover $\\{a_{i},b_{i},w_{i}\\}$
given estimates $\\{T_{0},T_{1},\dots,T_{4}\\}$.

###### Proof of Theorem 3.1.

We first set the parameters of Algorithm 6 as dictated by Theorem 5.1 (and its
proof) in the special case of $\ell=1$. Let $\varepsilon_{0}>0$ be chosen so
that

$\varepsilon_{0}<\frac{\min\\{\Phi(-c\sqrt{\log(1/\varepsilon
mdB)}),\varepsilon\\}\cdot s_{m}(W)}{8\sqrt{m}B^{2}},$

and $\eta_{0}$ to be the smaller of $\varepsilon_{0}/((1+B)\sqrt{m})$, and
whatever is specified Theorem 5.1 for $\varepsilon=\varepsilon_{0}$. Note that
$\Phi(-c\sqrt{\log(1/(\varepsilon mdB)})\geq\Omega\big{(}(\varepsilon
mdB)^{c^{2}/2}\min\\{1,\varepsilon mdB\\}\big{)}$.

We draw $N=\mathrm{poly}(d,m,B,1/s_{m}(W),1/\varepsilon_{0})$ i.i.d. samples
and run Algorithm 4 with the parameters $\eta_{1},\eta_{2},\eta_{3}$ as
described in the proof of Theorem 5.1. From the assumptions of Theorem 3.1, we
have that each $i\in[m]$ belongs to the “good set” $G$ as well. Hence, from
the guarantee of Theorem 5.1 we will obtain w.h.p. for each $i\in[m]$
estimates $\widetilde{a}_{i},\widetilde{b_{i}},\widetilde{w}_{i}$ (up to
relabeling the indices $[m]$) such that up an unknown sign
$\xi_{i}\in\\{1,-1\\}$ we have

$|a_{i}-\widetilde{a}_{i}|+|\widetilde{b_{i}}-\xi_{i}b_{i}|+\lVert\widetilde{w}_{i}-\xi_{i}w_{i}\rVert_{2}\leq\varepsilon_{0}<\varepsilon.$
(43)

Now consider the (ideal) linear system in the unknowns $\\{z_{i}:i\in[m]\\}$
given by $\widehat{f}_{1}=\sum_{i=1}^{m}z_{i}(a_{i}\xi_{i}w_{i})$; it has $d$
equations in $m\leq d$ unknowns. Let $M\coloneqq
W\text{diag}((\xi_{i}a_{i}:i\in[m]))$ be a $d\times m$ matrix representing the
above linear system as $Mz=\widehat{f}_{1}$. From Lemma 3.5,
$z^{*}_{i}=\xi_{i}\Phi(b_{i})$ is a solution. Moreover $M$ is well-
conditioned: since $|a_{i}|\in[1/B,B]$, we have $s_{1}(M)\leq Bs_{1}(W)\leq
B\sqrt{m}$, while $s_{m}(M)\geq s_{m}(W)/B$ (from the assumption on $W$).
Hence, this is a well-conditioned linear system with a unique solution
$z^{*}$.

Algorithm 3 solves the linear system $\widetilde{M}z=T_{1}$, where
$\widetilde{M}=\widetilde{W}\text{diag}(\widetilde{a})$; here each column of
$\widetilde{M}$ is close to its corresponding column of $M$, while the sample
estimate $T_{1}$ for $\widehat{f}_{1}$ satisfies $\lVert
T_{1}-\widehat{f}_{1}\rVert_{2}\leq\eta_{0}$. Let $\widetilde{z}$ be a
solution to the system $\widetilde{M}z=T_{1}$.

Observe that if $\lVert
z-\widetilde{z}\rVert_{\infty}\leq\|\widetilde{z}-z\|_{2}$ is at most
$\min_{i}|z^{*}_{i}|$, then Algorithm 3 recovers the signs correctly, since
$\widetilde{z}$ will not flip in sign. To calculate this perturbation first
observe that $i$th column of $E=\widetilde{M}-M$ has length at most

$\displaystyle\lVert\widetilde{a}_{i}\xi_{i}\widetilde{w}_{i}-a_{i}w_{i}\rVert_{2}$
$\displaystyle\leq|\widetilde{a}_{i}-a_{i}|\lVert\xi_{i}\widetilde{w}_{i}\rVert_{2}+a_{i}\lVert\xi_{i}\widetilde{w}_{i}-w_{i}\rVert_{2}\leq\varepsilon_{0}+B\varepsilon_{0}\leq\varepsilon_{0}(1+B).$

Hence $\lVert E\rVert_{2}=s_{1}(E)\leq\varepsilon_{0}(1+B)\sqrt{m}$. Moreover
by Weyl’s inequality $s_{m}(\widetilde{M})\geq s_{m}(M)-\lVert
E\rVert\geq\tfrac{1}{B}s_{m}(W)-\varepsilon_{0}(1+B)\geq s_{m}(W)/(2B)$ due to
our choice of parameter $\varepsilon_{0}$. From standard perturbation bounds
for linear systems, we have

$\displaystyle\lVert\widetilde{z}-z^{*}\rVert_{2}$
$\displaystyle\leq\Big{(}s_{m}(\widetilde{M})\Big{)}^{-1}\Big{(}\lVert
T_{1}-\widehat{f}_{1}\rVert_{2}+s_{1}(M-\widetilde{M})\lVert
z^{*}\rVert_{2}\Big{)}$
$\displaystyle\leq\frac{2B}{s_{m}(W)}\Big{(}\eta_{0}+\varepsilon_{0}(1+B)\sqrt{m}\Big{)}\leq\frac{4\varepsilon_{0}B^{2}m}{s_{m}(W)}$
$\displaystyle\leq\Phi(-c\sqrt{\log(1/(\varepsilon
mdB))})\leq\frac{1}{2}\min_{i\in[m]}|z^{*}_{i}|$

as required, due to our choice of $\varepsilon_{0}$. Hence the signs are also
recovered accurately. This along with (43) concludes the proof.

∎

##### Proof of Theorem 3.2

In order to establish Theorem 3.2 we draw
$N=\mathrm{poly}_{\ell}(d,m,B,1/s_{m}(W^{\odot\ell}),1/\varepsilon)$ i.i.d.
samples and run Algorithm 5 with the parameters
$\eta_{0},\eta_{1},\eta_{2},\eta_{3}$ as described in the proof of Theorem
5.1. From the guarantee of Theorem 5.1 we will obtain w.h.p., up to signs,
approximations for all units in $G$ up to an error of
$O(\frac{\varepsilon}{mdB})$. Furthermore, given these approximations the
guarantee of Lemma 5.8 tells us that w.h.p. the function $g(x)$ output by
Algorithm 5 will satisfy $\mathbb{E}_{x\sim\mathcal{N}(0,I_{d\times
d})}\big{(}f(x)-g(x)\big{)}^{2}\leq\varepsilon^{2}$.

##  6 Smoothed Analysis

We use the smoothed analysis framework of Spielman and Teng [ST04], which is a
beyond-worst-case-analysis paradigm that has been used to explain the
practical success of various algorithms. In smoothed analysis, the performance
of the algorithm is measured on a small random perturbation of the input
instance. We use the model studied in the context of parameter estimation and
tensor decomposition problems to obtain polynomial time guarantees under non-
degeneracy conditions [BCMV14, Vij20]. The smoothed analysis model for the
depth-2 neural RELU network setting is as follows:

  1. 1.

An adversary chooses set of parameters $a,b\in\mathbb{R}^{m}$ and
$W\in\mathbb{R}^{d\times m}$.

  2. 2.

The weight matrix $\widehat{W}\in\mathbb{R}^{d\times m}$ is obtained by a
small random i.i.d. perturbation as
$\widehat{W}_{ij}=W_{ij}+\xi_{i,j}~{}\forall i\in[d],j\in[m]$ where
$\xi_{i,j}\sim N(0,\tau^{2}/d)$. (Note that the average squared pertubation in
each column is $\tau^{2}$) 777Think of $\tau$ as a fairly small but inverse
polynomial quantity $1/\mathrm{poly}(n,d)$..

  3. 3.

Each sample $(x,f(x))$ is drawn i.i.d. with $x\sim N(0,I_{d\times d})$ and
$f(x)=a^{\top}\sigma(\widetilde{W}^{\top}x+b)$.

The goal is to design an algorithm that with high probability, estimates the
parameters $a,b,\widehat{W}$ up to some desired accuracy $\varepsilon$ in time
$\mathrm{poly}(m,d,1/\varepsilon,1/\tau)$. We now prove the following
corollary of Theorem 3.2.

##### Corollary 3.3

Suppose $\ell\in\mathbb{N}$ and $\varepsilon>0$ are constants in the smoothed
analysis model with smoothing parameter $\tau>0$, and also assume the ReLU
network $f(x)=a^{\top}\sigma(\widehat{W}^{\top}x+b)$ is $B$-bounded with
$m\leq 0.99\binom{d+\ell-1}{\ell}$. Then there is an algorithm that given
$N\geq\mathrm{poly}_{\ell}(m,d,1/\varepsilon,B,1/\tau)$ samples runs in
$\mathrm{poly}(N,m,d)$ time and with high probability finds a ReLU network
$g(x)={a^{\prime}}^{\top}\sigma({W^{\prime}}^{\top}x+{b^{\prime}})$ with at
most $m+2$ hidden units such that the $L_{2}$ error
$\operatorname*{\mathbb{E}}_{x\sim\mathcal{N}(0,I_{d\times
d})}[(f(x)-{g}(x))^{2}]\leq\varepsilon^{2}$. Furthermore there are constants
$c,c^{\prime}>0$ and signs $\xi_{i}\in\\{\pm 1\\}~{}\forall i\in[m]$, such
that in $\mathrm{poly}(N,m,d)$ time, for all $i\in[m]$ with
$|b_{i}|<c\sqrt{\log(1/(\varepsilon\cdot mdB))}$, we can recover
$(\widetilde{a}_{i},\widetilde{w}_{i},\widetilde{b}_{i})$, such that
$|a_{i}-\widetilde{a}_{i}|+\lVert
w_{i}-\xi_{i}\widetilde{w}_{i}\rVert_{2}+|b_{i}-\xi_{i}\widetilde{b}_{i}|<c^{\prime}\varepsilon/(mB)$.

###### Proof.

The proof of the corollary follows by combining Theorem 3.2 with existing
results on smoothed analysis [BCPV19] on the least singular value
$s_{m}(\widehat{W}^{\odot\ell})$. We apply Theorem 2.1 of [BCPV19] with
$\rho=\tau$, $U$ being the identity matrix to derive that for any $\delta>0$
and $m\leq(1-\delta)\binom{d+\ell-1}{\ell}$, we get with probability at least
$1-m\exp(-\Omega_{\ell}(\delta n))$ that

$s_{m}(\widehat{W}^{\odot\ell})\geq\frac{c_{\ell}}{\sqrt{m}}\Big{(}\frac{\tau}{d}\Big{)}^{\ell}.$

We then just apply Theorem 3.2 to conclude the proof.

∎

##  7 Conclusion

In this paper, we designed polynomial time algorithms for learning depth-2
neural networks with general ReLU activations (with non-zero bias terms), and
gave provable guarantees under mild non-degeneracy conditions. The results of
this work are theoretical in nature, in trying to understand whether efficient
algorithms exist for learning ReLU networks; hence we believe they do not have
any adverse societal impact. We addressed multiple challenges for learning
such ReLU network with non-zero bias terms throughout our analyses, that may
be more broadly useful in handling bias terms in the ReLU activations. We also
proved identifiability under minimal assumptions and adopted the framework of
smoothed analysis to establish beyond-worst-case guarantees. The major open
direction is to provide similar guarantees for networks of higher depth.

##  8 Acknowledgement

We thank Ainesh Bakshi for pointing us to an error in the previous version of
the paper.

## References

  * [ADH+19] Sanjeev Arora, Simon Du, Wei Hu, Zhiyuan Li, and Ruosong Wang.  Fine-grained analysis of optimization and generalization for overparameterized two-layer neural networks.  In International Conference on Machine Learning, pages 322–332. PMLR, 2019.
  * [APVZ14] Alexandr Andoni, Rina Panigrahy, Gregory Valiant, and Li Zhang.  Learning polynomials with neural networks.  In International conference on machine learning, pages 1908–1916. PMLR, 2014.
  * [AZLS19] Zeyuan Allen-Zhu, Yuanzhi Li, and Zhao Song.  A convergence theory for deep learning via over-parameterization.  In International Conference on Machine Learning, pages 242–252. PMLR, 2019.
  * [BBV04] Stephen Boyd, Stephen P Boyd, and Lieven Vandenberghe.  Convex optimization.  Cambridge university press, 2004.
  * [BCMV14] A. Bhaskara, M. Charikar, A. Moitra, and A. Vijayaraghavan.  Smoothed analysis of tensor decompositions.  In Proceedings of the Forty-Sixth Annual ACM Symposium on Theory of Computing, STOC ’14, page 594–603, New York, NY, USA, 2014. Association for Computing Machinery.
  * [BCPV19] Aditya Bhaskara, Aidao Chen, Aidan Perreault, and Aravindan Vijayaraghavan.  Smoothed analysis in unsupervised learning via decoupling.  In Proceedings of the 60th Annual IEEE Symposium on Foundations of Computer Science (FOCS). IEEE, 2019.
  * [BCV14] Aditya Bhaskara, Moses Charikar, and Aravindan Vijayaraghavan.  Uniqueness of tensor decompositions with applications to polynomial identifiability.  In Maria Florina Balcan, Vitaly Feldman, and Csaba Szepesvári, editors, Proceedings of The 27th Conference on Learning Theory, volume 35 of Proceedings of Machine Learning Research, pages 742–778, Barcelona, Spain, 2014. PMLR.
  * [Bha97] Rajendra Bhatia.  Matrix Analysis, volume 169.  Springer, 1997.
  * [BJW19] Ainesh Bakshi, Rajesh Jayaram, and David P Woodruff.  Learning two layer rectified neural networks in polynomial time.  In Conference on Learning Theory, pages 195–268. PMLR, 2019.
  * [CB18] Lénaïc Chizat and Francis Bach.  On the global convergence of gradient descent for over-parameterized models using optimal transport.  In Proceedings of the 32nd International Conference on Neural Information Processing Systems, NIPS’18, page 3040–3050, Red Hook, NY, USA, 2018. Curran Associates Inc.
  * [CKM20] Sitan Chen, Adam R Klivans, and Raghu Meka.  Learning deep relu networks is fixed-parameter tractable.  arXiv preprint arXiv:2009.13512, 2020.
  * [Dan17] Amit Daniely.  Sgd learns the conjugate kernel class of the network.  In Proceedings of the 31st International Conference on Neural Information Processing Systems, NIPS’17, page 2419–2427, Red Hook, NY, USA, 2017. Curran Associates Inc.
  * [DFS16] Amit Daniely, Roy Frostig, and Yoram Singer.  Toward deeper understanding of neural networks: The power of initialization and a dual view on expressivity.  In D. Lee, M. Sugiyama, U. Luxburg, I. Guyon, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 29. Curran Associates, Inc., 2016.
  * [DK20] Ilias Diakonikolas and Daniel M. Kane.  Small covers for near-zero sets of polynomials and learning latent variable models.  In 61st IEEE Annual Symposium on Foundations of Computer Science, FOCS 2020, Durham, NC, USA, November 16-19, 2020, pages 184–195. IEEE, 2020.
  * [DKKZ20] Ilias Diakonikolas, Daniel M Kane, Vasilis Kontonis, and Nikos Zarifis.  Algorithms and sq lower bounds for pac learning one-hidden-layer relu networks.  In Conference on Learning Theory, pages 1514–1539. PMLR, 2020.
  * [DZPS19] Simon S. Du, Xiyu Zhai, Barnabás Póczos, and Aarti Singh.  Gradient descent provably optimizes over-parameterized neural networks.  In 7th International Conference on Learning Representations, ICLR 2019, New Orleans, LA, USA, May 6-9, 2019. OpenReview.net, 2019.
  * [GGJ+20] Surbhi Goel, Aravind Gollakota, Zhihan Jin, Sushrut Karmalkar, and Adam Klivans.  Superpolynomial lower bounds for learning one-layer neural networks using gradient descent.  In International Conference on Machine Learning, pages 3587–3596. PMLR, 2020.
  * [GK19] Surbhi Goel and Adam R Klivans.  Learning neural networks with two nonlinear layers in polynomial time.  In Conference on Learning Theory, pages 1470–1499. PMLR, 2019.
  * [GKLW19] Rong Ge, Rohith Kuditipudi, Zhize Li, and Xiang Wang.  Learning two-layer neural networks with symmetric inputs.  In International Conference on Learning Representations, 2019.
  * [GLM18] Rong Ge, Jason D. Lee, and Tengyu Ma.  Learning one-hidden-layer neural networks with landscape design.  In International Conference on Learning Representations, 2018.
  * [GVX14] Navin Goyal, Santosh Vempala, and Ying Xiao.  Fourier pca and robust tensor decomposition.  STOC ’14, page 584–593, New York, NY, USA, 2014. Association for Computing Machinery.
  * [Har70] R. Harshman.  Foundations of the parafac procedure: Models and conditions for an "explanatory" multi-model factor analysis.  1970\.
  * [JGH21] Arthur Jacot, Franck Gabriel, and Clément Hongler.  Neural tangent kernel: Convergence and generalization in neural networks (invited paper).  In Proceedings of the 53rd Annual ACM SIGACT Symposium on Theory of Computing, STOC 2021, page 6, New York, NY, USA, 2021. Association for Computing Machinery.
  * [JGKA19] Majid Janzamin, Rong Ge, Jean Kossaifi, and Animashree Anandkumar.  Spectral learning on matrices and tensors.  Foundations and Trends in Machine Learning, 12, 11 2019.
  * [JSA15] Majid Janzamin, Hanie Sedghi, and Anima Anandkumar.  Beating the perils of non-convexity: Guaranteed training of neural networks using tensor methods.  arXiv preprint arXiv:1506.08473, 2015.
  * [LXS+19] Jaehoon Lee, Lechao Xiao, Samuel Schoenholz, Yasaman Bahri, Roman Novak, Jascha Sohl-Dickstein, and Jeffrey Pennington.  Wide neural networks of any depth evolve as linear models under gradient descent.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
  * [Moi18] Ankur Moitra.  Algorithmic Aspects of Machine Learning.  Cambridge University Press, 2018.
  * [MRT18] Mehryar Mohri, Afshin Rostamizadeh, and Ameet Talwalkar.  Foundations of machine learning.  MIT press, 2018.
  * [O’D14] Ryan O’Donnell.  Analysis of Boolean Functions.  Cambridge University Press, USA, 2014.
  * [P+15] Iosif Pinelis et al.  Exact rosenthal-type bounds.  Annals of Probability, 43(5):2511–2544, 2015.
  * [Ros70] Haskell P Rosenthal.  On the subspaces of l p (p> 2) spanned by sequences of independent random variables.  Israel Journal of Mathematics, 8(3):273–303, 1970.
  * [Rou20] Tim Roughgarden.  Beyond the Worst-Case Analysis of Algorithms.  Cambridge University Press, 2020.
  * [ST04] Daniel A. Spielman and Shang-Hua Teng.  Smoothed analysis of algorithms: Why the simplex algorithm usually takes polynomial time.  J. ACM, 51(3):385–463, 2004.
  * [Tur50] Paul Turán.  On the zeros of the polynomials of legendre.  Časopis pro pěstování matematiky a fysiky, 75(3):113–122, 1950.
  * [Vij20] Aravindan Vijayaraghavan.  Beyond the Worst-Case Analysis of Algorithms, chapter Efficient Tensor Decomposition.  In [Rou20], 2020.

##  Appendix A Expressions for the Hermite Coefficients

##### Lemma 3.5

The $k$’th Hermite expansion of $f(x)=a^{\mathsf{T}}\sigma(W^{\top}x+b)$,
$\hat{f}_{k}$, when $k=0,1$, is

$\hat{f}_{0}=\sum_{i=1}^{m}a_{i}[b_{i}\Phi(b_{i})+\frac{\exp(-\frac{b_{i}^{2}}{2})}{\sqrt{2\pi}}]~{},~{}\hat{f}_{1}=\sum_{i=1}^{m}a_{i}\Phi(b_{i})w_{i}$
(44)

when $k\geq 2$, the coefficients are

$\hat{f}_{k}=\sum_{i=1}^{m}(-1)^{k}\cdot a_{i}\cdot
He_{k-2}(b_{i})\cdot\frac{\exp(\frac{-b_{i}^{2}}{2})}{\sqrt{2\pi}}\cdot
w_{i}^{\otimes k}$ (45)

where the expectation is taken over $x\sim\mathcal{N}(0,I_{d})$ and
$\hat{f}_{k}$ is a $k$’th-order tensor.

###### Proof.

Note that Hermite polynomials can be written in terms of their generating
function [O’D14]

$He_{k}(x)=D^{(k)}_{t}\exp(t^{\mathsf{T}}x-\|t\|^{2}/2)|_{t=\mathbf{0}}$ (46)

Hence we can write $\hat{f}_{k}$ as

$\hat{f}_{k}=\mathbb{E}[f(x)He_{k}(x)]=\sum_{i=1}^{m}a_{i}\int_{w_{i}x+b_{i}\geq
0}(w_{i}^{\mathsf{T}}x+b_{i})\cdot
D^{(k)}_{t}\exp(t^{\mathsf{T}}x-\frac{\|t\|^{2}}{2})~{}d\mu|_{t=\mathbf{0}}$
(47)

where $d\mu=\frac{\exp(-\|x\|^{2}/2)}{(\sqrt{2\pi})^{d}}dx$ is the Gaussian
probability measure. Moving $d\mu$ into the exponential term, we get

$\hat{f}_{k}=\sum_{i=1}^{m}\frac{a_{i}}{(\sqrt{2\pi})^{d}}\int_{w_{i}x+b_{i}\geq
0}(w_{i}^{\mathsf{T}}x+b_{i})\cdot
D^{(k)}_{t}\exp(-\frac{\|x\|^{2}}{2}+t^{\mathsf{T}}x-\frac{\|t\|^{2}}{2})~{}dx|_{t=\mathbf{0}}$
(48) $=\sum_{i=1}^{m}\frac{a_{i}}{(\sqrt{2\pi})^{d}}\int_{w_{i}x+b_{i}\geq
0}(w_{i}^{\mathsf{T}}x+b_{i})\cdot
D^{(k)}_{t}\exp(-\frac{\|x-t\|^{2}}{2})~{}dx|_{t=\mathbf{0}}$ (49)
$=D^{(k)}_{t}[\sum_{i=1}^{m}\frac{a_{i}}{(\sqrt{2\pi})^{d}}\int_{w_{i}x+b_{i}\geq
0}(w_{i}^{\mathsf{T}}x+b_{i})\cdot\exp(-\frac{\|x-t\|^{2}}{2})~{}dx]_{t=\mathbf{0}}$
(50)

Denote

$I_{i}(t)=\frac{1}{(\sqrt{2\pi})^{d}}\int_{w_{i}x+b_{i}\geq
0}(w_{i}^{\mathsf{T}}x+b_{i})\cdot\exp(-\frac{\|x-t\|^{2}}{2})~{}dx$ (51)

then

$\hat{f}_{k}=D^{(k)}_{t}\sum_{i=1}^{m}a_{i}I_{i}(t)|_{t=\mathbf{0}}$ (52)

Now, let $y_{i}=w_{i}^{\mathsf{T}}x$ and $t_{w_{i}}=t^{\mathsf{T}}w_{i}$. To
evaluate $I_{i}(t)$ in terms of $y_{i}$, it suffices to only consider the
projection of $t$ onto $w_{i}$, $t_{w_{i}}$ with the remaining parts being
integrated out. Hence, we can rewrite $I_{i}(t)$ as

$I_{i}(t)=\frac{1}{(\sqrt{2\pi})^{d}}\int_{x^{\prime}\in\mathbb{R}^{d-1}}\exp(-\frac{\|x^{\prime}-t^{\prime}\|^{2}}{2})dx^{\prime}~{}\int_{y_{i}=-b_{i}}^{\infty}(y_{i}+b_{i})\cdot\exp(-\frac{\|y_{i}-t_{w_{i}}\|^{2}}{2})dy_{i}$
(53)
$=\frac{1}{\sqrt{2\pi}}\int_{y_{i}=-b_{i}}^{\infty}(y_{i}+b_{i})\cdot\exp(-\frac{\|y_{i}-t_{w_{i}}\|^{2}}{2})dy_{i}$
(54)
$=(t_{w_{i}}+b_{i})\Phi(t_{w_{i}}+b_{i})+\frac{\exp(-\frac{(t_{w_{i}}+b_{i})^{2}}{2})}{\sqrt{2\pi}}$
(55)

where $\Phi(z)$ is the standard Gaussian c.d.f. of $z$. We then have

$\hat{f}_{k}=\sum_{i=1}^{m}a_{i}\cdot
D^{(k)}_{t}[(t_{w_{i}}+b_{i})\Phi(t_{w_{i}}+b_{i})+\frac{\exp(-\frac{(t_{w_{i}}+b_{i})^{2}}{2})}{\sqrt{2\pi}}]_{t=\mathbf{0}}$
(56)

Therefore we have $\hat{f}_{0}$ and $\hat{f}_{1}$ as

$\hat{f}_{0}=\sum_{i=1}^{m}a_{i}[b_{i}\Phi(b_{i})+\frac{\exp(-\frac{b_{i}^{2}}{2})}{\sqrt{2\pi}}]~{},~{}\hat{f}_{1}=\sum_{i=1}^{m}a_{i}\Phi(b_{i})w_{i}$
(57)

Since we are taking the derivative with respect to $t$, for some function
$g(t_{w_{i}})$, by the chain rule we will have

$D^{(k)}_{t}g=\frac{d^{k}g}{dt_{w_{i}}^{k}}\cdot w_{i}^{\otimes k}$ (58)

Finally, recall Fact 1.3, the derivatives of a Gaussian p.d.f. can be
expressed in terms of Hermite polynomials, hence for $k\geq 2$

$\hat{f}_{k}=\sum_{i=1}^{m}a_{i}\cdot\frac{d^{k-2}}{dt_{w_{i}}^{k-2}}\frac{\exp(-\frac{(b_{i}+t_{w_{i}})^{2}}{2})}{\sqrt{2\pi}}\cdot
w_{i}^{\otimes k}|_{t_{w_{i}}=0}$ (59) $=\sum_{i=1}^{m}(-1)^{k-2}\cdot
a_{i}\cdot
He_{k-2}(b_{i}+t_{w_{i}})\cdot\frac{\exp(-\frac{(b_{i}+t_{w_{i}})^{2}}{2})}{\sqrt{2\pi}}\cdot
w_{i}^{\otimes k}|_{t_{w_{i}}=0}$ (60) $=\sum_{i=1}^{m}(-1)^{k}\cdot
a_{i}\cdot
He_{k-2}(b_{i})\cdot\frac{\exp(\frac{-b_{i}^{2}}{2})}{\sqrt{2\pi}}\cdot
w_{i}^{\otimes k}$ (61)

which proves the lemma. ∎

##### Proposition 3.6

Let
$\widetilde{f}(x)=\widetilde{a}^{\top}\sigma(\widetilde{W}^{\top}x+\widetilde{b})$
be the model trained using samples generated by the ground-truth ReLU network
$f(x)=a^{\top}\sigma(W^{\top}x+b)$. Then the statistical risk with respect to
the $\ell_{2}$ loss function can be expressed as follows

$\displaystyle L(\widetilde{a},\widetilde{b},\widetilde{W})$
$\displaystyle=\sum_{k\in\mathbb{N}}\frac{1}{k!}\Big{\|}T_{k}-\hat{f}_{k}\Big{\|}_{F}^{2}$
$\displaystyle~{}~{}\text{ where }~{}~{}T_{0}$
$\displaystyle=\sum_{i=1}^{m}\widetilde{a}_{i}(\widetilde{b}_{i}\Phi(\widetilde{b}_{i})+\frac{\exp(-\widetilde{b}_{i}^{2}/2)}{\sqrt{2\pi}}),~{}\text{
and
}~{}T_{1}=\sum_{i=1}^{m}\widetilde{a}_{i}\Phi(\widetilde{b}_{i})\widetilde{w}_{i}$
$\displaystyle\forall k\geq 2,~{}~{}T_{k}$
$\displaystyle=\sum_{i=1}^{m}(-1)^{k}\cdot\widetilde{a}_{i}\cdot
He_{k-2}(\widetilde{b_{i}})\cdot\frac{\exp(-\widetilde{b_{i}}^{2}/2)}{\sqrt{2\pi}}\cdot\widetilde{w}_{i}^{\otimes
k}$

###### Proof.

Let $\alpha\bowtie\alpha^{\prime}$ denote $\alpha$ being a permutation of
$\alpha^{\prime}$. Since $\bowtie$ is an equivalence relation, we can
partition $[d]^{*}$ into equivalence classes ($*$ is the Kleene star operator)
such that for some $\alpha\in[d]^{*}$,
$[\alpha]=\\{\alpha^{\prime}\in[d]^{*}~{}|~{}\alpha\bowtie\alpha^{\prime}\\}$.
Let $C$ be a subset of $[d]^{*}$ such that no pair of
$\alpha,\alpha^{\prime}\in C$ is in the same equivalence class. We can then
directly decompose the statistical risk as

$\mathbb{E}\Big{[}|\widetilde{f}(x)-f(x)|^{2}\Big{]}=\mathbb{E}\Big{[}\Big{|}\sum_{\alpha\in
C}\frac{T_{\alpha}He_{\alpha}(x)}{c_{\alpha}^{2}}-\sum_{\alpha\in
C}\frac{\hat{f}_{\alpha}He_{\alpha}(x)}{c_{\alpha}^{2}}\Big{|}^{2}\Big{]}$
(62)

where $c_{\alpha}^{2}=\mathbb{E}[He_{\alpha}(x)^{2}]$. Note that we omit $k$
here and directly write the Hermite polynomial obtained by differentiating
with respect to $x_{\alpha_{1}},...,x_{\alpha_{k}}$ as
$He_{\alpha}(x)\in\mathbb{R}$. The above equation can thus be further
simplified as

$\mathbb{E}\Big{[}\sum_{\alpha\in
C}\Big{(}\frac{(T_{\alpha}-\hat{f}_{\alpha})He_{\alpha}(x)}{c_{\alpha}^{2}}\Big{)}^{2}+\sum_{\begin{subarray}{c}\alpha\neq\alpha^{\prime}\\\
\alpha,\alpha^{\prime}\in
C\end{subarray}}\Big{(}\frac{(T_{\alpha}-\hat{f}_{\alpha})He_{\alpha}(x)}{c_{\alpha}^{2}}\Big{)}\Big{(}\frac{(T_{\alpha^{\prime}}-\hat{f}_{\alpha^{\prime}})He_{\alpha^{\prime}}(x)}{c_{\alpha^{\prime}}^{2}}\Big{)}\Big{]}$
(63) $=\sum_{\alpha\in
C}\frac{(T_{\alpha}-\hat{f}_{\alpha})^{2}}{c_{\alpha}^{2}}\mathbb{E}\Big{[}\frac{He_{\alpha}(x)^{2}}{c_{\alpha}^{2}}\Big{]}=\sum_{\alpha\in
C}\frac{(T_{\alpha}-\hat{f}_{\alpha})^{2}}{c_{\alpha}^{2}}$ (64)

since if both $\alpha,\alpha^{\prime}\in C$ and $\alpha\neq\alpha^{\prime}$,
$\mathbb{E}[He_{\alpha}(x)He_{\alpha^{\prime}}(x)]=0$. Next, we rewrite the
expression as

$\sum_{k\in\mathbb{N}}\sum_{\begin{subarray}{c}\alpha\in C\\\
|\alpha|=k\end{subarray}}\frac{(T_{\alpha}-\hat{f}_{\alpha})^{2}}{c_{\alpha}^{2}}=\sum_{k\in\mathbb{N}}\frac{1}{k!}\sum_{\begin{subarray}{c}\alpha\in
C\\\
|\alpha|=k\end{subarray}}\frac{k!}{c_{\alpha}^{2}}(T_{\alpha}-\hat{f}_{\alpha})^{2}=\sum_{k\in\mathbb{N}}\frac{1}{k!}\Big{\|}T_{k}-\hat{f}_{k}\Big{\|}_{F}^{2}$
(65)

The last equality is due to the fact that
$c_{\alpha}^{2}=\prod_{i=1}^{d}n_{i}!$, where $n_{i}$ is the number of times
that $i$ occurs in the multi-index $\alpha$, and therefore $k!/c_{\alpha}^{2}$
is the number of possible permutations of the elements in $\alpha$ with
$|\alpha|=k$ subjecting to 1 occurs $n_{1}$ times, 2 occurs $n_{2}$ times, …,
$d$ occurs $n_{d}$ times. Thus the proposition follows. ∎

##  Appendix B Robust Analysis for general $\ell$

In this section, we prove that Algorithm 4 and its algorithmic guarantee in
Theorem 5.1 (and Theorems 3.1 and 3.2).

We break down the proof into multiple parts.

###  B.1 Estimating the Hermite coefficients

To obtain the desired concentration bound, we first introduce an auxiliary
claim we will make use of in the following analysis.

######  Claim B.1.

For $a_{1},a_{2},...,a_{n}\in\mathbb{R}$ and $p\in\mathbb{N}$,
$|\sum_{i\in[n]}a_{i}|^{2p}\leq n^{2p}\max_{i\in[n]}|a_{i}|^{2p}\leq
n^{2p}\sum_{i\in[n]}|a_{i}|^{2p}$

###### Proof.

By triangle inequality,

$|\sum_{i\in[n]}a_{i}|\leq\sum_{i\in[n]}|a_{i}|\leq n\max_{i\in[n]}|a_{i}|\leq
n\sum_{i\in[n]}|a_{i}|\Rightarrow|\sum_{i\in[n]}a_{i}|^{2p}\leq
n^{2p}\max_{i\in[n]}|a_{i}|^{2p}\leq n^{2p}\sum_{i\in[n]}|a_{i}|^{2p}$

∎

Equipped with the essential claim, we are now ready to prove Lemma 5.2.

##### Lemma 5.2

For any $\eta>0$, if $T_{k}$ is estimated from $N\geq
c_{k}d^{k}m^{2}B^{4}\mathrm{poly}(\log(mdB/\eta))/\eta^{2}$ samples, then for
some constant $c_{k}>0$ that depends only on $k$, we have with probability at
least $1-(mdB)^{-\log(md)}$,

$\|T_{k}-\hat{f}_{k}\|_{F}\leq\eta.$ (66)

###### Proof.

Consider $p\in\mathbb{N}$, and a sum $S_{Y}=\sum_{j=1}^{N}Y_{j}$ of
independent zero-mean r.v.s with
$\tfrac{1}{N}\sum_{j=1}^{N}\operatorname*{\mathbb{E}}[Y_{j}^{2p}]\leq A_{2p}$
and $\tfrac{1}{N}\sum_{j=1}^{N}\operatorname*{\mathbb{E}}[Y_{j}^{2}]\leq
A_{2}$. Then by Rosenthal’s inequality (and Markov’s inequality)

$\displaystyle\operatorname*{\mathbb{E}}\Big{[}\Big{(}\sum_{j=1}^{N}Y_{j}\Big{)}^{2p}\Big{]}$
$\displaystyle\leq 2^{p\log(p)+2p+p^{2}}\cdot\max\\{NA_{2p},(NA_{2})^{p}\\}$
(67) $\displaystyle\text{And,
}\operatorname*{\mathbb{P}}\Big{[}\Big{|}\frac{1}{N}\sum_{j=1}^{N}Y_{j}\Big{|}>\eta\Big{]}$
$\displaystyle\leq
2^{p\log(p)+2p+p^{2}}\cdot\max\Big{\\{}\frac{A_{2p}}{N^{2p-1}\eta^{2p}},\Big{(}\frac{A_{2}}{N\eta^{2}}\Big{)}^{p}\Big{\\}}.$
(68)

Consider a fixed $\alpha\in[d]^{k}$ (an index of the tensor corresponding to
the $k$th Hermite coefficient); $|\alpha|=k$. Given samples
$\\{(x^{(j)},f(x^{(j)}):j\in[N]\\}$, the random variables of interest are
$Z_{j},Y_{j}$ are

$Z_{j}=\sum_{i=1}^{m}a_{i}\sigma(w_{i}^{\top}x^{(j)}+b_{i})He_{\alpha}(x^{(j)}),~{}~{}~{}~{}\text{
and }~{}~{}~{}Y_{j}\coloneqq Z_{j}-\operatorname*{\mathbb{E}}[Z_{j}].$

We will apply the above concentration inequality with the random variables
$Y_{j}$. We need bounds for $\operatorname*{\mathbb{E}}[Y_{j}^{2}]$ and
$\operatorname*{\mathbb{E}}[Y_{j}^{2p}]$. For convenience let
$Z\coloneqq\sum_{i=1}^{m}a_{i}\sigma(w_{i}^{\top}x+b_{i})He_{\alpha}(x)$, and
$Y\coloneqq Z-\operatorname*{\mathbb{E}}[Z]$. Note that by applying Claim B.1,
we can bound these quantities as

$\displaystyle\operatorname*{\mathbb{E}}[Y^{2p}]$
$\displaystyle=\operatorname*{\mathbb{E}}[(Z-\operatorname*{\mathbb{E}}[Z])^{2p}]\leq
2^{2p}(\operatorname*{\mathbb{E}}[Z^{2p}]+\operatorname*{\mathbb{E}}[Z]^{2p}),~{}~{}\text{where}$
$\displaystyle|\operatorname*{\mathbb{E}}[Z]|$
$\displaystyle=\Big{|}\sum_{i=1}^{m}a_{i}He_{k-2}(b_{i})\cdot\frac{\exp(-b_{i}^{2}/2)}{\sqrt{2\pi}}\cdot\prod_{t=1}^{k}w_{i}(\alpha(t))\Big{|}\leq
mB\sqrt{k!},$ $\displaystyle\operatorname*{\mathbb{E}}[Z^{2p}]$
$\displaystyle=\operatorname*{\mathbb{E}}\Big{[}\Big{(}\sum_{i=1}^{m}a_{i}\sigma(w_{i}^{\top}x+b_{i})\Big{)}^{2p}He_{\alpha}(x)^{2p}\Big{]},$

On the other hand, from Hölder’s inequality, we have
$(\sum_{i=1}^{m}|c_{i}||z_{i}|)^{2p}\leq(\lVert
c\rVert_{q^{*}}^{q^{*}})^{2p/q^{*}}\cdot\lVert z\rVert_{2p}^{2p}$ where
$q^{*}$ is the dual norm of $2p$ i.e., $2p/q^{*}=2p-1$. Hence, again combined
with Claim B.1, we have

$\displaystyle\operatorname*{\mathbb{E}}[Z^{2p}]$
$\displaystyle\leq\operatorname*{\mathbb{E}}\Big{[}\Big{(}\sum_{i=1}^{m}a_{i}^{q^{*}}\Big{)}^{2p/q^{*}}\Big{(}\sum_{i=1}^{m}\sigma(w_{i}^{\top}x+b_{i})^{2p}\Big{)}He_{\alpha}(x)^{2p}\Big{]}\leq(m^{2p-1}B^{2p})\sum_{i=1}^{m}\operatorname*{\mathbb{E}}\Big{[}(w_{i}^{\top}x+b_{i})^{2p}He_{\alpha}(x)^{2p}\Big{]}$
$\displaystyle\leq(2^{2p}m^{2p-1}B^{2p})\sum_{i=1}^{m}\operatorname*{\mathbb{E}}\Big{[}((w_{i}^{\top}x)^{2p}+b_{i}^{2p})He_{\alpha}(x)^{2p}\Big{]}$
$\displaystyle\leq(2mB)^{2p}\sum_{i=1}^{m}\Big{(}\operatorname*{\mathbb{E}}\Big{[}(w_{i}^{\top}x)^{2p}He_{\alpha}(x)^{2p}\Big{]}+B^{2p}\operatorname*{\mathbb{E}}\Big{[}He_{\alpha}(x)^{2p}\Big{]}\Big{)}$

We note that $w_{i}^{\top}x$ is a standard Gaussian since $\lVert
w_{i}\rVert_{2}=1$.

Now, let $He_{\alpha}(x)$ involve different indices of $x$ up to
$k_{1},k_{2},\dots,k_{d}$ times. Note that
$\sum_{t\in[d]}k_{t}\leq|\alpha|=k$. Using properties of Hermite polynomials,
we can bound $\operatorname*{\mathbb{E}}[He_{\alpha}(x)^{2p}]$ as

$\operatorname*{\mathbb{E}}[He_{\alpha}(x)^{2p}]=\operatorname*{\mathbb{E}}\Big{[}\Big{(}\sum_{\begin{subarray}{c}\sum_{t\in[d]}k_{t}\leq
k\end{subarray}}c_{k_{1}...k_{d}}\prod_{t\in[d]}x_{t}^{k_{t}}\Big{)}^{2p}\Big{]}$
$\leq\binom{d+k}{d}^{2p}\max_{\begin{subarray}{c}\sum_{t\in[d]}k_{t}\leq
k\end{subarray}}c_{k_{1}...k_{d}}^{2p}\operatorname*{\mathbb{E}}\Big{[}\prod_{t\in[d]}x_{t}^{2pk_{t}}\Big{]}\leq\binom{d+k}{d}^{2p}(2pk-1)!!(k!)^{2p}$
$\leq\Big{(}\binom{d+k}{d}\cdot
k!\Big{)}^{2p}(2pk)^{pk}=C_{1}^{2p}\cdot(2pk)^{pk}$

by setting $C_{1}=\binom{d+k}{d}\cdot k!$ and repetitively applying Claim B.1.
A similar argument also holds for
$\operatorname*{\mathbb{E}}[(w_{i}^{\top}x)^{2p}He_{\alpha}(x)^{2p}]$ by
Cauchy–Schwarz inequality

$\operatorname*{\mathbb{E}}[(w_{i}^{\top}x)^{2p}He_{\alpha}(x)^{2p}]\leq\sqrt{\operatorname*{\mathbb{E}}[(w_{i}^{\top}x)^{4p}]\operatorname*{\mathbb{E}}[He_{\alpha}(x)^{4p}]}\leq\sqrt{\operatorname*{\mathbb{E}}[(w_{i}^{\top}x)^{4p}]}\binom{d+k}{d}^{2p}(k!)^{2p}(4pk)^{pk}$
$\leq\Big{(}\sqrt{4p}\cdot\binom{d+k}{d}\cdot
k!\Big{)}^{2p}(4pk)^{pk}=(2^{k+2}p)^{p}\cdot C_{1}^{2p}\cdot(2pk)^{pk}$

since $w_{i}^{\top}x$ follows a standard Gaussian distribution. Hence we have

$\operatorname*{\mathbb{E}}[Z^{2p}]\leq
C_{1}^{2p}(2mB)^{2p}(2pk)^{pk}((2^{k+2}p)^{p}+B^{2p})m$

$A_{2p}=\operatorname*{\mathbb{E}}[Y^{2p}]\leq
2^{2p}\Big{(}(2mBC_{1})^{2p}(2pk)^{pk}((2^{k+2}p)^{p}+B^{2p})m+(mB\sqrt{k!})^{2p}\Big{)}$

$\Longrightarrow A_{2p}\leq(C_{2}(pk)^{k/2}mB^{2})^{2p}$

where $C_{2}=8C_{1}^{2k}$. Note that $p=1$ also gives the required bounds for
$\operatorname*{\mathbb{E}}[Y^{2}]$.

Now, setting $p:=\tfrac{1}{2}(\log(1/\eta)+\log(mdB))$, and applying
Rosenthal’s inequality (67) with
$N=\tfrac{c^{\prime}(k)}{\eta^{2}}m^{2}B^{4}\mathrm{poly}(\log(mdB/\eta))$, we
have for an appropriate constant $c^{\prime}(k)>0$

$\displaystyle\operatorname*{\mathbb{P}}\Big{[}\Big{|}\frac{1}{N}\sum_{j=1}^{N}Y_{j}\Big{|}>\eta\Big{]}$
$\displaystyle\leq
2^{p\log(p)+2p+p^{2}}\cdot\max\Big{\\{}\frac{A_{2p}}{N^{2p-1}\eta^{2p}},\Big{(}\frac{A_{2}}{N\eta^{2}}\Big{)}^{p}\Big{\\}}$
$\displaystyle\leq
2^{p\log(p)+2p+p^{2}}\cdot\max\Big{\\{}\frac{(C_{2}mB^{2})^{2p}(pk)^{pk}}{N^{2p-1}\eta^{2p}},\Big{(}\frac{(C_{2}mB^{2}k^{k/2})^{2}}{N\eta^{2}}\Big{)}^{p}\Big{\\}}$
$\displaystyle\leq\Big{(}\frac{1}{mdB}\Big{)}^{\log(mdB)+\log(1/\eta)},$

as required. Finally by setting $\eta=\eta^{\prime}/\sqrt{d^{k}}$ and a union
bound over all entries, we get that w.h.p., $\lVert
T_{k}-\hat{f}_{k}\rVert_{F}\leq\eta^{\prime}$, as required.

∎

###  B.2 Recovering the parameters under errors

####  B.2.1 Recovery of weight vectors $w_{i}$ for the terms in $G$

We first prove the following important lemma that shows that Jennrich’s
algorithm run with an appropriate choice of rank $k$ will recover the large
terms.

##### Lemma 5.3

Suppose $\varepsilon_{2}\in(0,\tfrac{1}{4})$, and
$\ell_{1},\ell_{2}\geq\ell,\ell_{3}>0$ be constants for some fixed $\ell$, and
$T=\mathsf{flatten}(\hat{f}_{\ell_{1}+\ell_{2}+\ell_{3}},\ell_{1},\ell_{2},\ell_{3})$
have decomposition $T=\sum_{i=1}^{m}\lambda_{i}(u_{i}\otimes v_{i}\otimes
z_{i})$ with $\lambda_{i}\in\mathbb{R}$ and unit vectors
$u_{i}=w_{i}^{\otimes\ell_{1}}\in\mathbb{R}^{d^{\ell_{1}}},v_{i}=w_{i}^{\otimes\ell_{2}}\in\mathbb{R}^{d^{\ell_{2}}},z_{i}=w_{i}^{\otimes\ell_{3}}\in\mathbb{R}^{d^{\ell_{3}}}$.
There exists
$\eta_{1}=\mathrm{poly}(\varepsilon_{2},s_{m}(W^{\odot\ell}))/\mathrm{poly}(m,d,B)>0$
and
$\varepsilon_{1}\coloneqq\max\\{2\varepsilon_{2},1/\mathrm{poly}(1/\varepsilon_{2},1/s_{m}(W^{\odot\ell}),d^{\ell_{1}+\ell_{2}+\ell_{3}},B)\\}$
such that if

$\displaystyle\|T-\widetilde{T}\|_{F}$
$\displaystyle\leq\eta^{\prime}_{1}\coloneqq\min\Big{\\{}\mathrm{poly}(\varepsilon_{2},s_{m}(W^{\odot\ell}))/\mathrm{poly}(m,d,B,1/\eta_{1}),\frac{\eta_{1}}{2}\Big{\\}},$
(69)

then Jennrich’s algorithm runs with rank
$k^{\prime}\coloneqq\operatorname*{argmax}_{r\leq
m}s_{r}(\text{flatten}(\widetilde{T},\ell_{1},\ell_{2}+\ell_{3},0))>\eta_{1}$
and w.h.p. outputs888 Note that one can also choose to pad the output with
zeros to output $m$ sets of parameters instead of $k^{\prime}$ if required.
$\\{\widetilde{\lambda}_{i},\widetilde{w}_{i}\\}_{i\in[k^{\prime}]}$ such that
there exists a permutation $\pi:[m]\to[m]$ and signs
$\xi_{i}\in\\{1,-1\\}~{}\forall i\in[m]$ satisfying:

$\displaystyle(i)~{}~{}\forall i\in[m],$
$\displaystyle~{}~{}|\lambda_{i}-\widetilde{\lambda}_{\pi(i)}|\leq\varepsilon_{2}^{2},\text{
and }$ (70) $\displaystyle(ii)~{}~{}\forall i\in[m],$ $\displaystyle\text{
s.t. }|\lambda_{i}|>\varepsilon_{1},~{}\text{we have}~{}\lVert w_{i}^{\otimes
t}-\xi_{\pi(i)}^{t}\widetilde{w}_{\pi(i)}^{\otimes
t}\rVert_{2}\leq\varepsilon_{2},~{}~{}~{}\forall t\in[2\ell].$ (71)

Before we proceed to the proof of this lemma, we first state and prove a
couple of simple claims. We use the following simple claim about the
assumptions of the theorem implying lower bounds on the least singular value
of the submatrices given by two columns of $W$.

######  Claim B.2.

Suppose the matrix $M_{j}\in\mathbb{R}^{d^{j}\times 2}$ formed by columns
$u^{\otimes j}$ and $v^{\otimes j}$ for $u,v\in\mathbb{S}^{d-1}$. Suppose
$s_{2}(M_{\ell})\geq\kappa$, then $s_{2}(M_{1})\geq\kappa/\sqrt{2\ell}$.

###### Proof.

Suppose $v=\alpha u+\sqrt{1-\alpha^{2}}u^{\perp}$ for some
$u^{\perp}\in\mathbb{S}^{d-1}$ that is perpendicular to $u$. It is easy to see
that

$\langle v^{\otimes\ell},u^{\otimes\ell}\rangle=\alpha^{\ell}.$

For two unit vectors $u,v$, the least singular value of the matrix given by
them as columns is

$\displaystyle\min_{\begin{subarray}{c}x,y\in\mathbb{R}\\\
x^{2}+y^{2}=1\end{subarray}}\lVert xu+yv\rVert_{2}$
$\displaystyle=\min_{\begin{subarray}{c}x,y\in\mathbb{R}\\\
x^{2}+y^{2}=1\end{subarray}}\sqrt{x^{2}+y^{2}+2xy\langle
u,v\rangle}=\min_{\begin{subarray}{c}x,y\in\mathbb{R}\\\
x^{2}+y^{2}=1\end{subarray}}\sqrt{1+2xy\langle u,v\rangle}=\sqrt{1-|\langle
u,v\rangle|}.$ $\displaystyle\text{Hence, }\kappa^{2}$
$\displaystyle=1-\alpha^{\ell}~{}~{}~{}\implies
s_{2}\big{(}[u\;v]\big{)}=\sqrt{1-\alpha}=\sqrt{1-(1-\kappa^{2})^{1/\ell}}\geq\frac{\kappa}{\sqrt{2\ell}}.$

∎

We use the following simple claim shows that if we obtain a rank-1 term which
is close, then the corresponding vectors are also close.

######  Claim B.3.

For any $\varepsilon>0,\ell\in\mathbb{N}$ with $\ell\geq 2$, suppose
$\alpha,\beta>0$, and $u,v\in\mathbb{S}^{d-1}$ satisfy $\lVert\alpha
u^{\otimes\ell}-\beta v^{\otimes\ell}\rVert_{F}\leq\varepsilon$, for some