LEARNING TWO-LAYER NEURAL NETWORKS WITH SYMMETRIC INPUTS
Year: 2019

Abstract:
We give a new algorithm for learning a two-layer neural network under a general class of input distributions. Assuming there is a ground-truth two-layer network

text: y = Aσ(Wx) + ξ,
whereA,W are weight matrices, ξ represents noise, and the number of neurons in the hidden layer is no larger than the input or output, our algorithm is guaranteed to recover the parameters A,W of the ground-truth network. The only requirement on the input x is that it is symmetric, which still allows highly complicated and structured input. Our algorithm is based on the method-of-moments framework and extends several results in tensor decompositions. We use spectral algorithms to avoid the complicated non-convex optimization in learning neural networks. Experiments show that our algorithm can robustly learn the ground-truth neural network with a small number of samples for many symmetric input distributions.

heading: 1 INTRODUCTION
text: Deep neural networks have been extremely successful in many tasks related to images, videos and reinforcement learning. However, the success of deep learning is still far from being understood in theory. In particular, learning a neural network is a complicated non-convex optimization problem, which is hard in the worst-case. The question of whether we can efficiently learn a neural network still remains generally open, even when the data is drawn from a neural network. Despite a lot of recent effort, the class of neural networks that we know how to provably learn in polynomial time is still very limited, and many results require strong assumptions on the input distribution.
In this paper we design a new algorithm that is capable of learning a two-layer1 neural network for a general class of input distributions. Following standard models for learning neural networks, we assume there is a ground truth neural network. The input data (x, y) is generated by first sampling the input x from an input distributionD, then computing y according to the ground truth network that is unknown to the learner. The learning algorithm will try to find a neural network f such that f(x) is as close to y as possible over the input distribution D. Learning a neural network is known to be a hard problem even in some simple settings (Goel et al., 2016; Brutzkus & Globerson, 2017), so we need to make assumptions on the network structure or the input distributionD, or both. Many works have worked with a simple input distribution (such as Gaussians) and try to learn more and more
1There are different ways to count the number of layers. Here by two-layer network we refer to a fullyconnected network with two layers of edges (two weight matrices). This is considered to be a three-layer network if one counts the number of layers for nodes (e.g. in Goel & Klivans (2017)) or a one-hidden layer network if one just counts the number of hidden layers.
complex networks (Tian, 2017; Brutzkus & Globerson, 2017; Li & Yuan, 2017; Soltanolkotabi, 2017; Zhong et al., 2017). However, the input distributions in real life are distributions of very complicated objects such as texts, images or videos. These inputs are highly structured, clearly not Gaussian and do not even have a simple generative model.
We consider a type of two-layer neural network, where the output y is generated as y = Aσ(Wx) + ξ. (1)
Here x ∈ Rd is the input, W ∈ Rk×d and A ∈ Rk×k are two weight matrices2. The function σ is the standard ReLU activation function σ(x) = max{x, 0} applied entry-wise to the vector Wx, and ξ is a noise vector that has E[ξ] = 0 and is independent of x. Although the network only has two layers, learning similar networks is far from trivial: even when the input distribution is Gaussian, Ge et al. (2017b) and Safran & Shamir (2018) showed that standard optimization objective can have bad local optimal solutions. Ge et al. (2017b) gave a new and more complicated objective function that does not have bad local minima.
For the input distributionD, our only requirement is thatD is symmetric. That is, for any x ∈ Rd, the probability of observing x ∼ D is the same as the probability of observing −x ∼ D. A symmetric distribution can still be very complicated and cannot be represented by a finite number of parameters. In practice, one can often think of the symmetry requirement as a “factor-2” approximation to an arbitrary input distribution: if we have arbitrary training samples, it is possible to augment the input data with their negations to make the input distribution symmetric, and it should take at most twice the effort in labeling both the original and augmented data. In many cases (such as images) the augmented data can be interpreted (for images it will just be negated colors) so reasonable labels can be obtained.

heading: 1.1 OUR RESULTS
text: When the input distribution is symmetric, we give the first algorithm that can learn a two-layer neural network. Our algorithm is based on the method-of-moments approach: first estimate some correlations between x and y, then use these information to recover the model parameters. More precisely we have Theorem 1 (informal). If the data is generated according to Equation (1), and the input distribution x ∼ D is symmetric. Given exact correlations between x, y of order at most 4, as long as A,W and input distribution are not degenerate, there is an algorithm that runs in poly(d) time and outputs a network Â, Ŵ of the same size that is effectively the same as the ground-truth network: for any input x, Âσ(Ŵx) = Aσ(Wx).
Of course, in practice we only have samples of (x, y) and cannot get the exact correlations. However, our algorithm is robust to perturbations, and in particular can work with polynomially many samples. Theorem 2 (informal). If the data is generated according to Equation (1), and the input distribution x ∼ D is symmetric. As long as the weight matricesA,W and input distributions are not degenerate, there is an algorithm that uses poly(d, 1/ ) time and number of samples and outputs a network Â, Ŵ of the same size that computes an -approximation function to the ground-truth network: for any input x, ‖Âσ(Ŵx)−Aσ(Wx)‖2 ≤ .
In fact, the algorithm recovers the original parameters A,W up to scaling and permutations. Here when we say weight matrices are not degenerate, we mean that the matrices A,W should be full rank, and in addition a certain distinguishing matrix that we define later in Section 2 is also full rank. We justify these assumptions using the smoothed analysis framework (Spielman & Teng, 2004).
In smoothed analysis, the input is not purely controlled by an adversary. Instead, the adversary can first generate an arbitrary instance (in our case, arbitrary weight matrices W,A and symmetric input distribution D), and the parameters for this instance will be randomly perturbed to yield a perturbed instance. The algorithm only needs to work with high probability on the perturbed instance. This limits the power of the adversary and prevents it from creating highly degenerate cases (e.g. choosing the weight matrices to be much lower rank than k). Roughly speaking, we show
2Here we assume A ∈ Rk×k for simplicity, our results can easily be generalized as long as the dimension of output is no smaller than the number of hidden units.
Theorem 3 (informal). There is a simple way to perturb the input distribution, W and A such that with high probability, the distance between the perturbed instance and original instance is at most λ, and our algorithm outputs an -approximation to the perturbed network with poly(d, 1/λ, 1/ ) time and number of samples.
In the rest of the paper, we will first review related works. Then in Section 2 we formally define the network and introduce some notations. Our algorithm is given in Section 3. Finally in Section 4 we run experiments to show that the algorithm can indeed learn the two-layer network efficiently and robustly. The experiments show that our algorithm works robustly with reasonable number of samples for different (symmetric) input distributions and weight matrices. Due to space constraints, the proof for polynomial number of samples (Theorem 2) and smoothed analysis (Theorem 3) are deferred to the appendix.

heading: 1.2 RELATED WORK
text: There are many works in learning neural networks, and they come in many different styles.
Non-standard Networks Some works focus on networks that do not use standard activation functions. Arora et al. (2014) gave an algorithm that learns a network with discrete variables. Livni et al. (2014) and follow-up works learn neural networks with polynomial activation functions. Oymak & Soltanolkotabi (2018) used the rank-1 tensor decomposition for learning a non-overlapping convolutional neural network with differentiable and smooth activation and Gaussian input.
ReLU network, Gaussian input When the input is Gaussian, Ge et al. (2017b) showed that for a two-layer neural network, although the standard objective does have bad local optimal solutions, one can construct a new objective whose local optima are all globally optimal. Several other works (Tian, 2017; Du et al., 2017b; Brutzkus & Globerson, 2017; Li & Yuan, 2017; Soltanolkotabi, 2017; Zhong et al., 2017; Zhang et al., 2018) extend this to different settings.
General input with score functions A closely related work (Janzamin et al., 2015) does not require the input distribution to be Gaussian, but still relies on knowing the score function of the input distribution (which in general cannot be estimated efficiently from samples). Recently, Gao et al. (2018) gave a way to design loss functions with desired properties for one-hidden-layer neural networks with general input distributions based on a new proposed local likelihood score function estimator. For general distributions (including symmetric ones) their estimator can still require number of samples that is exponential in dimension d (as in Assumption 1(d)).
General input distributions There are several lines of work that try to extend the learning results to more general distributions. Du et al. (2017a) showed how to learn a single neuron or a single convolutional filter under some conditions for the input distribution. Daniely et al. (2016); Zhang et al. (2016; 2017); Goel & Klivans (2017); Du & Goel (2018) used kernel methods to learn neural networks when the norm of the weights and input distributions are both bounded (and in general the running time and sample complexity in this line of work depend exponentially on the norms of weights/input). Recently, Du et al. (2018) showed that gradient descent minimizes the training error in an over-parameterized two-layer neural network. They only consider training error while our results also apply to testing error. The work that is most similar to our setting is Goel et al. (2018), where they showed how to learn a single neuron (or a single convolutional filter) for any symmetric input distribution. Our two-layer neural network model is much more complicated.
Method-of-Moments and Tensor Decomposition Our work uses method-of-moments, which has already been applied to learn many latent variable models (see Anandkumar et al. (2014) and references there). The particular algorithm that we use is inspired by an over-complete tensor decomposition algorithm FOOBI (De Lathauwer et al., 2007). Our smoothed analysis results are inspired by Bhaskara et al. (2014) and Ma et al. (2016), although our setting is more complicated and we need several new ideas.

heading: 2 PRELIMINARIES
text: In this section, we first describe the neural network model that we learn, and then introduce notations related to matrices and tensors. Finally we will define distinguishing matrix, which is a central object in our analysis.

heading: 2.1 NETWORK MODEL
text: We consider two-layer neural networks with d-dimensional input, k hidden units and k-dimensional output, as shown in Figure 1. We assume that k ≤ d. The input of the neural network is denoted by x ∈ Rd. Assume that the input x is i.i.d. drawn from a symmetric distribution D3. Let the two weight matrices in the neural network be W ∈ Rk×d and A ∈ Rk×k. The output y ∈ Rk is generated as follows: y = Aσ(Wx) + ξ, (2) where σ(·) is the element-wise ReLU function and ξ ∈ Rk is zero-mean random noise, which is independent with input x. Let the value of hidden units be h ∈ Rk, which is equal to σ(Wx). Denote i-th row of matrix W as w>i (i = 1, 2, ..., k). Also, let i-th column of matrix A be ai (i = 1, 2, ..., k). By property of ReLU activations, for any constant c > 0, scaling the i-th row of W by c while scaling the i-th column of A by 1/c does not change the function computed by the network. Therefore without loss of generality, we assume every row vector of W has unit norm.

heading: 2.2 NOTATIONS
text: We use [n] to denote the set {1, 2, · · · , n}. For two random variables X and Y , we say X d= Y if they come from the same distribution.
In the vector space Rn, we use 〈·, ·〉 to denote the inner product of two vectors, and use ‖·‖ to denote the Euclidean norm. We use ei to denote the i-th standard basis vector. For a matrix A ∈ Rm×n, let A[i,:] denote its i-th row vector, and let A[:,j] denote its j-th column vector. Let A’s singular values be σ1(A) ≥ σ2(A) ≥ · · · ≥ σmin(m,n)(A), and denote the smallest singular value be σmin(A) = σmin(m,n)(A). The condition number of matrix A is defined as κ(A) := σ1(A)/σmin(A). We use In to denote the identity matrix with dimension n× n. The spectral norm of a matrix is denoted as ‖ · ‖, and the Frobenius norm as ‖ · ‖F . We represent a d-dimensional linear subspace S by a matrix S ∈ Rn×d, whose columns form an orthonormal basis for subspace S. The projection matrix onto the subspace S is denoted by ProjS = SS
>, and the projection matrix onto the orthogonal subspace of S is denoted by ProjS⊥ = In − SS>. For matrix A ∈ Rm1×n1 , C ∈ Rm2×n2 , let the Kronecker product of A and C be A ⊗ C ∈ Rm1m2×n1n2 , which is defined as (A ⊗ C)(i1,i2),(j2,j2) = Ai1,i2Cj1,j2 . For a vector x ∈ Rd, the Kronecker product x ⊗ x has dimension d2. We denote the p-fold Kronecker product of x as x⊗p, which has dimension dp.
We often need to convert between vectors and matrices. For a matrixA ∈ Rm×n, let vec(A) ∈ Rmn be the vector obtained by stacking all the columns of A. For a vector a ∈ Rm2 , let mat(x) ∈ Rm×m denote the inverse mapping such that vec(mat(a)) = a. Let Rk×ksym be the space of all k × k
3Suppose the density function of distribution D is p(·), we assume p(x) = p(−x) for any x ∈ Rd
symmetric matrices, which has dimension ( k 2 ) + k. For convenience, we denote k2 = ( k 2 ) . For a symmetric matrix B ∈ Rk×ksym, we denote vec∗(B) ∈ Rk2+k as the vector obtained by stacking all the upper triangular entries (including diagonal entries) of B. Note that vec(B) still has dimension k2. For a vector b ∈ Rk2+k, let mat∗(b) ∈ Rk×ksym denote the inverse mapping of vec∗(·) such that vec∗(mat∗(B)) = b.

heading: 2.3 DISTINGUISHING MATRIX
text: A central object in our analysis is a large matrix whose columns are closely related to pairs of hidden variables. We call this the distinguishing matrix and define it below: Definition 1. Given a weight matrix W of the first layer, and the input distribution D, the distinguishing matrix ND ∈ Rd2×k2 is a matrix whose columns are indexed by ij where 1 ≤ i < j ≤ k, and
NDij = Ex∼D [ (w>i x)(w > j x)(x⊗ x)1{w>i xw>j x ≤ 0} ] .
Another related concept is the augmented distinguishing matrix M , which is a d2× (k2 + 1) matrix whose first k2 columns are exactly the same as distinguishing matrixN , and the last column (indexed by 0) is defined as
MD0 = Ex∼D [ x⊗ x ] .
For both matrices, when the input distribution is clear from context we use N or M and omit the superscript.
The exact reason for these definitions will only be clear after we explain the algorithm in Section 3. Our algorithm will require that these matrices are robustly full rank, in the sense that σmin(M) is lowerbounded. Intuitively, every column NDij looks at the expectation over samples that have opposite signs for weights wi, wj (w>i xw > j x ≤ 0, hence the name distinguishing matrix).
Requiring M and N to be full rank prevents several degenerate cases. For example, if two hidden units are perfectly correlated and always share the same sign for every input, this is very unnatural and requiring the distinguishing matrix to be full rank prevents such cases. Later in Section C we will also show that requiring a lowerbound on σmin(M) is not unreasonable: in the smoothed analysis setting where the nature can make a small perturbation on the input distribution D, we show that for any input distributionD, there exists simple perturbationsD′ that are arbitrarily close toD such that σmin(M D′) is lowerbounded.

heading: 3 OUR ALGORITHM
text: In this section, we describe our algorithm for learning the two-layer networks defined in Section 2.1. As a warm-up, we will first consider a single-layer neural network and recover the results in Goel et al. (2018) using method-of-moments. This will also be used as a crucial step in our algorithm. Due to space constraints we will only introduce algorithm and proof ideas, the detailed proof is deferred to Section A in appendix. Throughout this section, when we use E[·] without further specification the expectation is over the randomness x ∼ D and the noise ξ.

heading: 3.1 WARM-UP: LEARNING SINGLE-LAYER NETWORKS
text: We will first give a simple algorithm for learning a single-layer neural network. More precisely, suppose we are given samples (x1, y1), ..., (xn, yn) where xi ∼ D comes from a symmetric distribution, and the output yi is computed by
yi = σ(w >xi) + ξi. (3)
Here ξi’s are i.i.d. noises that satisfy E[ξi] = 0. Noise ξi is also assumed to be independent with input xi. The goal is to learn the weight vector w.
The idea of the algorithm is simple: we will estimate the correlations between x and y and the covariance of x, and then recover the hidden vectorw using these two estimates. The main challenge here is that y is not a linear function on x. Goel et al. (2018) gave a crucial observation that allows us to deal with the non-linearity:
Algorithm 1 Learning Single-layer Neural Networks Input: Samples (x1, y1), ..., (xn, yn) generated according to Equation (3). Output: Estimate of weight vector w. 1: Estimate v = 1n ∑n i=1 yixi.
2: Estimate C = 1n ∑n i=1 xix > i 3: return 2C−1v.
Lemma 1. Suppose x ∼ D comes from a symmetric distribution and y is computed as in (3), then
E[yx] = 1
2 E[xx>]w.
Importantly, the right hand side of Lemma 1 does not contain the ReLU function σ. This is true because if x comes from a symmetric distribution, averaging between x and −x can get rid of nonlinearities like ReLU or leaky-ReLU. Later we will prove a more general version of this lemma (Lemma 6).
Using this lemma, it is immediate to get a method-of-moments algorithm for learning w: we just need to estimate E[yx] and E[xx>], then we know w = 2(E[xx>])−1E[yx]. This is summarized in Algorithm 1.

heading: 3.2 LEARNING TWO-LAYER NETWORKS
text: In order to learn the weights of the network defined in Section 2.1, a crucial observation is that we have k outputs as well as k hidden-units. This gives a possible way to reduce the two-layer problem to the single-layer problem. For simplicity, we will consider the noiseless case in this section, where
y = Aσ(Wx). (4)
Let u ∈ Rk be a vector and consider u>y, it is clear that u>y = (u>A)σ(Wx). Let zi be the normalized version i-th row of A−1, then we know zi has the property that z>i A = λie > i where λi > 0 is a constant and ei is a basis vector.
The key observation here is that if u = zi, then u>A = λie>i . As a result, u >y = λie > i σ(Wx) = σ(λiw > i x) is the output of a single-layer neural network with weight equal to λiwi. If we know all the vectors {z1, ..., zk}, the input/output pairs (x, z>i y) correspond to single-layer networks with weight vectors {λiwi}. We can then apply the algorithm in Section 3.1 (or the algorithm in Goel et al. (2018)) to learn the weight vectors.
When u>A = λiei, we say that u>y is a pure neuron. Next we will design an algorithm that can find all vectors {zi}’s that generate pure neurons, and therefore reduce the problem of learning a two-layer network to learning a single-layer network.
Pure Neuron Detector In order to find the vector u that generates a pure neuron, we will try to find some property that is true if and only if the output can be represented by a single neuron.
Intuitively, using ideas similar to Lemma 1 we can get a property that holds for all pure neurons:
Lemma 2. Suppose ŷ = σ(w>x), then E[ŷ2] = 12w >E[xx>]w, and E[ŷx] = 12E[xx >]w. As a result we have
E[ŷ2] = 2E[ŷx>]E[xx>]−1E[ŷx].
As before, the ReLU activation does not appear because of the symmetric input distribution. For ŷ = u>y, we can estimate all of these moments (E[ŷ2],E[ŷx>],E[xx>]) using samples and check whether this condition is satisfied. However, the problem with this property is that even if z = u>y is not pure, it may still satisfy the property. More precisely, if ŷ = ∑k i=1 ciσ(w > i x), then we have
2E[ŷx>]E[xx>]−1E[ŷx]− E[ŷ2] = ∑
1≤i<j≤k
cicjE [ (w>i x)(w > j x)1{w>i xw>j x ≤ 0} ]
The additional terms may accidentally cancel each other which leads to a false positive. To address this problem, we consider a higher order moment:
Lemma 3. Suppose ŷ = σ(w>x), then E[ŷ2(x⊗ x)] = 2E [ ŷ · (E[ŷx>]E[xx>]−1x) · (x⊗ x) ] .
Moreover, if ŷ = ∑k i=1 ciσ(w > i x) where c ∈ Rk is a k-dimensional vector, we have
2E [ ŷ · (E[ŷx>]E[xx>]−1x) · (x⊗ x) ] − E[ŷ2(x⊗ x)] = ∑ 1≤i<j≤k cicjNij .
Here Nij’s are columns of the distinguishing matrix defined in Definition 1.
The important observation here is that there are k2 = ( k 2 ) extra terms in
2E [ ŷ · (E[ŷx>]E[xx>]−1x) · (x⊗ x) ] − E[ŷ2(x ⊗ x)] that are multiples of Nij , which are
d2 (or ( d+1
2
) considering their symmetry) dimensional objects. When the distinguishing matrix is
full rank, we know its columns Nij are linearly independent. In that case, if the sum of the extra terms is 0, then the coefficient in front of each Nij must also be 0. The coefficients are cicj which will be non-zero if and only if both ci, cj are non-zero, therefore to make all the coefficients 0 at most one of {ci} can be non-zero. This is summarized in the following Corollary: Corollary 1 (Pure Neuron Detector). Define f(u) := 2E [ (u>y) · (E[(u>y)x>]E[xx>]−1x) · (x⊗ x)
] − E[(u>y)2(x ⊗ x)]. Suppose the distinguishing matrix is full rank, if f(u) = 0 for unit vector u, then u must be equal to one of ±zi.
We will call the function f(u) a pure neuron detector, as u>y is a pure neuron if and only if f(u) = 0. Therefore, to finish the algorithm we just need to find all solutions for f(u) = 0.
Linearization The main obstacle in solving the system of equations f(u) = 0 is that every entry of f(u) is a quadratic function in u. The system of equations f(u) = 0 is therefore a system of quadratic equations. Solving a generic system of quadratic equations is NP-hard. However, in our case this can be solved by a technique that is very similar to the FOOBI algorithm for tensor decomposition (De Lathauwer et al., 2007). The key idea is to linearize the function by thinking of each degree 2 monomial uiuj as a separate variable. Now the number of variables is k2+k = ( k 2 ) +k and f is linear in this space. In other words, there exists a matrix T ∈ Rd2×(k2+k) such that Tvec∗(uu>) = f(u). Clearly, if u>y is a pure neuron, then Tvec∗(uu>) = f(u) = 0. That is, {vec∗(ziz>i )} are all in the nullspace of T . Later in Section A we will prove that the nullspace of T consists of exactly these vectors (and their combinations):
Lemma 4. Let T be the unique Rd2×(k2+k) matrix that satisfies T vec∗(uu>) = f(u) (where f(u) is defined as in Corollary 1), suppose the distinguishing matrix is full rank, then the nullspace of T is exactly the span of {vec∗(ziz>i )}.
Based on Lemma 4, we can just estimate the tensor T from the samples we are given, and its smallest singular directions would give us the span of {vec∗(ziz>i )}.
Finding zi’s from span of ziz>i ’s In order to reduce the problem to a single-layer problem, the final step is to find zi’s from span of ziz>i ’s. This is also a step that has appeared in FOOBI and more generally other tensor decomposition algorithms, and can be solved by a simultaneous diagonalization. Let Z be the matrix whose rows are zi’s, which means Z = diag(λ)A−1. Let X = Z>DXZ and Y = Z>DY Z be two random elements in the span of ziz>i , whereDX andDY are two random diagonal matrices. Both matrices X and Y can be diagonalized by matrix Z. In this case, if we compute XY −1 = Z>DXD−1Y (Z >)−1, since zi is a column of Z>, we know
XY −1zi = Z >DXD −1 Y (Z >)−1zi = Z >DXD −1 Y ei =
DX(i, i) DY (i, i) Z>ei = DX(i, i) DY (i, i) zi.
That is, zi is an eigenvector of XY −1! The matrix XY −1 can have at most k eigenvectors and there are k zi’s, therefore the zi’s are the only eigenvectors of XY −1.
Lemma 5. Given the span of ziz>i ’s, letX,Y be two random matrices in this span, with probability 1 the zi’s are the only eigenvectors of XY −1.
Using this procedure we can find all the zi’s (up to permutations and sign flip). Without loss of generality we assume z>i A = λie > i . The only remaining problem is that λi might be negative. However, this is easily fixable by checking E[z>i y] = λiE[σ(w>i x)]. Since σ(w>i x) is always nonnegative, E[z>i y] has the same sign as λi, and we can flip zi if E[z>i y] is negative.

