Learning Two-layer Neural Networks with Symmetric Inputs

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-layerThere are different ways to count the number of layers. Here by two-layer network we refer to a fully-connected 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 and Klivans, (2017)) or a one-hidden layer network if one just counts the number of hidden layers. 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 $\mathcal{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 $\mathcal{D}$. Learning a neural network is known to be a hard problem even in some simple settings (Goel et al.,, 2016; Brutzkus and Globerson,, 2017), so we need to make assumptions on the network structure or the input distribution $\mathcal{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 and Globerson,, 2017; Li and 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

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

If the data is generated according to Equation (1), and the input distribution $x\sim\mathcal{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 $\mbox{poly}(d)$ time and outputs a network $\hat{A},\hat{W}$ of the same size that is effectively the same as the ground-truth network: for any input $x$, $\hat{A}\sigma(\hat{W}x)=A\sigma(Wx)$.

Of course, in practice we only have samples of $(x,y)$ and cannot get the exact correlations. However, our algorithm is robust to perturbations, and in particular can work with polynomially many samples.

If the data is generated according to Equation (1), and the input distribution $x\sim\mathcal{D}$ is symmetric. As long as the weight matrices $A,W$ and input distributions are not degenerate, there is an algorithm that uses $\mbox{poly}(d,1/\epsilon)$ time and number of samples and outputs a network $\hat{A},\hat{W}$ of the same size that computes an $\epsilon$-approximation function to the ground-truth network: for any input $x$, $\|\hat{A}\sigma(\hat{W}x)-A\sigma(Wx)\|^{2}\leq\epsilon$.

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 and 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 $\mathcal{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

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 $\lambda$, and our algorithm outputs an $\epsilon$-approximation to the perturbed network with $\mbox{poly}(d,1/\lambda,1/\epsilon)$ 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.

2 Related Work

There are many works in learning neural networks, and they come in many different styles.

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 and Soltanolkotabi, (2018) used the rank-1 tensor decomposition for learning a non-overlapping convolutional neural network with differentiable and smooth activation and 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 and Globerson,, 2017; Li and Yuan,, 2017; Soltanolkotabi,, 2017; Zhong et al.,, 2017) extend this to different settings.

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

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 and Klivans, (2017); Du and 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.

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.

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 Notations

We use $[n]$ to denote the set $\{1,2,\cdots,n\}$. For two random variables $X$ and $Y$, we say $X\stackrel{{\scriptstyle d}}{{=}}Y$ if they come from the same distribution.

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:

Another related concept is the augmented distinguishing matrix $M$, which is a $d^{2}\times\left(k_{2}+1\right)$ matrix whose first $k_{2}$ columns are exactly the same as distinguishing matrix $N$, and the last column (indexed by ) is defined as

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 $\sigma_{min}(M)$ is lowerbounded. Intuitively, every column $N_{ij}^{\mathcal{D}}$ looks at the expectation over samples that have opposite signs for weights $w_{i},w_{j}$ ($w_{i}^{\top}xw_{j}^{\top}x\leq 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 $\sigma_{min}(M)$ is not unreasonable: in the smoothed analysis setting where the nature can make a small perturbation on the input distribution $\mathcal{D}$, we show that for any input distribution $\mathcal{D}$, there exists simple perturbations $\mathcal{D^{\prime}}$ that are arbitrarily close to $\mathcal{D}$ such that $\sigma_{min}(M^{D^{\prime}})$ is lowerbounded.

Our Algorithm

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}\sim\mathcal{D}$ comes from a symmetric distribution, and the output $y_{i}$ is computed by

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:

Suppose $x\sim\mathcal{D}$ comes from a symmetric distribution and $y$ is computed as in (3), then

Importantly, the right hand side of Lemma 1 does not contain the ReLU function $\sigma$. This is true because if $x$ comes from a symmetric distribution, averaging between $x$ and $-x$ can get rid of non-linearities like ReLU or leaky-ReLU. Later we will prove a more general version of this lemma (Lemma 6).

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

The key observation here is that if $u=z_{i}$, then $u^{\top}A=\lambda_{i}e_{i}^{\top}$. As a result, $u^{\top}y=\lambda_{i}e_{i}^{\top}\sigma(Wx)=\sigma(\lambda_{i}w_{i}^{\top}x)$ is the output of a single-layer neural network with weight equal to $\lambda_{i}w_{i}$. If we know all the vectors $\{z_{1},...,z_{k}\}$, the input/output pairs $(x,z_{i}^{\top}y)$ correspond to single-layer networks with weight vectors $\{\lambda_{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^{\top}A=\lambda_{i}e_{i}$, we say that $u^{\top}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.

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:

The additional terms may accidentally cancel each other which leads to a false positive. To address this problem, we consider a higher order moment:

Suppose $\hat{y}=\sigma(w^{\top}x)$, then

Here $N_{ij}$’s are columns of the distinguishing matrix defined in Definition 1.

We will call the function $f(u)$ a pure neuron detector, as $u^{\top}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$.

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 $\{\mbox{vec}^{*}(z_{i}z_{i}^{\top})\}$.

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}^{\top}$’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=\mbox{diag}(\lambda)A^{-1}$. Let $X=Z^{\top}D_{X}Z$ and $Y=Z^{\top}D_{Y}Z$ be two random elements in the span of $z_{i}z_{i}^{\top}$, 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^{\top}D_{X}D_{Y}^{-1}(Z^{\top})^{-1}$, since $z_{i}$ is a column of $Z^{\top}$, we know

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}$.

Given the span of $z_{i}z_{i}^{\top}$’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}$.

3 Detailed Algorithm and Guarantees

We can now give the full algorithm, see Algorithm 2. The main steps of this algorithm is as explained in the previous section. Steps 2 - 5 constructs the pure neuron detector and finds the span of $\mbox{vec}^{*}(z_{i}z_{i}^{\top})$ (as in Corollary 1); Steps 7 - 9 performs simultaneous diagonalization to get all the $z_{i}$’s; Steps 11, 12 calls Algorithm 1 to solve the single-layer problem and outputs the correct result.

We are now ready to state a formal version of Theorem 1:

It is easy to prove this theorem using the lemmas we have.

Now the output $z_{i}^{\top}y=\lambda_{i}\sigma(w_{i}^{\top}x)=\sigma(\lambda_{i}w_{i}^{\top}x)$ (again by property of ReLU function $\sigma$), by the design of Algorithm 1 we know $v_{i}=\lambda_{i}w_{i}$. We also know that $Z=\mbox{diag}(\lambda)A^{-1}$, therefore $Z^{-1}=A\mbox{diag}(\lambda)^{-1}$. Notice that $Z^{-1}\sigma(Vx)=A\mbox{diag}(\lambda)^{-1}\sigma(\mbox{diag}(\lambda)Wx)=A\sigma(Wx).$ These two scaling factors cancel each other, so the two networks compute the same function. ∎

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\times 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.

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.

2 Robustness to Noise

Figure 3 demonstrates the robustness of our algorithm to label noise $\xi$ for Gaussian and symmetric mixture of Gaussians input distributions. In this experiment, we fix the size of both $A$ and $W$ to be $10\times 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.

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\Sigma V^{\top}$ such that $U$ and $V$ are random orthonormal matrices and $\Sigma_{ii}=\lambda^{-i}$, where $\lambda$ 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.

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.

References

Appendix 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=\sigma(Wx)+\xi$) 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).

