Recovery Guarantees for One-hidden-layer Neural Networks

Introduction

Neural Networks (NNs) have achieved great practical success recently. Many theoretical contributions have been made very recently to understand the extraordinary performance of NNs. The remarkable results of NNs on complex tasks in computer vision and natural language processing inspired works on the expressive power of NNs [CSS16, CS16, RPK+16, DFS16, PLR+16, MPCB14, Tel16]. Indeed, several works found NNs are very powerful and the deeper the more powerful. However, due to the high non-convexity of NNs, knowing the expressivity of NNs doesn’t guarantee that the targeted functions will be learned. Therefore, several other works focused on the achievability of global optima. Many of them considered the over-parameterized setting, where the global optima or local minima close to the global optima will be achieved when the number of parameters is large enough, including [FB16, HV15, LSSS14, DPG+14, SS16, HM17]. This, however, leads to overfitting easily and can’t provide any generalization guarantees, which are actually the essential goal in most tasks.

A few works have considered generalization performance. For example, [XLS17] provide generalization bound under the Rademacher generalization analysis framework. Recently [ZBH+17] describe some experiments showing that NNs are complex enough that they actually memorize the training data but still generalize well. As they claim, this cannot be explained by applying generalization analysis techniques, like VC dimension and Rademacher complexity, to classification loss (although it does not rule out a margins analysis—see, for example, [Bar98]; their experiments involve the unbounded cross-entropy loss).

In this paper, we don’t develop a new generalization analysis. Instead we focus on parameter recovery setting, where we assume there are underlying ground-truth parameters and we provide recovery guarantees for the ground-truth parameters up to equivalent permutations. Since the parameters are exactly recovered, the generalization performance will also be guaranteed.

Several other techniques are also provided to recover the parameters or to guarantee generalization performance, such as tensor methods [JSA15] and kernel methods [AGMR17]. These methods require sample complexity $O(d^{3})$ or computational complexity $\widetilde{O}(n^{2})$, which can be intractable in practice. We propose an algorithm that has recovery guarantees for 1NN with sample complexity $\widetilde{O}(d)$ and computational time $\widetilde{O}(dn)$ under some mild assumptions.

Recently [Sha16] show that neither specific assumptions on the niceness of the input distribution or niceness of the target function alone is sufficient to guarantee learnability using gradient-based methods. In this paper, we assume data points are sampled from Gaussian distribution and the parameters of hidden neurons are linearly independent.

We distill some properties for activation functions, which are satisfied by a wide range of activations, including ReLU, squared ReLU, sigmoid and tanh. With these properties we show positive definiteness (PD) of the Hessian in the neighborhood of the ground-truth parameters given enough samples (Theorem 4.2). Further, for activations that are also smooth, we show local linear convergence is guaranteed using gradient descent.

We propose a tensor method to initialize the parameters such that the initialized parameters fall into the local positive definiteness area. Our contribution is that we reduce the sample/computational complexity from cubic dependency on dimension to linear dependency (Theorem 5.6).

Combining the above two results, we provide a globally converging algorithm (Algorithm 2) for smooth homogeneous activations satisfying the distilled properties. The whole procedure requires sample/computational complexity linear in dimension and logarithmic in precision (Theorem 6.1).

Related Work

The recent empirical success of NNs has boosted their theoretical analyses [FZK+16, Bal16, BMBY16, SBL16, APVZ14, AGMR17, GKKT17]. In this paper, we classify them into three main directions.

Expressive power is studied to understand the remarkable performance of neural networks on complex tasks. Although one-hidden-layer neural networks with sufficiently many hidden nodes can approximate any continuous function [Hor91], shallow networks can’t achieve the same performance in practice as deep networks. Theoretically, several recent works show the depth of NNs plays an essential role in the expressive power of neural networks [DFS16]. As shown in [CSS16, CS16, Tel16], functions that can be implemented by a deep network of polynomial size require exponential size in order to be implemented by a shallow network. [RPK+16, PLR+16, MPCB14, AGMR17] design some measures of expressivity that display an exponential dependence on the depth of the network. However, the increasing of the expressivity of NNs or its depth also increases the difficulty of the learning process to achieve a good enough model. In this paper, we focus on 1NNs and provide recovery guarantees using a finite number of samples.

2 Achievability of Global Optima

The global convergence is in general not guaranteed for NNs due to their non-convexity. It is widely believed that training deep models using gradient-based methods works so well because the error surface either has no local minima, or if they exist they need to be close in value to the global minima. [SCP16] present examples showing that for this to be true additional assumptions on the data, initialization schemes and/or the model classes have to be made. Indeed the achievability of global optima has been shown under many different types of assumptions.

In particular, [CHM+15] analyze the loss surface of a special random neural network through spin-glass theory and show that it has exponentially many local optima, whose loss is small and close to that of a global optimum. Later on, [Kaw16] eliminate some assumptions made by [CHM+15] but still require the independence of activations as [CHM+15], which is unrealistic. [SS16] study the geometric structure of the neural network objective function. They have shown that with high probability random initialization will fall into a basin with a small objective value when the network is over-parameterized. [LSSS14] consider polynomial networks where the activations are square functions, which are typically not used in practice. [HV15] show that when a local minimum has zero parameters related to a hidden node, a global optimum is achieved. [FB16] study the landscape of 1NN in terms of topology and geometry, and show that the level set becomes connected as the network is increasingly over-parameterized. [HM17] show that products of matrices don’t have spurious local minima and that deep residual networks can represent any function on a sample, as long as the number of parameters is larger than the sample size. [SC16] consider over-specified NNs, where the number of samples is smaller than the number of weights. [DPG+14] propose a new approach to second-order optimization that identifies and attacks the saddle point problem in high-dimensional non-convex optimization. They apply the approach to recurrent neural networks and show practical performance. [AGMR17] use results from tropical geometry to show global optimality of an algorithm, but it requires $(2n)^{k}\operatorname{poly}(n)$ computational complexity.

Almost all of these results require the number of parameters is larger than the number of points, which probably overfits the model and no generalization performance will be guaranteed. In this paper, we propose an efficient and provable algorithm for 1NNs that can achieve the underlying ground-truth parameters.

3 Generalization Bound / Recovery Guarantees

The achievability of global optima of the objective from the training data doesn’t guarantee the learned model to be able to generalize well on unseen testing data. In the literature, we find three main approaches to generalization guarantees.

1) Use generalization analysis frameworks, including VC dimension/Rademacher complexity, to bound the generalization error. A few works have studied the generalization performance for NNs. [XLS17] follow [SC16] but additionally provide generalization bounds using Rademacher complexity. They assume the obtained parameters are in a regularization set so that the generalization performance is guaranteed, but this assumption can’t be justified theoretically. [HRS16] apply stability analysis to the generalization analysis of SGD for convex and non-convex problems, arguing early stopping is important for generalization performance.

2) Assume an underlying model and try to recover this model. This direction is popular for many non-convex problems including matrix completion/sensing [JNS13, Har14, SL15, BLWZ17], mixed linear regression [ZJD16], subspace recovery [EV09] and other latent models [AGH+14].