heading: 3.3 DETAILED ALGORITHM AND GUARANTEES
text: We can now give the full algorithm, see Algorithm 2. The main steps of this algorithm is as explained in the previous section. Steps 2 - 5 constructs the pure neuron detector and finds the span of vec∗(ziz>i ) (as in Corollary 1); Steps 7 - 9 performs simultaneous diagonalization to get all the zi’s; Steps 11, 12 calls Algorithm 1 to solve the single-layer problem and outputs the correct result.
Algorithm 2 Learning Two-layer Neural Networks Input: Samples (x1, y1), ..., (xn, yn) generated according to Equation (4) Output: Weight matrices W and A. 1: {Finding span of vec(ziz>i )’s} 2: Estimate empirical moments Ê[xx>], Ê[yx>], Ê[y ⊗ x⊗3] and Ê[y ⊗ y ⊗ (x⊗ x)]. 3: Compute the vector f(u) = 2Ê [ (u>y)·(Ê[(u>y)x>]Ê[xx>]−1x)·(x⊗x) ] −Ê [ (u>y)2(x⊗x)
] where each entry is expressed as a degree-2 polynomial over u. 4: Construct matrix T ∈ Rd2×(k2+k) such that, Tvec∗(uu>) = f(u). 5: Compute the k-least right singular vectors of T , denoted by vec∗(U1), vec∗(U2), · · · , vec∗(Uk).
Let S be a k2 + k by k matrix, where the i-th column of S is vector vec∗(Ui). {Here span S is equal to span of {vec∗(ziz>i )}.} 6: {Finding zi’s from span} 7: Let X = mat∗(Sζ1), Y = mat∗(Sζ2), where ζ1 and ζ2 are two independently sampled k-
dimensional standard Gaussian vectors. 8: Let z1, ..., zk be eigenvectors of XY −1. 9: For each zi, if Ê[z>i y] < 0 let zi ← −zi. {Here zi’s are normalized rows of A−1.}
10: {Reduce to 1-Layer Problem} 11: For each zi, let vi be the output of Algorithm 1 with input (x1, z>i y1), ..., (xn, z > i yn). 12: Let Z be the matrix whose rows are z>i ’s, V be the matrix whose rows are v > i ’s. 13: return V , Z−1.
We are now ready to state a formal version of Theorem 1:
Theorem 4. SupposeA,W,E[xx>] and the distinguishing matrixN are all full rank, and Algorithm 2 has access to the exact moments, then the network returned by the algorithm computes exactly the same function as the original neural network.
It is easy to prove this theorem using the lemmas we have.
Proof. By Corollary 1, we know that after Step 5 of Algorithm 2, the span of columns of S is exactly equal to the span of {vec∗(ziz>i )}. By Lemma 5, we know the eigenvectors of XY −1 at Step 8 are exactly the normalized version of rows of A−1. Without loss of generality, we will fix the permutation and assume z>i A = λie > i . In Step 9, we use the fact that E[z>i y] = λiE[σ(w>i x)] where E[σ(w>i x)] is always positive because σ is the ReLU function. Therefore, after Step 9 we can assume all the λi’s are positive.
Now the output z>i y = λiσ(w > i x) = σ(λiw > i x) (again by property of ReLU function σ), by the design of Algorithm 1 we know vi = λiwi. We also know that Z = diag(λ)A−1, therefore Z−1 = Adiag(λ)−1. Notice that Z−1σ(V x) = Adiag(λ)−1σ(diag(λ)Wx) = Aσ(Wx). These two scaling factors cancel each other, so the two networks compute the same function.

heading: 4 EXPERIMENTS
text: In this section, we provide experimental results to validate the robustness of our algorithm for both Gaussian input distributions as well as more general symmetric distributions such as symmetric mixtures of Gaussians.
There are two important ways in which our implementation differs from our description in Section 3.3. First, our description of the simultaneous diagonalization step in our algorithm is mostly for simplicity of both stating and proving the algorithm. In practice we find it is more robust to draw 10k random samples from the subspace spanned by the last k right-singular vectors of T and compute the CP decomposition of all the samples (reshaped as matrices and stacked together as a tensor) via alternating least squares (Comon et al., 2009). As alternating least squares can also be unstable we repeat this step 10 times and select the best one. Second, once we have recovered and fixed A we use gradient descent to learn W , which compared to Algorithm 1 does a better job of ensuring the overall error will not explode even if there is significant error in recovering A. Crucially, these modifications are not necessary when the number of samples is large enough. For example, given 10,000 input samples drawn from a spherical Gaussian and A and W drawn as random 10 × 10 orthogonal matrices, our implementation of the original formulation of the algorithm was still able to recover both A and W with an average error of approximately 0.15 and achieve close to zero mean square error across 10 random trials.

heading: 4.1 SAMPLE EFFICIENCY
text: First we show that our algorithm does not require a large number of samples when the matrices are not degenerate. In particular, we generate random orthonormal matrices A and W as the ground truth, and use our algorithm to learn the neural network. As illustrated by Figure 2, regardless of the size of W and A our algorithm is able to recover both weight matrices with minimal error so long as the number of samples is a few times of the number of parameters. To measure the error in recovering A and W , we first normalize the columns of A and rows of W for both our learned parameters and the ground truth, pair corresponding columns and rows together, and then compute the squared distance between learned and ground truth parameters. Note in the rightmost plot of Figure 2, in order to compare the performance between different dimensions, we further normalize the recovering error by the dimension of W and A. It shows that the squared root of normalized error remains stable as the dimension of A and W grows from 10 to 32. In Figure 2, we also show the overall mean square error–averaged over all output units–achieved by our learned parameters.

heading: 4.2 ROBUSTNESS TO NOISE
text: Figure 3 demonstrates the robustness of our algorithm to label noise ξ for Gaussian and symmetric mixture of Gaussians input distributions. In this experiment, we fix the size of both A and W to be 10× 10 and again generate both parameters as random orthonormal matrices. The overall mean square error achieved by our algorithm grows almost perfectly in step with the amount of label noise,
indicating that our algorithm recovers the globally optimal solution regardless of the choice of input distribution.

heading: 4.3 ROBUSTNESS TO CONDITION NUMBER
text: We’ve already shown that our algorithm continues to perform well across a range of input distributions and even when A and W are high-dimensional. In all previous experiments however, we sampled A and W as random orthonormal matrices so as to control for their conditioning. In this experiment, we take the input distribution to be a random symmetric mixture of two Gaussians and vary the condition number of either A or W by sampling singular value decompositions UΣV > such that U and V are random orthonormal matrices and Σii = λ−i, where λ is chosen based on the desired condition number. Figure 4 respectively demonstrate that the performance of our algorithm remains steady so long as A and W are reasonably well-conditioned before eventually fluctuating. Moreover, even with these fluctuations the algorithm still recoversA andW with sufficient accuracy to keep the overall mean square error low.

heading: 5 CONCLUSION
text: Optimizing the parameters of a neural network is a difficult problem, especially since the objective function depends on the input distribution which is often unknown and can be very complicated. In this paper, we design a new algorithm using method-of-moments and spectral techniques to avoid the
complicated non-convex optimization for neural networks. Our algorithm can learn a network that is of similar complexity as the previous works, while allowing much more general input distributions.
There are still many open problems. The current result requires output to have the same (or higher) dimension than the hidden layer, and the hidden layer does not have a bias term. Removing these constraints are are immediate directions for future work. Besides the obvious ones of extending our results to more general distributions and more complicated networks, we are also interested in the relations to optimization landscape for neural networks. In particular, our algorithm shows there is a way to find the global optimal network in polynomial time, does that imply anything about the optimization landscape of the standard objective functions for learning such a neural network, or does it imply there exists an alternative objective function that does not have any local minima? We hope this work can lead to new insights for optimizing a neural network.

heading: ACKNOWLEDGEMENT
text: This work was supported by NSF CCF-1704656.

heading: A DETAILS OF EXACT ANALYSIS
text: In this section, we first provide the missing proofs for the lemmas appeared in Section 3. Then we discuss how to handle the noise case (i.e. y = σ(Wx) + ξ) and give the corresponding algorithm (Algorithm 3). At the end we also briefly discuss how to handle the case when the matrix A has more rows than columns (more outputs than hidden units).
Again, throughout the section when we write E[·], the expectation is taken over the randomness of x ∼ D and noise ξ.