where the expectation is taken over the input distribution.

There are two cases to consider: $p$ and $q$ are both even numbers or both odd numbers.

For the case where $p$ and $q$ are even numbers, we have

If $(a^{\top}x)\leq 0$, we know \big{(}\sigma(-a^{\top}x)\big{)}^{p}+\big{(}\sigma(a^{\top}x)\big{)}^{p}=(a^{\top}x)^{p}+0=(a^{\top}x)^{p}. Otherwise, we have \big{(}\sigma(-a^{\top}x)\big{)}^{p}+\big{(}\sigma(a^{\top}x)\big{)}^{p}=0+(a^{\top}x)^{p}=(a^{\top}x)^{p}. Thus,

For the other case where $p$ and $q$ are odd numbers, we have

Similarly, if $(a^{\top}x)\leq 0$, we know -\big{(}\sigma(-a^{\top}x)\big{)}^{p}+\big{(}\sigma(a^{\top}x)\big{)}^{p}=-(-a^{\top}x)^{p}+0=(a^{\top}x)^{p}. Otherwise, we have -\big{(}\sigma(-a^{\top}x)\big{)}^{p}+\big{(}\sigma(a^{\top}x)\big{)}^{p}=0+(a^{\top}x)^{p}=(a^{\top}x)^{p}. Thus,

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.