Without making any assumptions, those non-convex problems are intractable [AGKM12, GV15, SWZ17a, GG11, RSW16, SR11, HM13, AGM12, YCS14]. Recovery guarantees for NNs also need assumptions. Several different approaches under different assumptions are provided to have recovery guarantees on different NN settings.

Tensor methods [AGH+14, WTSA15, WA16, SWZ16] are a general tool for recovering models with latent factors by assuming the data distribution is known. Some existing recovery guarantees for NNs are provided by tensor methods [SA15, JSA15]. However, [SA15] only provide guarantees to recover the subspace spanned by the weight matrix and no sample complexity is given, while [JSA15] require $O(d^{3}/\epsilon^{2})$ sample complexity. In this paper, we use tensor methods as an initialization step so that we don’t need very accurate estimation of the moments, which enables us to reduce the total sample complexity from $1/\epsilon^{2}$ to $\log(1/\epsilon)$.

[ABGM14] provide polynomial sample complexity and computational complexity bounds for learning deep representations in unsupervised setting, and they need to assume the weights are sparse and randomly distributed in $$.

[Tia17] analyze 1NN by assuming Gaussian inputs in a supervised setting, in particular, regression and classification with a teacher. This paper also considers this setting. However, there are some key differences. a) [Tia17] require the second-layer parameters are all ones, while we can learn these parameters. b) In [Tia17], the ground-truth first-layer weight vectors are required to be orthogonal, while we only require linear independence. c) [Tia17] require a good initialization but doesn’t provide initialization methods, while we show the parameters can be efficiently initialized by tensor methods. d) In [Tia17], only the population case (infinite sample size) is considered, so there is no sample complexity analysis, while we show finite sample complexity.

Recovery guarantees for convolution neural network with Gaussian inputs are provided in [BG17], where they show a globally converging guarantee of gradient descent on a one-hidden-layer no-overlap convolution neural network. However, they consider population case, so no sample complexity is provided. Also their analysis depends on $\mathsf{ReLU}$ activations and the no-overlap case is very unlikely to be used in practice. In this paper, we consider a large range of activation functions, but for one-hidden-layer fully-connected NNs.

3) Improper Learning. In the improper learning setting for NNs, the learning algorithm is not restricted to output a NN, but only should output a prediction function whose error is not much larger than the error of the best NN among all the NNs considered. [ZLJ16, ZLW16] propose kernel methods to learn the prediction function which is guaranteed to have generalization performance close to that of the NN. However, the sample complexity and computational complexity are exponential. [AZS14] transform NNs to convex semi-definite programming. The works by [Bac14] and [BRV+05] are also in this direction. However, these methods are actually not learning the original NNs. Another work by [ZLWJ17] uses random initializations to achieve arbitrary small excess risk. However, their algorithm has exponential running time in $1/\epsilon$.

The paper is organized as follows. In Section 3, we present our problem setting and show three key properties of activations required for our guarantees. In Section 4, we introduce the formal theorem of local strong convexity and show local linear convergence for smooth activations. Section 5 presents a tensor method to initialize the parameters so that they fall into the basin of the local strong convexity region.

Problem Formulation

We consider the following regression problem. Given a set of $n$ samples

such that each sample $(x,y)\in S$ is sampled i.i.d. from this distribution, with

where $\phi(z)$ is the activation function, $k$ is the number of nodes in the hidden layer. The main question we want to answer is: How many samples are sufficient to recover the underlying parameters?

It is well-known that, training one hidden layer neural network is NP-complete [BR88]. Thus, without making any assumptions, learning deep neural network is intractable. Throughout the paper, we assume $x$ follows a standard normal distribution; the data is noiseless; the dimension of input data is at least the number of hidden nodes; and activation function $\phi(z)$ satisfies some reasonable properties.

Actually our results can be easily extended to multivariate Gaussian distribution with positive definite covariance and zero mean since we can estimate the covariance first and then transform the input to a standard normal distribution but with some loss of accuracy. Although this paper focuses on the regression problem, we can transform classification problems to regression problems if a good teacher is provided as described in [Tia17]. Our analysis requires $k$ to be no greater than $d$, since the first-layer parameters will be linearly dependent otherwise.

For activation function $\phi(z)$, we assume it is continuous and if it is non-smooth let its first derivative be left derivative. Furthermore, we assume it satisfies Property 3.1, 3.2, and 3.3. These properties are critical for the later analyses. We also observe that most activation functions actually satisfy these three properties.

The first derivative $\phi^{\prime}(z)$ is nonnegative and homogeneously bounded, i.e., $0\leq\phi^{\prime}(z)\leq L_{1}|z|^{p}$ for some constants $L_{1}>0$ and $p\geq 0$.

The second derivative $\phi^{\prime\prime}(z)$ is either (a) globally bounded $|\phi^{\prime\prime}(z)|\leq L_{2}$ for some constant $L_{2}$, i.e., $\phi(z)$ is $L_{2}$-smooth, or (b) $\phi^{\prime\prime}(z)=0$ except for $e$ ($e$ is a finite constant) points.

The first two properties are related to the first derivative $\phi^{\prime}(z)$ and the last one is about the second derivative $\phi^{\prime\prime}(z)$. At high level, Property 3.1 requires $\phi$ to be non-decreasing with homogeneously bounded derivative; Property 3.2 requires $\phi$ to be highly non-linear; Property 3.3 requires $\phi$ to be either smooth or piece-wise linear.

Positive Definiteness of Hessian

In this section, we study the Hessian of empirical risk near the ground truth. We consider the case when $v^{*}$ is already known. Note that for homogeneous activations, we can assume $v^{*}_{i}\in\{-1,1\}$ since $v\phi(z)=\frac{v}{|v|}\phi(|v|^{1/p}z)$, where $p$ is the degree of homogeneity. As $v_{i}^{*}$ only takes discrete values for homogeneous activations, in the next section, we show we can exactly recover $v^{*}$ using tensor methods with finite samples.

For a set of samples $S$, we define the Empirical Risk,

For a distribution ${\cal D}$, we define the Expected Risk,

Let’s calculate the gradient and the Hessian of $\widehat{f}_{S}(W)$ and $f_{\cal D}(W)$. For each $j\in[k]$, the partial gradient of $f_{\cal D}(W)$ with respect to $w_{j}$ can be represented as

For each $j,l\in[k]$ and $j\neq l$, the second partial derivative of $f_{\cal D}(W)$ for the $(j,l)$-th off-diagonal block is,

and for each $j\in[k]$, the second partial derivative of $f_{\cal D}(W)$ for the $j$-th diagonal block is

If Property 3.3(b) is satisfied, $\phi^{\prime\prime}(z)=0$ almost surely. So in this case the diagonal blocks of the empirical Hessian can be written as,

Now we show the Hessian of the objective near the global optimum is positive definite.