heading: A.1 MISSING PROOFS FOR SECTION 3
text: Single-layer: To get rid of the non-linearities like ReLU, we use the property of the symmetric distribution (similar to (Goel et al., 2018)). Here we provide a more general version (Lemma 6) instead of proving the specific Lemma 1. Note that Lemma 1 is the special case when a = w and p = q = 1 (here ξ does not affect the result since it has zero mean and is independent with x, thus E[yx] = E[σ(w>x)x]). Lemma 6. Suppose input x comes from a symmetric distribution, for any vector a ∈ Rd and any non-negative integers p and q satisfying that p+ q is an even number, we have
E [( σ(a>x) )p x⊗q ] = 1
2 E[(a>x)px⊗q],
where the expectation is taken over the input distribution. Proof. Since input x comes from a symmetric distribution, we know that E [( σ(a>x) )p x⊗q ] =
E [( σ(−a>x) )p (−x)⊗q ] . Thus, we have
E [( σ(a>x) )p x⊗q ] = 1
2
( E [( σ(−a>x) )p (−x)⊗q ] + E [( σ(a>x) )p x⊗q ]) .
There are two cases to consider: p and q are both even numbers or both odd numbers.
1. For the case where p and q are even numbers, we have 1
2
( E [( σ(−a>x) )p (−x)⊗q ] + E [( σ(a>x) )p x⊗q ]) = 1
2
( E [( σ(−a>x) )p x⊗q ] + E [( σ(a>x) )p x⊗q ]) = 1 2 E [(( σ(−a>x) )p + ( σ(a>x) )p) x⊗q ] .
If (a>x) ≤ 0, we know ( σ(−a>x) )p + ( σ(a>x) )p = (a>x)p + 0 = (a>x)p. Otherwise,
we have ( σ(−a>x) )p + ( σ(a>x) )p = 0 + (a>x)p = (a>x)p. Thus,
E [( σ(a>x) )p x⊗q ] = 1 2 E [(( σ(−a>x) )p + ( σ(a>x) )p) x⊗q ] = 1
2 E[(a>x)px⊗q].
2. For the other case where p and q are odd numbers, we have 1
2
( E [( σ(−a>x) )p (−x)⊗q ] + E [( σ(a>x) )p x⊗q ]) = 1 2 E [( − ( σ(−a>x) )p + ( σ(a>x) )p) x⊗q ] .
Similarly, if (a>x) ≤ 0, we know − ( σ(−a>x) )p + ( σ(a>x) )p = −(−a>x)p + 0 =
(a>x)p. Otherwise, we have − ( σ(−a>x) )p + ( σ(a>x) )p = 0 + (a>x)p = (a>x)p. Thus,
E [( σ(a>x) )p x⊗q ] = 1 2 E [( − ( σ(−a>x) )p + ( σ(a>x) )p) x⊗q ] = 1
2 E[(a>x)px⊗q].
Pure neuron detector: The first step in our algorithm is to construct a pure neuron detector based on Lemma 2 and Lemma 3. We will provide proofs for these two lemmas here.
Proof of Lemma 2. This proof easily follows from Lemma 6. Setting a = w, p = 2 and q = 0 in Lemma 6, we have E[ŷ2] = 12w >E[xx>]w. Similarly, we also know E[ŷx>] = 12w >E[xx>] and E[ŷx] = 12E[xx >]w. Thus, we have E[ŷ2] = 2E[ŷx>]E[xx>]−1E[ŷx].
Proof of Lemma 3. Here, we only prove the second equation, since the first equation is just a special case of the second equation. First, we rewrite ŷ = ∑k i=1 ciσ(w > i x) = u
>y by letting u>A = c>. Then we transform these two terms in the LHS as follows. Let’s look at E [ (u>y) ·
(E[(u>y)x>]E[xx>]−1x) · (x⊗ x) ] first. For E[(u>y)x>]E[xx>]−1, we have
E [ (u>y)x> ] E [ xx> ]−1 =E [ (u>Aσ(Wx))x> ] E [ xx>
]−1 = 1
2 u>AWE
[ xx> ] E [ xx> ]−1 = 1
2 u>AW
For any vector q ∈ Rd, consider 2E [ (u>y) · (q>x) · (x⊗ x) ] . We have
2E [ (u>y) · (q>x) · (x⊗ x) ] =2E[(u>Aσ(Wx))(q>x)(x⊗ x)] =E [ (u>AWx)(q>x)(x⊗ x)
] Let q> = E [ (u>y)x> ] E [ xx> ]−1 = 12u >AW , we have
2E [ (u>y) · (E[(u>y)x>]E[xx>]−1x) · (x⊗ x) ] = 1
2 E [ (u>AWx)2(x⊗ x) ] = 1
2 E [ 〈A>u,Wx〉2(x⊗ x) ] = 1
2 ∑ 1≤i≤k (A>u)2iE [ (w>i x) 2(x⊗ x) ] + ∑ 1≤i<j≤k (A>u)i(A >u)jE [ (w>i x)(w > j x)(x⊗ x) ] = 1
2 ∑ 1≤i≤k (A>uu>A)iiE [ (w>i x) 2(x⊗ x) ] + ∑ 1≤i<j≤k (A>uu>A)ijE [ (w>i x)(w > j x)(x⊗ x) ] .
(5) where the second equality holds due to Lemma 6.
Now, let’s look at the second term E [ (u>y)2(x⊗ x) ] .
E [ (u>y)2(x⊗ x) ] =E [ (u>Aσ(Wx))2(x⊗ x)
] =E [ 〈A>u, σ(Wx)〉2(x⊗ x)
] = ∑
1≤i≤k
(A>u)2iE [ σ(w>i x) 2(x⊗ x) ] + 2 ∑
1≤i<j≤k
(A>u)i(A >u)jE [ σ(w>i x)σ(w > j x)(x⊗ x) ] = ∑
1≤i≤k
(A>uu>A)iiE [ σ(w>i x) 2(x⊗ x) ] + 2 ∑
1≤i<j≤k
(A>uu>A)ijE [ σ(w>i x)σ(w > j x)(x⊗ x) ] = 1
2 ∑ 1≤i≤k (A>uu>A)iiE [ (w>i x) 2(x⊗ x) ] + 2 ∑ 1≤i<j≤k (A>uu>A)ijE [ σ(w>i x)σ(w > j x)(x⊗ x) ] .
(6)
Now, we subtract (5) by (6) to obtain 2E [ (u>y) · (E[(u>y)x>]E[xx>]−1x) · (x⊗ x) ] − E [ (u>y)2(x⊗ x) ] =
∑ 1≤i<j≤k (A>uu>A)ijE [ (w>i x)(w > j x)(x⊗ x) ] − 2 ∑ 1≤i<j≤k (A>uu>A)ijE [ σ(w>i x)σ(w > j x)(x⊗ x) ] =
∑ 1≤i<j≤k (A>uu>A)ijE[ [ (w>i x)(w > j x)(x⊗ x)1{w>i xw>j x ≤ 0} ] (7)
= ∑
1≤i<j≤k
(A>uu>A)ijNij , (8)
where (7) uses (9) of the following Lemma 7, and (8) uses the definition of distinguishing matrix N (Definition 1). Lemma 7. Given input x coming from a symmetric distribution, for any vector a, b ∈ Rd, we have
1 2 E [ (a>x)(b>x) ] − E [ σ(a>x)σ(b>x) ] = 1 2 E [ (a>x)(b>x)1{a>xb>x ≤ 0} ] and 1
2 E [ (a>x)(b>x)(x⊗x) ] −E [ σ(a>x)σ(b>x)(x⊗x) ] = 1 2 E [ (a>x)(b>x)(x⊗x)1{a>xb>x ≤ 0} ] , (9) where the expectation is taken over the input distribution.
Proof. Here we just prove the first identity, because the proof of the second one is almost identical. First, we rewrite 12E [ (a>x)(b>x) ] as follows
1 2 E [ (a>x)(b>x) ] = 1 2 E [ (a>x)(b>x)1{a>xb>x ≤ 0} ] + 1 2 E [ (a>x)(b>x)1{a>xb>x > 0} ] .
Thus, we only need to show that 12E [ (a>x)(b>x)1{a>xb>x > 0} ] = E [ σ(a>x)σ(b>x) ] . It’s
clear that E [ σ(a>x)σ(b>x) ] = E [ σ(a>x)σ(b>x)1{a>xb>x > 0} ] .
Since input x comes from a symmetric distribution, we have E [ σ(a>x)σ(b>x)1{a>xb>x > 0} ] =E [ σ(−a>x)σ(−b>x)1{a>xb>x > 0} ] = 1
2
( E [ σ(a>x)σ(b>x)1{a>xb>x > 0} ] + E [ σ(−a>x)σ(−b>x)1{a>xb>x > 0}
]) = 1 2 E [( σ(a>x)σ(b>x) + σ(−a>x)σ(−b>x) ) 1{a>xb>x > 0}
] = 1
2 E [ (a>x)(b>x)1{a>xb>x > 0} ] .
When a>xb>x > 0, we know a>x and b>x are both positive or both negative. In either case, we know that σ(a>x)σ(b>x) + σ(−a>x)σ(−b>x) = (a>x)(b>x). Thus, we have
E [ σ(a>x)σ(b>x) ] =E [ σ(a>x)σ(b>x)1{a>xb>x > 0} ] = 1
2 E [ (a>x)(b>x)1{a>xb>x > 0} ] ,
which finished our proof. Finding span: Now, we find the span of {vec∗(ziz>i )}. First, we recall that f(u) = 2E [ (u>y) ·
(E[(u>y)x>]E[xx>]−1x) · (x⊗ x) ] − E [ (u>y)2(x⊗ x) ] . Then, according to (8), we have
f(u) = 2E [ (u>y)·(E[(u>y)x>]E[xx>]−1x)·(x⊗x) ] −E [ (u>y)2(x⊗x) ] = ∑ 1≤i<j≤k (A>uu>A)ijNij .
It is not hard to verify that u>y is a pure neuron if and only if f(u) = 0. Note that f(u) = 0 is a system of quadratic equations. So we linearize it by increasing the dimension (i.e., consider uiuj as a single variable) similar to the FOOBI algorithm. Thus the number of variable is ( k 2 ) +k = k2 +k, i.e., ∃T ∈ Rd 2×(k2+k) such that Tvec∗(U) = f(U) = ∑
1≤i<j≤k
(A>UA)ijNij . (10)
Now, we prove the Lemma 4 which shows the null space of T is exactly the span of {vec∗(ziz>i )}.
Proof of Lemma 4. We divide the proof to the following two cases:
1. For any vector vec∗(U) belongs to the null space of T , we have Tvec∗(U) = 0. Note that the RHS of (10) equals to 0 if and only if A>UA is a diagonal matrix since the distinguishing matrix N is full column rank and A>UA is symmetric. Thus vec∗(U) belongs to the span of {vec∗(ziz>i )} since U = Z>DZ for some diagonal matrix D. 2. For any vector vec∗(U) belonging to the span of {vec∗(ziz>i )}, U is a linear combination of ziz>i ’s. Furthermore, Tvec ∗(U) is a linear combination of Tvec∗(ziz>i ). Note thatA >zi
only has one non-zero entry due to the definition of zi, for any i ∈ [k]. Thus all coefficients in the RHS of (10) are 0. We get Tvec∗(U) = 0.
Finding zi’s: Now, we prove the final Lemma 5 which finds all zi’s from the span of {vec∗(ziz>i )} by using simultaneous diagonalization. Given all zi’s, this two-layer network can be reduced to a single-layer one. Then one can use Algorithm 1 to recover the first layer parameters wi’s.
Proof of Lemma 5. As we discussed before this lemma, we have XY −1 = Z>DXD−1Y (Z >)−1. According to the following Lemma 8 (i.e., all diagonal elements ofDxD−1y are non-zero and distinct), the matrixXY −1 have k eigenvectors and there are k zi’s, therefore zi’s are the only eigenvectors of XY −1. Lemma 8. With probability 1, all diagonal elements of DX and DY are non-zero and all diagonal elements of DXD−1Y are distinct, where X = Z
>DXZ and Y = Z>DY Z are defined in Line 7 of Algorithm 2.
Proof. First, we know there exist diagonal matrices {Di : i ∈ [k]} such that mat∗(vec∗(Ui)) = Z>DiZ for all i ∈ [k], where {vec∗(Ui) : i ∈ [k]} are the k-least right singular vectors of T (see Line 5 of Algorithm 2). Then, let the vector di ∈ Rk be the diagonal elements of Di, for all i ∈ k. Let matrix Q ∈ Rk×k be a matrix where its i-th column is di. Then DX = diag(Qζ1) and DY = diag(Qζ2), where ζ1 and ζ2 are two random k-dimensional standard Gaussian vectors (see Line 7 of Algorithm 2).
Since {vec∗(Ui) : i ∈ [k]} are singular vectors of T , we know d1, d2, · · · , dk are independent. Thus, Q has full rank and none of its rows are zero vectors. Let i-th row of Q be q>i . Let’s consider DX first. In order for i-th diagonal element of DX to be zero, we need q>i ζ1 = 0. Since qi is not a zero vector, we know the solution space of qiζ1 = 0 is a lower-dimension manifold in Rk, which has zero measure. Since finite union of zero-measure sets still has measure zero, the event that zero valued elements exist in the diagonal of DX or DY happens with probability zero.
If i-th and j-th diagonal elements of DXD−1Y have same value, we have q > i ζ1(q > i ζ2) −1 = q>j ζ1(q > j ζ2)
−1. Again, we know the solution space is a lower-dimensional manifold in R2k space, with measure zero. Since finite union of zero-measure sets still has measure zero, the event that duplicated diagonal elements exist in DXD−1Y happens with probability zero.

heading: A.2 NOISY CASE
text: Now, we discuss how to handle the noisy case (i.e. y = σ(Wx) + ξ). The corresponding algorithm is described in Algorithm 3. Note that the noise ξ only affects the first two steps, i.e., pure neuron detector (Lemma 3) and finding span of vec∗(ziz>i ) (Lemma 4). It does not affect the last two steps, i.e., finding zi’s from the span (Lemma 5) and learning the reduced single-layer network. Because
Algorithm 3 Learning Two-layer Neural Networks with Noise Input: Samples (x1, y1), ..., (xn, yn) generated according to Equation (2) Output: Weight matrices W and A. 1: {Finding span of vec∗(ziz>i )} 2: Use first half of samples (i.e. {(xi, yi)}n/2i=1) to estimate empirical moments Ê[xx>], Ê[yx>],
Ê[yy>], Ê[y ⊗ x⊗3] and Ê[y ⊗ y ⊗ (x⊗ x)]. 3: Compute the vector f(u) = 2Ê [ (u>y) · (Ê[(u>y)x>]Ê[xx>]−1x) · (x⊗ x) ] − Ê [ (u>y)2(x⊗
x) ] + ( Ê [ (u>y)2 ] − 2Ê [ (u>y)x> ] Ê [ xx> ]−1Ê[(u>y)x])Ê[x ⊗ x] where each entry is expressed as a degree-2 polynomial over u. 4: Construct matrix T ∈ Rd2×(k2+k) such that, Tvec∗(uu>) = f(u). 5: Compute the k-least right singular vectors of T , denoted by vec∗(U1), vec∗(U2), · · · , vec∗(Uk).
Let S be a k2 + k by k matrix, where the i-th column of S is vector vec∗(Ui). {Here span S is equal to span of {vec∗(ziz>i )}.} 6: {Finding zi’s from span} 7: Let X = mat∗(Sζ1), Y = mat∗(Sζ2), where ζ1 and ζ2 are two independently sampled k-
dimensional standard Gaussian vectors. 8: Let z1, ..., zk be eigenvectors of XY −1. 9: For each zi, use the second half of samples {(xi, yi)}ni=n/2+1 to estimate Ê[z > i y]. If Ê[z>i y] < 0
let zi ← −zi. {Here zi’s are normalized rows of A−1.} 10: {Reduce to 1-Layer Problem} 11: For each zi, let vi be the output of Algorithm 1 with input {(xj , z>i yj)}nj=n/2+1. 12: Let Z be the matrix whose rows are z>i ’s, V be the matrix whose rows are v > i ’s. 13: return V , Z−1.
Lemma 5 is independent of the model and Lemma 1 is linear wrt. noise ξ, which has zero mean and is independent of input x.
Many of the steps in Algorithm 3 are designed with the robustness of the algorithm in mind. For example, in step 5 for the exact case we just need to compute the null space of T . However if we use the empirical moments the null space might be perturbed so that it has small singular values. The separation of the input samples into two halves is also to avoid correlations between the steps, and is not necessary if we have the exact moments.
Modification for pure neuron detector: Recall that in the noiseless case, our pure neuron detector contains a term E[(u>y)2(x ⊗ x)], which causes a noise square term in the noisy case. Here, we modify our pure neuron detector to cancel the extra noise square term. In the following lemma, we state our modified pure neuron detector in Equation 11, and give it a characterization.
Lemma 9. Suppose y = Aσ(Wx) + ξ, for any u ∈ Rk, we have f(u) := 2E [ (u>y) · (E[(u>y)x>]E[xx>]−1x) · (x⊗ x) ] − E [ (u>y)2(x⊗ x) ] + ( E [ (u>y)2 ] − 2E [ (u>y)x> ] E [ xx>
]−1E[(u>y)x])E[x⊗ x] (11) =
∑ 1≤i<j≤k (A>uu>A)ijNij − ( ∑ 1≤i<j≤k (A>uu>A)ijmij ) E[x⊗ x], (12)
where Nij’s are columns of the distinguishing matrix (Definition 1), and mij := E [ (w>i x)(w > j x)1{w>i xw>j x ≤ 0} ] .
We defer the proof of this lemma to the end of this section. Recall that the augmented distinguishing matrix M consists of the distinguishing matrix N plus column E[x ⊗ x]. Now, we need to assume the augmented distinguishing matrix M is full rank.
Modification for finding span: For Lemma 4, as we discussed above, here we assume the augmented distinguishing matrix M is full rank. The corresponding lemma is stated as follows (the proof is exactly the same as previous Lemma 4):
Lemma 10. Let T be the unique Rd2×(k2+k) matrix that satisfies T vec∗(uu>) = f(u) (defined in Equation 11), suppose the augmented distinguishing matrix is full rank, then the nullspace of T is exactly the span of {vec∗(ziz>i )}.
Similar to Theorem 4, we provide the following theorem for the noisy case. The proof is almost the same as Theorem 4 by using the noisy version lemmas (Lemmas 9 and 10).
Theorem 5. Suppose E[xx>], A,W and the augmented distinguishing matrix M are all full rank, and Algorithm 3 has access to the exact moments, then the network returned by the algorithm computes exactly the same function as the original neural network.
Now, we only need to prove Lemma 9 to finish this noise case.
Proof of Lemma 9. Similar to (5) and (6), we deduce these three terms in RHS of (11) one by one as follows. For the first term, it is exactly the same as (5) since the expectation is linear wrt. ξ. Thus, we have
2E [ (u>y) · (E[(u>y)x>]E[xx>]−1x) · (x⊗ x) ] = 1
2 ∑ 1≤i≤k (A>uu>A)iiE [ (w>i x) 2(x⊗ x) ] + ∑ 1≤i<j≤k (A>uu>A)ijE [ (w>i x)(w > j x)(x⊗ x) ] .
(13)
Now, let’s look at the second term E [ (u>y)2(x⊗ x) ] which is slightly different from (6) due to the noise ξ. Particularly, we add the third term to cancel this extra noise square term later.
E [ (u>y)2(x⊗ x) ] =E [ (u>(Aσ(Wx) + ξ))2(x⊗ x)
] =E [ (u>Aσ(Wx))2(x⊗ x) ] + E [ (u>ξ)2(x⊗ x)
] = 1
2 ∑ 1≤i≤k (A>uu>A)iiE [ (w>i x) 2(x⊗ x) ] + 2 ∑ 1≤i<j≤k (A>uu>A)ijE [ σ(w>i x)σ(w > j x)(x⊗ x) ] (14)
+ E [ (u>ξ)2(x⊗ x) ] , (15)
where (15) uses (6).
For the third term, we have
E [ (u>y)2 ] − 2E [ (u>y)x> ] E [ xx> ]−1E[(u>y)x] =E [ (u>Aσ(Wx))2 ] + E [ (u>ξ)2 ] − 1
2 E [ 〈A>u,Wx〉2 ] = ∑
1≤i≤k
(A>uu>A)iiE [ σ(w>i x) 2 ] − 1
2 ∑ 1≤i≤k (A>uu>A)iiE [ (w>i x) 2 ]
+ 2 ∑
1≤i<j≤k
(A>uu>A)ijE [ σ(w>i x)σ(w > j x) ] − ∑ 1≤i<j≤k (A>uu>A)ijE [ (w>i x)(w > j x) ] + E [ (u>ξ)2
] =−
∑ 1≤i<j≤k (A>uu>A)ijE [ (w>i x)(w > j x)1{w>i xw>j x ≤ 0} ] + E [ (u>ξ)2 ] =−
∑ 1≤i<j≤k (A>uu>A)ijmij + E [ (u>ξ)2 ] , (16)
where the third equality holds due to Lemma 6 and Lemma 7, and (16) uses the definition of mij .
Algorithm 4 Learning Two-layer Neural Networks with Non-square A Input: Samples (x1, y1), ..., (xn, yn) generated according to Equation (2). Output: Weight matrices W ∈ Rk×d and A ∈l×k. 1: Using half samples (i.e. {(xi, yi)}n/2i=1) to estimate empirical moments Ê[yx>]. 2: Let P be a l × k matrix, which columns are left singular vectors of Ê[yx>]. 3: Run Algorithm 3 on samples {(xi, P>yi)}ni=n/2. Let the output of Algorithm 3 be V,Z
−1. 4: return V , PZ−1.
Finally, we combine these three terms (13–16) as follows: f(u) = 2E [ (u>y) · (E[(u>y)x>]E[xx>]−1x) · (x⊗ x) ] − E [ (u>y)2(x⊗ x) ] + ( E [ (u>y)2 ] − 2E [ (u>y)x> ] E [ xx>
]−1E[(u>y)x])E[x⊗ x] =
∑ 1≤i<j≤k (A>uu>A)ijE [ (w>i x)(w > j x)(x⊗ x) ] − 2
∑ 1≤i<j≤k (A>uu>A)ijE [ σ(w>i x)σ(w > j x)(x⊗ x) ] − E [ (u>ξ)2 ] E [ x⊗ x ] + ( −
∑ 1≤i<j≤k (A>uu>A)ijmij + E [ (u>ξ)2 ]) E [ x⊗ x ] =
∑ 1≤i<j≤k (A>uu>A)ijNij − ( ∑ 1≤i<j≤k (A>uu>A)ijmij ) E[x⊗ x], (17)
where (17) uses (9) (same as (7)).

heading: A.3 EXTENSION TO NON-SQUARE A
text: In this paper, for simplicity, we have assumed that the dimension of output equals the number of hidden units and thus A is a k × k square matrix. Actually, our algorithm can be easily extended to the case where the dimension of output is at least the number of hidden units. In this section, we give an algorithm for this general case, by reducing it to the case where A is square. The pseudo-code is given in Algorithm 4.
Theorem 6. Suppose E[xx>],W,A and the augmented distinguishing matrix M are all full rank, and Algorithm 4 has access to the exact moments, then the network returned by the algorithm computes exactly the same function as the original neural network.
Proof. Let the ground truth parameters beA ∈ Rl×k andW ∈ Rk×d. The samples are generated by y = Aσ(Wx)+ξ, where the noise ξ is independent with input x. We have E[yx>] = 12AWE[xx
>]. Since bothW and E[xx>] are full-rank, we know the column span of E[yx>] are exactly the column span of A. Furthermore, we know the columns of P is a set of orthonormal basis for the column span of A.
For a ground truth neural network with weight matrices W and P>A, the generated sample will just be (x, P>y). According to Theorem 5, we know for any input x, we have Z−1σ(V x) = P>Aσ(Wx). Thus, we have
PZ−1σ(V x) = PP>Aσ(Wx) = Aσ(Wx),
where the second equality holds since PP> is just the projection matrix to the column span of A.

heading: B ROBUSTNESS OF MAIN ALGORITHM
text: In this section we will show that even if we do not have access to the exact moments, as long as the empirical moments are estimated with enough (polynomially many) samples, Algorithm 2 and
Algorithm 3 can still learn the parameters robustly. We will focus on Algorithm 3 as it is more general, the result for Algorithm 2 can be viewed as a corollary when the noise ξ = 0. Throughout this section, we will use V̂ , Ẑ−1 to denote the results of Algorithm 3 with empirical moments, and use V,Z−1 for the results when the algorithm has access to exact moments, similarly for other intermediate results. For the robustness of Algorithm 3, we prove the following theorem. Theorem 7. Assume that the norms of x, ξ, A are bounded by ‖x‖ ≤ Γ, ‖ξ‖ ≤ P2, ‖A‖ ≤ P1, the covariance matrix and the weight matrix are robustly full rank: σmin(E[xx>]) ≥ γ, σmin(A) ≥ β. Further assume that the augmented distinguishing matrix has smallest singular values σmin(M) ≥ α. For any small enough , for any δ < 1, given poly ( Γ, P1, P2, d, 1/ , 1/γ, 1/α, 1/β, 1/δ ) number of i.i.d. samples, let the output of Algorithm 3 be V̂ , Ẑ−1, we know with probability at least 1− δ,
‖Aσ(Wx)− Ẑ−1σ(V̂ x)‖ ≤ ,
for any input x.
In order to prove the above Theorem, we need to show that each step of Algorithm 3 is robust. We can divide Algorithm 3 into three steps: finding the span of vec∗(ziz>i )’s; finding zi’s from the span of vec∗(ziz>i )’s; recovering first layer using Algorithm 1. We will first state the key lemmas that prove every step is robust to noise, and finally combine them to show our main theorem.
First, we show that with polynomial number of samples, we can approximate the span of vec∗(ziz>i )’s in arbitrary accuracy. Let T̂ be the empirical estimate of T , which is the pure neuron detector matrix as defined in Algorithm 3. As shown in Lemma 10, the null space of T is exactly the span of vec∗(ziz>i )’s. We use standard matrix perturbation theory (see Section D.2) to show that the null space of T is robust to small perturbations. More precisely, in Lemma 11, we show that with polynomial number of samples, the span of k least singular vectors of T̂ is close to the null space of T .
Lemma 11. Under the same assumptions as in Theorem 7, let S ∈ R(k2+k)×k be the matrix whose k columns are the k least right singular vectors of T . Similarly define Ŝ ∈ R(k2+k)×k for empirical estimate T̂ . Then for any ≤ γ/2, for any δ < 1, givenO(d 3Γ14(ΓP1 √ k+P2) 6 log( dδ )
γ4 2 ) number of i.i.d. samples, we know with probability at least 1− δ,
‖SS> − ŜŜ>‖ ≤ √ 2
αβ2 .
The proof of the above lemma is in Section B.1. Basically, we need to lowerbound the spectral gap (k2-th singular value of T ) and to upperbound the Frobenius norm of T − T̂ . Standard matrix perturbation bound shows that if the perturbation is much smaller than the spectral gap, then the null space is preserved.
Next, we show that we can robustly find zi’s from the span of vec∗(ziz>i )’s. Since this step of the algorithm is the same as the simultaneous diagonalization algorithm for tensor decompositions, we use the robustness of simultaneous diagonalization (Bhaskara et al., 2014) to show that we can find zi’s robustly. The detailed proof is in Section B.2.
Lemma 12. Suppose that ‖SS> − ŜŜ>‖F ≤ , ‖A‖ ≤ P1, ‖ξ‖ ≤ P2, σmin(E[xx>]) ≥ γ, σmin(A) ≥ β. Let X̂ = mat∗(Ŝζ1), Ŷ = mat∗(Ŝζ2), where ζ1 and ζ2 are two independent standard Gaussian vectors. Let z1, · · · zk be the normalized row vectors of A−1. Let ẑ1, ..., ẑk be the eigenvectors of X̂Ŷ −1 (after sign flip). For any δ > 0 and small enough , with O( (ΓP1 √ k+P2) 2 log(d/δ) 2 ) number of i.i.d. samples in Ê[ẑ > i y], with probability at least 1− δ over the randomness of ζ1, ζ2 and i.i.d. samples, there exists a permutation π(i) ∈ [k] such that
‖ẑi − zπ[i]‖ ≤ poly(d, P1, 1/β, , 1/δ),
for any 1 ≤ i ≤ k.
Finally, given ẑi’s, the problem reduces to a one-layer problem. We will first give an analysis for Algorithm 1 as a warm-up. When we call Algorithm 1 from Algorithm 3, the situation is slightly different. Note we reserve fresh samples for this step, so that the samples used by Algorithm 1 are
still independent with the estimate ẑi (learned using the other set of samples). However, since ẑi is not equal to zi, this introduces an additional error term (ẑi − zi)>y which is not independent of x and cannot be captured by ξ. We modify the proof for Algorithm 1 to show that the algorithm is still robust as long as ‖ẑi − zi‖ is small enough. Lemma 13. Assume that ‖x‖ ≤ Γ, ‖A‖ ≤ P1, ‖ξ‖ ≤ P2 and σmin(E[xx>]) ≥ γ. Suppose that for each 1 ≤ i ≤ k, ‖ẑi − zi‖ ≤ τ . Then for any ≤ γ/2 and δ < 1, given O( (Γ2+P2Γ) 4 log( dδ )
γ4 2 )
number of samples for Algorithm 1, we know with probability at least 1− δ,
‖vi − v̂i‖ ≤ 2τ(ΓP1
√ k + P2)Γ
γ + 2 ,
for each 1 ≤ i ≤ k.
Combining the above three lemmas, we prove Theorem 7 in Section B.4.
B.1 ROBUST ANALYSIS FOR FINDING THE SPAN OF {VEC∗(ziz>j )}’S
We first prove that the step of finding the span of {vec∗(ziz>j )} is robust. The main idea is based on standard matrix perturbation bounds (see Section D.2). We first give a lowerbound on the k2-th singular value of T , giving a spectral gap between the smallest non-zero singular value and the null space. See the lemma below. The proof is given in Section B.1.1. Lemma 14. Suppose σmin(M) ≥ α, σmin(A) ≥ β, we know that matrix T has rank k2 and the k2-th singular value of T is lower bounded by αβ2.
Then we show that with enough samples the estimate T̂ is close enough to T , so Wedin’s Theorem (Lemma 25) implies the subspace found is also close to the true nullspace of T . The proof is deferred to Section B.1.2. Lemma 15. Assume that ‖x‖ ≤ Γ, ‖A‖ ≤ P1, ‖ξ‖ ≤ P2 and σmin(E[xx>]) ≥ γ > 0, then for any ≤ γ/2, for any 1 > δ > 0, given O(d 3Γ14(ΓP1 √ k+P2) 6 log( dδ )
γ4 2 ) number of i.i.d. samples, we know
‖T̂ − T‖F ≤ ,
with probability at least 1− δ.
Finally we combine the above two lemmas and show that the span of the least k right singular vectors of T̂ is close to the null space of T .
Lemma 11. Under the same assumptions as in Theorem 7, let S ∈ R(k2+k)×k be the matrix whose k columns are the k least right singular vectors of T . Similarly define Ŝ ∈ R(k2+k)×k for empirical estimate T̂ . Then for any ≤ γ/2, for any δ < 1, givenO(d 3Γ14(ΓP1 √ k+P2) 6 log( dδ )
γ4 2 ) number of i.i.d. samples, we know with probability at least 1− δ,
‖SS> − ŜŜ>‖ ≤ √ 2
αβ2 .
Proof. According to Lemma 15, given O(d 3Γ14(ΓP1
√ k+P2) 6 log( dδ )
γ4 2 ) number of i.i.d. samples, we know with probability at least 1− δ,
‖T̂ − T‖F ≤ . According to Lemma 14, we know σk2(T ) ≥ αβ2. Then, due to Lemma 27, we have
‖SS> − ŜŜ>‖ ≤ √
2‖T − T̂‖F σk2(T )
≤ √ 2
αβ2 .
B.1.1 LOWERBOUNDING k2-TH SINGULAR VALUE OF T
In order to lowerbound the k2-th singular value of T , we first express T as the product of four simpler matrices, T = MBCF , as illustrated in Figure 5. The definitions of these four matrices M ∈ Rd2×(k2+1), B ∈ R(k2+1)×k2 , C ∈ R(k2×k2), F ∈ Rk2×(k2+k) will be introduced later as we explain their effects. From Lemma 9, we know that
Tvec∗(U) = ∑
1≤i<j≤k
(A>UA)ijMij − ( ∑
1≤i<j≤k
(A>UA)ijmij ) E[x⊗ x],
for any symmetric k × k matrix U . Note that vec∗(U) is a (k2 + k)-dimensional vector. For convenience, we first use matrix F to transform vec∗(U) to vec(U), which has k2 dimensions. Matrix F is defined such that Fvec∗(U) = vec(U), for any k × k symmetric matrix U . Note that this is very easy as we just need to duplicate all the non-diagonal entries.
Second, we hope to get the coefficients (A>UA)ij’s. Notice that
vec(A>UA) = A> ⊗A>vec(U) = A> ⊗A>Fvec∗(U). Since we only care about the elements of A>UA at the ij-th position for 1 ≤ i < j ≤ k, we just pick corresponding rows of A> ⊗A> to construct our matrix C, which has dimension k2 × k2. The first matrix M is the augmented distinguishing matrix (see Definition 1). In order to better understand the reason that we need matrix B, let’s first re-write Tvec∗(U) in the following way:
Tvec∗(U) = ∑
1≤i<j≤k
(A>UA)ij ( Mij −mijE[x⊗ x] ) .
Thus, Tvec∗(U) is just a linear combination of (Mij − mijE[x ⊗ x])’s with coefficients equal to (A>UA)ij’s. We have already expressed coefficients (A>UA)ij’s using CFvec∗(U). Now, we just need to use matrix B to transform the augmented distinguishing matrix M to a d2× k2 matrix, with each column equal to (Mij −mijE[x ⊗ x]). In order to achieve this, the first k2 rows of B is just the identity matrix Ik2 , and the last row of B is [−m12,−m13, · · · ,−m1k,−m23,−m24, · · · ]>. With above characterization of T , we are ready to show that the k2-th singular value of T is lower bounded. Lemma 14. Suppose σmin(M) ≥ α, σmin(A) ≥ β, we know that matrix T has rank k2 and the k2-th singular value of T is lower bounded by αβ2.
Proof. Since matrix C has dimension k2 × k2, it’s clear that the rank of T is at most k2. We first prove that the rank of T is exactly k2.
Since the first k2 rows of B constitute the identity matrix Ik2 , we know B is a full-column rank matrix with rank equal to k2. We also know that matrix M is a full column rank matrix with rank k2 + 1. Thus, the product matrix MB is still a full-column rank matrix with rank k2. If we can prove that the product matrix CF has full-row rank equal to k2. It’s clear that T = MBCF also has rank k2. Next, we prove that CF has full-row rank.
Since σmin(A) ≥ β, we know A> ⊗ A> is full rank, and a subset of its rows C has full row rank. For the sake of contradiction, suppose that there exists non-zero vector a ∈ Rk2 , such that∑k2 l=1 al(CF )[l,:] = 0. Note that for any 1 ≤ l ≤ k2, (CF )[l,ii] = C[l,ii] for 1 ≤ i ≤ k and (CF )[l,ij] = C[l,ij] + C[l,ji] for 1 ≤ i < j ≤ k. Since C consists of a subset of rows of A> ⊗ A>, we know C[l,ij] = C[l,ji] for any l and any i < j. Thus,
∑k2 l=1 al(CF )[l,:] = 0 simply implies∑k2
l=1 alC[l,:] = 0, which breaks the fact that C is full-row rank. Thus, the assumption is false and CF has full-row rank.
Now, let’s prove that the k2-th singular value of T is lower bounded. We first show that in the product characterization of T , the smallest singular value of each individual matrix is lower bounded. According to the assumption, we know the smallest singular value of M is lower bounded by α. Since the first k2 rows of matrix B constitute a k2 × k2 identity matrix, we know
σmin(B) := min u:‖u‖≤1 ‖Bu‖ ≥ min u:‖u‖≤1 ‖Ik2×k2u‖ =: σmin(Ik2×k2) = 1,
where u is any k2-dimensional vector.
Since σmin(A) ≥ β, we know σmin(A>⊗A>) ≥ β2. According to the construction of C, we know C consists a subset of rows of A> ⊗A>. Denote the indices of the row not picked as S. We have
σmin(C) := min u:‖u‖≤1
‖u>C‖
= min v:‖v‖≤1 ∧ (vi=0,∀i∈S)
‖v>(A> ⊗A>)‖
≥ min v:‖v‖≤1
‖v>(A> ⊗A>)‖
=: σmin(A > ⊗A>)
≥β2,
where u has dimension k2 and v has dimension k2.
We lowerbound the smallest singular value of CF by showing that σmin(CF ) ≥ σmin(C). For any unit vector u ∈ Rk2 , we know [u>CF ]ii = [u>C]ii for any i and [u>CF ]ij = [u>C]ij + [u>C]ji for any i < j. We also know [u>C]ij = [u>C]ji for i < j. Thus, we know for any unit vector u, ‖u>CF‖ ≥ ‖u>C‖, which implies σmin(CF ) ≥ σmin(C). Finally, since in the beginning we have proved that matrix T has rank k2, the k2-th singular value is exactly the smallest non-zero singular value of T . Denote the smallest non-zero singular of T as σ+min(T ), we have
σk2(T ) =σ + min(T )
≥σmin(M)σmin(B)σmin(CF ) ≥αβ2,
where the first inequality holds because both M and B has full column rank.

heading: B.1.2 UPPERBOUNDING ‖T̂ − T‖F
text: In this section, we prove that given polynomial number of samples, ‖T̂ − T‖F is small with high probability. We do this by standard matrix concentration inequalities. Note that our requirements on the norm of x is just for convenience, and the same proof works as long as x has reasonable tail-behavior (e.g. sub-Gaussian).
Lemma 15. Assume that ‖x‖ ≤ Γ, ‖A‖ ≤ P1, ‖ξ‖ ≤ P2 and σmin(E[xx>]) ≥ γ > 0, then for any ≤ γ/2, for any 1 > δ > 0, given O(d 3Γ14(ΓP1 √ k+P2) 6 log( dδ )
γ4 2 ) number of i.i.d. samples, we know
‖T̂ − T‖F ≤ ,
with probability at least 1− δ.
Proof. In order to get an upper bound for ‖T̂ −T‖F , we first show that ‖T̂ −T‖2 is upper bounded. We know
‖T − T̂‖ = max v∈Rk2+k:‖v‖≤1 ‖(T − T̂ )v‖
≤ max U∈Rk×ksym:‖U‖F≤ √ 2 ‖(T − T̂ )vec∗(U)‖.
For any k × k symmetric matrix U with eigenvalue decomposition U = ∑k i=1 λiu
(i)(u(i))>, according to the definition of T , we know
max U∈Rk×ksym:‖U‖F≤ √ 2 ‖(T − T̂ )vec∗(U)‖ = max U∈Rk×ksym:‖U‖F≤ √ 2 ‖Tvec∗(U)− T̂vec∗(U)‖
≤ max u(i):‖u(i)‖≤ 4 √ 2 ‖ k∑ i=1 λi(f(u (i))− f̂(u(i)))‖
≤ max u(i):‖u(i)‖≤ 4 √ 2 k∑ i=1 |λi|‖f(u(i))− f̂(u(i)‖
≤ ( k∑ i=1 |λi| ) max u:‖u‖≤ 4 √ 2 ‖f(u)− f̂(u)‖
≤ √ k √√√√ k∑ i=1 λ2i max u:‖u‖≤ 4 √ 2 ‖f(u)− f̂(u)‖ = √ k‖U‖F max
u:‖u‖≤ 4 √ 2 ‖f(u)− f̂(u)‖
≤ √
2k max u:‖u‖≤2
‖f(u)− f̂(u)‖
where f̂(u) = T̂vec∗(uu>) and the fourth inequality uses the Cauchy-Schwarz inequality. Next, we only need to upper bound maxu:‖u‖≤2 ‖f(u)− f̂(u)‖. Recall that
f(u) =2E [ (u>y) · (E[(u>y)x>]E[xx>]−1x) · (x⊗ x) ] − E [ (u>y)2(x⊗ x) ] + ( E [ (u>y)2 ] − 2E [ (u>y)x> ] E [ xx>
]−1E[(u>y)x])E[x⊗ x], and
f̂(u) =2Ê [ (u>y) · (Ê[(u>y)x>]Ê[xx>]−1x) · (x⊗ x) ] − Ê [ (u>y)2(x⊗ x) ] + ( Ê [ (u>y)2 ] − 2Ê [ (u>y)x> ] Ê [ xx>
]−1Ê[(u>y)x])Ê[x⊗ x]. Notice that( E [ (u>y) ·(E[(u>y)x>]E[xx>]−1x) ·(x⊗x) ])> = E [ (u>y)x> ] E [ xx>
]−1E[(u>y)x(x⊗x)>]. We first show that given polynomial number of samples,∥∥∥2Ê[(u>y)x>]Ê[xx>]−1Ê[(u>y)x(x⊗ x)>]− 2E[(u>y)x>]E[xx>]−1E[(u>y)x(x⊗ x)>]∥∥∥ is upper bounded with high probability. Since each row of W has unit norm, we have ‖W‖ ≤ √ k. Due to the assumption that ‖x‖ ≤ Γ, ‖A‖ ≤ P1, ‖ξ‖ ≤ P2, we have
‖(u>y)x>‖ ≤‖u‖‖Aσ(Wx) + ξ‖‖x‖ ≤‖u‖(‖A‖‖W‖‖x‖+ ‖ξ‖)‖x‖
≤2Γ(ΓP1 √ k + P2).
According to Lemma 24, we know given O(Γ 2(ΓP1
√ k+P2) 2 log( dδ )
2 ) number of samples,∥∥∥Ê[(u>y)x>]− E[(u>y)x>]∥∥∥ ≤ , with probability at least 1− δ.
Similarly, we can show that given O(Γ 6(ΓP1
√ k+P2) 2 log( dδ )
2 ) number of samples,∥∥∥Ê[(u>y)x(x⊗ x)>]− E[(u>y)x(x⊗ x)>]∥∥∥ ≤ , with probability at least 1− δ.
Since ‖xx>‖ ≤ Γ2, we know that given O(Γ 4 log( dδ ) 2 ) number of samples,∥∥∥Ê[xx>]− E[xx>]∥∥∥ ≤ , with probability at least 1− δ. Suppose that ≤ γ/2 ≤ σmin(E[xx>])/2, we know Ê[xx>] has full rank. According to Lemma 29, we have∥∥∥Ê[xx>]−1 − E[xx>]−1∥∥∥ ≤ 2√2∥∥Ê[xx>]− E[xx>]∥∥
σ2min(E[xx>]) ≤ 2 √ 2 /γ2,
with probability at least 1− δ.
By union bound, we know for any < γ/2, given O(Γ 6(ΓP1
√ k+P2) 2 log( dδ )
2 ) number of samples, with probability at least 1− δ, we have∥∥∥Ê[(u>y)x>]− E[(u>y)x>]∥∥∥ ≤ ,∥∥∥Ê[xx>]−1 − E[xx>]−1∥∥∥ ≤ 2√2 /γ2,∥∥∥Ê[(u>y)x(x⊗ x)>]− E[(u>y)x(x⊗ x)>]∥∥∥ ≤ . Define
E1 := Ê [ (u>y)x> ] − E [ (u>y)x> ] ,
E2 := Ê [ xx> ]−1 − E[xx>]−1, E3 := Ê [ (u>y)x(x⊗ x)> ] − E [ (u>y)x(x⊗ x)> ] .
Then, we have∥∥∥2Ê[(u>y)x>]Ê[xx>]−1Ê[(u>y)x(x⊗ x)>]− 2E[(u>y)x>]E[xx>]−1E[(u>y)x(x⊗ x)>]∥∥∥ ≤2 ∥∥E1∥∥∥∥E2∥∥∥∥E3∥∥+ 2∥∥E1∥∥∥∥E[xx>]−1∥∥∥∥E[(u>y)x(x⊗ x)>]∥∥
+ 2 ∥∥E[(u>y)x>]∥∥∥∥E2∥∥∥∥E[(u>y)x(x⊗ x)>]∥∥+ 2∥∥E[(u>y)x>]∥∥∥∥E[xx>]−1∥∥∥∥E3∥∥
+ 2 ∥∥E1∥∥∥∥E2∥∥∥∥E[(u>y)x(x⊗ x)>]∥∥+ 2∥∥E[(u>y)x>]∥∥∥∥E2∥∥∥∥E3∥∥+ 2∥∥E1∥∥∥∥E[xx>]−1∥∥∥∥E3∥∥
≤2 (2√2 3
γ2 +
2Γ3(ΓP1 √ k + P2)
γ +
8 √ 2Γ4(ΓP1 √ k + P2) 2
γ2 +
2Γ(ΓP1 √ k + P2)
γ
+ 4 √ 2Γ3(ΓP1 √ k + P2) 2
γ2 +
4 √ 2Γ(ΓP1 √ k + P2) 2
γ2 + 2 γ ) =O( Γ4(ΓP1 √ k + P2) 2
γ2 ).
Thus, given O(Γ 14(ΓP1
√ k+P2) 6 log( dδ )
γ4 2 ) number of samples, we know∥∥∥2Ê[(u>y) · (Ê[(u>y)x>]Ê[xx>]−1x) · (x⊗ x)]− 2E[(u>y) · (E[(u>y)x>]E[xx>]−1x) · (x⊗ x)]∥∥∥ = ∥∥∥2Ê[(u>y)x>]Ê[xx>]−1Ê[(u>y)x(x⊗ x)>]− 2E[(u>y)x>]E[xx>]−1E[(u>y)x(x⊗ x)>]∥∥∥ ≤ ,
with probability at least 1− δ. Now, let’s consider the second term∥∥∥Ê[(u>y)2(x⊗ x)]− E[(u>y)2(x⊗ x)]∥∥∥. Since
∥∥(u>y)2(x ⊗ x)∥∥ ≤ 4Γ2(ΓP1√k + P2)2, according to Lemma 24, we know given O( Γ4(ΓP1 √ k+P2) 4 log( dδ )
2 ),∥∥∥Ê[(u>y)2(x⊗ x)]− E[(u>y)2(x⊗ x)]∥∥∥ ≤ , with probability at least 1− δ. Next, let’s look at the third term∥∥∥Ê[(u>y)2]Ê[x⊗ x]− E[(u>y)2]E[x⊗ x]∥∥∥. Again, using Lemma 24 and union bound, we know given O(Γ 2(ΓP1 √ k+P2) 2 log( dδ )
2 ) number of samples, we have ∥∥∥Ê[(u>y)2]− E[(u>y)2]∥∥∥ ≤ ,∥∥∥Ê[x⊗ x]− E[x⊗ x]∥∥∥ ≤ . Thus, Define
E4 := Ê [ (u>y)2 ] − E [ (u>y)2 ] E5 := Ê[x⊗ x]− E[x⊗ x].
Then, we have∥∥∥Ê[(u>y)2]Ê[x⊗ x]− E[(u>y)2]E[x⊗ x]∥∥∥ ≤∥∥E4∥∥∥∥E5∥∥+ ∥∥E4∥∥∥∥E[x⊗ x]∥∥+ ∥∥E[(u>y)2]∥∥∥∥E5∥∥ ≤ 2 + Γ2 + 4(ΓP1 √ k + P2) 2
=O((ΓP1 √ k + P2) 2 ).
Thus, we know that given O(Γ 2(ΓP1
√ k+P2) 6 log( dδ )
2 ) number of samples, we know∥∥∥Ê[(u>y)2]Ê[x⊗ x]− E[(u>y)2]E[x⊗ x]∥∥∥ ≤ , with probability at least 1− δ. Now, let’s bound the last term,∥∥∥2(Ê[(u>y)x>]Ê[xx>]−1Ê[(u>y)x])Ê[x⊗x]−2(E[(u>y)x>]E[xx>]−1E[(u>y)x])E[x⊗x]∥∥∥. Similar as the first term, we can show that given O(Γ 10(ΓP1 √ k+P2) 6 log( dδ )
γ4 2 ) number of samples, we have∥∥∥2(Ê[(u>y)x>]Ê[xx>]−1Ê[(u>y)x])Ê[x⊗x]−2(E[(u>y)x>]E[xx>]−1E[(u>y)x])E[x⊗x]∥∥∥ ≤ with probability at least 1− δ. Now, we are ready to combine our bound for each of four terms. By union bound, we know given O( Γ14(ΓP1 √ k+P2) 6 log( dδ )
γ4 2 ) number of samples,∥∥∥2Ê[(u>y) · (Ê[(u>y)x>]Ê[xx>]−1x) · (x⊗ x)]− 2E[(u>y) · (E[(u>y)x>]E[xx>]−1x) · (x⊗ x)]∥∥∥ ≤ ∥∥∥Ê[(u>y)2(x⊗ x)]− E[(u>y)2(x⊗ x)]∥∥∥ ≤ ,∥∥∥Ê[(u>y)2]Ê[x⊗ x]− E[(u>y)2]E[x⊗ x]∥∥∥ ≤ ,∥∥∥2(Ê[(u>y)x>]Ê[xx>]−1Ê[(u>y)x])Ê[x⊗ x]− 2(E[(u>y)x>]E[xx>]−1E[(u>y)x])E[x⊗ x]∥∥∥ ≤ ,
hold with probability at least 1− δ. Thus, we know
max u:‖u‖≤2
‖f(u)− f̂(u)‖ ≤ + + + = 4
with probability at least 1− δ. Recall that
‖T̂ − T‖F ≤ √ k2 + k‖T̂ − T‖
≤k √
2k max u:‖u‖≤2
‖f(u)− f̂(u)‖,
where the second inequality holds since ‖T̂ − T‖ ≤ √
2kmaxu:‖u‖≤2 ‖f(u)− f̂(u)‖.
Thus, we know given O(Γ 14(ΓP1
√ k+P2) 6 log( dδ )
γ4 2 ) number of samples, ‖T̂ − T‖F ≤ k √ 2k ≤ d √ 2d
with probability at least 1− δ. Thus, given O(d 3Γ14(ΓP1
√ k+P2) 6 log( dδ )
γ4 2 ) number of samples,
‖T̂ − T‖F ≤ with probability at least 1− δ.

heading: B.2 ROBUST ANALYSIS FOR SIMULTANEOUS DIAGONALIZATION
text: In this section, we will show that the simultaneous diagonalization step in our algorithm is robust. Let S and Ŝ be two (k2 + k) by k matrices, whose columns consist of the least k right singular vectors of T and T̂ respectively.
According to Lemma 11, we know with polynomial number of samples, the Frobenius norm of SS> − ŜŜ⊥ is well bounded. However, due to the rotation issue of subspace basis, we cannot conclude that ‖S − Ŝ‖F is small. Only after appropriate alignment, the difference between S and Ŝ becomes small. Lemma 16. Let S and Ŝ be two (k2 + k) by k matrices, whose columns consist of the least k right singular vectors of T and T̂ respectively. If ‖SS> − ŜŜ>‖F ≤ , there exists an rotation matrix R ∈ Rk×k satisfying RR> = R>R = Ik, such that
‖Ŝ − SR‖F ≤ 2 .
Proof. Since S has orthonormal columns, we have σk(SS>) = 1. Then, according to Lemma 35, we know there exists rotation matrix R such that
‖Ŝ − SR‖F ≤ ‖SS> − ŜŜ>‖F√
2( √ 2− 1) √ σk(SS>) ≤ 2 .
Let the k columns of S be vec∗(U1), vec∗(U2), · · · , vec∗(Uk). Note each Ui can be expressed as A−>DiA
−1, where Di is a diagonal matrix. LetQ be a k×k matrix, whose i-th column consists of the diagonal elements of Di, such that Qij equals the j-th diagonal element of Di. Let vec∗(X) = SRζ1, vec∗(Y ) = SRζ2, whereR is the rotation matrix in Lemma 16 and ζ1, ζ2 are two independent standard Gaussian vectors. Let DX = diag(QRζ1) and DY = diag(QRζ2). It’s not hard to check that X = A−>DXA−1 and Y = A−>DYA−1. Furthermore, we have XY −1 = A−>DXD−1Y A
>. Next, we show that the diagonal elements of DXD−1Y are well separated. Lemma 17. Assume that ‖A‖ ≤ P1, σmin(A) ≥ β. Then for any δ > 0 we know with probability at least 1− δ, we have sep(DXD−1Y ) ≥ poly(1/d, 1/P1, β, δ), where sep(DXD−1Y ) := mini 6=j |(DXD −1 Y )ii − (DXD −1 Y )jj |.
Proof. We first show that matrix Q is well-conditioned. Since Ui = A−>DiA−1, we have vec(Ui) = A−>⊗A−>vec(Di). Let U be a k2×k matrix whose columns consist of vec(Ui)’s. Also define Q̄ as a k2 × k matrix whose columns are vec(Di)’s. Note that matrix Q̄ only has k non-zero rows, which are exactly matrix Q. With the above definition, we have U = A−> ⊗ A−>Q̄. Since σmin(U) ≤ ‖A−> ⊗A−>‖σmin(Q̄), we have
σmin(Q̄) ≥ σmin(U)
‖A−> ⊗A−>‖ .
Notice that a subset of rows of U constitute matrix S, which is an orthonormal matrix. Thus, we have σmin(U) ≥ σmin(S) = 1. Since we assume σmin(A) ≥ β, we have
‖A−> ⊗A−>‖ = ‖A−>‖2 = 1 σmin(A)2 ≤ 1 β2 .
Thus, we have σmin(Q̄) ≥ β2, which implies σmin(Q) ≥ β2. We also know ‖U‖ ≥ σmin(A−> ⊗A−>)‖Q̄‖, thus
‖Q̄‖ ≤ ‖U‖ σmin(A−> ⊗A−>) .
Since ‖S‖ = 1, we know ‖U‖ ≤ √ 2. For the smallest singular value of A−> ⊗A−>, we have
σmin(A −> ⊗A−>) = σ2min(A−1) =
1 ‖A‖2 ≥ 1 P 21 .
Thus, we have ‖Q̄‖ ≤ √ 2P 21 , which implies ‖Q‖ ≤ √ 2P 21 .
Now, let’s prove that the diagonal elements of DXD−1Y are well-separated. Let q > i be the i-th row vector of Q. Then we know the i-th diagonal element of DXD−1Y is 〈qi,Rζ1〉 〈qi,Rζ2〉 . Since ‖Q‖ ≤ √ 2P 21 ,
we have ‖qi‖ ≤ √ 2P 21 for every row vector.
It’s not hard to show that with probability at least 1 − exp(−dΩ(1)), we have |〈qi, Rζ2〉| ≤ poly(d, P1) for each i. Now given ζ2 for which this happens, we have
〈qi,Rζ1〉 〈qi,Rζ2〉 − 〈qj ,Rζ1〉 〈qj ,Rζ2〉 =
ci〈qi, Rζ1〉 − cj〈qj , Rζ1〉, where ci, cj have magnitude as least poly(1/d, 1/P1). Since σmin(Q) ≥ β2, we know ‖Projq⊥j qi‖ ≥ β
2 (because otherwise there exists λ such that ‖(ei + λej)‖σmin(Q) ≤ ‖(ei + λej)>Q‖ = ‖Projq⊥j qi‖ < β
2, which is a contradiction). Let q⊥i,j = Projq⊥j qi = qi − λi,jqj , we can rewrite this as
〈qi, Rζ1〉 〈qi, Rζ2〉 − 〈qj , Rζ1〉 〈qj , Rζ2〉 = ci〈qi, Rζ1〉 − cj〈qj , Rζ1〉 = ci〈q⊥i,j , Rζ1〉 − (cj + λi,jci)〈qj , Rζ1〉.
By properties of Gaussians, we know 〈q⊥i,j , Rζ1〉 is independent of 〈qj , Rζ1〉, so we can first fix 〈qj , Rζ1〉 and apply anti-concentration of Gaussians (see Lemma 38) to 〈q⊥i,j , Rζ1〉. As a result we know with probability at least 1− δ/k2:
∣∣∣ 〈qi, Rζ1〉〈qi, Rζ2〉 − 〈qj , Rζ1〉〈qj , Rζ2〉 ∣∣∣ ≥ poly(1/d, 1/P1, β, δ).
By union bound, we know with probability at least 1− δ,
sep(DXD−1Y ) ≥ poly(1/d, 1/P1, β, δ).
Let X̂ = Ŝζ1 and Ŷ = Ŝζ2. Next, we prove that the eigenvectors of X̂Ŷ −1 are close to the eigenvectors of XY −1.
Lemma 12. Suppose that ‖SS> − ŜŜ>‖F ≤ , ‖A‖ ≤ P1, ‖ξ‖ ≤ P2, σmin(E[xx>]) ≥ γ, σmin(A) ≥ β. Let X̂ = mat∗(Ŝζ1), Ŷ = mat∗(Ŝζ2), where ζ1 and ζ2 are two independent standard Gaussian vectors. Let z1, · · · zk be the normalized row vectors of A−1. Let ẑ1, ..., ẑk be the eigenvectors of X̂Ŷ −1 (after sign flip). For any δ > 0 and small enough , with O( (ΓP1 √ k+P2) 2 log(d/δ) 2 ) number of i.i.d. samples in Ê[ẑ > i y], with probability at least 1− δ over the randomness of ζ1, ζ2 and i.i.d. samples, there exists a permutation π(i) ∈ [k] such that ‖ẑi − zπ[i]‖ ≤ poly(d, P1, 1/β, , 1/δ), for any 1 ≤ i ≤ k.
Proof. Let z′1, · · · , z′k be the eigenvectors of XY −1 (before sign flip step). Similarly define ẑ′1, · · · , ẑ′k for X̂Ŷ −1. We first prove that the eigenvectors of X̂Ŷ −1 are close to the eigenvectors of XY −1.
Let X̂ = X + EX and Ŷ = Y + EY . Then we have X̂Ŷ −1 = XY −1(I + F ) +G,
where F = −EY (I + Y −1EY )−1Y −1 and G = EX Ŷ −1. According to Lemma 34, we have ‖F‖ ≤ ‖EY ‖σmin(Y )−‖EY ‖ and ‖G‖ ≤ ‖EX‖ σmin(Ŷ )
. In order to bound the perturbation matrices ‖F‖ and ‖G‖, we need to first bound ‖EX‖, ‖EY ‖ and σmin(Y ), σmin(Ŷ ).
As we know, EX = X̂ − X = (Ŝ − SR)ζ1. According to Lemma 16, we have ‖Ŝ − SR‖ ≤ ‖Ŝ − SR‖F ≤ 2 . We also know with probability at least 1− exp(−dΩ(1)), ‖ζ1‖ ≤ poly(d). Thus, we have
‖EX‖ = ‖(Ŝ − SR)ζ1‖ ≤ ‖Ŝ − SR‖‖ζ1‖ ≤ poly( , d). Similarly, with probability at least 1− exp(−dΩ(1)), we also know ‖EY ‖ ≤ poly( , d). Now, we lower bound the smallest singular value of X and Y . Since X = A−>DXA−1, we have
σmin(X) ≥σ2min(A−1)σmin(DX)
= 1
‖A‖2 σmin(DX)
≥ 1 P 21 σmin(DX).
Since DX is a diagonal matrix, its smallest singular value equals the smallest absolute value of its diagonal element. Recall each diagonal element of DX is 〈qi, Rζ1〉, which follows a Gaussian distribution whose standard deviation is at least ‖qi‖ ≥ β2. By anti-concentration property of Gaussian (see Lemma 38), we know |〈qi, Rζ2〉| ≥ Ω(δβ2/k) for all i with probability 1 − δ/4. Thus, we have σmin(X) ≥ poly(1/d, 1/P1, β, δ). Similarly we have the same conclusion for Y . For small enough EY , we have σmin(Ŷ ) ≥ σmin(Y )− ‖EY ‖. Thus, for small enough , we have ‖F‖ ≤ poly(d, P1, 1/β, , 1/δ) and ‖G‖ ≤ poly(d, P1, 1/β, , 1/δ). In order to apply Lemma 33, we also need to bound κ(A−>) and ‖XY −1‖. Since σmin(A−>) ≥ 1P1 and ‖A
−>‖ ≤ 1/β, we have κ(A−>) ≤ P1/β. For the norm of XY −1, we have
‖XY −1‖ ≤ ‖X‖ σmin(Y ) ≤ ‖X‖poly(d, P1, 1/β, 1/δ),
where the second inequality holds because σmin(Y ) ≥ poly(1/d, 1/P1, β, δ). Recall that X = mat∗(SRζ1). It’s not hard to verify that with probability at least 1− exp(−dΩ(1)), we have ‖X‖ ≤ poly(d). Thus, we know ‖XY −1‖ ≤ poly(d, P1, 1/β, 1/δ). Similarly, we can also prove that ‖DXD−1Y ‖ ≤ poly(d, P1, 1/β, 1/δ) Overall, we have κ(A−>)(‖XY −1F‖+ ‖G‖) ≤κ(A−>)(‖XY −1‖‖F‖+ ‖G‖)
≤P1 β
( poly(d, P1, 1/β, 1/δ)poly(d, P1, 1/β, , 1/δ) + poly(d, P1, 1/β, , 1/δ) ) ≤poly(d, P1, 1/β, , 1/δ).
According to Lemma 17, we know with probability at least 1 − δ/4, sep(DXD−1Y ) ≥ poly(1/d, 1/P1, β, δ). Thus, by union bound, we know for small enough , with probability at least 1− δ, κ(A−>)(‖XY −1F‖+ ‖G‖) < sep(DXD−1Y )/(2k). According to Lemma 33, we know there exists a permutation π[i] ∈ [k], such that
‖ẑ′i − z′i‖ ≤3 ‖F‖‖DXD−1Y ‖+ ‖G‖ σmin(A−>)sep(DXD−1Y )
≤3poly(d, P1, 1/β, , 1/δ)poly(d, P1, 1/β, 1/δ) + poly(d, P1, 1/β, , 1/δ) 1/P1poly(1/d, 1/P1, β, 1/δ)
≤poly(d, P1, 1/β, , 1/δ).
with probability at least 1− δ. According to Lemma 5, the eigenvectors of XY −1 (after sign flip) are exactly the normalized rows of A−1 (up to permutation). Now the only issue is the sign of ẑ′i. By the robustness of sign flip step (see Lemma 18), we know for small enough , with O( (ΓP1 √ k+P2) 2 log(d/δ) 2 ) number of i.i.d. samples in Ê[ẑ>i y], with probability at least 1 − δ, the sign flip of ẑ′i is consistent with the sign flip of z′i.
In the following lemma, we show that the sign flip step of ẑi is robust. Lemma 18. Suppose that ‖x‖ ≤ Γ, ‖A‖ ≤ P1, ‖ξ‖ ≤ P2. Let z′1, · · · , z′k be the eigenvectors of XY −1 (before sign flip step). Similarly define ẑ′1, · · · , ẑ′k for X̂Ŷ −1. Suppose for each i, ‖z′i−ẑ′i‖ ≤ , where ≤ βγ
4Γ(1+ΓP1 √ k+P2)
. We know, for any δ < 1, with O( (ΓP1 √ k+P2) 2 log(d/δ) 2 ) number of
i.i.d. samples, Pr [ sign(E[〈z′i, y〉]) 6= sign(Ê[〈ẑ′i, y〉]) ] ≤ δ.
Proof. We first show that E[〈z′i, y〉] is bounded away from zero. Let Z ′ be a k × k matrix, whose rows are {z′i}. Without loss of generality, assume that Z ′ = diag(±λ)A−1. Since ‖A−1‖ ≤ 1/β, we know λi ≥ β, for each i. Thus, we have
|E[〈z′i, y〉]| = λiE[σ(w>i x)] ≥ βE[σ(w>i x)].
In order to lowerbound E[σ(w>i x)], let’s first look at E[σ(w>i x)x>].
‖E[σ(w>i x)x>]‖ = 1
2 ‖E[w>i xx>]‖
≥1 2 ‖wi‖σmin(E[xx>]) ≥γ 2 .
Now, let’s connect ‖E[σ(w>i x)x>]‖ with E[σ(w>i x)].
‖E[σ(w>i x)x>]‖ ≤E[‖σ(w>i x)x>‖] =E[σ(w>i x)‖x‖] ≤ΓE[σ(w>i x)].
Thus, we have E[σ(w>i x)] ≥ γ 2Γ , and further |E[〈z ′ i, y〉]| ≥ βγ 2Γ .
Now, let’s bound ‖E[〈z′i, y〉]− Ê[〈ẑ′i, y〉]‖.
‖E[〈z′i, y〉]− Ê[〈ẑ′i, y〉]‖ =‖E[〈z′i, y〉]− E[〈ẑ′i, y〉] + E[〈ẑ′i, y〉]− Ê[〈ẑ′i, y〉]‖ ≤‖E[〈z′i − ẑ′i, y〉]‖+ ‖E[〈ẑ′i, y〉]− Ê[〈ẑ′i, y〉]‖
≤(ΓP1 √ k + P2) + ‖E[〈ẑ′i, y〉]− Ê[〈ẑ′i, y〉]‖.
Note that we reserve fresh samples for this step, thus ẑ′i is independent with samples in Ê[〈ẑ′i, y〉]]. Since ‖〈ẑ′i, y〉‖ ≤ ΓP1 √ k + P2, we know with O( (ΓP1 √ k+P2) 2 log(d/δ) 2 ) number of samples, we have ‖E[〈ẑ′i, y〉]− Ê[〈ẑ′i, y〉]]‖ ≤ . Overall, we have ‖E[〈z′i, y〉] − Ê[〈ẑ′i, y〉]]‖ ≤ (1 + ΓP1 √ k + P2) . Combined with the fact that |E[〈z′i, y〉]| ≥ βγ 2Γ , we know as long as ≤ βγ 4Γ(1+ΓP1 √ k+P2) , with O( (ΓP1 √ k+P2) 2 log(d/δ) 2 ) number of samples, we have Pr[sign(E[〈z′i, y〉]) 6= sign(Ê[〈ẑ′i, y〉])] ≤ δ.

