Title: LEARNING TWO-LAYER NEURAL NETWORKS WITH SYMMETRIC INPUTS

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
where A, 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.
1 There 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.

1. INTRODUCTION:

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-layer 1 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 distribution D, 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 distribution D, or both. Many works have worked with a simple input distribution (such as Gaussians) and try to learn more and more 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σ(W x) + ξ.

(1)

Here x ∈ R d is the input, W ∈ R k×d and A ∈ R k×k are two weight matrices2 . The function σ is the standard ReLU activation function σ(x) = max{x, 0} applied entry-wise to the vector W x, 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 distribution D, our only requirement is that D is symmetric. That is, for any x ∈ R d , 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.

1.1. OUR RESULTS:

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 Â, Ŵ of the same size that is effectively the same as the ground-truth network: for any input x, Âσ( Ŵ x) = Aσ(W x).

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 matrices A, 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 Â, Ŵ of the same size that computes an -approximation function to the ground-truth network: for any input x, Âσ( Ŵ x) -Aσ(W x) 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 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.

1.2. RELATED WORK:

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.

2. PRELIMINARIES:

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.

2.1. NETWORK MODEL:

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 ∈ R d . 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 ∈ R k×d and A ∈ R k×k . The output y ∈ R k is generated as follows:

EQUATION

where σ(•) is the element-wise ReLU function and ξ ∈ R k is zero-mean random noise, which is independent with input x. Let the value of hidden units be h ∈ R k , which is equal to σ(W x). Denote i-th row of matrix W as w i (i = 1, 2, ..., k). Also, let i-th column of matrix A be a i (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.

σ σ σ h = σ(W x) ∈ R k x ∼ D y = Ah + ξ A ∈ R k×k W ∈ R k×d w 1 w 2 w k Figure 1: Network model.

2.2. NOTATIONS:

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 R n , we use •, • to denote the inner product of two vectors, and use • to denote the Euclidean norm. We use e i to denote the i-th standard basis vector. For a matrix A ∈ R m×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 I n 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 ∈ R n×d , whose columns form an orthonormal basis for subspace S. The projection matrix onto the subspace S is denoted by Proj S = SS , and the projection matrix onto the orthogonal subspace of S is denoted by Proj S ⊥ = I n -SS .

For matrix A ∈ R m1×n1 , C ∈ R m2×n2 , let the Kronecker product of A and C be A ⊗ C ∈ R m1m2×n1n2 , which is defined as (A ⊗ C) (i1,i2),(j2,j2) = A i1,i2 C j1,j2

. For a vector x ∈ R d , the Kronecker product x ⊗ x has dimension d 2 . We denote the p-fold Kronecker product of x as x ⊗p , which has dimension d p .

We often need to convert between vectors and matrices. For a matrix A ∈ R m×n , let vec(A) ∈ R mn be the vector obtained by stacking all the columns of A. For a vector a ∈ R m 2 , let mat(x) ∈ R m×m denote the inverse mapping such that vec(mat(a)) = a. Let R k×k sym be the space of all k × k symmetric matrices, which has dimension k 2 + k. For convenience, we denote k 2 = k 2 . For a symmetric matrix B ∈ R k×k sym , we denote vec * (B) ∈ R k2+k as the vector obtained by stacking all the upper triangular entries (including diagonal entries) of B. Note that vec(B) still has dimension k 2 . For a vector b ∈ R k2+k , let mat * (b) ∈ R k×k sym denote the inverse mapping of vec * (•) such that vec * (mat * (B)) = b.

2.3. DISTINGUISHING MATRIX:

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 N D ∈ R d 2 ×k2 is a matrix whose columns are indexed by ij where 1 ≤ i < j ≤ k, and

N D ij = E x∼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 d 2 × (k 2 + 1) matrix whose first k 2 columns are exactly the same as distinguishing matrix N , and the last column (indexed by 0) is defined as

M D 0 = E x∼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 N D ij looks at the expectation over samples that have opposite signs for weights w i , w j (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 distribution D, there exists simple perturbations D that are arbitrarily close to D such that σ min (M D ) is lowerbounded.

3. OUR ALGORITHM:

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 ξ.

3.1. WARM-UP: LEARNING SINGLE-LAYER NETWORKS:

We will first give a simple algorithm for learning a single-layer neural network. More precisely, suppose we are given samples (x 1 , y 1 ), ..., (x n , y n ) where x i ∼ D comes from a symmetric distribution, and the output y i is computed by

EQUATION

Here ξ i 's are i.i.d. noises that satisfy E[ξ i ] = 0. Noise ξ i is also assumed to be independent with input x i . 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 vector w 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 (x 1 , y 1 ), ..., (x n , y n ) generated according to Equation (3).

Output: Estimate of weight vector w.

1: Estimate v = 1 n n i=1 y i x i . 2: Estimate C = 1 n n i=1 x i x i 3: return 2C -1 v.

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 ]) -1 E[yx]

. This is summarized in Algorithm 1.

3.2. LEARNING TWO-LAYER NETWORKS:

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

EQUATION

Let u ∈ R k be a vector and consider u y, it is clear that u y = (u A)σ(W x). Let z i be the normalized version i-th row of A -1 , then we know z i has the property that z i A = λ i e i where λ i > 0 is a constant and e i is a basis vector.

The key observation here is that if u = z i , then u A = λ i e i . As a result, u y = λ i e i σ(W x) = σ(λ i w i x) is the output of a single-layer neural network with weight equal to λ i w i . If we know all the vectors {z 1 , ..., z k }, the input/output pairs (x, z i y) correspond to single-layer networks with weight vectors {λ i w i }. 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 = λ i e i , we say that u y is a pure neuron. Next we will design an algorithm that can find all vectors {z i }'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 ŷ = σ(w x), then E[ŷ 2 ] = 1 2 w E[xx ]w, and E[ŷx] = 1 2 E[xx ]w. As a result we have E[ŷ 2 ] = 2E[ŷx ]E[xx ] -1 E[ŷx].