As we can see from Theorem 4.2, $\rho(\sigma_{k})$ from Property 3.2 plays an important role for positive definite (PD) property. Interestingly, many popular activations, like ReLU, sigmoid and tanh, have $\rho(\sigma_{k})>0$, while some simple functions like linear ($\phi(z)=z$) and square ($\phi(z)=z^{2}$) functions have $\rho(\sigma_{k})=0$ and their Hessians are rank-deficient. Another important numbers are $\kappa$ and $\lambda$, two different condition numbers of the weight matrix, which directly influences the positive definiteness. If $W^{*}$ is rank deficient, $\lambda\to\infty$, $\kappa\to\infty$ and we don’t have PD property. In the best case when $W^{*}$ is orthogonal, $\lambda=\kappa=1$. In the worse case, $\lambda$ can be exponential in $k$. Also $W$ should be close enough to $W^{*}$. In the next section, we provide tensor methods to initialize $w_{i}^{*}$ and $v_{i}^{*}$ such that they satisfy the conditions in Theorem 4.2.

For the PD property to hold, we need the samples to be independent of the current parameters. Therefore, we need to do resampling at each iteration to guarantee the convergence in iterative algorithms like gradient descent. The following theorem provides the linear convergence guarantee of gradient descent for smooth activations.

Let $W$ be the current iterate satisfying $\|W-W^{*}\|\leq\operatorname{poly}(1/\nu,1/k,1/\lambda,\rho/\sigma_{1}^{2p})\|W^{*}\|$. Let $S$ denote a set of i.i.d. samples from distribution ${\cal D}$ (defined in (1)) with $|S|\geq d\cdot\operatorname{poly}(\log d,t,k,\nu,\tau,\lambda,\sigma_{1}^{2p}/\rho)$ and let the activation function satisfy Property 3.1,3.2 and 3.3(a). Define $m_{0}:=\Theta(v_{\min}^{2}\rho(\sigma_{k})/(\kappa^{2}\lambda))$ and $M_{0}:=\Theta(kv_{\max}^{2}\sigma_{1}^{2p})$. If we perform gradient descent with step size $1/M_{0}$ on $\widehat{f}_{S}(W)$ and obtain the next iterate,

then with probability at least $1-d^{-\Omega(t)}$,

We provide the proofs in the Appendix D.1

Tensor Methods for Initialization

In this section, we show that Tensor methods can recover the parameters $W^{*}$ to some precision and exactly recover $v^{*}$ for homogeneous activations.

It is known that most tensor problems are NP-hard [Hås90, HL13] or even hard to approximate [SWZ17b]. However, by making some assumptions, tensor decomposition method becomes efficient [AGH+14, WTSA15, WA16, SWZ16]. Here we utilize the noiseless assumption and Gaussian inputs assumption to show a provable and efficient tensor methods.

where $A_{1,i,j}=v_{i}\otimes v_{i}\otimes e_{j}\otimes e_{j}$, $A_{2,i,j}=v_{i}\otimes e_{j}\otimes v_{i}\otimes e_{j}$, $A_{3,i,j}=e_{j}\otimes v_{i}\otimes v_{i}\otimes e_{j}$, $A_{4,i,j}=v_{i}\otimes e_{j}\otimes e_{j}\otimes v_{i}$, $A_{5,i,j}=e_{j}\otimes v_{i}\otimes e_{j}\otimes v_{i}$ and $A_{6,i,j}=e_{j}\otimes e_{j}\otimes v_{i}\otimes v_{i}$.

Denote $\overline{w}=w/\|w\|$. Now let’s calculate some moments.

According to Definition 5.1, we have the following results,

For each $j\in$, $M_{j}=\sum_{i=1}^{k}v_{i}^{*}m_{j,i}\overline{w}_{i}^{*\otimes j}$.

Note that some $m_{j,i}$’s will be zero for specific activations. For example, for activations with symmetric first derivatives, i.e., $\phi^{\prime}(z)=\phi^{\prime}(-z)$, like sigmoid and erf, we have $\phi(z)+\phi(-z)$ being a constant and $M_{2}=0$ since $\gamma_{0}(\sigma)=\gamma_{2}(\sigma)$. Another example is $\mathsf{ReLU}$. $\mathsf{ReLU}$ functions have vanishing $M_{3}$, i.e., $M_{3}=0$, as $\gamma_{3}(\sigma)=3\gamma_{1}(\sigma)$. To make tensor methods work, we make the following assumption.

Assume the activation function $\phi(z)$ satisfies the following conditions: 1. If $M_{j}\neq 0$, then $m_{j,i}\neq 0$ for all $i\in[k]$. 2. At least one of $M_{3}$ and $M_{4}$ is non-zero. 3. If $M_{1}=M_{3}=0$, then $\phi(z)$ is an even function, i.e., $\phi(z)=\phi(-z)$. 4. If $M_{2}=M_{4}=0$, then $\phi(z)$ is an odd function, i.e., $\phi(z)=-\phi(-z)$.

If $\phi(z)$ is an odd function then $\phi(z)=-\phi(-z)$ and $v\phi(w^{\top}x)=-v\phi(-w^{\top}x)$. Hence we can always assume $v>0$. If $\phi(z)$ is an even function, then $v\phi(w^{\top}x)=v\phi(-w^{\top}x)$. So if $w$ recovers $w^{*}$ then $-w$ also recovers $w^{*}$. Note that $\mathsf{ReLU}$, leaky $\mathsf{ReLU}$ and squared $\mathsf{ReLU}$ satisfy Assumption 5.3. We further define the following non-zero moments.

According to Definition 5.1 and 5.4, we have,

$P_{2}=\sum_{i=1}^{k}v_{i}^{*}m_{j_{2},i}(\alpha^{\top}\overline{w}_{i}^{*})^{{j_{2}-2}}\overline{w}_{i}^{*\otimes 2}$ and $P_{3}=\sum_{i=1}^{k}v_{i}^{*}m_{j_{3},i}(\alpha^{\top}\overline{w}_{i}^{*})^{{j_{3}-3}}\overline{w}_{i}^{*\otimes 3}$.

2 Algorithm

Now we briefly introduce how to use a set of samples with size linear in dimension to recover the ground truth parameters to some precision. As shown in the previous section, we have a rank-$k$ 3rd-order moment $P_{3}$ that has tensor decomposition formed by $\{\overline{w}_{1}^{*},\overline{w}_{2}^{*},\cdots,\overline{w}_{k}^{*}\}$. Therefore, we can use the non-orthogonal decomposition method [KCL15] to decompose the corresponding estimated tensor $\widehat{P}_{3}$ and obtain an approximation of the parameters. The precision of the obtained parameters depends on the estimation error of $P_{3}$, which requires $\Omega(d^{3}/\epsilon^{2})$ samples to achieve $\epsilon$ error. Also, the time complexity for tensor decomposition on a $d\times d\times d$ tensor is $\Omega(d^{3})$.

The detailed algorithm is shown in Algorithm 1. First, we randomly partition the dataset into three subsets each with size $\widetilde{O}(d)$. Then apply the power method on $\widehat{P}_{2}$, which is the estimation of $P_{2}$ from $S_{2}$, to estimate $V$. After that, the non-orthogonal tensor decomposition (KCL)[KCL15] on $\widehat{R}_{3}$ outputs $\widehat{u}_{i}$ which estimates $s_{i}V^{\top}\overline{w}_{i}^{*}$ for $i\in[k]$ with unknown sign $s_{i}\in\{-1,1\}$. Hence $\overline{w}_{i}^{*}$ can be estimated by $s_{i}V\widehat{u}_{i}$. Finally we estimate the magnitude of $w_{i}^{*}$ and the signs $s_{i},v_{i}^{*}$ in the RecMagSign function for homogeneous activations. We discuss the details of each procedure and provide PowerMethod and RecMagSign algorithms in Appendix E.