heading: B.3 ROBUST ANALYSIS FOR RECOVERING FIRST LAYER WEIGHTS
text: We will first show that Algorithm 1 is robust.
Theorem 8. Assume that ‖x‖ ≤ Γ, ‖w‖ ≤ 1, |ξ| ≤ P2 and σmin(E[xx>]) ≥ γ. Then for any ≤ γ/2, for any 1 > δ > 0, given O( (Γ 2+P2Γ) 4 log( dδ )
γ4 2 ) number of i.i.d. samples, we know with probability at least 1− δ,
‖ŵ − w‖ ≤ ,
where ŵ is the learned weight vector.
Proof. We first show that given polynomial number of i.i.d. samples, ‖Ê[yx] − E[yx]‖ is upper bounded with high probability, where Ê[yx] is the empirical estimate of E[yx]. Due to the assumption that ‖x‖ ≤ Γ, ‖w‖ ≤ 1, ‖ξ‖ ≤ P2, we have
‖yx‖ =‖(σ(w>x) + ξ)x‖ ≤‖w>xx‖+ ‖ξx‖ ≤‖w‖‖x‖2 + |ξ|‖x‖ ≤Γ2 + P2Γ
According to Lemma 24, we know given O( (Γ 2+P2Γ) 2 log( dδ )
2 ) number of samples, with probability at least 1− δ, we have ‖Ê[yx]− E[yx]‖ ≤ .
Since ‖xx>‖ ≤ Γ2, we know that given O(Γ 4 log( dδ ) 2 ) number of samples, ‖Ê[xx >]− E[xx>]‖ ≤ with probability at least 1− δ. Suppose that ≤ γ/2 ≤ σmin(E[xx>])/2, we know Ê[xx>] has full rank. According to Lemma 29, we have∥∥∥Ê[xx>]−1 − E[xx>]−1∥∥∥ ≤ 2√2∥∥Ê[xx>]− E[xx>]∥∥
σ2min(E[xx>]) ≤ 2 √ 2 /γ2,
with probability at least 1− δ.
By union bound, we know for any ≤ γ/2, given O( (Γ 2+P2Γ) 2 log( dδ )
2 ) number of samples, with probability at least 1− δ, we have
‖Ê[xx>]−1 − E[xx>]−1‖ ≤ 2 √ 2 /γ2
‖Ê[yx]− E[yx]‖ ≤
Denote
E1 :=Ê[xx>]−1 − E[xx>]−1
E2 :=Ê[yx]− E[yx].
Then, we have
‖ŵ − w‖ =2‖Ê[xx>]−1Ê[yx]− E[xx>]−1E[yx]‖ ≤2‖E1‖‖E2‖+ 2‖E1‖‖E[yx]‖+ 2‖E[xx>]−1‖‖E2‖
≤4 √ 2 2
γ2 +
4 √ 2(Γ2 + P2Γ)
γ2 +
2
γ
≤O( (Γ 2 + P2Γ)
γ2 ).
Thus, given O( (Γ 2+P2Γ) 4 log( dδ )
γ4 2 ) number of samples, with probability at least 1− δ, we have
‖ŵ − w‖ ≤ .
Now let’s go back to the call to Algorithm 1 in Algorithm 3. Let zi’s be the normalized rows of A−1, and let ẑi’s be the eigenvectors of X̂Ŷ −1 (with correct sign). From Lemma 12, we know {ẑi} are close to {zi} with permutation. Without loss of generality, we assume the permutation here is just an identity mapping, which means ‖zi − ẑi‖ is small for each i. For each zi, let vi be the output of Algorithm 1 given infinite number of inputs (x, z>i y). For each ẑi, let v̂i be the output of Algorithm 1 given only finite number of samples (x, ẑ>i y). In this section, we show that suppose ‖zi − ẑi‖ is bounded, with polynomial number of samples, ‖vi − v̂i‖ is also bounded.
The input for Algorithm 1 is (x, ẑ>i y). We view ẑ > i y as the summation of z > i y and a noise term (ẑi − zi)>y. Here, the issue is that the noise term (ẑi − zi)>y is not independent with the sample (x, ẑ>i y), which makes the robust analysis in Theorem 8 not applicable. On the other hand, since we reserve a separate set of samples for Algorithm 1, the estimate ẑi is independent with the samples (x, y)’s used by Algorithm 1. Thus, the samples (x, ẑi>y)’s here are still i.i.d., which enables us to use matrix concentration bounds to show the robustness here. Lemma 13. Assume that ‖x‖ ≤ Γ, ‖A‖ ≤ P1, ‖ξ‖ ≤ P2 and σmin(E[xx>]) ≥ γ. Suppose that for each 1 ≤ i ≤ k, ‖ẑi − zi‖ ≤ τ . Then for any ≤ γ/2 and δ < 1, given O( (Γ2+P2Γ) 4 log( dδ )
γ4 2 )
number of samples for Algorithm 1, we know with probability at least 1− δ,
‖vi − v̂i‖ ≤ 2τ(ΓP1
√ k + P2)Γ
γ + 2 ,
for each 1 ≤ i ≤ k.
Proof. For any i, let’s bound ‖vi − v̂i‖. As we know, vi = 2E[xx>]−1E[z>i yx] and v̂i = 2Ê[xx>]−1Ê[ẑ>i yx]. Thus, in order to bound ‖vi − v̂i‖, we only need to show∥∥E[xx>]−1E[ẑ>i yx]− E[xx>]−1E[z>i yx]∥∥ and ∥∥E[xx>]−1E[ẑ>i yx]− Ê[xx>]−1Ê[ẑ>i yx]∥∥ are both bounded.
The first term can be bounded as follows.∥∥E[xx>]−1E[ẑ>i yx]− E[xx>]−1E[z>i yx]∥∥ =∥∥E[xx>]−1E[(ẑi − zi)>yx]∥∥ ≤ ∥∥E[xx>]−1∥∥∥∥ẑi − zi∥∥∥∥Aσ(Wx) + ξ∥∥∥∥x∥∥ ≤τ(ΓP1 √ k + P2)Γ
γ
We can use standard matrix concentration bounds to upper bound the second term. By similar analysis of Theorem 1, we know given O( (Γ 2+P2Γ) 4 log( dδ )
γ4 2 ) number of i.i.d. samples, with probability at least 1− δ, ∥∥E[xx>]−1E[ẑ>i yx]− Ê[xx>]−1Ê[ẑ>i yx]∥∥ ≤ .
Overall, we have ‖vi − v̂i‖ = ∥∥2E[xx>]−1E[z>i yx]− 2Ê[xx>]−1Ê[ẑ>i yx]∥∥ ≤2 ∥∥E[xx>]−1E[ẑ>i yx]− E[xx>]−1E[z>i yx]∥∥+ 2∥∥E[xx>]−1E[ẑ>i yx]− Ê[xx>]−1Ê[ẑ>i yx]∥∥
≤2τ(ΓP1 √ k + P2)Γ
γ + 2 .
By union bound, we know given O( (Γ 2+P2Γ) 4 log( dδ )
γ4 2 ) number of i.i.d. samples, with probability at least 1− δ,
‖vi − v̂i‖ ≤ 2τ(ΓP1
√ k + P2)Γ
γ + 2 .
for any 1 ≤ i ≤ k.