where the second equality holds due to Lemma 6.

where (7) uses (9) of the following Lemma 7, and (8) uses the definition of distinguishing matrix $N$ (Definition 1). $\Box$

where the expectation is taken over the input distribution.

Since input $x$ comes from a symmetric distribution, we have

When $a^{\top}xb^{\top}x>0$, we know $a^{\top}x$ and $b^{\top}x$ are both positive or both negative. In either case, we know that $\sigma(a^{\top}x)\sigma(b^{\top}x)+\sigma(-a^{\top}x)\sigma(-b^{\top}x)=(a^{\top}x)(b^{\top}x)$. Thus, we have

It is not hard to verify that $u^{\top}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 $\binom{k}{2}+k=k_{2}+k$, i.e.,

Now, we prove the Lemma 4 which shows the null space of $T$ is exactly the span of $\{\mbox{vec}^{*}(z_{i}z_{i}^{\top})\}$.

Proof of Lemma 4. We divide the proof to the following two cases:

For any vector $\mbox{vec}^{*}(U)$ belongs to the null space of $T$, we have $T\mbox{vec}^{*}(U)=0$. Note that the RHS of (10) equals to 0 if and only if $A^{\top}UA$ is a diagonal matrix since the distinguishing matrix $N$ is full column rank and $A^{\top}UA$ is symmetric. Thus $\mbox{vec}^{*}(U)$ belongs to the span of $\{\mbox{vec}^{*}(z_{i}z_{i}^{\top})\}$ since $U=Z^{\top}DZ$ for some diagonal matrix $D$.

For any vector $\mbox{vec}^{*}(U)$ belonging to the span of $\{\mbox{vec}^{*}(z_{i}z_{i}^{\top})\}$, $U$ is a linear combination of $z_{i}z_{i}^{\top}$’s. Furthermore, $T\mbox{vec}^{*}(U)$ is a linear combination of $T\mbox{vec}^{*}(z_{i}z_{i}^{\top})$. Note that $A^{\top}z_{i}$ only has one non-zero entry due to the definition of $z_{i}$, for any $i\in[k]$. Thus all coefficients in the RHS of (10) are 0. We get $T\mbox{vec}^{*}(U)=0$.

Finding $z_{i}$’s: Now, we prove the final Lemma 5 which finds all $z_{i}$’s from the span of $\{\mbox{vec}^{*}(z_{i}z_{i}^{\top})\}$ 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^{\top}D_{X}D_{Y}^{-1}(Z^{\top})^{-1}$. According to the following Lemma 8 (i.e., all diagonal elements of $D_{x}D_{y}^{-1}$ 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}$. $\Box$

With probability $1$, all diagonal elements of $D_{X}$ and $D_{Y}$ are non-zero and all diagonal elements of $D_{X}D_{Y}^{-1}$ are distinct, where $X=Z^{\top}D_{X}Z$ and $Y=Z^{\top}D_{Y}Z$ are defined in Line 7 of Algorithm 2.

A.2 Noisy Case

Now, we discuss how to handle the noisy case (i.e. $y=\sigma(Wx)+\xi$). The corresponding algorithm is described in Algorithm 3. Note that the noise $\xi$ only affects the first two steps, i.e., pure neuron detector (Lemma 3) and finding span of $\mbox{vec}^{*}(z_{i}z_{i}^{\top})$ (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 Lemma 5 is independent of the model and Lemma 1 is linear wrt. noise $\xi$, 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 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):

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

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. $\xi$. Thus, we have

where the third equality holds due to Lemma 6 and Lemma 7, and (16) uses the definition of $m_{ij}$.

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

where (17) uses (9) (same as (7)). $\Box$

A.3 Extension to Non-square A𝐴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\times 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.

For a ground truth neural network with weight matrices $W$ and $P^{\top}A$, the generated sample will just be $(x,P^{\top}y)$. According to Theorem 5, we know for any input $x$, we have $Z^{-1}\sigma(Vx)=P^{\top}A\sigma(Wx)$. Thus, we have

where the second equality holds since $PP^{\top}$ is just the projection matrix to the column span of $A$. ∎