3 Theoretical Analysis

We formally present our theorem for Algorithm 1, and provide the proof in the Appendix E.2.

Global Convergence

Combining the positive definiteness of the Hessian near the global optimal in Section 4 and the tensor initialization methods in Section 5, we come up with the overall globally converging algorithm Algorithm 2 and its guarantee Theorem 6.1.

with probability at least $1-d^{-\Omega(t)}$.

This follows by combining Theorem 4.4 and Theorem 5.6.

Numerical Experiments

First we show the number of samples required to recover the parameters for different dimensions. We fix $k=5$, change $d$ for $d=10,20,\cdots,100$ and $n$ for $n=1000,2000,\cdots,10000$. For each pair of $d$ and $n$, we run $10$ trials. We say a trial successfully recovers the parameters if there exists a permutation $\pi:[k]\rightarrow[k]$, such that the returned parameters $W$ and $v$ satisfy

We record the recovery rates and represent them as grey scale in Fig. 1(a). As we can see from Fig. 1(a), the least number of samples required to have 100% recovery rate is about proportional to the dimension.

Next we test the tensor initialization. We show the error between the output of the tensor method and the ground truth parameters against the number of samples under different dimensions in Fig 1(b). The pure dark blocks indicate, in at least one of the 10 trials, $\sum_{i=1}^{k}v^{{(0)}}_{i}\neq\sum_{i=1}^{k}v^{*}_{i}$, which means $v^{(0)}_{i}$ is not correctly initialized. Let $\Pi(k)$ denote the set of all possible permutations $\pi:[k]\rightarrow[k]$. The grey scale represents the averaged error,

over 10 trials. As we can see, with a fixed dimension, the more samples we have the better initialization we obtain. We can also see that to achieve the same initialization error, the sample complexity required is about proportional to the dimension.

We also compare different initialization methods for gradient descent in Fig. 1(c). We fix $d=10,k=5,n=10000$ and compare three different initialization approaches, (\@slowromancapi@) Let both $v$ and $W$ be initialized from tensor methods, and then do gradient descent for $W$ while $v$ is fixed; (\@slowromancapii@) Let both $v$ and $W$ be initialized from random Gaussian, and then do gradient descent for both $W$ and $v$; (\@slowromancapiii@) Let $v=v^{*}$ and $W$ be initialized from random Gaussian, and then do gradient descent for $W$ while $v$ is fixed. As we can see from Fig 1(c), Approach (\@slowromancapi@) is the fastest and Approach (\@slowromancapii@) doesn’t converge even if more iterations are allowed. Both Approach (\@slowromancapi@) and (\@slowromancapiii@) have linear convergence rate when the objective value is small enough, which verifies our local linear convergence claim.

Conclusion

As shown in Theorem 6.1, the tensor initialization followed by gradient descent will provide a globally converging algorithm with linear time/sample complexity in dimension, logarithmic in precision and polynomial in other factors for smooth homogeneous activation functions. Our distilled properties for activation functions include a wide range of non-linear functions and hopefully provide an intuition to understand the role of non-linear activations played in optimization. Deeper neural networks and convergence for SGD will be considered in the future.

References

Appendix A Notation

We use $\|T\|$ to denote the operator norm of the tensor $T$, i.e.,

For any function $f$, we define $\widetilde{O}(f)$ to be $f\cdot\log^{O(1)}(f)$. In addition to $O(\cdot)$ notation, for two functions $f,g$, we use the shorthand $f\lesssim g$ (resp. $\gtrsim$) to indicate that $f\leq Cg$ (resp. $\geq$) for an absolute constant $C$. We use $f\eqsim g$ to mean $cf\leq g\leq Cf$ for constants $c,C$.

Appendix B Preliminaries

In this section, we introduce some lemmata and corollaries that will be used in the proofs.

We provide some facts that will be used in the later proofs.

Let $z$ denote a fixed $d$-dimensional vector, then for any $C\geq 1$ and $n\geq 1$, we have

This follows by Proposition 1.1 in [HKZ12]. ∎

This follows by Proposition 1.1 in [HKZ12]. ∎

Part (\@slowromancapii@). We first show how to lower bound $\sigma_{k}(\overline{W})$,

It remains to upper bound $\sigma_{1}(\overline{W})$,

We show how to simplify $\|V_{\perp}^{\top}U\|$,

Similarly we can simplify $\|U_{\perp}^{\top}V\|$,

For each $i\in[d_{3}]$, let $b_{i}$ denote the $i$-th column of $B$. We can upper bound $\|CB\|$,

Since all the three components $|u^{\top}x|$, $|v^{\top}x|$, $|w^{\top}x|$ are positive and related to a common random vector $x$, we can show a lower bound,

B.2 Matrix Bernstein

In many proofs we need to bound the difference between some population matrices/tensors and their empirical versions. Typically, the classic matrix Bernstein inequality requires the norm of the random matrix be bounded almost surely (e.g., Theorem 6.1 in [Tro12]) or the random matrix satisfies subexponential property (Theorem 6.2 in [Tro12]) . However, in our cases, most of the random matrices don’t satisfy these conditions. So we derive the following lemmata that can deal with random matrices that are not bounded almost surely or follow subexponential distribution, but are bounded with high probability.

Then we have for any $0<\epsilon<1$ and $t\geq 1$, if

with probability at least $1-1/d^{2t}-n\gamma$,

In the next a few paragraphs, we will upper bound the above three terms.

The first term in Eq. (4). Denote $\xi^{c}$ as the complementary set of $\xi$, thus $\Pr[\xi^{c}_{i}]\leq\gamma$. By a union bound over $i\in[n]$, with probability $1-n\gamma$, $\|B_{i}\|\leq m$ for all $i\in[n]$. Thus $\widehat{M}=\widehat{B}$.

The second term in Eq. (4). For a matrix $B$ sampled from ${\cal B}$, we use $\xi$ to denote the event that $\xi=\{\|B\|\leq m\}$. Then, we can upper bound $\|\overline{M}-\overline{B}\|$ in the following way,

Since $\epsilon<1$, we also have $\|\overline{M}-\overline{B}\|\leq\frac{1}{2}\|\overline{B}\|$ and $\frac{3}{2}\|\overline{B}\|\geq\|\overline{M}\|\geq\frac{1}{2}\|\overline{B}\|$.

The third term in Eq. (4). We can bound $\|\widehat{M}-\overline{M}\|$ by Matrix Bernstein’s inequality [Tro12].

for $t\geq 1$, we have with probability at least $1-1/d^{2t}$,

Putting it all together, we have for $0<\epsilon<1$, if

with probability at least $1-1/d^{2t}-n\gamma$,