heading: B.4 PROOF OF THEOREM 7
text: Proof of Theorem 7. Combining Lemma 11, Lemma 12 and Lemma 13, we know given poly ( Γ, P1, P2, d, 1/ , 1/γ, 1/α, 1/β, 1/δ ) number of i.i.d. samples, with probability at least 1− δ,
‖zi − ẑi‖ ≤ , ‖vi − v̂i‖ ≤ for each 1 ≤ i ≤ k. Let V be a k×dmatrix whose rows are vi’s. Similarly define matrix V̂ for v̂i’s. Since ‖vi− v̂i‖ ≤ for any i, we know every row vector of V − V̂ has norm at most , which implies ‖V − V̂ ‖ ≤ √ k .
Let Z be a k × k matrix whose rows are zi’s. Similarly define matrix Ẑ for ẑi’s. Again, we have ‖Z−Ẑ‖ ≤ √ k . In order to show ‖Z−1−Ẑ−1‖ is small using standard matrix perturbation bounds (Lemma 29), we need to lower bound σmin(Z). Notice that Z is just matrix A−1 with normalized row vectors. As we know, σmin(A−1) ≥ 1/P1, and ‖A−1‖ ≤ 1/β, which implies that every row vector of A−1 has norm at most 1/β. Let Dz be the diagonal matrix whose i, i-th entry is the norm of i-th row ofA−1, then Z = D−1z A −1, and we know σmin(Z) ≥ σmin(D−1z )σmin(A−1) ≥ β/P1.
Then, according to Lemma 29, as long as ≤ β2P1 , we have
‖Z−1 − Ẑ−1‖ ≤ 2 √ 2 ‖Z − Ẑ‖ σ2min(Z) ≤ 2 √ 2P 21 √ k /β2.
Define E1 :=Z −1 − Ẑ−1
E2 :=V − V̂ . We know ‖E1‖ ≤ 2 √ 2P 21 √ k /β2 and ‖E2‖ ≤ √ k . In order to bound ‖V ‖, we can bound the norm of its row vectors. We have, ‖vi‖ =‖2E[xx>]−1E[z>i yx]‖
≤2Γ(ΓP1 √ k + P2)
γ ,
which implies ‖V ‖ ≤ 2 √ kΓ(ΓP1 √ k+P2)
γ .Now we can bound ‖Z −1σ(V x)−Ẑ−1σ(V̂ x)‖ as follows.
‖Z−1σ(V x)− Ẑ−1σ(V̂ x)‖ ≤‖E1‖‖E2x‖+ ‖E1‖‖V x‖+ ‖Z−1‖‖E2x‖
≤2 √ 2P 21 Γk 2
β2 +
4 √
2P 21 Γ 2(ΓP1 √ k + P2)k β2γ + P1Γ √ k β ,
where the first inequality holds since ‖σ(V x)‖ ≤ ‖V x‖ and ‖σ(V̂ x)− σ(V x)‖ ≤ ‖V̂ x− V x‖. Thus, we know given poly ( Γ, P1, P2, d, 1/ , 1/γ, 1/α, 1/β, 1/δ ) number of i.i.d. samples, with probability at least 1− δ, ‖Aσ(Wx)− Ẑ−1σ(V̂ x)‖ = ‖Z−1σ(V x)− Ẑ−1σ(V̂ x)‖ ≤ ,
where the first equality holds because Aσ(Wx) = Z−1σ(V x), as shown in Theorem 5.