Appendix 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 $\xi=0$. Throughout this section, we will use $\hat{V},\hat{Z}^{-1}$ to denote the results of Algorithm 3 with empirical moments, and use $V,Z^{-1}$ for the results when the algorithm has access to exact moments, similarly for other intermediate results. For the robustness of Algorithm 3, we prove the following theorem.

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 $\mbox{vec}^{*}(z_{i}z_{i}^{\top})$’s; finding $z_{i}$’s from the span of $\mbox{vec}^{*}(z_{i}z_{i}^{\top})$’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 $\mbox{vec}^{*}(z_{i}z_{i}^{\top})$’s in arbitrary accuracy. Let $\hat{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 $\mbox{vec}^{*}(z_{i}z_{i}^{\top})$’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 $\hat{T}$ is close to the null space of $T$.

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-\hat{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 $\mbox{vec}^{*}(z_{i}z_{i}^{\top})$’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.

Finally, given $\hat{z}_{i}$’s, the problem reduces to a one-layer problem. We will first give an analysis for Algorithm 1 as a warm-up. When we call Algorithm 1 from Algorithm 3, the situation is slightly different. Note we reserve fresh samples for this step, so that the samples used by Algorithm 1 are still independent with the estimate $\hat{z}_{i}$ (learned using the other set of samples). However, since $\hat{z}_{i}$ is not equal to $z_{i}$, this introduces an additional error term $(\hat{z}_{i}-z_{i})^{\top}y$ which is not independent of $x$ and cannot be captured by $\xi$. We modify the proof for Algorithm 1 to show that the algorithm is still robust as long as $\|\hat{z}_{i}-z_{i}\|$ is small enough.

Combining the above three lemmas, we prove Theorem 7 in Section B.4.

We first prove that the step of finding the span of $\{\mbox{vec}^{*}(z_{i}z_{j}^{\top})\}$ 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.

Suppose $\sigma_{\min}(M)\geq\alpha,\sigma_{\min}(A)\geq\beta$, we know that matrix $T$ has rank $k_{2}$ and the $k_{2}$-th singular value of $T$ is lower bounded by $\alpha\beta^{2}$.

Then we show that with enough samples the estimate $\hat{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.

Finally we combine the above two lemmas and show that the span of the least $k$ right singular vectors of $\hat{T}$ is close to the null space of $T$.

According to Lemma 15, given $O(\frac{d^{3}\Gamma^{14}(\Gamma P_{1}\sqrt{k}+P_{2})^{6}\log(\frac{d}{\delta})}{\gamma^{4}\epsilon^{2}})$ number of i.i.d. samples, we know with probability at least $1-\delta$,

According to Lemma 14, we know $\sigma_{k_{2}}(T)\geq\alpha\beta^{2}.$ Then, due to Lemma 27, we have

for any symmetric $k\times k$ matrix $U$.

Note that $\mbox{vec}^{*}(U)$ is a $(k_{2}+k)$-dimensional vector. For convenience, we first use matrix $F$ to transform $\mbox{vec}^{*}(U)$ to $\mbox{vec}(U)$, which has $k^{2}$ dimensions. Matrix $F$ is defined such that $F\mbox{vec}^{*}(U)=\mbox{vec}(U)$, for any $k\times 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^{\top}UA)_{ij}$’s. Notice that

Since we only care about the elements of $A^{\top}UA$ at the $ij$-th position for $1\leq i<j\leq k$, we just pick corresponding rows of $A^{\top}\otimes A^{\top}$ to construct our matrix $C$, which has dimension $k_{2}\times 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\mbox{vec}^{*}(U)$ in the following way:

With above characterization of $T$, we are ready to show that the $k_{2}$-th singular value of $T$ is lower bounded.

Since matrix $C$ has dimension $k_{2}\times 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_{k_{2}}$, 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 $MB$ 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=MBCF$ also has rank $k_{2}$. Next, we prove that $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 $\alpha$. Since the first $k_{2}$ rows of matrix $B$ constitute a $k_{2}\times k_{2}$ identity matrix, we know

where $u$ is any $k_{2}$-dimensional vector.

Since $\sigma_{\min}(A)\geq\beta$, we know $\sigma_{\min}(A^{\top}\otimes A^{\top})\geq\beta^{2}$. According to the construction of $C$, we know $C$ consists a subset of rows of $A^{\top}\otimes A^{\top}$. Denote the indices of the row not picked as $S$. We have

where $u$ has dimension $k_{2}$ and $v$ has dimension $k^{2}$.

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 $\sigma_{\min}^{+}(T)$, we have

where the first inequality holds because both $M$ and $B$ has full column rank. ∎

In this section, we prove that given polynomial number of samples, $\|\hat{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).

In order to get an upper bound for $\|\hat{T}-T\|_{F}$, we first show that $\|\hat{T}-T\|_{2}$ is upper bounded. We know

For any $k\times k$ symmetric matrix $U$ with eigenvalue decomposition $U=\sum_{i=1}^{k}\lambda_{i}u^{(i)}(u^{(i)})^{\top}$, according to the definition of $T$, we know

where $\hat{f}(u)=\hat{T}\mbox{vec}^{*}(uu^{\top})$ and the fourth inequality uses the Cauchy-Schwarz inequality. Next, we only need to upper bound $\max_{u:\|u\|\leq 2}\|f(u)-\hat{f}(u)\|$. Recall that

We first show that given polynomial number of samples,

Since each row of $W$ has unit norm, we have $\|W\|\leq\sqrt{k}$. Due to the assumption that $\|x\|\leq\Gamma,\|A\|\leq P_{1},\|\xi\|\leq P_{2}$, we have

According to Lemma 24, we know given $O(\frac{\Gamma^{2}(\Gamma P_{1}\sqrt{k}+P_{2})^{2}\log(\frac{d}{\delta})}{\epsilon^{2}})$ number of samples,

Similarly, we can show that given $O(\frac{\Gamma^{6}(\Gamma P_{1}\sqrt{k}+P_{2})^{2}\log(\frac{d}{\delta})}{\epsilon^{2}})$ number of samples,

Since $\|xx^{\top}\|\leq\Gamma^{2}$, we know that given $O(\frac{\Gamma^{4}\log(\frac{d}{\delta})}{\epsilon^{2}})$ number of samples,

By union bound, we know for any $\epsilon<\gamma/2,$ given $O(\frac{\Gamma^{6}(\Gamma P_{1}\sqrt{k}+P_{2})^{2}\log(\frac{d}{\delta})}{\epsilon^{2}})$ number of samples, with probability at least $1-\delta$, we have

Thus, given $O(\frac{\Gamma^{14}(\Gamma P_{1}\sqrt{k}+P_{2})^{6}\log(\frac{d}{\delta})}{\gamma^{4}\epsilon^{2}})$ number of samples, we know

Since \big{\|}(u^{\top}y)^{2}(x\otimes x)\big{\|}\leq 4\Gamma^{2}(\Gamma P_{1}\sqrt{k}+P_{2})^{2}, according to Lemma 24, we know given $O(\frac{\Gamma^{4}(\Gamma P_{1}\sqrt{k}+P_{2})^{4}\log(\frac{d}{\delta})}{\epsilon^{2}}),$

Again, using Lemma 24 and union bound, we know given $O(\frac{\Gamma^{2}(\Gamma P_{1}\sqrt{k}+P_{2})^{2}\log(\frac{d}{\delta})}{\epsilon^{2}})$ number of samples, we have

Thus, we know that given $O(\frac{\Gamma^{2}(\Gamma P_{1}\sqrt{k}+P_{2})^{6}\log(\frac{d}{\delta})}{\epsilon^{2}})$ number of samples, we know

Similar as the first term, we can show that given $O(\frac{\Gamma^{10}(\Gamma P_{1}\sqrt{k}+P_{2})^{6}\log(\frac{d}{\delta})}{\gamma^{4}\epsilon^{2}})$ number of samples, we have

Now, we are ready to combine our bound for each of four terms. By union bound, we know given $O(\frac{\Gamma^{14}(\Gamma P_{1}\sqrt{k}+P_{2})^{6}\log(\frac{d}{\delta})}{\gamma^{4}\epsilon^{2}})$ number of samples,

hold with probability at least $1-\delta$. Thus, we know

where the second inequality holds since $\|\hat{T}-T\|\leq\sqrt{2k}\max_{u:\|u\|\leq 2}\|f(u)-\hat{f}(u)\|.$

Thus, we know given $O(\frac{\Gamma^{14}(\Gamma P_{1}\sqrt{k}+P_{2})^{6}\log(\frac{d}{\delta})}{\gamma^{4}\epsilon^{2}})$ number of samples,

with probability at least $1-\delta$. Thus, given $O(\frac{d^{3}\Gamma^{14}(\Gamma P_{1}\sqrt{k}+P_{2})^{6}\log(\frac{d}{\delta})}{\gamma^{4}\epsilon^{2}})$ number of samples,

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 $\hat{S}$ be two $(k_{2}+k)$ by $k$ matrices, whose columns consist of the least $k$ right singular vectors of $T$ and $\hat{T}$ respectively.

According to Lemma 11, we know with polynomial number of samples, the Frobenius norm of $SS^{\top}-\hat{S}\hat{S}^{\perp}$ is well bounded. However, due to the rotation issue of subspace basis, we cannot conclude that $\|S-\hat{S}\|_{F}$ is small. Only after appropriate alignment, the difference between $S$ and $\hat{S}$ becomes small.

Since $S$ has orthonormal columns, we have $\sigma_{k}(SS^{\top})=1$. Then, according to Lemma 35, we know there exists rotation matrix $R$ such that

Let the $k$ columns of $S$ be $\mbox{vec}^{*}(U_{1}),\mbox{vec}^{*}(U_{2}),\cdots,\mbox{vec}^{*}(U_{k})$. Note each $U_{i}$ can be expressed as $A^{-\top}D_{i}A^{-1}$, where $D_{i}$ is a diagonal matrix. Let $Q$ be a $k\times 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 $\mbox{vec}^{*}(X)=SR\zeta_{1},\mbox{vec}^{*}(Y)=SR\zeta_{2}$, where $R$ is the rotation matrix in Lemma 16 and $\zeta_{1},\zeta_{2}$ are two independent standard Gaussian vectors. Let $D_{X}=\mbox{diag}(QR\zeta_{1})$ and $D_{Y}=\mbox{diag}(QR\zeta_{2})$. It’s not hard to check that $X=A^{-\top}D_{X}A^{-1}$ and $Y=A^{-\top}D_{Y}A^{-1}$. Furthermore, we have $XY^{-1}=A^{-\top}D_{X}D_{Y}^{-1}A^{\top}$. Next, we show that the diagonal elements of $D_{X}D_{Y}^{-1}$ are well separated.

Assume that $\|A\|\leq P_{1},\sigma_{\min}(A)\geq\beta$. Then for any $\delta>0$ we know with probability at least $1-\delta$, we have

where $\mbox{sep}(D_{X}D_{Y}^{-1}):=\min_{i\neq j}|(D_{X}D_{Y}^{-1})_{ii}-(D_{X}D_{Y}^{-1})_{jj}|$.

We first show that matrix $Q$ is well-conditioned. Since $U_{i}=A^{-\top}D_{i}A^{-1}$, we have $\mbox{vec}(U_{i})=A^{-\top}\otimes A^{-\top}\mbox{vec}(D_{i}).$ Let $U$ be a $k^{2}\times k$ matrix whose columns consist of $\mbox{vec}(U_{i})$’s. Also define $\bar{Q}$ as a $k^{2}\times k$ matrix whose columns are $\mbox{vec}(D_{i})$’s. Note that matrix $\bar{Q}$ only has $k$ non-zero rows, which are exactly matrix $Q$. With the above definition, we have $U=A^{-\top}\otimes A^{-\top}\bar{Q}$. Since $\sigma_{\min}(U)\leq\|A^{-\top}\otimes A^{-\top}\|\sigma_{\min}(\bar{Q})$, we have

Notice that a subset of rows of $U$ constitute matrix $S$, which is an orthonormal matrix. Thus, we have $\sigma_{\min}(U)\geq\sigma_{\min}(S)=1$. Since we assume $\sigma_{\min}(A)\geq\beta$, we have

Thus, we have $\sigma_{\min}(\bar{Q})\geq\beta^{2}$, which implies $\sigma_{\min}(Q)\geq\beta^{2}$.

We also know $\|U\|\geq\sigma_{\min}(A^{-\top}\otimes A^{-\top})\|\bar{Q}\|$, thus

Since $\|S\|=1$, we know $\|U\|\leq\sqrt{2}.$ For the smallest singular value of $A^{-\top}\otimes A^{-\top}$, we have

Thus, we have $\|\bar{Q}\|\leq\sqrt{2}P_{1}^{2}$, which implies $\|Q\|\leq\sqrt{2}P_{1}^{2}$.

Now, let’s prove that the diagonal elements of $D_{X}D_{Y}^{-1}$ are well-separated. Let $q_{i}^{\top}$ be the $i$-th row vector of $Q$. Then we know the $i$-th diagonal element of $D_{X}D_{Y}^{-1}$ is $\frac{\langle q_{i},R\zeta_{1}\rangle}{\langle q_{i},R\zeta_{2}\rangle}$. Since $\|Q\|\leq\sqrt{2}P_{1}^{2}$, we have $\|q_{i}\|\leq\sqrt{2}P_{1}^{2}$ for every row vector.

It’s not hard to show that with probability at least $1-\exp(-d^{\Omega(1)})$, we have $|{\langle q_{i},R\zeta_{2}\rangle}|\leq\mbox{poly}(d,P_{1})$ for each $i$. Now given $\zeta_{2}$ for which this happens, we have $\frac{\langle q_{i},R\zeta_{1}\rangle}{\langle q_{i},R\zeta_{2}\rangle}-\frac{\langle q_{j},R\zeta_{1}\rangle}{\langle q_{j},R\zeta_{2}\rangle}=c_{i}\langle q_{i},R\zeta_{1}\rangle-c_{j}\langle q_{j},R\zeta_{1}\rangle$, where $c_{i},c_{j}$ have magnitude as least $\mbox{poly}(1/d,1/P_{1})$. Since $\sigma_{min}(Q)\geq\beta^{2}$, we know $\|\mbox{Proj}_{q_{j}^{\perp}}q_{i}\|\geq\beta^{2}$ (because otherwise there exists $\lambda$ such that $\|(e_{i}+\lambda e_{j})\|\sigma_{min}(Q)\leq\|(e_{i}+\lambda e_{j})^{\top}Q\|=\|\mbox{Proj}_{q_{j}^{\perp}}q_{i}\|<\beta^{2}$, which is a contradiction). Let $q_{i,j}^{\perp}=\mbox{Proj}_{q_{j}^{\perp}}q_{i}=q_{i}-\lambda_{i,j}q_{j}$, we can rewrite this as

By properties of Gaussians, we know $\langle q_{i,j}^{\perp},R\zeta_{1}\rangle$ is independent of $\langle q_{j},R\zeta_{1}\rangle$, so we can first fix $\langle q_{j},R\zeta_{1}\rangle$ and apply anti-concentration of Gaussians (see Lemma 38) to $\langle q_{i,j}^{\perp},R\zeta_{1}\rangle$. As a result we know with probability at least $1-\delta/k^{2}$:

By union bound, we know with probability at least $1-\delta,$

Let $\hat{X}=\hat{S}\zeta_{1}$ and $\hat{Y}=\hat{S}\zeta_{2}$. Next, we prove that the eigenvectors of $\hat{X}\hat{Y}^{-1}$ are close to the eigenvectors of $XY^{-1}$.

Let $z_{1}^{\prime},\cdots,z_{k}^{\prime}$ be the eigenvectors of $XY^{-1}$ (before sign flip step). Similarly define $\hat{z}_{1}^{\prime},\cdots,\hat{z}_{k}^{\prime}$ for $\hat{X}\hat{Y}^{-1}$. We first prove that the eigenvectors of $\hat{X}\hat{Y}^{-1}$ are close to the eigenvectors of $XY^{-1}$.

Let $\hat{X}=X+E_{X}$ and $\hat{Y}=Y+E_{Y}$. Then we have

where $F=-E_{Y}(I+Y^{-1}E_{Y})^{-1}Y^{-1}$ and $G=E_{X}\hat{Y}^{-1}.$ According to Lemma 34, we have $\|F\|\leq\frac{\|E_{Y}\|}{\sigma_{\min}(Y)-\|E_{Y}\|}$ and $\|G\|\leq\frac{\|E_{X}\|}{\sigma_{\min}(\hat{Y})}$. In order to bound the perturbation matrices $\|F\|$ and $\|G\|$, we need to first bound $\|E_{X}\|,\|E_{Y}\|$ and $\sigma_{\min}(Y),\sigma_{\min}(\hat{Y})$.

As we know, $E_{X}=\hat{X}-X=(\hat{S}-SR)\zeta_{1}$. According to Lemma 16, we have $\|\hat{S}-SR\|\leq\|\hat{S}-SR\|_{F}\leq 2\epsilon.$ We also know with probability at least $1-\exp(-d^{\Omega(1)}),$ $\|\zeta_{1}\|\leq\mbox{poly}(d).$ Thus, we have

Similarly, with probability at least $1-\exp(-d^{\Omega(1)})$, we also know $\|E_{Y}\|\leq\mbox{poly}(\epsilon,d).$

Now, we lower bound the smallest singular value of $X$ and $Y$. Since $X=A^{-\top}D_{X}A^{-1}$, we have

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 $\langle q_{i},R\zeta_{1}\rangle$, which follows a Gaussian distribution whose standard deviation is at least $\|q_{i}\|\geq\beta^{2}$. By anti-concentration property of Gaussian (see Lemma 38), we know $|{\langle q_{i},R\zeta_{2}\rangle}|\geq\Omega(\delta\beta^{2}/k)$ for all $i$ with probability $1-\delta/4$. Thus, we have $\sigma_{\min}(X)\geq\mbox{poly}(1/d,1/P_{1},\beta,\delta)$. Similarly we have the same conclusion for $Y$. For small enough $E_{Y}$, we have $\sigma_{\min}(\hat{Y})\geq\sigma_{\min}(Y)-\|E_{Y}\|.$

Thus, for small enough $\epsilon$, we have $\|F\|\leq\mbox{poly}(d,P_{1},1/\beta,\epsilon,1/\delta)$ and $\|G\|\leq\mbox{poly}(d,P_{1},1/\beta,\epsilon,1/\delta)$. In order to apply Lemma 33, we also need to bound $\kappa(A^{-\top})$ and $\|XY^{-1}\|$. Since $\sigma_{\min}(A^{-\top})\geq\frac{1}{P_{1}}$ and $\|A^{-\top}\|\leq 1/\beta$, we have $\kappa(A^{-\top})\leq P_{1}/\beta.$ For the norm of $XY^{-1}$, we have

where the second inequality holds because $\sigma_{\min}(Y)\geq\mbox{poly}(1/d,1/P_{1},\beta,\delta).$ Recall that $X=\mbox{mat}^{*}(SR\zeta_{1})$. It’s not hard to verify that with probability at least $1-\exp(-d^{\Omega(1)})$, we have $\|X\|\leq\mbox{poly}(d)$. Thus, we know $\|XY^{-1}\|\leq\mbox{poly}(d,P_{1},1/\beta,1/\delta)$. Similarly, we can also prove that $\|D_{X}D_{Y}^{-1}\|\leq\mbox{poly}(d,P_{1},1/\beta,1/\delta)$

According to Lemma 17, we know with probability at least $1-\delta/4$, $\text{sep}(D_{X}D_{Y}^{-1})\geq\mbox{poly}(1/d,1/P_{1},\beta,\delta).$ Thus, by union bound, we know for small enough $\epsilon$, with probability at least $1-\delta$,

According to Lemma 33, we know there exists a permutation $\pi[i]\in[k]$, such that

In the following lemma, we show that the sign flip step of $\hat{z}_{i}$ is robust.

Suppose that $\|x\|\leq\Gamma,\|A\|\leq P_{1},\|\xi\|\leq P_{2}.$ Let $z_{1}^{\prime},\cdots,z_{k}^{\prime}$ be the eigenvectors of $XY^{-1}$ (before sign flip step). Similarly define $\hat{z}_{1}^{\prime},\cdots,\hat{z}_{k}^{\prime}$ for $\hat{X}\hat{Y}^{-1}$. Suppose for each $i$, $\|z_{i}^{\prime}-\hat{z}_{i}^{\prime}\|\leq\epsilon$, where $\epsilon\leq\frac{\beta\gamma}{4\Gamma(1+\Gamma P_{1}\sqrt{k}+P_{2})}$. We know, for any $\delta<1$, with $O(\frac{(\Gamma P_{1}\sqrt{k}+P_{2})^{2}\log(d/\delta)}{\epsilon^{2}})$ number of i.i.d. samples,