We show that the four Properties in Lemma B.7 are satisfied. Define function $B(x)=h(x)xx^{\top}$.

(\@slowromancapi@) $\|B(x)\|=\|h(x)xx^{\top}\|=|h(x)|\|x\|^{2}$.

Applying Lemma B.7, we obtain, for any $0<\epsilon<1$ and $t\geq 1$, if

Appendix C Properties of Activation Functions

We can easily verify that $\mathsf{ReLU}$ , leaky $\mathsf{ReLU}$ and squared $\mathsf{ReLU}$ satisfy Property 3.2 by calculating $\rho(\sigma)$ in Property 3.2, which is shown in Table 1. Property 3.1 for $\mathsf{ReLU}$ , leaky $\mathsf{ReLU}$ and squared $\mathsf{ReLU}$ can be verified since they are non-decreasing with bounded first derivative. $\mathsf{ReLU}$ and leaky $\mathsf{ReLU}$ are piece-wise linear, so they satisfy Property 3.3(b). Squared $\mathsf{ReLU}$ is smooth so it satisfies Property 3.3(a).

The equality in the first inequality happens when $\phi^{\prime}(\sigma\cdot z)$ is a constant a.e.. The equality in the second inequality happens when $|\phi^{\prime}(\sigma\cdot z)|$ is a constant a.e., which is invalidated by the non-linearity and smoothness condition. The equality in the third inequality holds only when $\phi^{\prime}(z)=0$ a.e., which leads to a constant function under non-decreasing condition. Therefore, $\rho(\sigma)>0$ for any smooth non-decreasing non-linear activations with bounded symmetric first derivatives. The statements about linear activations and quadratic activation follow direct calculations. ∎

Appendix D Local Positive Definiteness of Hessian

The main idea of the proof follows the following inequalities,

The proof sketch is first to bound the range of the eigenvalues of $\nabla^{2}f_{\cal D}(W^{*})$ (Lemma D.3) and then bound the spectral norm of the remaining error, $\|\nabla^{2}\widehat{f}_{S}(W)-\nabla^{2}f_{\cal D}(W^{*})\|$. $\|\nabla^{2}\widehat{f}_{S}(W)-\nabla^{2}f_{\cal D}(W^{*})\|$ can be further decomposed into two parts, $\|\nabla^{2}\widehat{f}_{S}(W)-H\|$ and $\|H-\nabla^{2}f_{\cal D}(W^{*})\|$, where $H$ is $\nabla^{2}f_{\cal D}(W)$ if $\phi$ is smooth, otherwise $H$ is a specially designed matrix . We can upper bound them when $W$ is close enough to $W^{*}$ and there are enough samples. In particular, if the activation satisfies Property 3.3(a), see Lemma D.10 for bounding $\|H-\nabla^{2}f_{\cal D}(W^{*})\|$ and Lemma D.11 for bounding $\|H-\nabla^{2}\widehat{f}_{S}(W)\|$. If the activation satisfies Property 3.3(b), see Lemma D.15.

Finally we can complete the proof by setting $\delta=O(v_{\min}^{2}\rho(\sigma_{1})/(kv_{\max}^{2}\kappa^{2}\lambda\sigma_{1}^{2p}))$ in Lemma D.11 and Lemma D.15, setting $\|W-W^{*}\|\lesssim v_{\min}^{2}\rho(\sigma_{k})/(k\kappa^{2}\lambda v_{\max}^{2}\sigma_{1}^{p})$ in Lemma D.10 and setting $\|W-W^{*}\|\leq v_{\min}^{4}\rho^{2}(\sigma_{k})\sigma_{k}/(k^{2}\kappa^{4}{\lambda}^{2}v_{\max}^{4}\sigma_{1}^{4p})$ in Lemma D.15. ∎

D.1.2 Local Linear Convergence of Gradient Descent

Although Theorem D.1 gives upper and lower bounds for the spectrum of the Hessian w.h.p., it only holds when the current set of parameters $W$ are independent of samples. When we use iterative methods, like gradient descent, to optimize this objective, the next iterate calculated from the current set of samples will depend on this set of samples. Therefore, we need to do resampling at each iteration. Here we show that for activations that satisfies Properties 3.1, 3.2 and 3.3(a), linear convergence of gradient descent is guaranteed. To the best of our knowledge, there is no linear convergence guarantees for general non-smooth objective. So the following proposition also applies to smooth objectives only, which excludes $\mathsf{ReLU}$.

Let $S$ denote a set of i.i.d. samples from distribution ${\cal D}$ (defined in (1)) Let the activation function satisfy Property 3.1,3.2 and 3.3(a). Define

and perform gradient descent with step size $1/M_{0}$ on $\widehat{f}_{S}(W^{c})$ and obtain the next iterate,

then with probability at least $1-d^{-\Omega(t)}$,

To prove Theorem D.2, we need to show the positive definite properties on the entire line between the current iterate and the optimum by constructing a set of anchor points, which are independent of the samples. Then we apply traditional analysis for the linear convergence of gradient descent.

In particular, given a current iterate $W^{c}$, we set $d^{(p+1)/2}$ anchor points $\{W^{a}\}_{a=1,2,\cdots,d^{(p+1)/2}}$ equally along the line $\xi W^{*}+(1-\xi)W^{c}$ for $\xi\in$.

According to Theorem D.1, by setting $t\leftarrow t+(p+1)/2$, we have with probability at least $1-d^{-(t+(p+1)/2)}$ for each anchor point $\{W^{a}\}$,

Then given an anchor point $W^{a}$, according to Lemma D.16, we have with probability $1-2d^{-(t+(p+1)/2)}$, for any points $W$ between $(W^{a-1}+W^{a})/2$ and $(W^{a}+W^{a+1})/2$,

Finally by applying union bound over these $d^{(p+1)/2}$ small intervals, we have with probability at least $1-d^{-t}$ for any points $W$ on the line between $W^{c}$ and $W^{*}$,

Now we can apply traditional analysis for linear convergence of gradient descent.

We can upper bound $\|\nabla\widehat{f}_{S}(W^{c})\|_{F}^{2}$,

where the third equality holds by setting $\eta=\frac{1}{M_{0}}$. ∎

D.2 Positive Definiteness of Population Hessian at the Ground Truth

The goal of this Section is to prove Lemma D.3.

If $\phi(z)$ satisfies Property 3.1,3.2 and 3.3, we have the following property for the second derivative of function $f_{\cal D}(W)$ at $W^{*}$,

The proof directly follows Lemma D.6 (Section D.2.2) and Lemma D.7(Section D.2.3). ∎

The main idea is to explicitly calculate the LHS of Eq (8), then reformulate the equation and find a lower bound represented by $\alpha_{0},\alpha_{1},\alpha_{2},\beta_{0},\beta_{2}$.

Further, we can rewrite the diagonal term in the following way,

We can rewrite the off-diagonal term in the following way,

$A+B_{1}+B_{2}+B_{3}=C_{1}+C_{2}+C_{3}+C_{4}+C_{5}$.