heading: C SMOOTHED ANALYSIS FOR DISTINGUISHING MATRICES
text: In smoothed analysis, it’s clear that after adding small Gaussian perturbations, matrix A and W will become robustly full rank with reasonable probability (Lemma 36). In this section, we will focus on the tricky part, using smoothed analysis framework to show that it is natural to assume the distinguishing matrix is robustly full rank. We will consider two settings. In the first case, the input distribution is the Gaussian distribution N (0, Id), and the weights for the first layer matrix W is perturbed by a small Gaussian noise. In this case we show that the augmented distinguishing matrix M has smallest singular value σmin(M) that depends polynomially on the dimension and the amount of perturbation. This shows that for the Gaussian input distribution, σmin(M) is lower bounded as long as W is in general position. In the second case, we will fix a full rank weight matrix W , and consider an arbitrary symmetric input distribution D. There is no standard way of perturbing a symmetric distribution, we give a simple perturbation D′ that can be arbitrarily close to D, and prove that σmin(MD ′ ) is lowerbounded.
Perturbing W for Gaussian Input We first consider the case when the input follows standard Gaussian distribution N (0, Id). The weight matrix W is perturbed to W̃ where
W̃ = W + ρE. (18)
Here E ∈ Rk×d is a random matrix whose entries are i.i.d. standard Gaussians. We will use M̃ to denote the perturbed version of the augmented distinguishing matrix M . Recall that the columns of M̃ has the form:
M̃ij = E [ (w̃>i x)(w̃ > j x)(x⊗ x)1{w̃>i xw̃>j x ≤ 0} ] ,
where w̃>i is the i-th row of W̃ . Also, since M̃ is the augmented distinguishing matrix it has a final column M̃0 = vec(Id). We show that the smallest singular value of M̃ is lower bounded with high probability. Theorem 9. Suppose that k ≤ d/5, and the input follows standard Gaussian distributionN (0, Id). Given any weight matrix W with ‖wi‖ ≤ τ for each row vector, let W̃ be a perturbed version of W according to Equation (18) and M̃ be the perturbed augmented distinguishing matrix. With probability at least 1− exp(−dΩ(1)), we have
σmin(M̃) ≥ poly(1/τ, 1/d, ρ).
We will prove this Theorem in Section C.1.
Perturbing the Input Distribution Our algorithm works for a general symmetric input distribution D. However, we cannot hope to get a result like Theorem 9 for every symmetric input distribution D. As a simple example, if D is just concentrated on 0, then we do not get any information about weights and the problem is highly degenerate. Therefore, we must specify a way to perturb the input distribution.
We define a perturbation that is parametrized by a random Gaussian matrix Q and a parameter λ ∈ (0, 1). The random matrix Q is used to generate a Gaussian distribution DQ with a random covariance matrix. To sample a point in DQ, first sample n ∼ N (0, Id), and then output Qn. The (Q,λ) perturbation of a distribution D, which we denote by DQ,λ is a mixture between the distributionD and the distributionDQ. More precisely, to sample x fromDQ,λ, pick z as a Bernoulli random variable where Pr[z = 1] = λ and Pr[z = 0] = 1 − λ; pick x′ according to D and pick x′′ = Qn according to distribution DQ, then let
x = { x′ z = 0 x′′ z = 1
Intuitively, the (Q,λ) perturbation of a distribution D mixes the distribution D with a Gaussian distribution DQ with covariance matrix QQ>. Since both D and DQ are symmetric, their mixture is also symmetric. Also, the TV-distance between D and DQ,λ is bounded by λ. Throughout this section we will use D′ to denote the perturbed distribution DQ,λ
We show that given any input distribution, after applying (Q,λ)-perturbation with a random Gaussian matrix Q, the smallest singular value of the augmented distinguishing matrix MD ′ is lower bounded. Recall that MD ′ is defined as
MD ′ ij = Ex∼D′ [ (w>i x)(w > j x)(x⊗ x)1{w>i xw>j x ≤ 0} ] ,
as the first ( k 2 ) columns and has Ex∼D′ [x⊗ x] as the last column. Theorem 10. Given weight matrix W with ‖wi‖ ≤ τ for each row vector and symmetric input distribution D. Suppose that k ≤ d/7 and σmin(W ) ≥ ρ, after applying (Q,λ)-perturbations to yield perturbed input distribution D′, where Q is a d× d matrix whose entries are i.i.d. Gaussians, we have with probability at least 1− exp(−dΩ(1)) over the randomness of Q,
σmin(M D′) ≥ poly(1/τ, 1/d, ρ, λ).
We will prove this later in Section C.3.

heading: C.1 SMOOTHED ANALYSIS FOR GAUSSIAN INPUTS
text: In this section, we will prove Theorem 9, as restated below: Theorem 9. Suppose that k ≤ d/5, and the input follows standard Gaussian distributionN (0, Id). Given any weight matrix W with ‖wi‖ ≤ τ for each row vector, let W̃ be a perturbed version of W according to Equation (18) and M̃ be the perturbed augmented distinguishing matrix. With probability at least 1− exp(−dΩ(1)), we have
σmin(M̃) ≥ poly(1/τ, 1/d, ρ).
To prove this theorem, recall the definition of Mij : mat(Mij) = E [ (w>i x)(w > j x)(xx >)1{w>i xw>j x ≤ 0} ] .
Since Gaussian distribution is highly symmetric, for every direction u that is orthogonal to both wi and wj , we have u>mat(Mij)u be a constant. We can compute this constant as
mij := E [ (w>i x)(w > j x)1{w>i xw>j x ≤ 0} ] .
This implies that if we consider mat(Mij) − mijId, it is going to be a matrix whose rows and columns are in span of wi and wj . In fact we can compute the matrix explicitly as the following lemma: Lemma 19. Suppose input x follows standard Gaussian distributionN (0, Id), and suppose weight matrix W has full-row rank, then for any 1 ≤ i < j ≤ k, we have
mat(Mij) = 1
π
( φij cos(φij)− sin(φij) ) ‖wi‖‖wj‖Id +
φij π (wiw > j + wjw > i )
− sin(φij) π ( ‖wi‖ ‖wj‖ wjw > j + ‖wj‖ ‖wi‖ wiw > i ),
where 0 < φij < π is the angle between weight vectors wi and wj .
Of course, the same lemma would be applicable to W̃ , so we have an explicit formula for M̃ij . We will bound the smallest singular value using the idea of leave-one-out distance (as previously used in Rudelson & Vershynin (2009)).
Leave-one-out Distance Leave-one-out distance is a metric that is closely related to the smallest singular value but often much easier to estimate. Definition 2. For a matrix A ∈ Rd×n(d ≥ n), the leave-one-out distance d(A) is defined to be the smallest distance between a column of A to the span of other columns. More precisely, let Ai be the i-th column of A and S−i be the span of all the columns except for Ai, then
d(A) := min i∈[n] ‖(Id − ProjS−i)Ai‖.
Rudelson & Vershynin (2009) showed that one can lowerbound the smallest singular value of a matrix by its leave-one-out distance.
Lemma 20 (Rudelson & Vershynin (2009)). For matrix A ∈ Rd×n(d ≥ n), we always have d(A) ≥ σmin(A) ≥ 1√nd(A).
Therefore, to bound σmin(M̃) we just need to lowerbound d(M̃). We use the ideas similar to Bhaskara et al. (2014) and Ma et al. (2016). Since every column of M̃ (except for M̃0) is random, we will try to show that even if we condition on all the other columns, because of the randomness in M̃ij , the distance between M̃ij to the span of other columns is large. However, there are several obstacles in this approach:
1. The augmented distinguishing matrix M̃ has a special column M̃0 = vec(Id) that does not have any randomness.
2. The closed form expression for M̃ (as in Lemma 19) has complicated coefficients that are not linear in the vectors w̃i and w̃j .
3. The columns of M̃ij are not independent with each other, so if we condition on all the other columns, M̃ij is no longer random.
To address the first obstacle, we will prove a stronger version of Lemma 20 that allows a special column.
Lemma 21. Let A ∈ Rd×(n+1) (d ≥ n + 1) be an arbitrary matrix whose columns are A0, A1, ..., An. For any i = 1, 2, ..., n, let S−i be the subspace spanned by all the other columns (includingA0) except forAi, and let d′(A) := mini=1,...,n ‖(Id−ProjS−i)Ai‖. Suppose the column A0 has norm √ d and A1, ..., An has norm at most C, then
σmin(A) ≥ min (√ nC2d
4nC2 + d ,
√ d
4n2C2 + nd d′(A)
) .
This lemma shows that if we can bound the leave-one-out distance for all but one column, then the smallest singular value of the matrix is still lowerbounded as long as the columns do not have very different norms. We defer the proof to Section C.2.
For the second obstacle, we show that these coefficients are lowerbounded with high probability. Therefore we can condition on the event that all the coefficients are large enough.
Lemma 22. Given weight vectors wi and wj with norm ‖wi‖, ‖wj‖ ≤ τ , let w̃i = wi + ρεi, w̃j = wj+ρεj where εi, εj are i.i.d. Gaussian random vectors. With probability at least 1−exp(−dΩ(1)), we know ‖w̃i‖ ≤ τ + √ 3ρ2d/2, ‖w̃j‖ ≤ τ + √ 3ρ2d/2 and
φ̃ij ≥ √ ρ2(d− 2)
√ 2τ + √ 3ρ2d ,
where φ̃ij is the angle between w̃i and w̃j . In particular, if W̃ = W + ρE where E is an i.i.d. Gaussian random matrix, with probability at least 1−exp(−dΩ(1)), for all i, ‖w̃i‖ ≤ τ+ √ 3ρ2d/2,
and for all i < j, the coefficient φ̃ij/π in front of the term w̃iw̃>j + w̃jw̃ > i is at least
√ ρ2(d−2)
( √ 2τ+ √ 3ρ2d)π .
This lemma intuitively says that after the perturbation w̃i and w̃j cannot be close to co-linear. We defer the detailed proof to Section C.2.
For the final obstacle, we use ideas very similar to Ma et al. (2016) which decouples the randomness of the columns.
Proof of Theorem 9. Let E1 be the event that Lemma 22 does not hold. Event E1 will be one of the bad events (but note that we do not condition on E1 not happening, we use a union bound at the end).
We partition [d] into two disjoint subsets L1, L2 of size d/2. Let M̃ ′ be the set of rows of M̃ indexed by L1 × L2. That is, the columns of M̃ ′ are
M̃ ′ij = φ̃ij π (w̃i,L1⊗ w̃j,L2 + w̃j,L1⊗ w̃i,L2)− sin(φ̃ij) π ( ‖w̃i‖ ‖w̃j‖ w̃j,L1⊗ w̃j,L2 + ‖w̃j‖ ‖w̃i‖ w̃i,L1⊗ w̃i,L2),
for i < j, where w̃i,L denotes the restriction of vector w̃i to the subset L. Note that the restriction of vec(Id) to the rows indexed by L1 × L2 is just an all zero vector.
We will focus on a column M̃ ′ij with i < j and try to prove it has a large distance to the span of all the other columns. Let Vij be the span of all other columns, which is equal to Vij = span{M̃ ′kl : k < l ∧ (k, l) 6= (i, j)} (note that we do not need to consider M̃0 because that column is 0 when restricted to L1 × L2.
It’s clear that Vij is correlated with M̃ ′ij , which is bad for the proof. To get around this problem, we follow the idea of Ma et al. (2016) and define the following subspace that contains Vij ,
V̂ij = span { w̃k,L1 ⊗ x, x⊗ w̃k,L2 , w̃j,L1 ⊗ x, x⊗ w̃i,L2 ∣∣∣k /∈ {i, j}, x ∈ Rd/2}. By definition Vij ⊂ V̂ij , and thus V̂ ⊥ij ⊂ V ⊥ij , where V ⊥ij denotes the orthogonal subspace of Vij . Observe that w̃j,L1 ⊗ w̃i,L2 , w̃j,L1 ⊗ w̃j,L2 , w̃i,L1 ⊗ w̃i,L2 ∈ V̂ij , thus
ProjV̂ ⊥ij M̃ ′ ij = φ̃ij π ProjV̂ ⊥ij ( w̃i,L1 ⊗ w̃j,L2 ) .
Note that w̃i,L1 ⊗ w̃j,L2 is independent with V̂ij . Moreover, subspace V̂ij has dimension at most (k− 2)d/2 + (k− 2)d/2 + d/2 + d/2 = (k− 1)d < 45 · d2
4 . Then by Lemma 31, we know that with probability at least 1− exp(−dΩ(1)),
ProjV̂ ⊥ij ( w̃i,L1 ⊗ w̃j,L2 ) ≥ poly(1/d, ρ).
Let E2 be the event that this inequality does not hold for some i, j.
Let Sij = span{M̃0, M̃kl : k < l ∧ (k, l) 6= (i, j)}. Now we know when neither bad events E1 or E2 happens, for every pair i < j,
ProjS⊥ijM̃ij ≥ProjV ⊥ij M̃ ′ ij
≥ProjV̂ ⊥ij M̃ ′ ij
= φ̃ij π ProjV̂ ⊥ij ( w̃i,L1 ⊗ w̃j,L2 ) ≥poly(1/τ, 1/d, ρ).
Currently, we have proved that for any i < j, the distance between column M̃ij and the span of other columns is at least inverse polynomial. To use Lemma 21 we just need to give a bound on the norms of these columns. By Lemma 22, we know when E1 does not happen
∀i, ‖w̃i‖ ≤ τ + √ 3ρ2d
2 ,
where τ is the uniform upper bound of the norm of every row vector of W . Let τ̃ = τ + √
3ρ2d 2 , we
know τ̃ = poly(τ, d, ρ).
Thus, we have
‖M̃ij‖ ≤ 1
π
( φ̃ij | cos(φ̃ij)|+ sin(φ̃ij) ) ‖w̃i‖‖w̃j‖ √ d
+ φ̃ij π (‖w̃i‖‖w̃j‖+ ‖w̃j‖‖w̃i‖)
+ sin(φ̃ij)
π (‖w̃i‖‖w̃j‖+ ‖w̃j‖‖w̃i‖)
≤ τ̃ 2
π (π + 1)
√ d+ 2τ̃2 + 2τ̃2
π .
Thus, there exists C = poly(τ, d, ρ), such that ‖M̃ij‖ ≤ C for every i < j. Now applying Lemma 21 immediately gives the result.

heading: C.2 PROOF OF AUXILIARY LEMMAS FOR SECTION C.1
text: We will first prove the characterization for columns in the augmented distinguishing matrix.
Proof of Lemma 19. For simplicity, we start by assuming that every weight vectorwi has unit norm. At the end of the proof we will discuss how to incorporate the norms of wi, wj . Also throughout the proof we will abuse notation to use Mij as its matrix form mat(Mij).
Let Sij be the subspace spanned by wi and wj . Let S⊥ij be the orthogonal subspace of Sij . Let {e(i,j)1 , e (i,j) 2 } be a set of orthonormal basis for Sij such that e (i,j) 1 = wi and 〈e (i,j) 2 , wj〉 > 0. We use matrix Sij ∈ Rd×2 to represent subspace Sij , which matrix has e(i,j)1 and e (i,j) 2 as two columns. Also, let S⊥ij be a d× (d− 2) matrix, whose columns constitute an orthonormal basis of S⊥ij .
Let ProjSij = SijS > ij , and ProjS⊥ij = Id − SijS > ij . Then, we have
Mij = ProjS⊥ijMij + ProjSijMij
= (ProjS⊥ijMij +mijProjSijId) + (ProjSijMij −mijProjSijId) = (ProjS⊥ijMij +mijProjSijId) + (ProjSijMij −mijSijS > ij )
First, we show that ProjS⊥ijMij +mijProjSijId = mijId,
which is equivalent to proving that ProjS⊥ijMij = mijProjS⊥ijId. It’s obvious that the column span of ProjS⊥ijMij belongs to the subspace S ⊥ ij . Actually, the row span of ProjS⊥ijMij also belongs to the subspace S⊥ij . To show this, let’s consider u>(ProjS⊥ijMij)v, where u ∈ S ⊥ ij and v ∈ Sij .
u>(ProjS⊥ijMij)v = u >(Id − SijS>ij )Mijv
= (u> − u>SijS>ij )Mijv = u>Mijv,
where the last equality holds since u ∈ S⊥ij is orthogonal to e (i,j) 1 and e (i,j) 2 . We also know that u>Mijv = u >E [ (w>i x)(w > j x)(xx >)1{w>i xw>j x ≤ 0} ] v
= E [ (w>i x)(w > j x)(u >x)(v>x)1{w>i xw>j x ≤ 0} ]
= E [ (w>i x)(w > j x)(v >x)1{w>i xw>j x ≤ 0} ] E[(u>x)]
= 0,
where the third equality holds because u>x is independent with w>i x,w > j x and v >. Note since u is orthogonal with wi, wj , v, we know for standard Gaussian vector x, random variable u>x is independent with w>i x,w > j x, v >x.
Since the column span and row span of ProjS⊥ijMij both belong to the subspace S ⊥ ij , there must exist a (d − 2) × (d − 2) matrix C, such that ProjS⊥ijMij = S ⊥ ijC(S ⊥ ij ) >. We only need to show this matrix C must be mijId−2. In order to show this, we prove for any u, v ∈ S⊥ij , u>(ProjS⊥ijMij)v = miju >v.
u>(ProjS⊥ijMij)v = u >Mijv = u>E [ (w>i x)(w > j x)(xx >)1{w>i xw>j x ≤ 0} ] v
= E [ (w>i x)(w > j x)(u >x)(v>x)1{w>i xw>j x ≤ 0} ]
= E [ (w>i x)(w > j x)1{w>i xw>j x ≤ 0} ] E[u>xv>x] = miju >E[xx>]v
= miju >v,
where the fourth equality holds because u>x, v>x are independent with w>i x,w > j x.
Thus, we know Mij =(ProjS⊥ijMij +mijProjSijId) + (ProjSijMij −mijSijS > ij )
=mijId + (ProjSijMij −mijSijS > ij ).
Let’s now compute the closed form for mij . Recall that mij := E [ (w>i x)(w > j x)1{w>i xw>j x ≤ 0} ] .
Note, we only need to consider input x within subspace Sij , which subspace has dimension two. Using the polar representation of two-dimensional Gaussian random variables (r is the radius and θ is the angle), we have
mij = 1
2π ∫ ∞ 0 r3 exp(−r 2 2 )dr ∫ π 2
π 2−φij
2 cos(θ) cos(θ + φij)dθ = 1
π (φij cos(φij)− sin(φij)).
Next, we compute the closed form of ProjSijMij . Note ProjSijMij is symmetric, because ProjSijMij = Mij−mijId+mijSijS > ij , andMij , Id and SijS > ij are all symmetric. It’s obvious that the column span of ProjSijMij belongs to subspace Sij . Combined with the fact that ProjSijMij is symmetric, we know the row span of ProjSijMij also belongs to the subspace Sij . Thus, matrix ProjSijMij can be represented as a linear combination of (e (i,j) 1 )(e (i,j) 1 ) >, (e (i,j) 1 )(e (i,j) 2 ) >, (e (i,j) 2 )(e (i,j) 1 ) > and (e(i,j)2 )(e (i,j) 2 ) >, which means
ProjSijMij = c (i,j) 11 (e (i,j) 1 )(e (i,j) 1 ) >+c (i,j) 12 (e (i,j) 1 )(e (i,j) 2 ) >+c (i,j) 21 (e (i,j) 2 )(e (i,j) 1 ) >+c (i,j) 22 (e (i,j) 2 )(e (i,j) 2 ) >,
where c(i,j)11 , c (i,j) 12 , c (i,j) 21 and c (i,j) 22 are four coefficients. Now, we only need to figure out the four coefficients of this linear combination. Similar as the computation for mij , we use polar integration to show that,
c (i,j) 11 =〈ProjSijMij , (e (i,j) 1 )(e (i,j) 1 ) T 〉
= 1
2π ∫ ∞ 0 r5 exp(−r 2 2 )dr ∫ π 2
π 2−φij
( cos3(θ) cos(θ + φij) + cos(θ) cos 3(θ + φij) ) dθ
= 1
4π (12φij cos(φij)− 9 sin(φij)− sin(3φij)),
where the first equality holds because e(i,j)1 is orthogonal with e (i,j) 2 . Similarly, we can show that
c (i,j) 22 =〈ProjSijMij , (e (i,j) 2 )(e (i,j) 2 ) T 〉
= 1
2π ∫ ∞ 0 r5 exp(−r 2 2 )dr ∫ π 2
π 2−φij
( cos(θ) cos(θ + φij) sin 2(θ) + cos(θ) cos(θ + φij) sin 2(θ + φij) ) dθ
= 1
4π (4φij cos(φij)− 7 sin(φij) + sin(3φij)),
and
c (i,j) 21 =〈ProjSijMij , (e (i,j) 2 )(e (i,j) 1 ) T 〉
= 1
2π ∫ ∞ 0 r5 exp(−r 2 2 )dr ∫ π 2
π 2−φij
( − cos2(θ) cos(θ + φij) sin(θ) + cos(θ) cos2(θ + φij) sin(θ + φij) ) dθ
= 1
π (φij sin(φij)− cos(φij) sin2(φij)).
It’s easy to check that c(i,j)12 = c (i,j) 21 . Let M ′ ij be ProjSijMij−mijSijSij . Then, according to above computation, we know
M ′ij =ProjSijMij −mijSijSij
=(c (i,j) 11 −mij)(e (i,j) 1 )(e (i,j) 1 ) > + c (i,j) 12 (e (i,j) 1 )(e (i,j) 2 ) > + c (i,j) 21 (e (i,j) 2 )(e (i,j) 1 ) > + (c (i,j) 22 −mij)(e (i,j) 2 )(e (i,j) 2 ) >
= 1
4π
( 8φij cos(φij)− 5 sin(φij)− sin(3φij) ) (e (i,j) 1 )(e (i,j) 1 ) >
+ 1
4π
( − 3 sin(φij) + sin(3φij) ) (e (i,j) 2 )(e (i,j) 2 ) >
+ 1
π
( φij sin(φij)− cos(φij) sin2(φij) )( (e (i,j) 1 )(e (i,j) 2 ) > + (e (i,j) 2 )(e (i,j) 1 ) >). Since e(i,j)1 = wi and e (i,j) 2 = 1 sin(φij)
wj − cot(φij)wi, we can also express M ′ij as a linear combination of wiw>i , wjw > j , wiw > j and wjw > i .
M ′ij = 1
4π
( 8φij cos(φij)− 5 sin(φij)− sin(3φij) ) wiw > i
+ 1
4π
( − 3 sin(φij) + sin(3φij) ) (
wj sin(φij) − cot(φij)wi)( wj sin(φij) − cot(φij)wi)>
+ 1
π
( φij sin(φij)− cos(φij) sin2(φij) )( wi(
wj sin(φij) − cot(φij)wi)> + ( wj sin(φij) − cot(φij)wi)w>i ) = φij π (wiw > j + wjw > i )− sin(φij) π (wjw > j + wiw > i ).
Thus,
Mij =mijId +M ′ ij
= 1
π (φij cos(φij)− sin(φij))Id + φij π (wiw > j + wjw > i )− sin(φij) π (wjw > j + wiw > i ).
Finally, if the rows wi, wj do not have unit norm, let w̄i = wi/‖wi‖, w̄j = wj/‖wj‖, we know mij = E [ (w>i x)(w > j x)1{w>i xw>j x ≤ 0} ] = ‖wi‖‖wj‖E [ (w̄>i x)(w̄ > j x)1{w>i xw>j x ≤ 0}
] = 1
π (φij cos(φij)− sin(φij))‖wi‖‖wj‖.
Here we used the fact that the indicator variable does not change whether we use wi, wj or w̄i, w̄j . Similarly,
Mij = E [ (w>i x)(w > j x)xx > 1{w>i xw>j x ≤ 0} ] = ‖wi‖‖wj‖E [ (w̄>i x)(w̄ > j x)xx > 1{w>i xw>j x ≤ 0}
] = 1
π (φij cos(φij)− sin(φij))‖wi‖‖wj‖Id + φij π (wiw > j + wjw > i )− sin(φij) π ( ‖wi‖ ‖wj‖ wjw > j + ‖wj‖ ‖wi‖ wiw > i ).
Now we can prove the lemmas used to handle the two obstacles. First we give the stronger leaveone-out distance bound.
Proof of Lemma 21. The smallest singular value of A can be defined as follows:
σmin(A) := min u:‖u‖=1
‖Au‖.
Suppose u∗ ∈ argminu:‖u‖=1‖Au‖. Let u∗i be the coordinate corresponding to the column Ai, for 0 ≤ i ≤ n. We consider two cases here. If |u∗0| ≥ √ 4nC2 4nC2+d , then we have
σmin(A) =‖Au∗‖ = ∥∥u∗0A0 + ∑
1≤i≤n
u∗iAi ∥∥
≥ ∥∥u∗0A0∥∥− ∥∥ ∑
1≤i≤n
u∗iAi ∥∥
≥ √ 4nC2
4nC2 + d
√ d− ( ∑ 1≤i≤n |u∗i |)C
≥ √ 4nC2d
4nC2 + d − √ n
√ d
4nC2 + d C
=
√ nC2d
4nC2 + d ,
where the third inequality uses Cauchy-Schwarz inequality. If |u∗0| < √ 4nC2 4nC2+d , we know ∑ 1≤i≤n |u∗i |2 ≥ d 4nC2+d . Let k ∈ argmax1≤i≤n|u ∗ i |. We know that
|u∗k| ≥ √ d 4n2C2+nd . Thus,
σmin(A) ≥ ∥∥u∗kAk + ∑
i:i 6=k
u∗iAi ∥∥
=|u∗k| ∥∥Ak + ∑
i:i 6=k
u∗i u∗k Ai ∥∥
≥|u∗k|‖(Id − ProjS−k)Ak‖ ≥|u∗k|d′(A)
≥ √
d
4n2C2 + nd d′(A).
Above all, we know that the smallest singular value of A is lower bounded as follows, σmin(A) ≥ min (√ nC2d
4nC2 + d ,
√ d
4n2C2 + nd d′(A)
) .
Next we give the bound on the angle between two perturbed vectors w̃i and w̃j .
Proof of Lemma 22. According to the definition of ρ-perturbation, we know w̃i = wi + ρεi, w̃j = wj+ρεj , where εi, εj are i.i.d. standard Gaussian vectors. First, we show that with high probability, the projection of w̃i on the orthogonal subspace of w̃j is lower bounded. Denote the subspace spanned by w̃j as Sw̃j , and denote the subspace spanned by {w̃j , wi} as Sw̃j∪wi . Thus, we have∥∥ProjS⊥
w̃j
w̃i ∥∥ ≥∥∥ProjS⊥
w̃j∪wi (wi + ρεi) ∥∥ =ρ ∥∥ProjS⊥ w̃j∪wi εi ∥∥
where S⊥w̃j is the orthogonal subspace of Sw̃j .
Fix εj , then S⊥w̃j∪wi is a fixed subspace of R d with dimension d − 2. Let U be a d × (d − 2) matrix, whose columns constitute a set of orthonormal basis for the subspace S⊥w̃j∪wi . Thus, it’s
not hard to check that ProjS⊥ w̃j∪wi
εi d = Uε, where ε ∈ Rd−2 is a standard Gaussian vector. Denote
Y := ‖ProjS⊥ w̃j∪wi
εi‖2 d = ‖Uε‖2 = ‖ε‖2, which is a chi-squared random variable with (d − 2)
degrees of freedom. According to the tail bound for chi-squared random variable, we have
Pr[| 1 d− 2 Y − 1| ≥ t] ≤ 2e −(d−2)t2 8 , ∀t ∈ (0, 1).
Let t = 12 , we know that with probability at least 1− 2 exp( −(d−2) 32 ),
Y ≥ d− 2 2 .
Thus, we have ‖ProjS⊥ w̃j w̃i‖ ≥ ρ‖ProjS⊥ w̃j∪wi
εi‖ ≥ √ ρ2(d−2) 2 . Recall that
‖ProjS⊥ w̃j w̃i‖ = sin(φ̃ij)‖w̃i‖.
We also know
‖w̃i‖ =‖wi + ρεi‖ ≤‖wi‖+ ρ‖εi‖ =τ + ρ‖εi‖,
where the last equality holds since ‖wi‖ ≤ τ . Note ‖εi‖2 is another chi-squared random variable with d degrees of freedom. Similar as above, we can show that with probability at least 1− 2 exp(−d32 ),
‖εi‖2 ≤ 3d
2 .
By union bound, we know with probability at least 1− 2 exp(−d32 )− 2 exp( −(d−2) 32 ), sin(φ̃ij)‖w̃i‖ ≥ √ ρ2(d− 2)
2 , ‖w̃i‖ ≤τ + √ 3ρ2d
2 .
Combined with the fact that φ̃ij ≥ sin(φ̃ij) when φ̃ij ∈ [0, π], we know with probability at least 1− 2 exp(−d32 )− 2 exp( −(d−2) 32 ),
φ̃ij ≥ sin(φ̃ij)
= sin(φ̃ij)‖w̃i‖ ‖w̃i‖
≥
√ ρ2(d−2)
2 τ + √
3ρ2d 2
=
√ ρ2(d− 2)
√ 2τ + √ 3ρ2d
Given W̃ = W + ρE, where E is an i.i.d. Gaussian matrix, by union bound, we know with probability at least 1− exp(−dΩ(1)),
∀i, ‖w̃i‖ ≤ τ + √ 3ρ2d/2
∀i < j, φ̃ij π ≥
√ ρ2(d− 2)
( √ 2τ + √ 3ρ2d)π

heading: C.3 SMOOTHED ANALYSIS FOR GENERAL INPUTS
text: In this section, we show that starting from any well-conditioned weight matrix W , and any symmetric input distribution D, how to perturb the distribution locally to D′ so that the smallest singular value of MD ′ is at least inverse polynomial.
Recall the definition of (Q,λ)-perturbation: we mix the original distribution D with a distribution DQ which is just a Gaussian N (0, QQ>). To create a sample x in D′, with probability 1 − λ we draw a sample from D; otherwise we draw a standard Gaussian n ∼ N (0, Id) and let x = Qn. We will prove Theorem 10 which we restate below:
Theorem 10. Given weight matrix W with ‖wi‖ ≤ τ for each row vector and symmetric input distribution D. Suppose that k ≤ d/7 and σmin(W ) ≥ ρ, after applying (Q,λ)-perturbations to yield perturbed input distribution D′, where Q is a d× d matrix whose entries are i.i.d. Gaussians, we have with probability at least 1− exp(−dΩ(1)) over the randomness of Q,
σmin(M D′) ≥ poly(1/τ, 1/d, ρ, λ).
To prove this, let us first take a look at the structure of augmented distinguishing matrix for these distributions. Let MD, MDQ , MD ′ be the augmented distinguishing matrices for distributions D, DQ and D′ respectively. Since D′ is a mixture of D and DQ, and the augmented distinguishing matrix is defined as expectations over samples, we immediately have
MD ′ = (1− λ)MD + λMDQ .
Our proof will go in two steps. First we will show that σmin(MDQ) is large. Then we will show that even mixing with MD will not significantly reduce the smallest singular value, so σmin(MD ′ ) is also large. In addition to the techniques that we developed in Section C.1, we need two ideas that we call noise domination and subspace decoupling to solve the new challenges here.
Noise Domination First let us focus on σmin(MDQ). This instance has weight W and input distribution N (0, QQ>). Let MWQ be the augmented distinguishing matrix for an instance with weight WQ and input distribution N (0, Id). Our first observation shows that MDQ and MWQ are closely related, and we only need to analyze the smallest singular value of MWQ. The problem now is very similar to what we did in Theorem 9, except that the weight WQ is not an i.i.d. Gaussian matrix. However, we will still be able to use Theorem 9 as a black-box because the amount of noise in WQ is in some sense dominating the noise in a standard Gaussian. More precisely, we use the following simple claim:
Claim 1. Suppose property P holds for N (µ, Id) for any µ, and the property P is convex (in the sense that if P holds for two distributions it also holds for their mixture), then for any covariance matrix Σ Id, we know P also holds for N (µ,Σ).
Intuitively the claim says that if the property holds for a Gaussian distribution with smaller variance regardless of the mean, then it will also hold for a Gaussian distribution with larger variance. The proof is quite simple:
Proof. Let Σ′ = Σ − Id, by assumption we know Σ′ is still a positive semidefinite matrix. Let x ∼ N (µ,Σ), x′ ∼ N (µ,Σ′) and δ ∼ N (0, Id), by property of Gaussians it is easy to see that x d = x′ + δ. Let dx, dx′ and dδ be the density function for x, x′, δ respectively, then we know for any point u dx(u) = Ex′∼N (µ,Σ′)[dδ(u− x′)].
That is, N (µ,Σ) is a mixture of N (x′, I). Since property P is true for all N (x′, I), it is also true for N (µ,Σ).
With this claim we can immediately use the result of Theorem 9 to show σmin(MDQ) is large.
Subspace Decoupling Next we need to consider the mixtureMD ′ . The worry here is that although σmin(M DQ) is large, mixing with D might introduce some cancellations and make σmin(MD ′ ) much smaller. To prove that this cannot happen with high probability, the key observation is that in the first step, to prove σmin(MWQ) is large we have only used the property of WQ. If we let Q̄ be the projection of Q to the orthogonal space of row span of W , then Q̄ is still a Gaussian random matrix even if we condition on the value of WQ! Therefore in the second step we will use the additional randomness in Q̄ to show that the cancellation cannot happen. The idea of partitioning the randomness of Gaussian matrices has been widely used in analysis of approximate message passing algorithms. The actual proof is more involved and we will need to partition the Gaussian matrix Q into more parts in order to handle the special column in the augmented distinguishing matrix .
Now we are ready to give the full proof of Theorem 10
Proof of Theorem 10. Let us first recall the definition of augmented distinguishing matrix: MD ′
is a d2 by (k2 + 1) matrix, where the first k2 columns consist of
MD ′ ij := Ex∼D′ [ (w>i x)(w > j x)(x⊗ x)1{w>i xw>j x ≤ 0} ] ,
and the last column is Ex∼D′ [x ⊗ x]. According to the definition of (Q,λ)-perturbation, if we let DQ be N (0, QQ>), then we have
MD ′ = (1− λ)MD + λMDQ .
In the first step, we will try to analyze MDQ . The first k2 columns of this matrix MDQ can be written as:
M DQ ij = Ex∼DQ [ (w>i x)(w > j x)(x⊗ x)1{w>i xw>j x ≤ 0} ] = En∼N (0,Id) [ (w>i Qn)(w > j Qn)(Qn⊗Qn)1{w>i Qnw>j Qn ≤ 0}
] = Q⊗QEn∼N (0,Id) [ (w>i Qn)(w > j Qn)(n⊗ n)1{w>i Qnw>j Qn ≤ 0}
] for any i < j, and the last column is
Ex∼DQ [x⊗ x] = En∼N (0,Id)[Qn⊗Qn] = Q⊗QEn∼N (0,Id)[n⊗ n].
Except for the factor Q⊗Q, the remainder of these columns are exactly the same as the augmented distinguishing matrix of a network whose first layer weight matrix is WQ and input distribution is N (0, Id). We use MWQ to denote the augmented distinguishing matrix of such a network, then we have
MDQ = Q⊗QMWQ.
Therefore we can first analyze the smallest singular value of MWQ. Let W̃ = WQ. Note that Q is a Gaussian matrix, and W is fixed, so WQ is also a Gaussian random matrix except its entries are not i.i.d. More precisely, there are only correlations within columns of WQ, and for any column of WQ, the covariance matrix is WW>. Since the smallest singular value of W is at least ρ, we know σmin(WW
>) ≥ ρ2. Let the covariance matrix of WQ be ΣWQ ∈ Rkd×kd, which has smallest singular value at least ρ2. Therefore we know ΣWQ ρ2Ikd. It’s not hard to verify that with probability at least 1− exp(−dΩ(1)), the norm of every row of WQ is upper bounded by poly(τ, d). By Claim 1, any convex property that holds for any N (0, ρ2Ikd) perturbation must also hold for ΣWQ
4. Thus, we know with probability at least 1− exp(−dΩ(1)), σmin(M WQ) ≥ poly(1/τ, 1/d, ρ).
To prepare for the next step, we will rewrite MWQ as the product of two matrices. According to the closed form of MWQij in Lemma 19, we know each column of M
WQ can be expressed as a linear combination of w̃i ⊗ w̃j’s and vec(Id). Therefore:
MWQ = [ W̃> ⊗ W̃>, vec(Id) ] R,
4The conclusion of Theorem 9 is clearly convex because it is a probability, and probabilities are linear in terms of mixture of distributions.
where matrix R has dimension (k2 + 1)× (k2 + 1). It’s not hard to verify that
σmin(M WQ) ≤ ∥∥∥[W̃> ⊗ W̃>, vec(Id)]∥∥∥σmin(R). Thus,
σmin(R) ≥ σmin(M WQ)∥∥∥[W̃> ⊗ W̃>, vec(Id)]∥∥∥ . Note that W is a k × d matrix with ‖wi‖ ≤ τ for every row vector, and W̃ = WQ, where Q is an standard Gaussian matrix. Thus, similar as the proof in Lemma 22, we can show that with probability at least 1− exp(−dΩ(1)),∥∥∥[W̃> ⊗ W̃>, vec(Id)]∥∥∥ ≤ ∥∥∥[W̃> ⊗ W̃>, vec(Id)]∥∥∥
F ≤ poly(τ, d).
Thus, we know σmin(R) ≥ poly(1/τ, 1/d, ρ).
Now we will try to perform the second step using the idea of subspace decoupling. Let ProjW be the projection matrix to the row span of W , and let ProjW⊥ = Id − ProjW . Let Q̄ = ProjW⊥Q. Let the columns of U ∈ Rd×(d−k) be a set of orthonormal basis for the orthogonal subspace W⊥. By symmetry of Gaussian, ProjWQ is independent with ProjW⊥Q. Thus, from now on we will condition on ProjWQ, and still treat ProjW⊥Q as a Gaussian random matrix. More precisely, ProjW⊥Q has the same distribution as UP , where P ∈ R(d−k)×d is a standard Gaussian matrix. We further decouple the P part into two subspaces (this is done mostly to handle the special column in augmented distinguishing matrix). Let the rows of V ∈ Rk×d be a set of orthonormal basis for the row span of W̃ = WQ. And let the rows of V ⊥ ∈ R(d−k)×d be a a set of orthonormal basis for the orthogonal subspace W̃⊥. We can then decompose UP into the row span of V and V ⊥ as follows,
UP d = UP1V ⊥ + UP2V,
where P1 ∈ R(d−k)×(d−k) and P2 ∈ R(d−k)×k are two independent standard Gaussian matrices. After this decomposition, we have
Q̄⊗ Q̄ [ W̃> ⊗ W̃>, vec(Id) ] R
d =UP ⊗ UP [ W̃> ⊗ W̃>, vec(Id) ] R
d =(UP1V ⊥ + UP2V )⊗ (UP1V ⊥ + UP2V ) [ W̃> ⊗ W̃>, vec(Id) ] R
=U ⊗ U(P1V ⊥ + P2V )⊗ (P1V ⊥ + P2V ) [ W̃> ⊗ W̃>, vec(Id) ] R
=U ⊗ U [ (P1V ⊥ + P2V )W̃ > ⊗ (P1V ⊥ + P2V )W̃>, vec((P1V ⊥ + P2V )(P1V ⊥ + P2V )>) ] R
=U ⊗ U [ P2V W̃ > ⊗ P2V W̃>, vec(P1P>1 + P2P>2 ) ] R,
where the last equality holds because the row span of V ⊥ is orthogonal to the column span of W̃>.
Now, we go back to matrix MD ′ . Let ProjW⊥⊗W⊥ be ProjW⊥ ⊗ ProjW⊥ . We have,
σmin(M D′) =σmin ( (1− λ)MD + λMDQ ) =σmin ( (1− λ)MD + λQ⊗Q [ W̃> ⊗ W̃>, vec(Id) ] R )
≥σmin ( (1− λ)ProjW⊥⊗W⊥MD + λProjW⊥⊗W⊥Q⊗Q [ W̃> ⊗ W̃>, vec(Id) ] R )
=σmin ( (1− λ)ProjW⊥⊗W⊥MD + λQ̄⊗ Q̄ [ W̃> ⊗ W̃>, vec(Id) ] R )
=σmin ( (1− λ)ProjW⊥⊗W⊥MD + λU ⊗ U [ P2V W̃ > ⊗ P2V W̃>, vec(P1P>1 + P2P>2 ) ] R )
Since R has full row rank, we know that the row span of ProjW⊥⊗W⊥M D belongs to the row span of R. According to the definition of U , it’s also clear that the column span of ProjW⊥⊗W⊥M D belongs to the column span of U ⊗ U . Thus, there exists matrix C ∈ R(d−k)2×(k2+1) such that
ProjW⊥⊗W⊥M D = U ⊗ UCR.
Thus,
σmin(M D′) =σmin ( (1− λ)MD + λMDQ ) ≥σmin ( (1− λ)ProjW⊥⊗W⊥MD + λU ⊗ U [ P2V W̃ > ⊗ P2V W̃>, vec(P1P>1 + P2P>2 ) ] R )
=σmin
( λU ⊗ U (1− λ λ C + [ P2V W̃ > ⊗ P2V W̃>, vec(P1P>1 + P2P>2 ) ]) R ) .
Note that C only depends on U and R, U only depends on W , and R only depends on WQ. With WQ fixed, C is also fixed. Clearly, C is independent with P1 and P2. For convenience, denote
H := 1− λ λ
C + [ P2V W̃ > ⊗ P2V W̃>, vec(P1P>1 + P2P>2 ) ] .
Now, let’s prove that the smallest singular value of matrix H ∈ R(d−k)2×(k2+1) is lower bounded using leave-one-out distance. Let’s first consider its submatrix Ĥ which consists of the first k2 columns of H . Note that within random matrix P2V W̃>, every row are independent with each other. Within each row, the covariance matrix is W̃W̃>. Recall that W̃ is a random matrix whose covariance ΣWQ ρ2Ikd, we can again apply Claim 1 with the property proved in Lemma 37. As a result, with probability at least 1− exp(−dΩ(1)),
σmin(W̃ ) ≥ poly(1/d, ρ).
Thus the covariance matrix of each row of P2V W̃> has smallest singular value at least γ := poly(1/d, ρ).
We can view P2V W̃> as the summation of two independent Gaussian matrix, one of which has covariance matrix γI(d−k)k. For this matrix, we will do something very similar to Theorem 9 in order to lowerbound its smallest singular value.
Claim 2. For a random matrix K ∈ R(d−k)×k that is equal to Ko + E where E is a Gaussian random matrix whose entries have variance γ. If d ≥ 7k, for any subspace SC that is independent of K and has dimension at most k2 + 1, the leave-one-out distance d(ProjS⊥CK ⊗ K) is at least poly(γ, 1/d).
The proof idea is similar as Theorem 9, and we try to apply Lemma 31 to K ⊗K. In the proof we should think of K := P2V W̃>, and denote i-th column of K as Ki. We also think of the space SC as the column span of C.
As we did in Theorem 9, we partition [d − k] into 2 disjoint subsets L1 and L2 of size (d − k)/2. Let Ĥ ′ be the set of rows of Ĥ indexed by L1 × L2.
We fix a column Ĥ ′ij , i 6= j ∈ [k]. Let S = span{Ĥ ′kl : (k, l) 6= (i, j)}. It’s clear that S is correlated with Ĥ ′ij . Let C
′ be the set of rows of C indexed by L1 × L2. Let SC′ be the column span of C ′, which has dimension at most k2 + 1. We define the following subspace that contains S,
Ŝ = SC′ ∪ span{Kj,L1 ⊗ x, x⊗Ki,L2 ,Kl,L1 ⊗ x, x⊗Kl,L2 ∣∣∣x ∈ R(d−k)/2, l /∈ {i, j}}
Therefor by definition S ⊂ Ŝ, and thus Ŝ⊥ ⊂ S⊥, where S⊥ denotes the orthogonal subspace of S. Notice that Ĥ ′ij = Ki,L1 ⊗Kj,L2 + λ1−λC ′ ij is independent with Ŝ, assuming C is fixed. Moreover, Ŝ has dimension at most
(d− k)/2 + (d− k)/2 + (k − 2)(d− k)/2 + (k − 2)(d− k)/2 + k2 + 1 ≤ 4 5
(d− k)2
4 ,
if k ≤ d/7. Then, according to Lemma 31, we know with probability at least 1− exp(−dΩ(1)),∥∥ProjŜ⊥Ĥ ′ij∥∥ ≥ poly(1/d, ρ). For the column Ĥ ′ii, i ∈ [k], we define subspace Ŝ slightly different,
Ŝ = SC′ ∪ span{Kl,L1 ⊗ x, x⊗Kl,L2 ∣∣∣x ∈ R(d−k)/2, l 6= i}.
Here the dimension of Ŝ is also smaller than (d− k)2/5, assuming that k ≤ d/7. We can similarly show that with probability at least 1− exp(−dΩ(1)),∥∥ProjŜ⊥Ĥ ′ii∥∥ ≥ poly(1/d, ρ). Thus, by union bound, we know that the leave-one-out distance of matrix Ĥ ′ is lower bounded by poly(1/d, ρ).
Now, let’s add the additional column vec(P1P>1 + P2P > 2 ) into consideration. For convenience we denote this column by b. We will first prove that the vector b has large norm when projected to the orthogonal subspace of columns in Ĥ , then we will combine this with the fact that σmin(Ĥ) is large to show that σmin(H) is also large (this last step is very similar to Lemma 21).
We know matrix Ĥ only depends on the randomness of P2. Thus, with P2 fixed, all columns in H are fixed except for b. Now, we rely on the randomness in P1 to show that the distance between b and the span of other columns in H is lower bounded. In order to get ride of the correlation within column b, we also need to consider a subset of its rows indexed by L1 × L2, denoted by b′. Let the first column of P1 be p ∈ Rd−k, and the submatrix consisting of other columns be P̂1. Let SĤ′ be the column span of Ĥ ′. Let
ŜĤ′ = SĤ′ ∪ SC′ ∪ span(vec(P̂1P̂ > 1 ) ′, vec(P2P>2 ) ′),
where vec(P̂1P̂>1 ) ′ is the restriction of vec(P̂1P̂>1 ) to the rows indexed by L1 ×L2. The dimension of ŜĤ′ is at most k 2 + k2 + 1 + 2 ≤ (d− k)2/12, assuming that k ≤ d/7. Clearly,
ProjS⊥ Ĥ′ b′ ≥ProjŜ⊥ Ĥ′ b′
=ProjŜ⊥ Ĥ′ pL1 ⊗ pL2 ,
where pL is the restriction of p to rows indexed by L. Note that pL1 and pL2 are two independent standard Gaussian vectors. Thus, according to Lemma 31, we know with probability at least 1 − exp(−dΩ(1)), the distance between b′ and the column span of Ĥ ′ is at least poly(1/d).
Claim 3. For any matrix A ∈ R (d−k)2 4 ×k 2 and a vector v ∈ R (d−k)2 4 , if the leave-one-out distance d(A) ≥ δ, ‖ProjA⊥b‖ ≥ ζ, and ‖v‖ ≤ C1, let B ∈ R (d−k)2 4 ×(k 2+1) be the matrix that is the concatenation of A and v, then the leave-one-out distance d(B) ≥ poly(ζ, δ, 1/C1).
The proof idea is similar as the proof in Lemma 21. In the proof, we should think ofA := Ĥ ′, v := b′ and B := H ′, where H ′ is the subset of rows of H indexed by L1 × L2. We know that d(Ĥ ′) ≥ δ = poly(1/d, ρ), ‖ProjA⊥b‖ ≥ ζ = poly(1/d). It’s not hard to show that with probability at least 1− exp(−dΩ(1)),
‖b′‖ ≤poly(d).
Thus, there exists C1 = poly(d), such that ‖b′‖ ≤ C1. We already proved that the leave-one-out distance of b′ in H ′ is lower bounded. We only need to show the leave-one-out distance for the first k2 columns, which are Ĥ ′ij , i, j ∈ [k].
For any i, j ∈ [k], the leave-one-out distance for Ĥ ′ij within matrix H ′ can be expressed as follows
min ckl,cb
‖Ĥ ′ij + ∑
(k,l) 6=(i,j)
cklĤ ′ kl + cbb ′‖
Let {c∗kl, c∗b} be one set of the optimal solutions to minckl,cb ‖Ĥ ′ij + ∑ (k,l)6=(i,j) Ĥ ′ kl + cbb
′‖. If c∗b = 0, we immediately have
‖Ĥ ′ij + ∑
(k,l) 6=(i,j)
c∗klĤ ′ kl + c ∗ bb ′‖
=‖Ĥ ′ij + ∑
(k,l) 6=(i,j)
c∗klĤ ′ kl‖
≥min ckl ‖Ĥ ′ij + ∑ (k,l) 6=(i,j) cklĤ ′ kl‖
≥δ,
where the last inequality holds because the leave-one-out distance of matrix Ĥ ′ is lower bounded by δ.
If c∗b 6= 0, we need to be more careful. In this case, we have,
‖Ĥ ′ij + ∑
(k,l)6=(i,j)
c∗klĤ ′ kl + c ∗ bb ′‖
=|c∗b | ∥∥∥ 1 c∗b Ĥ ′ij + ∑ (k,l)6=(i,j) c∗kl c∗b Ĥ ′kl + b ′ ∥∥∥
≥|c∗b |ζ,
where the last inequality holds because the distance of b′ to the column span of Ĥ ′ is lower bounded. If |c∗b | ≥ δ2C1 , we have ‖Ĥ ′ ij + ∑ (k,l)6=(i,j) c ∗ klĤ ′ kl + c ∗ bb ′‖ ≥ δζ2C1 .
If |c∗b | < δ2C1 , we have
‖Ĥ ′ij + ∑
(k,l)6=(i,j)
c∗klĤ ′ kl + c ∗ bb ′‖
≥‖Ĥ ′ij + ∑
(k,l)6=(i,j)
c∗klĤ ′ kl‖ − ‖c∗bb′‖
≥min ckl ‖Ĥ ′ij + ∑ (k,l)6=(i,j) cklĤ ′ kl‖ − |c∗b |‖b′‖
≥δ − δ 2C1 C1
=δ/2.
Thus, the leave-one-out distance of H ′ is lower bounded by poly(ζ, δ, 1/C1). Recall that with probability at least 1 − exp(−dΩ(1)), we have δ = poly(1/d, ρ), ζ = poly(1/d), C1 = poly(d). Thus, we have d(H ′) ≥ poly(1/d, ρ). According to Lemma 20, we have
σmin(H ′) ≥ 1√
k2 + 1 d(H ′)
≥poly(1/d, ρ).
Finally, we put everything together. SinceH ′ is a full column rank matrix, we know that σmin(H) ≥ σmin(H ′) ≥ poly(1/d, ρ). By union bound, we know with probability at least 1− exp(−dΩ(1)),
σmin(H) ≥poly(1/d, ρ) σmin(R) ≥poly(1/τ, 1/d, ρ).
Since U is an orthonormal matrix, we know σmin(U ⊗U) = 1. According to Eq. 19, we know with probability at least 1− exp(−dΩ(1)),
σmin(M D′) ≥σmin(λU ⊗ UHR) ≥λσmin(U ⊗ U)σmin(H)σmin(R) ≥poly(1/τ, 1/d, ρ, λ)
where the second inequality holds since all of U ⊗ U , H and R have full column rank.

heading: D TOOLS
text: In this section, we collect some known results on matrix perturbations and concentration bounds. Basically, we used matrix concentration bounds to do the robust analysis and used matrix perturbation bounds to do the smoothed analysis. We also proved several corollaries that are useful in our setting.

heading: D.1 MATRIX CONCENTRATION BOUNDS
text: Matrix concentration bounds tell us that with enough number of independent samples, the empirical mean of a random matrix can converge to the mean of this matrix. Lemma 23 (Matrix Bernstein; Theorem 1.6 in Tropp (2012)). Consider a finite sequence {Zk} of independent, random matrices with dimension d1 × d2. Assume that each random matrix satisfies
E[Zk] = 0 and ‖Zk‖ ≤ R almost surely. Define
σ2 := max {∥∥∑
k
E[ZkZ∗k ] ∥∥,∥∥∑
k
E[Z∗kZk] ∥∥}.

heading: Then, for all t ≥ 0,
text: Pr {∥∥∑
k
Zk ∥∥ ≥ t} ≤ (d1 + d2) exp( −t2/2
σ2 +Rt/3
) .
As a corollary, we have: Lemma 24. Consider a finite sequence {Z1, Z2 · · ·Zm} of independent, random matrices with dimension d1 × d2. Assume that each random matrix satisfies
‖Zk‖ ≤ R, 1 ≤ k ≤ m.

heading: Then, for all t ≥ 0,
text: Pr {∥∥ m∑
k=1
(Zk − E[Zk]) ∥∥ ≥ t} ≤ (d1 + d2) exp( −t2/2
4mR2 + (2Rt)/3
) .
Proof. For each k, let Z ′k = Zk − E[Zk] be the new random matrices. It’s clear that E[Z ′k] = 0 and ‖Z ′k‖ ≤ 2R. For the variance,
σ2 ≤ m∑ k=1 ‖E[Z ′kZ ′∗k ]‖ (19)
≤ m∑ k=1 4R2 (20)
=4mR2. (21) (22)
Thus, according to Lemma 23, we have Pr {∥∥ m∑
k=1
(Zk − E[Zk]) ∥∥ ≥ t} = Pr{∥∥ m∑
k=1
(Z ′k) ∥∥ ≥ t} ≤ (d1 + d2) exp( −t2/2
4mR2 + (2Rt)/3
) .

heading: D.2 MATRIX PERTURBATION BOUNDS
text: Perturbation Bound for Singular Vectors For singular vectors, the perturbation is bounded by Wedin’s Theorem. Lemma 25 (Wedin’s theorem; Theorem 4.1, p.260 in Stewart & Sun (1990).). Given matrices A,E ∈ Rm×n with m ≥ n. Let A have the singular value decomposition
A = [U1, U2, U3] [ Σ1 0 0 Σ2 0 0 ] [V1, V2] >
Let Â = A + E, with analogous singular value decomposition. Let Φ be the matrix of canonical angles between the column span of U1 and that of Û1, and Θ be the matrix of canonical angles between the column span of V1 and that of V̂1. Suppose that there exists a δ such that
min i,j |[Σ1]i,i − [Σ2]j,j | > δ and min i |[Σ1]i,i| > δ,
then
‖ sin(Φ)‖2 + ‖ sin(Θ)‖2 ≤ 2‖E‖ 2
δ2 .
In order to show the robustness of least k right singular vectors of T , we combine Wedin’s theorem with the following Lemma. Lemma 26 (Theorem 4.5, p.92 in Stewart & Sun (1990).). Let Φ be the matrix of canonical angles between the column span of U and that of Û , then
‖ProjÛ − ProjU‖ = ‖ sin(Φ)‖.
The exact lemma used in our proof is the following corollary in Ge et al. (2015). Lemma 27 (Lemma G.5 in Ge et al. (2015)). Given matrix A,E ∈ Rm×n with m ≥ n. Suppose that the A has rank k. Let S and Ŝ be the subspaces spanned by the first k right singular vectors of A and Â = A+ E, respectively. Then, we have:
‖ProjŜ⊥ − ProjS⊥‖ ≤ √
2‖E‖F σk(A) .
Perturbation Bound for pseudo-inverse With a lowerbound on σmin(A), we can get bounds for the perturbation of pseudo-inverse. Lemma 28 (Theorem 3.4 in Stewart (1977)). Consider the perturbation of a matrix A ∈ Rm×n : B = A+ E. Assume that rank(A) = rank(B) = n, then
‖B† −A†‖ ≤ √ 2‖A†‖‖B†‖‖E‖.
The following corollary is particularly useful for us. Lemma 29 (Lemma G.8 in Ge et al. (2015)). Consider the perturbation of a matrix A ∈ Rm×n : B = A+ E where ‖E‖ ≤ σmin(A)/2. Assume that rank(A) = rank(B) = n, then
‖B† −A†‖ ≤ 2 √ 2‖E‖/σmin(A)2.
Perturbation Bound for Tensor To lowerbound the leave-one-out distance in augmented distinguishing matrix , we use the following Lemma as the main tool. Lemma 30 (Theorem 3.6 in Bhaskara et al. (2014)). For any constant δ ∈ (0, 1), given any subspace V of dimension δdl in Rdl , there exist vectors v1, v2, · · · , vr in V with unit norm, such that for random (ρ/ √ d)-perturbations x̃(1), x̃(2), · · · , x̃(l) ∈ Rd of any vector x(1), x(2), · · · , x(l) ∈ Rd, we know with probability at least 1− exp(−δd1/(2l)l),
∃j ∈ [r], 〈vj , x̃(1) ⊗ x̃(2) ⊗ · · · ⊗ x̃(l)〉 ≥ ρl( 1 d )3 l .
For second-order tensor, we have the following corollary.
Lemma 31. For any constant δ ∈ (0, 1), given any subspace V of dimension δd2 in Rd2 , there exist vectors v1, v2, · · · , vr in V with unit norm, such that for random (ρ/ √ d)-perturbations x̃(1), x̃(2) ∈ Rd of any vector x(1), x(2) ∈ Rd, we know with probability at least 1− exp(−δd1/16),
∃j ∈ [r], 〈vj , x̃(1) ⊗ x̃(2)〉 ≥ ρ2( 1
d )9.
Perturbation Bound for Eigendecomposition Here, We restate some generic results from Bhaskara et al. (2014) on the stability of a matrix’s eigendecomposition under perturbation. Let M and M̂ be two n× n mtrices such that M = UDU−1 and M̂ = M(I + E) + F . Definition 3 (Definition A.1 in Bhaskara et al. (2014)). Let sep(D) = mini6=j |Dii −Djj |.
The following Lemma guarantees that the eigenvalues of M̂ are distinct if the perturbation are not too large. Lemma 32 (Lemma A.2 in Bhaskara et al. (2014)). If κ(U)(‖ME‖ + ‖F‖) < sep(D)/2n, then the eigenvalues of M̂ are distinct and diagonalizable.
The following Lemma further upperbound the difference between corresponding eigenvectors. Lemma 33 (Lemma A.3 in Bhaskara et al. (2014)). Let u1, ..., un and û1, ..., ûn respectively be the eigenvectors of M and M̂ , ordered by their corresponding eigenvalues. If κ(U)(‖ME‖+ ‖F‖) < sep(D)/2n, then for all i we have ‖ûi − ui‖ ≤ 3σmax(E)σmax(D)+σmax(F )σmin(U)sep(D) .
In the setting of simultaneous diagonalization, let Na = Ta + Ea and Nb = Tb + Eb, we have NaN −1 b = TaT −1 b (I + F ) +G, where F = −Eb(I +T−1b Eb)−1T −1 b and G = EaN −1 b . The following lemma bound the maximum singular value of perturbation matrix F and G.
Lemma 34 (Claim A.5 in Bhaskara et al. (2014)). σmax(F ) ≤ σmax(Eb)σmin(Tb)−σmax(Eb) and σmax(G) ≤ σmax(Ea) σmin(Nb)
Alignment of Subspace Basis. Due to the rotation issue, we cannot conclude that ‖S − Ŝ‖ is small even we know ‖SS>− ŜŜ>‖ is bounded. The following Lemma shows that after appropriate alignment, S is indeed close to Ŝ.
Lemma 35 (Lemma 6 in Ge et al. (2017a)). Given matrices S, Ŝ ∈ Rd×r, we have
min Z>Z=ZZ>=Ir
‖Ŝ − SZ‖2F ≤ ‖SS> − ŜŜ>‖2F
2( √ 2− 1)σr(SS>)

heading: D.3 SMALLEST SINGULAR VALUE OF RANDOM MATRICES
text: For a random rectangular matrix where each element is an inpdependent Gaussian variable, Rudelson & Vershynin (2009) gives the following result: Lemma 36 (Theorem 1.1 in Rudelson & Vershynin (2009)). Let A ∈ Rm×n and suppose that m ≥ n. Assume that the entries of A are independent standard Gaussian variable, then for every > 0, with probability at least 1 − (C )m−n+1 + e−C′n, where C,C ′ are two absolute constants, we have:
σn(A) ≥ ( √ m− √ n− 1).
However, in our setting, we are more interested in fixed matrices perturbed by Gaussian variables. The smallest singular value of these “perturbed rectangular matrices” can be bounded as follows. Lemma 37 (Lemma G.16 in Ge et al. (2015)). Let A ∈ Rm×n and suppose that m ≥ 3n. If all the entries of A are independently ρ-perturbed to yield Ã, then for any > 0, with probability at least 1− (C )0.25m, for some absolute constant C, the smallest singular value of Ã is bounded below by:
σn(Ã) ≥ ρ √ m.

heading: D.4 ANTI-CONCENTRATION
text: We use the anti-concentration property for Gaussian random variables in our proof of Lemma 17. Lemma 38 (Anti-concentration in Carbery & Wright (2001)). Let x ∈ Rn be a Gaussian variable x ∈ N(0, I), for any polynomial p(x) of degree d, there exists a constant κ such that
Pr [ |p(x)| ≤ √ V ar[p(x)] ] ≤ κ 1/d.

Review summaries:
This paper pushes forward our understanding of learning neural networks. The authors show that they can learn a two-layer (one hidden layer) NN, under the assumption that the input distribution is symmetric. The authors convincingly argue that this is not an excessive limitation, particularly in view of the fact that this is intended to be a theoretical contribution. Specifically, the main result of the paper relies on the concept of smoothed analysis. It states that give data generated from a network, the input distribution can be perturbed so that their algorithm then returns an epsilon solution. 

This paper studies the problem of learning the parameters of a two-layer (or one-hidden layer) ReLU network $y=A\sigma(Wx)$, under the assumption that the distribution of $x$ is symmetric. The main technique here is the "pure neuron detector", which is a high-order moment function of a vector. It can be proved that the pure neuron detector is zero if and only if the vector is equal to the row vector of A^{-1}. Hence, we can "purify" the two layer neural network into independent one layer neural networks, and solve the problem easily.