As before, the ReLU activation does not appear because of the symmetric input distribution. For ŷ = u y, we can estimate all of these moments (E[ŷ 2 ], E[ŷ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 ŷ = k i=1 c i σ(w i x), then we have

2E[ŷx ]E[xx ] -1 E[ŷx] -E[ŷ 2 ] = 1≤i<j≤k c i c j E (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 ŷ = σ(w x), then

E[ŷ 2 (x ⊗ x)] = 2E ŷ • (E[ŷx ]E[xx ] -1 x) • (x ⊗ x) . Moreover, if ŷ = k i=1 c i σ(w i x) where c ∈ R k is a k-dimensional vector, we have 2E ŷ • (E[ŷx ]E[xx ] -1 x) • (x ⊗ x) -E[ŷ 2 (x ⊗ x)] = 1≤i<j≤k c i c j N ij .

Here N ij 's are columns of the distinguishing matrix defined in Definition 1.

The important observation here is that there are

k 2 = k 2 extra terms in 2E ŷ • (E[ŷx ]E[xx ] -1 x) • (x ⊗ x) -E[ŷ 2 (x ⊗ x)

] that are multiples of N ij , which are d 2 (or d+1 2 considering their symmetry) dimensional objects. When the distinguishing matrix is full rank, we know its columns N ij are linearly independent. In that case, if the sum of the extra terms is 0, then the coefficient in front of each N ij must also be 0. The coefficients are c i c j which will be non-zero if and only if both c i , c j are non-zero, therefore to make all the coefficients 0 at most one of {c i } 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 ] -1 x) • (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 ±z i .

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 u i u j as a separate variable. Now the number of variables is k 2 +k = k 2 +k and f is linear in this space. In other words, there exists a matrix T ∈ R d 2 ×(k2+k) such that T vec * (uu ) = f (u). Clearly, if u y is a pure neuron, then T vec * (uu ) = f (u) = 0. That is, {vec * (z i z 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 R d 2 ×(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 * (z i z 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 * (z i z i )}.

Finding z i 's from span of z i z i 's In order to reduce the problem to a single-layer problem, the final step is to find z i 's from span of z i z 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 z i 's, which means Z = diag(λ)A -1 . Let X = Z D X Z and Y = Z D Y Z be two random elements in the span of z i z i , where D X and D Y 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 D X D -1 Y (Z ) -1 , since z i is a column of Z , we know XY -1 z i = Z D X D -1 Y (Z ) -1 z i = Z D X D -1 Y e i = D X (i, i) D Y (i, i) Z e i = D X (i, i) D Y (i, i) z i .

That is, z i is an eigenvector of XY -1 ! The matrix XY -1 can have at most k eigenvectors and there are k z i 's, therefore the z i 's are the only eigenvectors of XY -1 .

Lemma 5. Given the span of z i z i 's, let X, Y be two random matrices in this span, with probability 1 the z i 's are the only eigenvectors of XY -1 .

Using this procedure we can find all the z i 's (up to permutations and sign flip). Without loss of generality we assume z i A = λ i e i . The only remaining problem is that λ i might be negative. 

(u) = 2 Ê (u y)•( Ê[(u y)x ] Ê[xx ] -1 x)•(x⊗x) -Ê (u y) 2 (x⊗x)

where each entry is expressed as a degree-2 polynomial over u.

4: Construct matrix T ∈ R d 2 ×(k2+k) such that, T vec * (uu ) = f (u). 5: Compute the k-least right singular vectors of T , denoted by vec * (U 1 ), vec * (U 2 ), • • • , vec * (U k ).

Let S be a k 2 + k by k matrix, where the i-th column of S is vector vec * (U i ). {Here span S is equal to span of {vec * (z i z i )}.} 6: {Finding z i 's from span} 7: Let X = mat * (Sζ 1 ), Y = mat * (Sζ 2 ), where ζ 1 and ζ 2 are two independently sampled kdimensional standard Gaussian vectors. 8: Let z 1 , ..., z k be eigenvectors of XY -1 . 9: For each z i , if Ê[z i y] < 0 let z i ← -z i . {Here z i 's are normalized rows of A -1 .} 10: {Reduce to 1-Layer Problem} 11: For each z i , let v i be the output of Algorithm 1 with input (x 1 , z i y 1 ), ..., (x n , z i y n ). 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. Suppose A, W, E[xx ] and the distinguishing matrix N 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 * (z i z 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 = λ i e i . In Step 9, we use the fact that E

[z i y] = λ i E[σ(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) = σ(λ i w i x) (again by property of ReLU function σ), by the design of Algorithm 1 we know v i = λ i w i . We also know that Z = diag(λ)A -1 , therefore

Z -1 = Adiag(λ) -1 . Notice that Z -1 σ(V x) = Adiag(λ) -1 σ(diag(λ)W x) = Aσ(W x

). These two scaling factors cancel each other, so the two networks compute the same function. 

4. EXPERIMENTS:

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.

4.1. SAMPLE EFFICIENCY:

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.

4.2. ROBUSTNESS TO NOISE:

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.

4.3. ROBUSTNESS TO CONDITION NUMBER:

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 recovers A and W with sufficient accuracy to keep the overall mean square error low.

5. CONCLUSION:

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.

A DETAILS OF EXACT ANALYSIS:

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 = σ(W x) + ξ) 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 ξ.

A.1 MISSING PROOFS FOR SECTION 3:

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 ∈ R d 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) p x ⊗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) p x ⊗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) p x ⊗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[ŷ 2 ] = 1 2 w E[xx ]w. Similarly, we also know E[ŷx ] = 1 2 w E[xx ] and E[ŷx] = 1 2 E[xx ]w. Thus, we have E[ŷ 2 ] = 2E[ŷx ]E[xx ] -1 E[ŷ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 ŷ = k i=1 c i σ(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 ] -1 x) • (x ⊗ x) first. For E[(u y)x ]E[xx ] -1 , we have E (u y)x E xx -1 =E (u Aσ(W x))x E xx -1 = 1 2 u AW E xx E xx -1 = 1 2 u AW For any vector q ∈ R d , consider 2E (u y) • (q x) • (x ⊗ x) . We have 2E (u y) • (q x) • (x ⊗ x) =2E[(u Aσ(W x))(q x)(x ⊗ x)] =E (u AW x)(q x)(x ⊗ x) Let q = E (u y)x E xx -1 = 1 2 u AW , we have 2E (u y) • (E[(u y)x ]E[xx ] -1 x) • (x ⊗ x) = 1 2 E (u AW x) 2 (x ⊗ x) = 1 2 E A u, W x 2 (x ⊗ x) = 1 2 1≤i≤k (A u) 2 i E (w i x) 2 (x ⊗ x) + 1≤i<j≤k (A u) i (A u) j E (w i x)(w j x)(x ⊗ x) = 1 2 1≤i≤k (A uu A) ii E (w i x) 2 (x ⊗ x) + 1≤i<j≤k (A uu A) ij E (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

EQUATION

Now, we subtract ( 5) by ( 6) to obtain

EQUATION

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 ∈ R d , we have

EQUATION

) 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

1 2 E (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

1 2 E (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 * (z

i z i )}. First, we recall that f (u) = 2E (u y) • (E[(u y)x ]E[xx ] -1 x) • (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 ] -1 x)•(x⊗x) -E (u y) 2 (x⊗x) = 1≤i<j≤k (A uu A) ij N ij .

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 u i u j as a single variable) similar to the FOOBI algorithm. Thus the number of variable is k

EQUATION

Now, we prove the Lemma 4 which shows the null space of T is exactly the span of {vec * (z i z 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 T vec * (U ) = 0. Note that the RHS of (10) equals to 0 if and only if A U A is a diagonal matrix since the distinguishing matrix N is full column rank and A U A is symmetric. Thus vec * (U ) belongs to the span of {vec * (z i z i )} since U = Z DZ for some diagonal matrix D.

2. For any vector vec * (U ) belonging to the span of {vec * (z

i z i )}, U is a linear combination of z i z i 's. Furthermore, T vec * (U ) is a linear combination of T vec * (z i z i ).

Note that A z i only has one non-zero entry due to the definition of z i , for any i ∈ [k]. Thus all coefficients in the RHS of (10) are 0. We get T vec * (U ) = 0.

Finding z i 's: Now, we prove the final Lemma 5 which finds all z i 's from the span of {vec * (z i z i )} by using simultaneous diagonalization. Given all z i '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 w i 's.

Proof of Lemma 5. As we discussed before this lemma, we have

XY -1 = Z D X D -1 Y (Z ) -1 .

According to the following Lemma 8 (i.e., all diagonal elements of D x D -1 y are non-zero and distinct), the matrix XY -1 have k eigenvectors and there are k z i 's, therefore z i 's are the only eigenvectors of XY -1 . Lemma 8. With probability 1, all diagonal elements of D X and D Y are non-zero and all diagonal elements of D X D -1 Y are distinct, where X = Z D X Z and Y = Z D Y Z are defined in Line 7 of Algorithm 2.

Proof. First, we know there exist diagonal matrices

{D i : i ∈ [k]} such that mat * (vec * (U i )) = Z D i Z for all i ∈ [k], where {vec * (U i ) : i ∈ [k]

} are the k-least right singular vectors of T (see Line 5 of Algorithm 2). Then, let the vector d i ∈ R k be the diagonal elements of D i , for all i ∈ k. Let matrix Q ∈ R k×k be a matrix where its i-th column is d i . Then D X = diag(Qζ 1 ) and D Y = diag(Qζ 2 ), where ζ 1 and ζ 2 are two random k-dimensional standard Gaussian vectors (see Line 7 of Algorithm 2).

Since {vec * (U i ) : i ∈ [k]} are singular vectors of T , we know d 1 , d 2 , • • • , d k 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 D X first. In order for i-th diagonal element of D X to be zero, we need q i ζ 1 = 0. Since q i is not a zero vector, we know the solution space of q i ζ 1 = 0 is a lower-dimension manifold in R k , 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 D X or D Y happens with probability zero.

If i-th and j-th diagonal elements of

D X D -1 Y 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 R 2k space, with measure zero. Since finite union of zero-measure sets still has measure zero, the event that duplicated diagonal elements exist in D X D -1 Y happens with probability zero.

A.2 NOISY CASE Now, we discuss how to handle the noisy case (i.e. y = σ(W x) + ξ). 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 * (z i z i ) (Lemma 4). It does not affect the last two steps, i.e., finding z i 's from the span (Lemma 5) and learning the reduced single-layer network. Because 

i )} n/2 i=1 ) to estimate empirical moments Ê[xx ], Ê[yx ], Ê[yy ], Ê[y ⊗ x ⊗3 ] and Ê[y ⊗ y ⊗ (x ⊗ x)]. 3: Compute the vector f (u) = 2 Ê (u y) • ( Ê[(u y)x ] Ê[xx ] -1 x) • (x ⊗ x) -Ê (u y) 2 (x ⊗ x) + Ê (u y) 2 -2 Ê (u y)x Ê xx -1 Ê (u y)x Ê[x ⊗ x]

where each entry is expressed as a degree-2 polynomial over u.

4: Construct matrix T ∈ R d 2 ×(k2+k) such that, T vec * (uu ) = f (u). 5: Compute the k-least right singular vectors of T , denoted by vec * (U 1 ), vec * (U 2 ), • • • , vec * (U k ).

Let S be a k 2 + k by k matrix, where the i-th column of S is vector vec * (U i ). {Here span S is equal to span of {vec * (z i z i )}.} 6: {Finding z i 's from span} 7: Let X = mat * (Sζ 1 ), Y = mat * (Sζ 2 ), where ζ 1 and ζ 2 are two independently sampled kdimensional standard Gaussian vectors. 8: Let z 1 , ..., z k be eigenvectors of XY -1 . 9: For each z i , use the second half of samples

{(x i , y i )} n i=n/2+1 to estimate Ê[z i y]. If Ê[z i y] < 0 let z i ← -z i .

{Here z i 's are normalized rows of A -1 .} 10: {Reduce to 1-Layer Problem} 11: For each z i , let v i be the output of Algorithm 1 with input {(x j , z i y j )} n j=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σ(W x) + ξ, for any u ∈ R k , we have

EQUATION

where N ij 's are columns of the distinguishing matrix (Definition 1), and

m ij := 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 R d 2 ×(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 * (z i z 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

EQUATION

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.

EQUATION

where ( 15) uses (6).

For the third term, we have

EQUATION

)

where the third equality holds due to Lemma 6 and Lemma 7, and ( 16) uses the definition of m ij . 3: Run Algorithm 3 on samples {(x i , P y i )} n i=n/2 . Let the output of Algorithm 3 be V, Z -1 . 4: return V , P Z -1 .

Finally, we combine these three terms (13-16) as follows:

EQUATION

where ( 17) uses (9) (same as ( 7)).

A.3 EXTENSION TO NON-SQUARE A:

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 be A ∈ R l×k and W ∈ R k×d . The samples are generated by y = Aσ(W x)+ξ, where the noise ξ is independent with input x. We have

E[yx ] = 1 2 AW E[xx ].

Since both W 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σ(W x). Thus, we have

P Z -1 σ(V x) = P P Aσ(W x) = Aσ(W x),

where the second equality holds since P P is just the projection matrix to the column span of A.

B ROBUSTNESS OF MAIN ALGORITHM:

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 , Ẑ-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 ≤ Γ, ξ ≤ P 2 , A ≤ P 1 , 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 Γ, P 1 , P 2 , d, 1/ , 1/γ, 1/α, 1/β, 1/δ number of i.i.d. samples, let the output of Algorithm 3 be V , Ẑ-1 , we know with probability at least 1 -δ,

Aσ(W x) -Ẑ-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 * (z i z i )'s; finding z i 's from the span of vec * (z i z 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 * (z i z 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 * (z i z 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 Ŝ ∈ R (k2+k)×k for empirical estimate T . Then for any ≤ γ/2, for any δ < 1, 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 -δ,

SS -Ŝ Ŝ ≤ √ 2 αβ 2 .

The proof of the above lemma is in Section B.1. Basically, we need to lowerbound the spectral gap (k 2 -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 z i 's from the span of vec * (z i z 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 z i 's robustly. The detailed proof is in Section B.2. 

Lemma 12. Suppose that SS -Ŝ Ŝ F ≤ , A ≤ P 1 , ξ ≤ P 2 , σ min (E[xx ]) ≥ γ, σ min (A) ≥ β. Let X = mat * ( Ŝζ 1 ), Ŷ = mat * ( Ŝζ 2 ),

(i) ∈ [k] such that ẑi -z π[i] ≤ poly(d, P 1 , 1/β, , 1/δ), for any 1 ≤ i ≤ k.

Finally, given ẑ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 ẑi (learned using the other set of samples). However, since ẑi is not equal to z i , this introduces an additional error term (ẑ i -z i ) 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 ẑi -z i is small enough. Lemma 13. Assume that x ≤ Γ, A ≤ P 1 , ξ ≤ P 2 and σ min (E[xx ]) ≥ γ. Suppose that for each 1 ≤ i ≤ k, ẑi -z i ≤ τ . 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 -δ,

v i -vi ≤ 2τ (ΓP 1 √ k + P 2 )Γ γ + 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 * (z i z j )}'S

We first prove that the step of finding the span of {vec * (z i z j )} is robust. The main idea is based on standard matrix perturbation bounds (see Section D.2). We first give a lowerbound on the k 2 -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 k 2 and the k 2 -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 ≤ P 1 , ξ ≤ P 2 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 Ŝ ∈ R (k2+k)×k for empirical estimate T . Then for any ≤ γ/2, for any δ < 1, 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 -δ,

SS -Ŝ Ŝ ≤ √ 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 In order to lowerbound the k 2 -th singular value of T , we first express T as the product of four simpler matrices, T = M BCF , as illustrated in Figure 5 . The definitions of these four matrices

SS -Ŝ Ŝ ≤ √ 2 T -T F σ k2 (T ) ≤ √ 2 αβ 2 . M ∈ R d 2 ×(k 2 +1) B ∈ R (k 2 +1)×k 2 C ∈ R k 2 ×k 2 F ∈ R k 2 ×(k 2 +k) T = M BCF ∈ R d 2 ×(k2+k) d 2 k 2 + 1 k 2 k 2 k 2 + k recall that k 2 := k 2

M ∈ R d 2 ×(k2+1) , B ∈ R (k2+1)×k2 , C ∈ R (k2×k 2 ) , F ∈ R k 2 ×(k2+k)

will be introduced later as we explain their effects. From Lemma 9, we know that

T vec * (U ) = 1≤i<j≤k (A U A) ij M ij - 1≤i<j≤k (A U A) ij m ij E[x ⊗ x],

for any symmetric k × k matrix U .

Note that vec * (U ) is a (k 2 + k)-dimensional vector. For convenience, we first use matrix F to transform vec * (U ) to vec(U ), which has k 2 dimensions. Matrix F is defined such that F vec * (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 U A) ij 's. Notice that vec(A U A) = A ⊗ A vec(U ) = A ⊗ A F vec * (U ).

Since we only care about the elements of A U A 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 k 2 × k 2 .

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 T vec * (U ) in the following way:

T vec * (U ) = 1≤i<j≤k (A U A) ij M ij -m ij E[x ⊗ x] . Thus, T vec * (U ) is just a linear combination of (M ij -m ij E[x ⊗ x])'s with coefficients equal to (A U A) ij 's.

We have already expressed coefficients (A U A) ij 's using CF vec * (U ). Now, we just need to use matrix B to transform the augmented distinguishing matrix M to a d 2 × k 2 matrix, with each column equal to

(M ij -m ij E[x ⊗ x]).

In order to achieve this, the first k 2 rows of B is just the identity matrix I k2 , and the last row of

B is [-m 12 , -m 13 , • • • , -m 1k , -m 23 , -m 24 , • • • ] .

With above characterization of T , we are ready to show that the k 2 -th singular value of T is lower bounded. Lemma 14. Suppose σ min (M ) ≥ α, σ min (A) ≥ β, we know that matrix T has rank k 2 and the k 2 -th singular value of T is lower bounded by αβ 2 .

Proof. Since matrix C has dimension k 2 × k 2 , it's clear that the rank of T is at most k 2 . We first prove that the rank of T is exactly k 2 .

Since the first k 2 rows of B constitute the identity matrix I k2 , we know B is a full-column rank matrix with rank equal to k 2 . We also know that matrix M is a full column rank matrix with rank k 2 + 1. Thus, the product matrix M B is still a full-column rank matrix with rank k 2 . If we can prove that the product matrix CF has full-row rank equal to k 2 . It's clear that T = M BCF also has rank k 2 . 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 ∈ R k2 , such that k2 l=1 a l (CF ) [l,:] = 0. Note that for any

1 ≤ l ≤ k 2 , (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 a l (CF ) [l,:] = 0 simply implies k2 l=1 a l C [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 k 2 -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 k 2 rows of matrix B constitute a k 2 × k 2 identity matrix, we know

σ min (B) := min u: u ≤1 Bu ≥ min u: u ≤1 I k2×k2 u =: σ min (I k2×k2 ) = 1,

where u is any k 2 -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 k 2 and v has dimension k 2 .

We lowerbound the smallest singular value of CF by showing that σ min (CF ) ≥ σ min (C). For any unit vector u ∈ R k2 , 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 k 2 , the k 2 -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.

B.1.2 UPPERBOUNDING T -T F:

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 ≤ P 1 , ξ ≤ P 2 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∈R k 2 +k : v ≤1 (T -T )v ≤ max U ∈R k×k sym : U F ≤ √ 2 (T -T )vec * (U ) .

For any k × k symmetric matrix U with eigenvalue decomposition U = k i=1 λ i u (i) (u (i) ) , according to the definition of T , we know max

U ∈R k×k sym : U F ≤ √ 2 (T -T )vec * (U ) = max U ∈R k×k sym : U F ≤ √ 2 T vec * (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 λ 2 i 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 max u: u ≤2 f (u) -f (u) . Recall that

f (u) =2E (u y) • (E[(u y)x ]E[xx ] -1 x) • (x ⊗ x) -E (u y) 2 (x ⊗ x) + E (u y) 2 -2E (u y)x E xx -1 E (u y)x E[x ⊗ x],

and

f (u) =2 Ê (u y) • ( Ê[(u y)x ] Ê[xx ] -1 x) • (x ⊗ x) -Ê (u y) 2 (x ⊗ x) + Ê (u y) 2 -2 Ê (u y)x Ê xx -1 Ê (u y)x Ê[x ⊗ x].

Notice that

E (u y)•(E[(u y)x ]E[xx ] -1 x)•(x⊗x) = E (u y)x E xx -1 E (u y)x(x⊗x) .

We first show that given polynomial number of samples,

2 Ê (u y)x Ê xx -1 Ê (u y)x(x ⊗ x) -2E (u y)x E xx -1 E (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 ≤ P 1 , ξ ≤ P 2 , we have

(u y)x ≤ u Aσ(W x) + ξ x ≤ u ( A W x + ξ ) x ≤2Γ(ΓP 1 √ k + P 2 ).

According to Lemma 24, we know given O(

Γ 2 (ΓP1 √ k+P2) 2 log( d δ )

2

) number of samples,

Ê (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,

Ê (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,

Ê xx -E xx ≤ , with probability at least 1 -δ. Suppose that ≤ γ/2 ≤ σ min (E[xx ])/2, we know Ê[xx ] has full rank.

According to Lemma 29, we have

Ê xx -1 -E xx -1 ≤ 2 √ 2 Ê xx -E xx σ 2 min (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

Ê (u y)x -E (u y)x ≤ , Ê xx -1 -E xx -1 ≤ 2 √ 2 /γ 2 , Ê (u y)x(x ⊗ x) -E (u y)x(x ⊗ x) ≤ . Define E 1 := Ê (u y)x -E (u y)x , E 2 := Ê xx -1 -E xx -1 , E 3 := Ê (u y)x(x ⊗ x) -E (u y)x(x ⊗ x) .

Then, we have

2 Ê (u y)x Ê xx -1 Ê (u y)x(x ⊗ x) -2E (u y)x E xx -1 E (u y)x(x ⊗ x) ≤2 E 1 E 2 E 3 + 2 E 1 E xx -1 E (u y)x(x ⊗ x) + 2 E (u y)x E 2 E (u y)x(x ⊗ x) + 2 E (u y)x E xx -1 E 3 + 2 E 1 E 2 E (u y)x(x ⊗ x) + 2 E (u y)x E 2 E 3 + 2 E 1 E xx -1 E 3 ≤2 2 √ 2 3 γ 2 + 2Γ 3 (ΓP 1 √ k + P 2 ) γ + 8 √ 2Γ 4 (ΓP 1 √ k + P 2 ) 2 γ 2 + 2Γ(ΓP 1 √ k + P 2 ) γ + 4 √ 2Γ 3 (ΓP 1 √ k + P 2 ) 2 γ 2 + 4 √ 2Γ(ΓP 1 √ k + P 2 ) 2 γ 2 + 2 γ =O( Γ 4 (ΓP 1 √ k + P 2 ) 2 γ 2 ).

Thus, given O(

Γ 14 (ΓP1 √ k+P2) 6 log( d δ ) γ 4 2

) number of samples, we know

2 Ê (u y) • ( Ê[(u y)x ] Ê[xx ] -1 x) • (x ⊗ x) -2E (u y) • (E[(u y)x ]E[xx ] -1 x) • (x ⊗ x) = 2 Ê (u y)x Ê xx -1 Ê (u y)x(x ⊗ x) -2E (u y)x E xx -1 E (u y)x(x ⊗ x) ≤ ,

with probability at least 1 -δ.

Now, let's consider the second term

Ê (u y) 2 (x ⊗ x) -E (u y) 2 (x ⊗ x) . Since (u y) 2 (x ⊗ x) ≤ 4Γ 2 (ΓP 1 √ k + P 2 ) 2

, according to Lemma 24, we know given

O( Γ 4 (ΓP1 √ k+P2) 4 log( d δ ) 2 ), Ê (u y) 2 (x ⊗ x) -E (u y) 2 (x ⊗ x) ≤ ,

with probability at least 1 -δ.

Next, let's look at the third term

Ê (u y) 2 Ê[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

Ê (u y) 2 -E (u y) 2 ≤ , Ê[x ⊗ x] -E[x ⊗ x] ≤ .

Thus, Define:

E 4 := Ê (u y) 2 -E (u y) 2 E 5 := Ê[x ⊗ x] -E[x ⊗ x].

Then, we have

Ê (u y) 2 Ê[x ⊗ x] -E (u y) 2 E[x ⊗ x] ≤ E 4 E 5 + E 4 E[x ⊗ x] + E (u y) 2 E 5 ≤ 2 + Γ 2 + 4(ΓP 1 √ k + P 2 ) 2 =O((ΓP 1 √ k + P 2 ) 2 ).

Thus, we know that given O(

Γ 2 (ΓP1 √ k+P2) 6 log( d δ )

2

) number of samples, we know

Ê (u y) 2 Ê[x ⊗ x] -E (u y) 2 E[x ⊗ x] ≤ ,

with probability at least 1 -δ.

Now, let's bound the last term,

2 Ê (u y)x Ê xx -1 Ê (u y)x Ê[x⊗x]-2 E (u y)x E xx -1 E (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 Ê (u y)x Ê xx -1 Ê (u y)x Ê[x⊗x]-2 E (u y)x E xx -1 E (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, 2 Ê (u y) • ( Ê[(u y)x ] Ê[xx ] -1 x) • (x ⊗ x) -2E (u y) • (E[(u y)x ]E[xx ] -1 x) • (x ⊗ x) ≤ Ê (u y) 2 (x ⊗ x) -E (u y) 2 (x ⊗ x) ≤ , Ê (u y) 2 Ê[x ⊗ x] -E (u y) 2 E[x ⊗ x] ≤ , 2 Ê (u y)x Ê xx -1 Ê (u y)x Ê[x ⊗ x] -2 E (u y)x E xx -1 E (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 ≤ k 2 + k T -T ≤k √ 2k max u: u ≤2 f (u) -f (u) ,

where the second inequality holds since

T -T ≤ √ 2k max u: 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 -δ.

B.2 ROBUST ANALYSIS FOR SIMULTANEOUS DIAGONALIZATION:

In this section, we will show that the simultaneous diagonalization step in our algorithm is robust. Let S and Ŝ be two (k 2 + 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 -Ŝ Ŝ⊥ is well bounded. However, due to the rotation issue of subspace basis, we cannot conclude that S -Ŝ F is small. Only after appropriate alignment, the difference between S and Ŝ becomes small. Lemma 16. Let S and Ŝ be two (k 2 + k) by k matrices, whose columns consist of the least k right singular vectors of T and T respectively. If SS -Ŝ Ŝ F ≤ , there exists an rotation matrix

R ∈ R k×k satisfying RR = R R = I k , such that Ŝ -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 

Ŝ -SR F ≤ SS -Ŝ Ŝ F 2( √ 2 -1) σ k (SS ) ≤ 2 . Let the k columns of S be vec * (U 1 ), vec * (U 2 ), • • • , vec * (U k ). Note each U i can be expressed as A -D i A -1 , where D i is a diagonal matrix. Let Q be a k × k matrix, whose i-th column consists of the diagonal elements of D i , such that Q ij equals the j-th diagonal element of D i . Let vec * (X) = SRζ 1 , vec * (Y ) = SRζ 2 ,

Y = diag(QRζ 2 ). It's not hard to check that X = A -D X A -1 and Y = A -D Y A -1 . Furthermore, we have XY -1 = A -D X D -1 Y A . Next, we show that the diagonal elements of D X D -1

Y are well separated. Lemma 17. Assume that A ≤ P 1 , σ min (A) ≥ β. Then for any δ > 0 we know with probability at least 1 -δ, we have sep(

D X D -1 Y ) ≥ poly(1/d, 1/P 1 , β, δ), where sep(D X D -1 Y ) := min i =j |(D X D -1 Y ) ii -(D X D -1 Y ) jj |.

Proof. We first show that matrix Q is well-conditioned. Since

U i = A -D i A -1 , we have vec(U i ) = A -⊗A -vec(D i ).

Let U be a k 2 ×k matrix whose columns consist of vec(U i )'s. Also define Q as a k 2 × k matrix whose columns are vec(D i )'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 -) = σ 2 min (A -1 ) = 1 A 2 ≥ 1 P 2 1 . Thus, we have Q ≤ √ 2P 2 1 , which implies Q ≤ √ 2P 2 1 . Now, let's prove that the diagonal elements of D X D -1

Y are well-separated. Let q i be the i-th row vector of Q. Then we know the i-th diagonal element of

EQUATION

1 for every row vector. It's not hard to show that with probability at least 1 -exp(-d Ω(1) ), we have | q i , Rζ 2 | ≤ poly(d, P 1 ) for each i. Now given ζ 2 for which this happens, we have qi,Rζ1 qi,Rζ2 -qj ,Rζ1 qj ,Rζ2 = c i q i , Rζ 1 -c j q j , Rζ 1 , where c i , c j have magnitude as least poly(1/d, 1/P 1 ). Since σ min (Q) ≥ β 2 , we know Proj q ⊥ j q i ≥ β 2 (because otherwise there exists λ such that (e i + λe j ) σ min (Q) ≤ (e i + λe j ) Q = Proj q ⊥ j q i < β 2 , which is a contradiction). Let q ⊥ i,j = Proj q ⊥ j q i = q i -λ i,j q j , we can rewrite this as

q i , Rζ 1 q i , Rζ 2 - q j , Rζ 1 q j , Rζ 2 = c i q i , Rζ 1 -c j q j , Rζ 1 = c i q ⊥ i,j , Rζ 1 -(c j + λ i,j c i ) q j , Rζ 1 .

By properties of Gaussians, we know q ⊥ i,j , Rζ 1 is independent of q j , Rζ 1 , so we can first fix q j , 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 -δ/k 2 :

q i , Rζ 1 q i , Rζ 2 - q j , Rζ 1 q j , Rζ 2 ≥ poly(1/d, 1/P 1 , β, δ).

By union bound, we know with probability at least 1 -δ,

sep(D X D -1 Y ) ≥ poly(1/d, 1/P 1 , β, δ).

Let X = Ŝζ 1 and Ŷ = Ŝζ 2 . Next, we prove that the eigenvectors of X Ŷ -1 are close to the eigenvectors of XY -1 . 

Lemma 12. Suppose that SS -Ŝ Ŝ F ≤ , A ≤ P 1 , ξ ≤ P 2 , σ min (E[xx ]) ≥ γ, σ min (A) ≥ β. Let X = mat * ( Ŝζ 1 ), Ŷ = mat * ( Ŝζ 2 ),

(i) ∈ [k] such that ẑi -z π[i] ≤ poly(d, P 1 , 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 ẑ 1 , • • • , ẑ k for X Ŷ -1 . We first prove that the eigenvectors of X Ŷ -1 are close to the eigenvectors of XY -1 . Let X = X + E X and Ŷ = Y + E Y . Then we have X Ŷ -1 = XY -1 (I + F ) + G,

where

F = -E Y (I + Y -1 E Y ) -1 Y -1 and G = E X Ŷ -1 .

According to Lemma 34, we have

F ≤ E Y σmin(Y )-E Y and G ≤ E X σmin( Ŷ )

. In order to bound the perturbation matrices F and G , we need to first bound E X , E Y and σ min (Y ), σ min ( Ŷ ).

As we know, E:

X = X -X = ( Ŝ -SR)ζ 1 . According to Lemma 16, we have Ŝ -SR ≤ Ŝ -SR F ≤ 2 .

We also know with probability at least 1 -exp(-d Ω(1) ), ζ 1 ≤ poly(d). Thus, we have

E X = ( Ŝ -SR)ζ 1 ≤ Ŝ -SR ζ 1 ≤ poly( , d).

Similarly, with probability at least 1 -exp(-d Ω( 1) ), we also know E Y ≤ poly( , d).

Now, we lower bound the smallest singular value of X and

Y . Since X = A -D X A -1 , we have σ min (X) ≥σ 2 min (A -1 )σ min (D X ) = 1 A 2 σ min (D X ) ≥ 1 P 2 1 σ min (D X ).

Since D X is a diagonal matrix, its smallest singular value equals the smallest absolute value of its diagonal element. Recall each diagonal element of D X is q i , Rζ 1 , which follows a Gaussian distribution whose standard deviation is at least q i ≥ β 2 . By anti-concentration property of Gaussian (see Lemma 38), we know | q i , Rζ 2 | ≥ Ω(δβ 2 /k) for all i with probability 1 -δ/4. Thus, we have σ min (X) ≥ poly(1/d, 1/P 1 , β, δ). Similarly we have the same conclusion for Y .

For small enough E Y , we have

σ min ( Ŷ ) ≥ σ min (Y ) -E Y .

Thus, for small enough , we have F ≤ poly(d, P 1 , 1/β, , 1/δ) and G ≤ poly(d, P 1 , 1/β, , 1/δ). In order to apply Lemma 33, we also need to bound κ(A -) and XY -1 . Since σ min (A -) ≥ 1 P1 and A - ≤ 1/β, we have κ(A -) ≤ P 1 /β. For the norm of XY -1 , we have

XY -1 ≤ X σ min (Y ) ≤ X poly(d, P 1 , 1/β, 1/δ),

where the second inequality holds because σ min (Y ) ≥ poly(1/d, 1/P 1 , β, δ). 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, P 1 , 1/β, 1/δ). Similarly, we can also prove that

D X D -1 Y ≤ poly(d, P 1 , 1/β, 1/δ) Overall, we have κ(A -)( XY -1 F + G ) ≤κ(A -)( XY -1 F + G ) ≤ P 1 β poly(d, P 1 , 1/β, 1/δ)poly(d, P 1 , 1/β, , 1/δ) + poly(d, P 1 , 1/β, , 1/δ) ≤poly(d, P 1 , 1/β, , 1/δ).

According to Lemma 17, we know with probability at least 1 -δ/4, sep(D X D -1 Y ) ≥ poly(1/d, 1/P 1 , β, δ). Thus, by union bound, we know for small enough , with probability at least

1 -δ, κ(A -)( XY -1 F + G ) < sep(D X D -1 Y )/(2k).

According to Lemma 33, we know there exists a permutation π

[i] ∈ [k], such that ẑ i -z i ≤3 F D X D -1 Y + G σ min (A -)sep(D X D -1 Y ) ≤3 poly(d, P 1 , 1/β, , 1/δ)poly(d, P 1 , 1/β, 1/δ) + poly(d, P 1 , 1/β, , 1/δ) 1/P 1 poly(1/d, 1/P 1 , β, 1/δ) ≤poly(d, P 1 , 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 ẑ 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 Ê[ẑ i y]

, with probability at least 1 -δ, the sign flip of ẑ i is consistent with the sign flip of z i .

In the following lemma, we show that the sign flip step of ẑi is robust.

Lemma 18. Suppose that x ≤ Γ, A ≤ P 1 , ξ ≤ P 2 . Let z 1 , • • • , z k be the eigenvectors of XY -1 (before sign flip step). Similarly define ẑ 1 , • • • , ẑ k for X Ŷ -1 . Suppose for each i, z i -ẑ 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 ]) = sign( Ê[ ẑ 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 ]| = λ i E[σ(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 w i σ 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 ] -Ê[ ẑ i , y ] . E[ z i , y ] -Ê[ ẑ i , y ] = E[ z i , y ] -E[ ẑ i , y ] + E[ ẑ i , y ] -Ê[ ẑ i , y ] ≤ E[ z i -ẑ i , y ] + E[ ẑ i , y ] -Ê[ ẑ i , y ] ≤(ΓP 1 √ k + P 2 ) + E[ ẑ i , y ] -Ê[ ẑ i , y ] .

Note that we reserve fresh samples for this step, thus ẑ i is independent with samples in Ê[ ẑ i , y ]].

Since ẑ i , y ≤ ΓP 1 √ k + P 2 , we know with O( (ΓP1 √ k+P2) 2 log(d/δ) 2

) number of samples, we have

E[ ẑ i , y ] -Ê[ ẑ i , y ]] ≤ . Overall, we have E[ z i , y ] -Ê[ ẑ i , y ]] ≤ (1 + ΓP 1 √ k + P 2 ) . 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 ]) = sign( Ê[ ẑ i , y ])] ≤ δ.

B.3 ROBUST ANALYSIS FOR RECOVERING FIRST LAYER WEIGHTS:

We will first show that Algorithm 1 is robust.

Theorem 8. Assume that x ≤ Γ, w ≤ 1, |ξ| ≤ P 2 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 ≤ ,

where ŵ is the learned weight vector.

Proof. We first show that given polynomial number of i.i.d. samples,

Ê[yx] -E[yx] is upper bounded with high probability, where Ê[yx] is the empirical estimate of E[yx]. Due to the assump- tion that x ≤ Γ, w ≤ 1, ξ ≤ P 2 , we have yx = (σ(w x) + ξ)x ≤ w xx + ξx ≤ w x 2 + |ξ| x ≤Γ 2 + P 2 Γ

According to Lemma 24, we know given O(

(Γ 2 +P2Γ) 2 log( d δ )

2

) number of samples, with probability at least 1 -δ, we have

Ê[yx] -E[yx] ≤ . Since xx ≤ Γ 2 , we know that given O( Γ 4 log( d δ ) 2 ) number of samples, Ê[xx ] -E[xx ] ≤ with probability at least 1 -δ. Suppose that ≤ γ/2 ≤ σ min (E[xx ])/2, we know Ê[xx ] has full rank. According to Lemma 29, we have Ê xx -1 -E xx -1 ≤ 2 √ 2 Ê xx -E xx σ 2 min (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

Ê[xx ] -1 -E[xx ] -1 ≤ 2 √ 2 /γ 2 Ê[yx] -E[yx] ≤ Denote E 1 := Ê[xx ] -1 -E[xx ] -1 E 2 := Ê[yx] -E[yx].

Then, we have

ŵ -w =2 Ê[xx ] -1 Ê[yx] -E[xx ] -1 E[yx] ≤2 E 1 E 2 + 2 E 1 E[yx] + 2 E[xx ] -1 E 2 ≤ 4 √ 2 2 γ 2 + 4 √ 2(Γ 2 + P 2 Γ) γ 2 + 2 γ ≤O( (Γ 2 + P 2 Γ) γ 2 ).

Thus, given O(

(Γ 2 +P2Γ) 4 log( d δ ) γ 4 2

) number of samples, with probability at least 1 -δ, we have ŵ -w ≤ . Now let's go back to the call to Algorithm 1 in Algorithm 3. Let z i 's be the normalized rows of A -1 , and let ẑi 's be the eigenvectors of X Ŷ -1 (with correct sign). From Lemma 12, we know {ẑ i } are close to {z i } with permutation. Without loss of generality, we assume the permutation here is just an identity mapping, which means z i -ẑi is small for each i.

For each z i , let v i be the output of Algorithm 1 given infinite number of inputs (x, z i y). For each ẑi , let vi be the output of Algorithm 1 given only finite number of samples (x, ẑ i y). In this section, we show that suppose z i -ẑi is bounded, with polynomial number of samples, v i -vi is also bounded.

The input for Algorithm 1 is (x, ẑ i y). We view ẑ i y as the summation of z i y and a noise term

(ẑ i -z i ) y.

Here, the issue is that the noise term (ẑ i -z i ) y is not independent with the sample (x, ẑ 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 ẑi is independent with the samples (x, y)'s used by Algorithm 1. Thus, the samples (x, ẑ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 ≤ P 1 , ξ ≤ P 2 and σ min (E[xx ]) ≥ γ. Suppose that for each 1 ≤ i ≤ k, ẑi -z i ≤ τ . 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 -δ,

v i -vi ≤ 2τ (ΓP 1 √ k + P 2 )Γ γ + 2 , for each 1 ≤ i ≤ k.

Proof. For any i, let's bound v i -vi . As we know,

v i = 2E[xx ] -1 E[z i yx] and vi = 2 Ê[xx ] -1 Ê[ẑ i yx].

Thus, in order to bound v i -vi , we only need to show

E[xx ] -1 E[ẑ i yx] -E[xx ] -1 E[z i yx] and E[xx ] -1 E[ẑ i yx] -Ê[xx ] -1 Ê[ẑ i yx]

are both bounded.

The first term can be bounded as follows.

E[xx ] -1 E[ẑ i yx] -E[xx ] -1 E[z i yx] = E[xx ] -1 E[(ẑ i -z i ) yx] ≤ E[xx ] -1 ẑi -z i Aσ(W x) + ξ x ≤ τ (ΓP 1 √ k + P 2 )Γ γ

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 ] -1 E[ẑ i yx] -Ê[xx ] -1 Ê[ẑ i yx] ≤ . Overall, we have v i -vi = 2E[xx ] -1 E[z i yx] -2 Ê[xx ] -1 Ê[ẑ i yx] ≤2 E[xx ] -1 E[ẑ i yx] -E[xx ] -1 E[z i yx] + 2 E[xx ] -1 E[ẑ i yx] -Ê[xx ] -1 Ê[ẑ i yx] ≤ 2τ (ΓP 1 √ k + P 2 )Γ γ + 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 -δ,

v i -vi ≤ 2τ (ΓP 1 √ k + P 2 )Γ γ + 2 .

for any 1 ≤ i ≤ k.

B.4 PROOF OF THEOREM 7:

Proof of Theorem 7. Combining Lemma 11, Lemma 12 and Lemma 13, we know given poly Γ, P 1 , P 2 , d, 1/ , 1/γ, 1/α, 1/β, 1/δ number of i.i.d. samples, with probability at least

1 -δ, z i -ẑi ≤ , v i -vi ≤ for each 1 ≤ i ≤ k.

Let V be a k × d matrix whose rows are v i 's. Similarly define matrix V for vi 's. Since v i -vi ≤ 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 z i 's. Similarly define matrix Ẑ for ẑi 's. Again, we have Z -Ẑ ≤ √ k . In order to show Z -1 -Ẑ-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/P 1 , and A -1 ≤ 1/β, which implies that every row vector of A -1 has norm at most 1/β. Let D z be the diagonal matrix whose i, i-th entry is the norm of i-th row of A -1 , then Z = D -1 z A -1 , and we know σ min (Z) ≥ σ min (D -1 z )σ min (A -1 ) ≥ β/P 1 . Then, according to Lemma 29, as long as ≤ β 2P1 , we have

Z -1 -Ẑ-1 ≤ 2 √ 2 Z -Ẑ σ 2 min (Z) ≤ 2 √ 2P 2 1 √ k /β 2 . Define E 1 :=Z -1 -Ẑ-1 E 2 :=V -V . We know E 1 ≤ 2 √ 2P 2 1 √ k /β 2 and E 2 ≤ √ k .

In order to bound V , we can bound the norm of its row vectors. We have,

v i = 2E[xx ] -1 E[z i yx] ≤ 2Γ(ΓP 1 √ k + P 2 ) γ , which implies V ≤ 2 √ kΓ(ΓP1 √ k+P2) γ . Now we can bound Z -1 σ(V x)-Ẑ-1 σ( V x) as follows. Z -1 σ(V x) -Ẑ-1 σ( V x) ≤ E 1 E 2 x + E 1 V x + Z -1 E 2 x ≤ 2 √ 2P 2 1 Γk 2 β 2 + 4 √ 2P 2 1 Γ 2 (ΓP 1 √ k + P 2 )k β 2 γ + P 1 Γ √ 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 Γ, P 1 , P 2 , d, 1/ , 1/γ, 1/α, 1/β, 1/δ number of i.i.d. samples, with probability at least 1 -δ,

Aσ(W x) -Ẑ-1 σ( V x) = Z -1 σ(V x) -Ẑ-1 σ( V x) ≤ ,

where the first equality holds because Aσ(W x) = Z -1 σ(V x), as shown in Theorem 5.

C SMOOTHED ANALYSIS FOR DISTINGUISHING MATRICES:

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, I d ), 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 (M D ) is lowerbounded.

Perturbing W for Gaussian Input We first consider the case when the input follows standard Gaussian distribution N (0, I d ). The weight matrix W is perturbed to W where

EQUATION

Here E ∈ R k×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 x w 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(I d ). 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 distribution N (0, I d ).

Given any weight matrix W with w i ≤ τ 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 D Q with a random covariance matrix. To sample a point in D Q , first sample n ∼ N (0, I d ), and then output Qn. The (Q, λ) perturbation of a distribution D, which we denote by D Q,λ is a mixture between the distribution D and the distribution D Q . More precisely, to sample x from D Q,λ , 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 D Q , 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 D Q with covariance matrix QQ . Since both D and D Q are symmetric, their mixture is also symmetric. Also, the TV-distance between D and D Q,λ is bounded by λ. Throughout this section we will use D to denote the perturbed distribution D Q,λ

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 M D is lower bounded. Recall that M D is defined as

M D ij = E x∼D (w i x)(w j x)(x ⊗ x)1{w i xw j x ≤ 0} ,

as the first k 2 columns and has E x∼D [x ⊗ x] as the last column. Theorem 10. Given weight matrix W with w i ≤ τ 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.

C.1 SMOOTHED ANALYSIS FOR GAUSSIAN INPUTS:

In this section, we will prove Theorem 9, as restated below: Theorem 9. Suppose that k ≤ d/5, and the input follows standard Gaussian distribution N (0, I d ).

Given any weight matrix W with w i ≤ τ 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

M ij : mat(M ij ) = 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 w i and w j , we have u mat(M ij )u be a constant. We can compute this constant as

m ij := E (w i x)(w j x)1{w i xw j x ≤ 0} .

This implies that if we consider mat(M ij ) -m ij I d , it is going to be a matrix whose rows and columns are in span of w i and w j . In fact we can compute the matrix explicitly as the following lemma: Lemma 19. Suppose input x follows standard Gaussian distribution N (0, I d ), and suppose weight matrix W has full-row rank, then for any

1 ≤ i < j ≤ k, we have mat(M ij ) = 1 π φ ij cos(φ ij ) -sin(φ ij ) w i w j I d + φ ij π (w i w j + w j w i ) - sin(φ ij ) π ( w i w j w j w j + w j w i w i w i ),

where 0 < φ ij < π is the angle between weight vectors w i and w j .

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 ∈ R d×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 A i be the i-th column of A and S -i be the span of all the columns except for A i , then Rudelson & Vershynin (2009) showed that one can lowerbound the smallest singular value of a matrix by its leave-one-out distance.

d(A) := min i∈[n] (I d -Proj S-i )A i .

Lemma 20 (Rudelson & Vershynin (2009) ). For matrix A ∈ R d×n (d ≥ n), we always have

d(A) ≥ σ min (A) ≥ 1 √ n d(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(I d ) 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 ∈ R d×(n+1) (d ≥ n + 1) be an arbitrary matrix whose columns are A 0 , A 1 , ..., A n . For any i = 1, 2, ..., n, let S -i be the subspace spanned by all the other columns (including A 0 ) except for A i , and let d (A) := min i=1,...,n (I d -Proj S-i )A i . Suppose the column A 0 has norm √ d and A 1 , ..., A n has norm at most C, then

σ min (A) ≥ min nC 2 d 4nC 2 + d , d 4n 2 C 2 + 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 w i and w j with norm w i , w j ≤ τ , let w i = w i + ρε i , w j = w j +ρε j where ε i , ε j are i.i.d. Gaussian random vectors. With probability at least 1-exp(-d Ω(1) ), we know w i ≤ τ + 3ρ 2 d/2, w j ≤ τ + 3ρ 2 d/2 and

φ ij ≥ ρ 2 (d -2) √ 2τ + 3ρ 2 d ,

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ρ 2 d/2, and for all i < j, the coefficient φ ij /π in front of the term w i w j + w j w i is at least

√ ρ 2 (d-2) ( √ 2τ + √ 3ρ 2 d)π

.

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 E 1 be the event that Lemma 22 does not hold. Event E 1 will be one of the bad events (but note that we do not condition on E 1 not happening, we use a union bound at the end).

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 π .

Since the column span and row span of Proj S ⊥ ij M ij both belong to the subspace S ⊥ ij , there must exist a

(d -2) × (d -2) matrix C, such that Proj S ⊥ ij M ij = S ⊥ ij C(S ⊥ ij )

. We only need to show this matrix C must be m ij I d-2 . In order to show this, we prove for any u

, v ∈ S ⊥ ij , u (Proj S ⊥ ij M ij )v = m ij u v. u (Proj S ⊥ ij M ij )v = u M ij v = 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] = m ij u E[xx ]v = m ij u v,

where the fourth equality holds because u x, v x are independent with w i x, w j x.

Thus, we know

M ij =(Proj S ⊥ ij M ij + m ij Proj Sij I d ) + (Proj Sij M ij -m ij S ij S ij ) =m ij I d + (Proj Sij M ij -m ij S ij S ij ).

Let's now compute the closed form for m ij . Recall that

m ij := E (w i x)(w j x)1{w i xw j x ≤ 0} .

Note, we only need to consider input x within subspace S ij , 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 )(e (i,j) 1

m ij = 1 2π ∞ 0 r 3 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 Proj:

Sij M ij . Note Proj Sij M ij is symmetric, because Proj Sij M ij = M ij -m ij I d +m ij S ij S ij ,

) , (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

Proj Sij M ij = 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 )(e (i,j) 1

) T = 1 2π ∞ 0 r 5 exp(- r 2 2 )dr π 2 π 2 -φij cos 3 (θ) 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 = P roj Sij M ij , (e (i,j) 2

)(e (i,j) 2 )(e (i,j) 1

) T = 1 2π ∞ 0 r 5 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 )),

) T = 1 2π ∞ 0 r 5 exp(- r 2 2 )dr π 2 π 2 -φij -cos 2 (θ) cos(θ + φ ij ) sin(θ) + cos(θ) cos 2 (θ + φ ij ) sin(θ + φ ij ) dθ = 1 π (φ ij sin(φ ij ) -cos(φ ij ) sin 2 (φ ij )).

It's easy to check that c

(i,j) 12 = c (i,j) 21 . Let M ij be Proj Sij M ij -m ij S ij S ij .

Then, according to above computation, we know

M ij =Proj Sij M ij -m ij S ij S ij =(c (i,j) 11 -m ij )(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 -m ij )(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 ) sin 2 (φ ij ) (e (i,j) 1

)(e (i,j) 2

) + (e (i,j) 2

)(e (i,j) 1

) .

Since e (i,j) 1

= w i and e (i,j) 2 = 1 sin(φij ) w j -cot(φ ij )w i , we can also express M ij as a linear combination of w i w i , w j w j , w i w j and w j w i .

M ij = 1 4π 8φ ij cos(φ ij ) -5 sin(φ ij ) -sin(3φ ij ) w i w i + 1 4π -3 sin(φ ij ) + sin(3φ ij ) ( w j sin(φ ij ) -cot(φ ij )w i )( w j sin(φ ij ) -cot(φ ij )w i ) + 1 π φ ij sin(φ ij ) -cos(φ ij ) sin 2 (φ ij ) w i ( w j sin(φ ij ) -cot(φ ij )w i ) + ( w j sin(φ ij ) -cot(φ ij )w i )w i = φ ij π (w i w j + w j w i ) - sin(φ ij ) π (w j w j + w i w i ).

Thus,

M ij =m ij I d + M ij = 1 π (φ ij cos(φ ij ) -sin(φ ij ))I d + φ ij π (w i w j + w j w i ) - sin(φ ij ) π (w j w j + w i w i ).

Finally, if the rows w i , w j do not have unit norm, let wi = w i / w i , wj = w j / w j , we know

m ij = E (w i x)(w j x)1{w i xw j x ≤ 0} = w i w j E ( w i x)( w j x)1{w i xw j x ≤ 0} = 1 π (φ ij cos(φ ij ) -sin(φ ij )) w i w j .

Here we used the fact that the indicator variable does not change whether we use w i , w j or wi , wj . Similarly,

M ij = E (w i x)(w j x)xx 1{w i xw j x ≤ 0} = w i w j E ( w i x)( w j x)xx 1{w i xw j x ≤ 0} = 1 π (φ ij cos(φ ij ) -sin(φ ij )) w i w j I d + φ ij π (w i w j + w j w i ) - sin(φ ij ) π ( w i w j w j w j + w j w i w i w 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 * ∈ argmin u: u =1 Au . Let u * i be the coordinate corresponding to the column A i , for 0 ≤ i ≤ n. We consider two cases here.

If |u * 0 | ≥ 4nC 2

4nC 2 +d , then we have

σ min (A) = Au * = u * 0 A 0 + 1≤i≤n u * i A i ≥ u * 0 A 0 - 1≤i≤n u * i A i ≥ 4nC 2 4nC 2 + d √ d -( 1≤i≤n |u * i |)C ≥ 4nC 2 d 4nC 2 + d - √ n d 4nC 2 + d C = nC 2 d 4nC 2 + d ,

where the third inequality uses Cauchy-Schwarz inequality.

If |u * 0 | < 4nC 2 4nC 2 +d , we know 1≤i≤n |u * i | 2 ≥ d 4nC 2 +d . Let k ∈ argmax 1≤i≤n |u * i |. We know that |u * k | ≥ d 4n 2 C 2 +nd . Thus, σ min (A) ≥ u * k A k + i:i =k u * i A i =|u * k | A k + i:i =k u * i u * k A i ≥|u * k | (I d -Proj S -k )A k ≥|u * k |d (A) ≥ d 4n 2 C 2 + nd d (A).

Above all, we know that the smallest singular value of A is lower bounded as follows,

σ min (A) ≥ min nC 2 d 4nC 2 + d , d 4n 2 C 2 + 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 = w i + ρε i , w j = w j + ρε 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 S wj , and denote the subspace spanned by { w j , w i } as S wj ∪wi . Thus, we have

Proj S ⊥ w j w i ≥ Proj S ⊥ w j ∪w i (w i + ρε i ) =ρ Proj S ⊥ w j ∪w i ε i

where S ⊥ wj is the orthogonal subspace of S wj .

Fix ε j , then S ⊥ wj ∪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 ⊥ wj ∪wi . Thus, it's not hard to check that Proj S ⊥ w j ∪w i

ε i d = U ε, where ε ∈ R d-2 is a standard Gaussian vector. Denote Y := Proj S ⊥ w j ∪w i ε 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)t 2 8

, ∀t ∈ (0, 1).

Let t = 1 2 , we know that with probability at least 1 -2 exp( -(d-2)

32:

),

Y ≥ d -2 2 .

Thus, we have Proj S ⊥ w j

w i ≥ ρ Proj S ⊥ w j ∪w i ε i ≥ ρ 2 (d-2) 2 . Recall that Proj S ⊥ w j w i = sin( φ ij ) w i .

We also know

w i = w i + ρε i ≤ w i + ρ ε i =τ + ρ ε i ,

where the last equality holds since w i ≤ τ . 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( -d 32 ), ε i 2 ≤ 3d 2 .

By union bound, we know with probability at least 1 -2 exp( -d 32 ) -2 exp( -(d-2)

),

sin( φ ij ) w i ≥ ρ 2 (d -2) 2 , w i ≤τ + 3ρ 2 d 2 .

Combined with the fact that φ ij ≥ sin( φ ij ) when φ ij ∈ [0, π], we know with probability at least 1 -2 exp( -d 32 ) -2 exp( -(d-2)

),

φ ij ≥ sin( φ ij ) = sin( φ ij ) w i w i ≥ ρ 2 (d-2) 2 τ + 3ρ 2 d 2 = ρ 2 (d -2) √ 2τ + 3ρ 2 d

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ρ 2 d/2 ∀i < j, φ ij π ≥ ρ 2 (d -2) ( √ 2τ + 3ρ 2 d)π C.3 SMOOTHED ANALYSIS FOR GENERAL INPUTS

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 M D is at least inverse polynomial.

Recall the definition of (Q, λ)-perturbation: we mix the original distribution D with a distribution D Q 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, I d ) and let x = Qn. We will prove Theorem 10 which we restate below:

Theorem 10. Given weight matrix W with w i ≤ τ 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 M D , M D Q , M D be the augmented distinguishing matrices for distributions D, D Q and D respectively. Since D is a mixture of D and D Q , and the augmented distinguishing matrix is defined as expectations over samples, we immediately have

M D = (1 -λ)M D + λM D Q .

Our proof will go in two steps. First we will show that σ min (M D Q ) is large. Then we will show that even mixing with M D will not significantly reduce the smallest singular value, so σ min (M D ) 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 (M D Q ). This instance has weight W and input distribution N (0, QQ ). Let M W Q be the augmented distinguishing matrix for an instance with weight W Q and input distribution N (0, I d ). Our first observation shows that M D Q and M W Q are closely related, and we only need to analyze the smallest singular value of M W Q . The problem now is very similar to what we did in Theorem 9, except that the weight W Q 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 W Q 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 (µ, I d ) 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 Σ I d , 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 Σ = Σ -I d , by assumption we know Σ is still a positive semidefinite matrix. Let x ∼ N (µ, Σ), x ∼ N (µ, Σ ) and δ ∼ N (0, I d ), by property of Gaussians it is easy to see that d x and d δ be the density function for x, x , δ respectively, then we know for any point u d

x d = x + δ. Let d x ,

x (u) = E x ∼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 (M D Q ) is large.

Subspace Decoupling Next we need to consider the mixture M D . The worry here is that although σ min (M D Q ) is large, mixing with D might introduce some cancellations and make σ min (M D ) much smaller. To prove that this cannot happen with high probability, the key observation is that in the first step, to prove σ min (M W Q ) is large we have only used the property of W Q. 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 W Q! 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: M D is a d 2 by (k 2 + 1) matrix, where the first k 2 columns consist of

M D ij := E x∼D (w i x)(w j x)(x ⊗ x)1{w i xw j x ≤ 0} , and the last column is E x∼D [x ⊗ x].

According to the definition of (Q, λ)-perturbation, if we let D Q be N (0, QQ ), then we have

M D = (1 -λ)M D + λM D Q .

In the first step, we will try to analyze M D Q . The first k 2 columns of this matrix M D Q can be written as:

M D Q ij = E x∼D Q (w i x)(w j x)(x ⊗ x)1{w i xw j x ≤ 0} = E n∼N (0,I d ) (w i Qn)(w j Qn)(Qn ⊗ Qn)1{w i Qnw j Qn ≤ 0} = Q ⊗ QE n∼N (0,I d ) (w i Qn)(w j Qn)(n ⊗ n)1{w i Qnw j Qn ≤ 0}

for any i < j, and the last column is

E x∼D Q [x ⊗ x] = E n∼N (0,I d ) [Qn ⊗ Qn] = Q ⊗ QE n∼N (0,I d ) [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 W Q and input distribution is N (0, I d ). We use M W Q to denote the augmented distinguishing matrix of such a network, then we have

M D Q = Q ⊗ QM W Q .

Therefore we can first analyze the smallest singular value of M W Q . Let W = W Q. Note that Q is a Gaussian matrix, and W is fixed, so W Q is also a Gaussian random matrix except its entries are not i.i.d. More precisely, there are only correlations within columns of W Q, and for any column of W Q, the covariance matrix is W W . Since the smallest singular value of W is at least ρ, we know σ min (W W ) ≥ ρ 2 . Let the covariance matrix of W Q be Σ W Q ∈ R kd×kd , which has smallest singular value at least ρ 2 . Therefore we know Σ W Q ρ 2 I kd . It's not hard to verify that with probability at least 1 -exp(-d Ω(1) ), the norm of every row of W Q is upper bounded by poly(τ, d). By Claim 1, any convex property that holds for any N (0, ρ 2 I kd ) perturbation must also hold for Σ W Q4 . Thus, we know with probability at least 1 -exp(-d Ω(1) ),

σ min (M W Q ) ≥ poly(1/τ, 1/d, ρ).

To prepare for the next step, we will rewrite M W Q as the product of two matrices. According to the closed form of M W Q ij in Lemma 19, we know each column of M W Q can be expressed as a linear combination of w i ⊗ w j 's and vec(I d ). Therefore:

M W Q = W ⊗ W , vec(I d ) R,

where matrix R has dimension (k 2 + 1) × (k 2 + 1). It's not hard to verify that σ min (M W Q ) ≤ W ⊗ W , vec(I d ) σ min (R).

Thus,

σ min (R) ≥ σ min (M W Q ) W ⊗ W , vec(I d )

.

Note that W is a k × d matrix with w i ≤ τ for every row vector, and W = W Q, 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(I d ) ≤ W ⊗ W , vec(I d ) 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 Proj W be the projection matrix to the row span of W , and let Proj d×(d-k) be a set of orthonormal basis for the orthogonal subspace W ⊥ . By symmetry of Gaussian, Proj W Q is independent with Proj W ⊥ Q. Thus, from now on we will condition on Proj W Q, and still treat Proj W ⊥ Q as a Gaussian random matrix. More precisely, Proj W ⊥ Q has the same distribution as U P , where P ∈ R (d-k)×d is a standard Gaussian matrix.

W ⊥ = I d -Proj W . Let Q = Proj W ⊥ Q. Let the columns of U ∈ R

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 ∈ R k×d be a set of orthonormal basis for the row span of W = W Q. 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 U P into the row span of V and V ⊥ as follows, U P d = U P 1 V ⊥ + U P 2 V,

where P 1 ∈ R (d-k)×(d-k) and P 2 ∈ R (d-k)×k are two independent standard Gaussian matrices. After this decomposition, we have

Q ⊗ Q W ⊗ W , vec(I d ) R d =U P ⊗ U P W ⊗ W , vec(I d ) R d =(U P 1 V ⊥ + U P 2 V ) ⊗ (U P 1 V ⊥ + U P 2 V ) W ⊗ W , vec(I d ) R =U ⊗ U (P 1 V ⊥ + P 2 V ) ⊗ (P 1 V ⊥ + P 2 V ) W ⊗ W , vec(I d ) R =U ⊗ U (P 1 V ⊥ + P 2 V ) W ⊗ (P 1 V ⊥ + P 2 V ) W , vec((P 1 V ⊥ + P 2 V )(P 1 V ⊥ + P 2 V ) ) R =U ⊗ U P 2 V W ⊗ P 2 V W , vec(P 1 P 1 + P 2 P 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 M D . Let Proj W ⊥ ⊗W ⊥ be Proj W ⊥ ⊗ Proj W ⊥ . We have,

σ min (M D ) =σ min (1 -λ)M D + λM D Q =σ min (1 -λ)M D + λQ ⊗ Q W ⊗ W , vec(I d ) R ≥σ min (1 -λ)Proj W ⊥ ⊗W ⊥ M D + λProj W ⊥ ⊗W ⊥ Q ⊗ Q W ⊗ W , vec(I d ) R =σ min (1 -λ)Proj W ⊥ ⊗W ⊥ M D + λ Q ⊗ Q W ⊗ W , vec(I d ) R

=σ min (1 -λ)Proj W ⊥ ⊗W ⊥ M D + λU ⊗ U P 2 V W ⊗ P 2 V W , vec(P 1 P 1 + P 2 P 2 ) R Since R has full row rank, we know that the row span of Proj W ⊥ ⊗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 Proj W ⊥ ⊗W ⊥ M D belongs to the column span of U ⊗ U . Thus, there exists matrix C ∈ R (d-k) 2 ×(k 2 +1) such that Proj W ⊥ ⊗W ⊥ M D = U ⊗ U CR.

Thus,:

σ min (M D ) =σ min (1 -λ)M D + λM D Q ≥σ min (1 -λ)Proj W ⊥ ⊗W ⊥ M D + λU ⊗ U P 2 V W ⊗ P 2 V W , vec(P 1 P 1 + P 2 P 2 ) R =σ min λU ⊗ U 1 -λ λ C + P 2 V W ⊗ P 2 V W , vec(P 1 P 1 + P 2 P 2 ) R .

Note that C only depends on U and R, U only depends on W , and R only depends on W Q. With W Q fixed, C is also fixed. Clearly, C is independent with P 1 and P 2 . For convenience, denote H := 1 -λ λ C + P 2 V W ⊗ P 2 V W , vec(P 1 P 1 + P 2 P 2 ) . Now, let's prove that the smallest singular value of matrix H ∈ R (d-k) 2 ×(k 2 +1) is lower bounded using leave-one-out distance. Let's first consider its submatrix Ĥ which consists of the first k 2 columns of H. Note that within random matrix P 2 V 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 Σ W Q ρ 2 I kd , 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 P 2 V W has smallest singular value at least γ := poly(1/d, ρ).

We can view P 2 V 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 K o + E where E is a Gaussian random matrix whose entries have variance γ. If d ≥ 7k, for any subspace S C that is independent of K and has dimension at most k 2 + 1, the leave-one-out distance d(Proj S ⊥ C K ⊗ 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 := P 2 V W , and denote i-th column of K as K i . We also think of the space S C as the column span of C.

As we did in Theorem 9, we partition [d -k] into 2 disjoint subsets L 1 and L 2 of size (d -k)/2. Let Ĥ be the set of rows of Ĥ indexed by L 1 × L 2 .

We fix a column Ĥ ij , i = j ∈ [k]. Let S = span{ Ĥ kl : (k, l) = (i, j)}. It's clear that S is correlated with Ĥ ij . Let C be the set of rows of C indexed by L 1 × L 2 . Let S C be the column span of C , which has dimension at most k 2 + 1. We define the following subspace that contains S, Ŝ = S C ∪ span{K j,L1 ⊗ x, x ⊗ K i,L2 , K l,L1 ⊗ x, x ⊗ K l,L2 x ∈ R (d-k)/2 , l / ∈ {i, j}}

Therefor by definition S ⊂ Ŝ, and thus Ŝ⊥ ⊂ S ⊥ , where S ⊥ denotes the orthogonal subspace of S.

Notice that Ĥ ij = K i,L1 ⊗ K j,L2 + λ 1-λ C ij is independent with Ŝ, assuming C is fixed. Moreover, Ŝ has dimension at most According to Lemma 20, we have

σ min (H ) ≥ 1 √ k 2 + 1 d(H )

≥poly(1/d, ρ).

Finally, we put everything together. Since H 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 ⊗ U HR) ≥λσ 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.

D TOOLS:

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.

D.1 MATRIX CONCENTRATION BOUNDS:

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 {Z k } of independent, random matrices with dimension d 1 × d 2 . Assume that each random matrix satisfies

E[Z k ] = 0 and Z k ≤ R almost surely. Define σ 2 := max k E[Z k Z * k ] , k E[Z * k Z k ] .

Then, for all t ≥ 0, 

EQUATION

:

Here we assume A ∈ R k×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.

Suppose the density function of distribution D is p(•), we assume p(x) = p(-x) for any x ∈ R d

The conclusion of Theorem 9 is clearly convex because it is a probability, and probabilities are linear in terms of mixture of distributions.

Figures and Tables:

Figure 2: Error in recovering W , A and outputs ("MSE") for different numbers of training samples and different dimensions of W and A. Each point is the result of averaging across five trials, where on the left W and A are both drawn as random 10 × 10 orthonormal matrices and in the center as 32 × 32 orthonormal matrices. On the right, given 10, 000 training samples we plot the square root of the algorithm's error normalized by the dimension of W and A, which are again drawn as random orthonormal matrices. The input distribution is a spherical Gaussian.

Figure 3: Error in recovering W , A and outputs ("MSE") for different amounts of label noise. Each point is the result of averaging across five trials with 10,000 training samples, where for each trial W and A are both drawn as 10 × 10 orthonormal matrices. The input distribution on the left is a spherical Gaussian and on the right a mixture of two Gaussians with one component based at the all-ones vector and the other component at its reflection.

Figure4: Error in recovering W , A and outputs ("MSE"), on the left for different levels of conditioning of W and on the right for A. Each point is the result of averaging across five trials with 20,000 training samples, where for each trial one parameter is drawn as a random orthonormal matrix while the other as described in Section 4.3. The input distribution is a mixture of Gaussians with two components, one based at the all-ones vector and the other at its reflection.

where ζ 1 and ζ 2 are two independent standard Gaussian vectors. Let z 1 , • • • z k be the normalized row vectors of A -1 . Let ẑ1 , ..., ẑk be the eigenvectors of X Ŷ -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 Ê[ẑ i y], with probability at least 1 -δ over the randomness of ζ 1 , ζ 2 and i.i.d. samples, there exists a permutation π

Figure 5: Characterize T as the product of four matrices.

where R is the rotation matrix in Lemma 16 and ζ 1 , ζ 2 are two independent standard Gaussian vectors. Let D X = diag(QRζ 1 ) and D

. Now, we only need to figure out the four coefficients of this linear combination. Similar as the computation for m ij , we use polar integration to show that, c (i,j) 11 = P roj Sij M ij , (e (i,j) 1

(d -k)/2 + (d -k)/2 + (k -2)(d -k)/2 + (k -2)(d -k)/2 + k 2 + 1 ≤ 4 5 (d -k) 2 4 , Let {c * kl , c * b } be one set of the optimal solutions to min c kl ,c b Ĥ ij + (k,l) =(i,j) Ĥ kl + c b b . If c *b = 0, we immediately have Ĥ inequality holds because the leave-one-out distance of matrix Ĥ is lower bounded by δ.If c * b = 0, we need to be more careful. In this case, we have,Ĥ ij + (k,l) =(i,j) inequality holds because the distance of b to the column span of Ĥ is lower bounded. If |c * b | ≥ δ 2C1 , we have Ĥ ij + (k,l) =(i,j) c * kl Ĥ kl + c * b b ≥ δζ 2C1 . leave-one-out distanceof H is lower bounded by poly(ζ, δ, 1/C 1 ). Recall that with probability at least 1 -exp(-d Ω(1) ), we have δ = poly(1/d, ρ), ζ = poly(1/d), C 1 = poly(d). Thus, we have d(H ) ≥ poly(1/d, ρ).

Pr k Z k ≥ t ≤ (d 1 + d 2 ) exp -t 2 /2 σ 2 + Rt/3 .As a corollary, we have:Lemma 24. Consider a finite sequence {Z 1 , Z 2 • • • Z m } of independent, random matrices with dimension d 1 × d 2 . Assume that each random matrix satisfies Z k ≤ R, 1 ≤ k ≤ m. Then, for all t ≥ 0, Pr m k=1 (Z k -E[Z k ]) ≥ t ≤ (d 1 + d 2 ) exp -t 2 /2 4mR 2 + (2Rt)/3 . Proof. For each k, let Z k = Z k -E[Z k ] be the new random matrices. It's clear that E[Z k ] = 0 and Z k ≤ 2R.For the variance,

Thus, according to Lemma 23, we havePr m k=1 (Z k -E[Z k ]) ≥ t = Pr m k=1 (Z k ) ≥ t ≤ (d 1 + d 2 ) exp -t 2 /2 4mR 2 + (2Rt)/3.

However, this is easily fixable by checkingE[z i y] = λ i E[σ(w i x)]. Since σ(w i x) is always nonnegative, E[z i y]has the same sign as λ i , and we can flip z i if E[z i y] is negative.

Table content: <table><tr><td>3.3 DETAILED ALGORITHM AND GUARANTEES</td></tr><tr><td>We can now give the full algorithm, see Algorithm 2. The main steps of this algorithm is as ex-</td></tr><tr><td>plained in the previous section. Steps 2 -5 constructs the pure neuron detector and finds the span of</td></tr><tr><td>vec  *  (z i z i ) (as in Corollary 1); Steps 7 -9 performs simultaneous diagonalization to get all the z i 's;</td></tr><tr><td>Steps 11, 12 calls Algorithm 1 to solve the single-layer problem and outputs the correct result.</td></tr><tr><td>Algorithm 2 Learning Two-layer Neural Networks</td></tr><tr><td>Input: Samples (x 1 , y 1 ), ..., (x n , y n ) generated according to Equation (4)</td></tr><tr><td>Output: Weight matrices W and A.</td></tr><tr><td>1: {Finding span of vec(z i z i )'s} 2: Estimate empirical moments Ê[xx ], Ê[yx ], Ê[y ⊗ x ⊗3 ] and Ê[y ⊗ y ⊗ (x ⊗ x)].</td></tr></table>

Algorithm 3 Learning Two-layer Neural Networks with Noise Input: Samples (x 1 , y 1 ), ..., (x n , y n ) generated according to Equation (2) Output: Weight matrices W and A.

Table content: <table/>

Algorithm 4 Learning Two-layer Neural Networks with Non-square A Input: Samples (x 1 , y 1 ), ..., (x n , y n ) generated according to Equation (2). Output: Weight matrices W ∈ R k×d and A ∈ l×k . 1: Using half samples (i.e. {(x i , y i )} : Let P be a l × k matrix, which columns are left singular vectors of Ê[yx ].

Table content: <table><tr><td>n/2</td></tr></table>

where ζ 1 and ζ 2 are two independent standard Gaussian vectors. Let z 1 , • • • z k be the normalized row vectors of A -1 . Let ẑ1 , ..., ẑk be the eigenvectors of X Ŷ -1 (after sign flip). For any δ > 0 and small enough , with i.i.d. samples in Ê[ẑ i y], with probability at least 1 -δ over the randomness of ζ 1 , ζ 2 and i.i.d. samples, there exists a permutation π

Table content: <table><tr><td>O( (ΓP1</td><td>√</td><td>k+P2) 2 log(d/δ)</td></tr></table>

and M ij , I d and S ij S ij are all symmetric. It's obvious that the column span of Proj Sij M ij belongs to subspace S ij . Combined with the fact that Proj Sij M ij is symmetric, we know the row span of Proj Sij M ij also belongs to the subspace S ij . Thus, matrix Proj Sij M ij can be represented as a linear combination of (e

Table content: <table><tr><td>(i,j)</td></tr><tr><td>1</td></tr></table>

Peer review (iclr2019):
 This is a strong theory paper and I recommend to accept. Paper Summary: This paper studies the problem of learning a two-layer fully connected neural network where both the output layer and the first layer are unknown. In contrast to previous papers in this line which require the input distribution being standard Gaussian, this paper only requires the input distribution is symmetric. This paper proposes an algorithm which only uses polynomial samples and runs in polynomial time. The algorithm proposed in this paper is based on the method-of-moments framework and several new techniques that are specially designed to exploit this two-layer architecture and the symmetric input assumption. This paper also presents experiments to illustrate the effectiveness of the proposed approach (though in experiments, the algorithm is slightly modified). Novelty: 1. This paper extends the key observation by Goel et al. 2018 to higher orders (Lemma 6). I believe this is an important generalization as it is very useful in studying multi-neuron neural networks. 2. This paper proposes the notation, distinguishing matrix, which is a natural concept to study multi-neuron neural networks in the population level. 3. The “Pure Neuron Detector” procedure is very interesting, as it reduces the problem of learning a group of weights to a much easier problem, learning a single weight vector. Clarity: This paper is well written. Major comments: My major concern is on the requirement of the output dimension. In the main text, this paper assumes the output dimension is the same as the number of neurons and in the appendix, the authors show this condition can be relaxed to the output dimension being larger than the number of neurons. This is a strong assumption, as in practice, the output dimension is usually 1 for many regression problems or the number of classes for classification problems. Furthermore, this assumption is actually crucial for the algorithm proposed in this paper. If the output dimension is small, then the “Pure Neuron Detection” step does work. Please clarify if I understand incorrectly. If this is indeed the case, I suggest discussing this strong assumption in the main text and listing the problem of relaxing it as an open problem. Minor comments: 1. I suggest adding the following papers to the related work section in the final version: https://arxiv.org/abs/1805.06523 https://arxiv.org/abs/1810.02054 https://arxiv.org/abs/1810.04133 https://arxiv.org/abs/1712.00779 These paper are relatively new but very relevant. 2. There are many typos in the references. For example, “relu” should be ReLU.