The key properties we need are, for two vectors $a,b$, $\|a+b\|^{2}=\|a\|^{2}+2\langle a,b\rangle+\|b\|^{2}$; for two matrices $A,B$, $\|A+B\|_{F}^{2}=\|A\|_{F}^{2}+2\langle A,B\rangle+\|B\|_{F}^{2}$. Then, we have

D.2.2 Lower Bound on the Eigenvalues of Hessian for Non-orthogonal Case

First we show the lower bound of the eigenvalues. The main idea is to reduce the problem to a $k$-by-$k$ problem and then lower bound the eigenvalues using orthogonal weight matrices.

If $\phi(z)$ satisfies Property 3.1,3.2 and 3.3, we have the following property for the second derivative of function $f_{\cal D}(W)$ at $W^{*}$,

We use ${\cal D}_{d}$ to denote Gaussian distribution ${\cal N}(0,I_{d})$, ${\cal D}_{d-k}$ to denote Gaussian distribution ${\cal N}(0,I_{d-k})$, and ${\cal D}_{k}$ to denote Gaussian distribution ${\cal N}(0,I_{k})$. Then we can rewrite formulation (10) (removing $v_{\min}^{2}$) as

We calculate $A,B,C$ separately. First, we can show

Third, we have $C=0$ since $U_{\perp}^{\top}x$ is independent of $w_{i}^{*\top}x$ and $U^{\top}x$. Thus, putting them all together,

For $B$, similar to the proof of Lemma D.4, we have,

Note that $1=\|a\|^{2}=\|b\|^{2}+\|c\|^{2}$. Thus, we finish the proof for the lower bound. ∎

D.2.3 Upper Bound on the Eigenvalues of Hessian for Non-orthogonal Case

If $\phi(z)$ satisfies Property 3.1,3.2 and 3.3, we have the following property for the second derivative of function $f_{\cal D}(W)$ at $W^{*}$,

Similarly, we can calculate the upper bound of the eigenvalues by

where the first inequality follows Hölder’s inequality, the second inequality follows by Property 3.1, the third inequality follows by $\|w_{i}^{*}\|\leq\sigma_{1}(W^{*})$, and the last inequality by Cauchy-Schwarz inequality.

D.3 Error Bound of Hessians near the Ground Truth for Smooth Activations

The goal of this Section is to prove Lemma D.8

Let $\phi(z)$ satisfy Property 3.1, Property 3.2 and Property 3.3(a). Let $W$ satisfy $\|W-W^{*}\|\leq\sigma_{k}/2$. Let $S$ denote a set of i.i.d. samples from the distribution defined in (1). Then for any $t\geq 1$ and $0<\epsilon<1/2$, if

then we have, with probability at least $1-1/d^{\Omega(t)}$,

This follows by combining Lemma D.10 and Lemma D.11 directly. ∎

The goal of this Section is to prove Lemma D.10.

Note that if $\|W-W^{*}\|\leq\sigma_{k}(W^{*})/2$, we have $\sigma_{k}(W^{*})/2\leq\sigma_{i}(W)\leq\frac{3}{2}\sigma_{1}(W^{*})$ for all $i\in[k]$ by Weyl’s inequality. By definition of singular value, we have $\sigma_{k}(W^{*})\leq\|w_{i}^{*}\|\leq\sigma_{1}(W^{*})$. By definition of spectral norm, we have $\|w_{i}-w_{i}^{*}\|\leq\|W-W^{*}\|$. Thus, we can lower bound $\|w_{i}\|$,

Similarly, we have $\|w_{i}\|\geq\frac{1}{2}\|w_{i}^{*}\|$. ∎

If $\phi(z)$ satisfies Property 3.1, Property 3.2 and Property 3.3(a), then for any $W$ with $\|W-W^{*}\|\leq\sigma_{k}/2$, we have

In the next a few paragraphs, we first show how to bound the off-diagonal term, and then show how to bound the diagonal term.

where the first step follows by definition of $\Delta_{i,l}$, the second step follows by definition of spectral norm and $v_{i}^{*}v_{l}^{*}\leq v_{\max}^{2}$, the third step follows by triangle inequality, the fourth step follows by linearity of expectation, the fifth step follows by Property 3.3(a), i.e., $|\phi^{\prime}(w_{i}^{\top}x)-\phi^{\prime}(w_{i}^{*\top}x)|\leq L_{2}|(w_{i}-w_{i}^{*})^{\top}x|$, the sixth step follows by Property 3.1, i.e., $\phi^{\prime}(z)\leq L_{1}|z|^{p}$, the seventh step follows by Fact B.6, and the last step follows by $\|a\|^{2}=1$, $\|w_{i}-w_{i}^{*}\|\leq\|W-W^{*}\|$, $\|w_{i}\|\leq\frac{3}{2}\|w_{i}^{*}\|$, and $\|w_{i}^{*}\|\leq\sigma_{1}(W^{*})$.

Note that the proof for the off-diagonal terms also applies to bounding the second-term in the diagonal block Eq. (D.3.1). Thus we only need to show how to bound the first term in the diagonal block Eq. (D.3.1).

where the first step follows by $y=\sum_{r=1}^{k}v_{r}^{*}\phi(w_{r}^{*\top}x)$, the second step follows by definition of spectral norm and $v_{r}^{*}v_{i}^{*}\leq|v_{\max}|^{2}$, the third step follows by Property 3.3(a), i.e., $|\phi^{\prime\prime}(w_{i}^{\top}x)|\leq L_{2}$, the fourth step follows by $|\phi(w_{r}^{\top}x)-\phi(w_{r}^{*\top}x)\leq\max_{z\in[w_{r}^{\top}x,w_{r}^{*\top}x]}|\phi^{\prime}(z)||(w_{r}-w_{r}^{*})^{\top}x|$, the fifth step follows Property 3.1, i.e., $|\phi^{\prime}(z)|\leq L_{1}|z|^{p}$, the sixth step follows by $\max_{z\in[w_{r}^{\top}x,w_{r}^{*\top}x]}|z|^{p}\leq(|w_{r}^{\top}x|^{p}+|w_{r}^{*\top}x|^{p})$, the seventh step follows by Fact B.6.

Putting it all together, we can bound the error by

D.3.2 Empirical and Population Difference for Smooth Activations

Note that Bernstein inequality requires the spectral norm of each random matrix to be bounded almost surely. However, since we assume Gaussian distribution for $x$, $\|x\|^{2}$ is not bounded almost surely. The main idea is to do truncation and then use Matrix Bernstein inequality. Details can be found in Lemma B.7 and Corollary B.8.

Let $\phi(z)$ satisfy Property 3.1,3.2 and 3.3(a). Let $W$ satisfy $\|W-W^{*}\|\leq\sigma_{k}/2$. Let $S$ denote a set of i.i.d. samples from distribution ${\cal D}$ (defined in (1)). Then for any $t\geq 1$ and $0<\epsilon<1/2$, if

then we have, with probability at least $1-d^{-\Omega(t)}$,

Define $\Delta=\nabla^{2}f_{\cal D}(W)-\nabla^{2}\widehat{f}_{S}(W)$. Let’s first consider the diagonal blocks. Define

Let’s further decompose $\Delta_{i,i}$ into $\Delta_{i,i}=\Delta_{i,i}^{(1)}+\Delta_{i,i}^{(2)}$, where

Combining Claims. D.12, D.13 and D.14, and taking a union bound over $k^{2}$ different $\Delta_{i,j}$, we obtain if $n\geq\epsilon^{-2}\kappa^{2}(W^{*})\tau d\cdot\operatorname{poly}(t,\log d)$ for any $\epsilon\in(0,1/2)$, with probability at least $1-1/d^{t}$,

For each $i\in[k]$, if $n\geq d\operatorname{poly}(\log d,t)$

According to Properties 3.1,3.2 and 3.3(a), we have for each $x\in S$,

(\@slowromancapi@) Bounding $|h_{r}(x)|$.

According to Fact B.1, we have for any constant $t\geq 1$, with probability $1-1/(nd^{4t})$,

(\@slowromancapii@) Bounding $\|\overline{B}_{r}\|$.

where the first step follows by definition of spectral norm, and last step follows by Fact B.6. Using Fact B.6, we can also prove an upper bound $\|\overline{B}_{r}\|$, $\|\overline{B}_{r}\|\lesssim L_{1}L_{2}\|w_{r}\|^{p}\|w_{r}-w_{r}^{*}\|$.

By applying Corollary B.8, if $n\geq\epsilon^{-2}d\operatorname{poly}(\log d,t)$, then with probability $1-1/d^{4t}$,

For each $i\in[k]$, if $n\geq\epsilon^{-2}d\tau\operatorname{poly}(\log d,t)$ , then

Recall the definition of $\Delta_{i,i}^{(2)}$.

(\@slowromancapi@) Bounding $|h_{i}(x)|$.

For any constant $t\geq 1$, $(\phi^{\prime}(w_{i}^{\top}x))^{2}\leq L_{1}^{2}|w_{i}^{\top}x|^{2p}\lesssim L_{1}^{2}\|w_{i}\|^{2p}t^{p}\log^{p}(n)$ with probability $1-1/(nd^{4t})$ according to Fact B.1.

(\@slowromancapii@) Bounding $\|\overline{B}_{i}\|$.

where $\overline{w}_{i}=w_{i}/\|w_{i}\|$ and $v$ is a unit vector orthogonal to $w_{i}$ such that $a=\alpha\overline{w}_{i}+\beta v$. Now from Property 3.2, we have

By applying Corollary B.8, we have, for any $0<\epsilon<1$, if $n\geq\epsilon^{-2}d\frac{\|w_{i}\|^{4p}}{\rho^{2}(\|w_{i}\|)}\operatorname{poly}(\log d,t)$ the following bound holds

Therefore, if $n\geq\epsilon^{-2}d\tau\operatorname{poly}(\log d,t)$, where $\tau=\frac{(3\sigma_{1}/2)^{4p}}{\min_{\sigma\in[\sigma_{k}/2,3\sigma_{1}/2]}\rho^{2}(\sigma)}$, we have with probability $1-1/d^{4t}$

For each $i\in[k],j\in[k],i\neq j$, if $n\geq\epsilon^{-2}\kappa^{2}\tau d\operatorname{poly}(\log d,t)$, then

Recall the definition of off-diagonal blocks $\Delta_{i,l}$,

(\@slowromancapi@) Bounding $|h_{i,l}(x)|$.

For any constant $t\geq 1$, we have with probability $1-1/(nd^{4t})$

where the third step follows by Fact B.1.

(\@slowromancapii@) Bounding $\|\overline{B}_{i,l}\|$.

For simplicity, we just set $\|b\|=1$ and $\|c\|=0$ to lower bound,

where the second step is from Property 3.2 and the fact that $\sigma_{k}/2\leq\sigma_{2}(V)\leq\sigma_{1}(V)\leq 3\sigma_{1}/2$.

where the first step follows by Property 3.1, the second step follows by Fact B.6, and the last step follows by $\|v_{1}\|\leq\sigma_{1}$, $\|v_{2}\|\leq\sigma_{1}$ and $\|b\|\leq 1$. Similarly, we can upper bound,

Therefore, applying Corollary B.8, we have, if $n\geq\epsilon^{-2}\kappa^{2}(W^{*})\tau d\operatorname{poly}(\log d,t)$, then

holds with probability at least $1-1/d^{4t}$. ∎

D.4 Error Bound of Hessians near the Ground Truth for Non-smooth Activations

The goal of this Section is to prove Lemma D.15,

Let $\phi(z)$ satisfy Property 3.1,3.2 and 3.3(b). Let $W$ satisfy $\|W-W^{*}\|\leq\sigma_{k}/2$. Let $S$ denote a set of i.i.d. samples from the distribution defined in (1). Then for any $t\geq 1$ and $0<\epsilon<1/2$, if

with probability at least $1-d^{-\Omega(t)}$,

As we noted previously, when Property 3.3(b) holds, the diagonal blocks of the empirical Hessian can be written as, with probability $1$,

Note that $H\neq\nabla^{2}f_{\operatorname*{{\mathcal{D}}}}(W)$. However we can still bound $\|H-\nabla^{2}\widehat{f}_{S}(W)\|$ and $\|H-\nabla^{2}f_{\operatorname*{{\mathcal{D}}}}(W^{*})\|$ when we have enough samples and $\|W-W^{*}\|$ is small enough. The proof for $\|H-\nabla^{2}\widehat{f}_{S}(W)\|$ basically follows the proof of Lemma D.11 as $\Delta_{ii}^{(2)}$ in Eq. (14) and $\Delta_{il}$ in Eq. (16) forms the blocks of $H-\nabla^{2}\widehat{f}_{S}(W)$ and we can bound them without smoothness of $\phi(\cdot)$.

Now we focus on $H-\nabla^{2}f_{\operatorname*{{\mathcal{D}}}}(W^{*})$. We again consider each block.

We used the boundedness of $\phi^{\prime\prime}(z)$ to prove Lemma D.10. Here we can’t use this condition. Without smoothness, we will stop at the following position.

where the first step follows by definition of spectral norm, the second step follows by triangle inequality, and the last step follows by linearity of expectation.

where the first step follows by $a=Ub+U_{\bot}c$, the last step follows by $(a+b)^{2}\leq 2a^{2}+2b^{2}$.

By Property 3.3(b), we have $e$ exceptional points which have $\phi^{\prime\prime}(z)\neq 0$. Let these $e$ points be $p_{1},p_{2},\cdots,p_{e}$. Note that if $v_{i}^{\top}z$ and $v_{i}^{*\top}z$ are not separated by any of these exceptional points, i.e., there exists no $j\in[e]$ such that $v_{i}^{\top}z\leq p_{j}\leq v_{i}^{*\top}z$ or $v_{i}^{*\top}z\leq p_{j}\leq v_{i}^{\top}z$, then we have $\phi^{\prime}(v_{i}^{\top}z)=\phi^{\prime}(v_{i}^{*\top}z)$ since $\phi^{\prime\prime}(s)$ are zeros except for $\{p_{j}\}_{j=1,2,\cdots,e}$. So we consider the probability that $v_{i}^{\top}z,v_{i}^{*\top}z$ are separated by any exception point. We use $\xi_{j}$ to denote the event that $v_{i}^{\top}z,v_{i}^{*\top}z$ are separated by an exceptional point $p_{j}$. By union bound, $1-\sum_{j=1}^{e}\Pr[\xi_{j}]$ is the probability that $v_{i}^{\top}z,v_{i}^{*\top}z$ are not separated by any exceptional point. The first term of Equation (D.4) can be bounded as,

where the first step follows by if $v_{i}^{\top}z,v_{i}^{*\top}z$ are not separated by any exceptional point then $\phi^{\prime}(v_{i}^{\top}z)=\phi^{\prime}(v_{i}^{*\top}z)$ and the last step follows by Hölder’s inequality and Property 3.1.

It remains to upper bound $\Pr_{z\sim\operatorname*{{\mathcal{D}}}_{3}}[\xi_{j}]$. First note that if $v_{i}^{\top}z,v_{i}^{*\top}z$ are separated by an exceptional point, $p_{j}$, then $|v_{i}^{*\top}z-p_{j}|\leq|v_{i}^{\top}z-v_{i}^{*\top}z|\leq\|v_{i}-v_{i}^{*}\|\|z\|$. Therefore,

Note that $(\frac{v_{i}^{*\top}z}{\|z\|\|v_{i}^{*}\|}+1)/2$ follows Beta(1,1) distribution which is uniform distribution on $$.

where the first step is because we can view $\frac{v_{i}^{*\top}z}{\|z\|}$ and $\frac{p_{j}}{\|z\|}$ as two independent random variables: the former is about the direction of $z$ and the later is related to the magnitude of $z$. Thus, we have

Finally combining Eq. (17), Eq. (19) and Eq. (20), we have

D.5 Positive Definiteness for a Small Region

Here we introduce a lemma which shows that the Hessian of any $W$, which may be dependent on the samples but is very close to an anchor point, is close to the Hessian of this anchor point.

with probability at least $1-d^{-t}$, for any $W$ (which is not necessarily to be independent of samples) satisfying $\|W^{a}-W\|\leq\sigma_{k}/4$, we have

We further decompose $\Delta_{i,i}$ into $\Delta_{i,i}=\Delta_{i,i}^{(1)}+\Delta_{i,i}^{(2)}$, where

We can further decompose $\Delta_{i,i}^{(1)}$ into $\Delta_{i,i}^{(1,1)}$ and $\Delta_{i,i}^{(1,2)}$,

Combining Claim D.17 and Claim D.18 , we have if

Therefore, combining Eq. (22), Claim D.19 and Claim D.20, we complete the proof. ∎

For each $i\in[k]$, if $n\geq d\operatorname{poly}(\log d,t)$, then

Define function $h_{1}(x)=\|x\|^{p+1}$ and $h_{2}(x)=|w_{q}^{*\top}x|^{p}|(w_{q}^{*}-w_{q}^{a})^{\top}x|$. Note that $h_{1}$ and $h_{2}$ don’t contain $W$ which maybe depend on the samples. Therefore, we can use the modified matrix Bernstein inequality (Corollary B.8) to bound $\Delta_{i,i}^{(1)}$.

(\@slowromancapi@) Bounding $|h_{1}(x)|$.

By Fact B.2, we have $h_{1}(x)\lesssim(td\log n)^{(p+1)/2}$ with probability at least $1-1/(nd^{4t})$.

Therefore we have with probability at least $1-1/d^{t}$,

For each $i\in[k]$, if $n\geq d\operatorname{poly}(\log d,t)$, then

holds with probability at least $1-1/d^{4t}$.

Recall the definition of $\Delta_{i,i}^{(1,2)}$,

In order to upper bound the $\|\Delta_{i,i}^{(1,2)}\|$, it suffices to upper bound the spectral norm of this quantity,

By Property 3.3, we have $|\phi^{\prime\prime}(w_{i}^{\top}x)-\phi^{\prime\prime}(w_{i}^{a\top}x)|\leq 2L_{2}$.

For the second part $|w_{q}^{*\top}x|^{p}|(w_{q}^{*}-w_{q}^{a})^{\top}x|xx^{\top}$, according to Eq. (D.3.2), we have with probability $1-d^{-t}$, if $n\geq d\operatorname{poly}(\log d,t)$,

For each $i\in[k]$, if $n\geq d\operatorname{poly}(\log d,t)$, then

Applying Corollary B.8 finishes the proof.

For each $i\in[k],j\in[k],i\neq j$, if $n\geq d\operatorname{poly}(\log d,t)$, then

Applying Corollary B.8 completes the proof. ∎

Appendix E Tensor Methods

We describe the details of each procedure in Algorithm 1 in this section.

a) Compute the subspace estimation from $\widehat{P}_{2}$ (Algorithm 3). Note that the eigenvalues of $P_{2}$ and $\widehat{P}_{2}$ are not necessarily nonnegative. However, only $k$ of the eigenvalues will have large magnitude. So we can first compute the top-k eigenvectors/eigenvalues of both $C\cdot I+\widehat{P}_{2}$ and $C\cdot I-\widehat{P}_{2}$, where $C$ is large enough such that $C\geq 2\|P_{2}\|$. Then from the $2k$ eigen-pairs, we pick the top-$k$ eigenvectors with the largest eigenvalues in magnitude, which is executed in TopK in Algorithm 3. For the outputs of TopK, $k_{1},k_{2}$ are the numbers of picked largest eigenvalues from $C\cdot I+\widehat{P}_{2}$ and $C\cdot I-\widehat{P}_{2}$ respectively and $\pi_{1}(i)$ returns the original index of $i$-th largest eigenvalue from $C\cdot I+\widehat{P}_{2}$ and similarly $\pi_{2}$ is for $C\cdot I-\widehat{P}_{2}$. Finally orthogonalizing the picked eigenvectors leads to an estimation of the subspace spanned by $\{w_{1}^{*}\;w_{2}^{*}\;\cdots\;w_{k}^{*}\}$. Also note that forming $\widehat{P}_{2}$ takes $O(n\cdot d^{2})$ time and each step of the power method doing a multiplication between a $d\times d$ matrix and a $d\times k$ matrix takes $k\cdot d^{2}$ time by a naive implementation. Here we reduce this complexity from $O((k+n)d^{2})$ to $O(knd)$. The idea is to compute each step of the power method without explicitly forming $\widehat{P}_{2}$. We take $P_{2}=M_{2}$ as an example; other cases are similar. In Algorithm 3, let the step $\widehat{P}_{2}V$ be calculated as $\widehat{P}_{2}V=\frac{1}{|S|}\sum_{(x,y)\in S}y(x(x^{\top}V)-V)$. Now each iteration only needs $O(knd)$ time. Furthermore, the number of iterations required will be a small number, since the power method has a linear convergence rate and as an initialization method we don’t need a very accurate solution. The detailed algorithm is shown in Algorithm 3. The approximation error bound of $\widehat{P}_{2}$ to $P_{2}$ is provided in Lemma E.5. Lemma E.6 provides the theoretical bound for Algorithm 3.