Provable ICA with Unknown Gaussian Noise, and Implications for Gaussian Mixtures and Autoencoders

Introduction

We present an algorithm (with rigorous performance guarantees) for a basic statistical problem. Suppose $\eta$ is an independent $n$-dimensional Gaussian random variable with an unknown covariance matrix $\Sigma$ and $A$ is an unknown $n\times n$ matrix. We are given samples of the form $y=Ax+\eta$ where $x$ is a random variable whose components are independent and have a fourth moment strictly less than that of a standard Gaussian random variable. The most natural case is when $x$ is chosen uniformly at random from $\{+1,-1\}^{n}$, although our algorithms in even the more general case above. Our goal is to reconstruct an additive approximation to the matrix $A$ and the covariance matrix $\Sigma$ running in time and using a number of samples that is polynomial in $n$ and $\frac{1}{\epsilon}$, where $\epsilon$ is the target precision (see Theorem 1.1) This problem arises in several research directions within machine learning: Independent Component Analysis (ICA), Deep Learning, Gaussian Mixture Models (GMM), etc. We describe these connections next, and known results (focusing on algorithms with provable performance guarantees, since that is our goal).

Most obviously, the above problem can be seen as an instance of Independent Component Analysis (ICA) with unknown Gaussian noise. ICA has an illustrious history with applications ranging from econometrics, to signal processing, to image segmentation. The goal generally involves finding a linear transformation of the data so that the coordinates are as independent as possible . This is often accomplished by finding directions in which the projection is “non-Gaussian” . Clearly, if the datapoint $y$ is generated as $Ax$ (i.e., with no noise $\eta$ added) then applying linear transformation $A^{-1}$ to the data results in samples $A^{-1}y$ whose coordinates are independent. This restricted case was considered by Comon and Frieze, Jerrum and Kannan , and their goal was to recover an additive approximation to $A$ efficiently and using a polynomial number of samples. (We will later note a gap in their reasoning, albeit fixable by our methods. See also recent papers by Anandkumar et al., Hsu and Kakade, that do not use local search and avoids this issue.) To the best of our knowledge, there are currently no known algorithms with provable guarantees for the more general case of ICA with Gaussian noise (this is especially true if the covariance matrix is unknown, as in our problem), although many empirical approaches are known. (eg. , the issue of “empirical” vs “rigorous” is elaborated upon after Theorem 1.1.)

The second view of our problem is as a concisely described Gaussian Mixture Model. Our data is generated as a mixture of $2^{n}$ identical Gaussian components (with an unknown covariance matrix) whose centers are the points $\left\{Ax:x\in\left\{-1,1\right\}^{n}\right\}$, and all mixing weights are equal. Notice, this mixture of $2^{n}$ Gaussians admits a concise description using $O(n^{2})$ parameters. The problem of learning Gaussian mixtures has a long history, and the popular approach in practice is to use the EM algorithm , though it has no worst-case guarantees (the method may take a very long time to converge, and worse, may not always converge to the correct solution). An influential paper of Dasgupta initiated the program of designing algorithms with provable guarantees, which was improved in a sequence of papers . But in the current setting, it is unclear how to apply any of the above algorithms (including $EM$) since the trivial application would keep track of exponentially many parameters – one for each component. Thus, new ideas seem necessary to achieve polynomial running time.

The third view of our problem is as a simple form of autoencoding . This is a central notion in Deep Learning, where the goal is to obtain a compact representation of a target distribution using a multilayered architecture, where a complicated function (the target) can be built up by composing layers of a simple function (called the autoencoder ). The main tenet is that there are interesting functions which can be represented concisely using many layers, but would need a very large representation if a “shallow” architecture is used instead). This is most useful for functions that are “highly varying” (i.e. cannot be compactly described by piecewise linear functions or other “simple” local representations). Formally, it is possible to represent using just (say) $n^{2}$ parameters, some distributions with $2^{n}$ “varying parts” or “interesting regions.” The Restricted Boltzmann Machine (RBM) is an especially popular autoencoder in Deep Learning, though many others have been proposed. However, to the best of our knowledge, there has been no successful attempt to give a rigorous analysis of Deep Learning. Concretely, if the data is indeed generated using the distribution represented by an RBM, then do the popular algorithms for Deep Learning learn the model parameters correctly and in polynomial time? Clearly, if the running time were actually found to be exponential in the number of parameters, then this would erode some of the advantages of the compact representation.

We give a provable algorithm for ICA with unknown Gaussian noise. We have not made an attempt to optimize the quoted running time of this model, but we emphasize that this is in fact the first algorithm with provable guarantees for this problem and moreover we believe that in practice our algorithm will run almost as fast as the usual ICA algorithms, which are its close relatives.

There is an algorithm that recovers the unknown $A$ and $\Sigma$ up to additive error $\epsilon$ in each entry in time that is polynomial in $n,\|A\|_{2},\|\Sigma\|_{2},\nicefrac{{1}}{{\epsilon}},\nicefrac{{1}}{{\lambda_{\min}(A)}}$ where $\|\cdot\|_{2}$ denotes the operator norm and $\lambda_{\min}(\cdot)$ denotes the smallest eigenvalue.

The classical approach for ICA initiated in Comon and Frieze, Jerrum and Kannan ) is for the noiseless case in which $y=Ax$. The first step is whitening, which applies a suitable linear transformation that makes the variance the same in all directions, thus reducing to the case where $A$ is a rotation matrix. Given samples $y=Rx$ where $R$ is a rotation matrix, the rows of $R$ can be found in principle by computing the vectors $u$ that are local minima of $E[(u\cdot y)^{4}]$. Subsequently, a number of works (see e.g. ) have focused on giving algorithms that are robust to noise. A popular approach is to use the fourth order cumulant (as an alternative to the fourth order moment) as a method for “denoising,” or any one of a number of other functionals whose local optima reveal interesting directions. However, theoretical guarantees of these algorithms are not well understood.

The above procedures in the noise-free model can almost be made rigorous (i.e., provably polynomial running time and number of samples), except for one subtlety: it is unclear how to use local search to find all optima in polynomial time. In practice, one finds a single local optimum, projects to the subspace orthogonal to it and continues recursively on a lower-dimensional problem. However, a naive implementation of this idea is unstable since approximation errors can accumulate badly, and to the best of our knowledge no rigorous analysis has been given prior to our work. (This is not a technicality: in some similar settings the errors are known to blow up exponentially .) One of our contributions is a modified local search that avoids this potential instability and finds all local optima in this setting. (Section 4.2.)

Our major new contribution however is dealing with noise that is an unknown Gaussian. This is an important generalization, since many methods used in ICA are quite unstable to noise (and a wrong estimate for the covariance could lead to bad results). Here, we no longer need to assume we know even rough estimates for the covariance. Moreover, in the context of Gaussian Mixture Models this generalization corresponds to learning a mixture of many Gaussians where the covariance of the components is not known in advance.

We design new tools for denoising and especially whitening in this setting. Denoising uses the fourth order cumulant instead of the fourth moment used in and whitening involves a novel use of the Hessian of the cumulant. Even then, we cannot reduce to the simple case $y=Rx$ as above, and are left with a more complicated functional form (see “quasi-whitening” in Section 2.) Nevertheless, we can reduce to an optimization problem that can be solved via local search, and which remains amenable to a rigorous analysis. The results of the local optimization step can be then used to simplify the complicated functional form and recover $A$ as well as the noise $\Sigma$. We defer many of our proofs to the supplementary material section, due to space constraints.

In order to avoid cluttered notation, we have focused on the case in which $x$ is chosen uniformly at random from $\{-1,+1\}^{n}$, although our algorithm and analysis work under the more general conditions that the coordinates of $x$ are (i) independent and (ii) have a fourth moment that is less than three (the fourth moment of a Gaussian random variable). In this case, the functional $P(u)$ (see Lemma 2.2) will take the same form but with weights depending on the exact value of the fourth moment for each coordinate. Since we already carry through an unknown diagonal matrix $D$ throughout our analysis, this generalization only changes the entries on the diagonal and the same algorithm and proof apply.

Denoising and quasi-whitening

As mentioned, our approach is based on the fourth order cumulant. The cumulants of a random variable are the coefficients of the Taylor expansion of the logarithm of the characteristic function . Let $\kappa_{r}(X)$ be the $r^{th}$ cumulant of a random variable $X$. We make use of:

(i) If $X$ has mean zero, then $\kappa_{4}(X)=\mathop{\bf E\/}[X^{4}]-3\mathop{\bf E\/}[X^{2}]^{2}$. (ii) If $X$ is Gaussian with mean $\mu$ and variance $\sigma^{2}$, then $\kappa_{1}(X)=\mu$, $\kappa_{2}(X)=\sigma^{2}$ and $\kappa_{r}(X)=0$ for all $r>2$. (iii) If $X$ and $Y$ are independent, then $\kappa_{r}(X+Y)=\kappa_{r}(X)+\kappa_{r}(Y)$.

Proof: The crucial observation is that $u^{T}y=u^{T}Ax+u^{T}\eta$ is the sum of two independent random variables, $Ax$ and $\eta$ and that $P(u)=-\kappa_{4}(u^{T}Ax+u^{T}\eta)=-\kappa_{4}(u^{T}Ax)-\kappa_{4}(u^{T}\eta)=-\kappa_{4}(u^{T}Ax)$. So in fact, the functional $P(u)$ is invariant under additive Gaussian noise independent of the variance matrix $\Sigma$. This vastly simplifies our computation:

Furthermore $\mathop{\bf E\/}[(u^{T}Ax)^{2}]^{2}=(u^{T}AA^{T}u)^{2}$ and we conclude that

In prior works on ICA, whitening refers to reducing to the case where $y=Rx$ for some some rotation matrix $R$. Here we give a technique to reduce to the case where $y=RDx+\eta^{\prime}$ where $\eta^{\prime}$ is some other Gaussian noise (still unknown), $R$ is a rotation matrix and $D$ is a diagonal matrix that depends upon $A$. We call this quasi-whitening. Quasi-whitening suffices for us since local search using the objective function $\kappa_{4}(u^{T}y)$ will give us (approximations to) the rows of $RD$, from which we will be able to recover $A$.

Quasi-whitening involves computing the Hessian of $P(u)$, which recall is the matrix of all 2nd order partial derivatives of $P(u)$. Throughout this section, we will denote the Hessian operator by $\mathcal{H}$. In matrix form, the Hessian of $P(u)$ is

Let $D_{A}(u)$ be a diagonal matrix in which the $k^{th}$ entry is $24(A_{k}\cdot u)^{2}$.

Of course, the exact Hessian of $P(u)$ is unavailable and we will instead compute an empirical approximation $\widehat{P}(u)$ to $P(u)$ (given many samples from the distribution), and we will show that the Hessian of $\widehat{P}(u)$ is a good approximation to the Hessian of $P(u)$.

Given $2N$ samples $y_{1},y^{\prime}_{1},y_{2},y^{\prime}_{2}...,y_{N},y^{\prime}_{N}$ of the random variable $y$, let

Our first step is to show that the expectation of the Hessian of $\widehat{P}(u)$ is exactly the Hessian of $P(u)$. In fact, since the expectation of $\widehat{P}(u)$ is exactly $P(u)$ (and since $\widehat{P}(u)$ is an analytic function of the samples and of the vector $u$), we can interchange the Hessian operator and the expectation operator. Roughly, one can imagine the expectation operator as an integral over the possible values of the random samples, and as is well-known in analysis, one can differentiate under the integral provided that all functions are suitably smooth over the domain of integration.

$\mathop{\bf E\/}_{y,y^{\prime}}[-(u^{T}y)^{4}+3(u^{T}y)^{2}(u^{T}y^{\prime})^{2}]=P(u)$

This claim follows immediately from the definition of $P(u)$, and since $y$ and $y^{\prime}$ are independent.

$\mathcal{H}(P(u))=\mathop{\bf E\/}_{y,y^{\prime}}[\mathcal{H}(-(u^{T}y)^{4}+3(u^{T}y)^{2}(u^{T}y^{\prime})^{2})]$

Next, we compute the two terms inside the expectation:

$\mathcal{H}((u^{T}y)^{4})=12(u^{T}y)^{2}yy^{T}$

$\mathcal{H}((u^{T}y)^{2}(u^{T}y^{\prime})^{2})=2(u^{T}y^{\prime})^{2}yy^{T}+2(u^{T}y)^{2}y^{\prime}(y^{\prime})^{T}+4(u^{T}y)(u^{T}y^{\prime})(y(y^{\prime})^{T}+(y^{\prime})y^{T})$

Let $\lambda_{min}(A)$ denote the smallest eigenvalue of $A$. Our analysis also requires bounds on the entries of $D_{A}(u_{0})$:

If $u_{0}$ is chosen uniformly at random then with high probability for all $i$,

Proof: We can bound $\max_{i=1}^{n}|A_{i}\cdot u|$ by $\max_{i=1}^{n}\|A_{i}\|_{2}\frac{\log n}{\sqrt{n}}$ thus the bound for $\max_{i=1}^{n}(D_{A}(u_{0}))_{i,i}$ follows. Note that with high probability the minimum absolute value of $n$ Gaussian random variables is at least $1/n^{2}$, hence $\min_{i=1}^{n}(D_{A}(u_{0}))_{i,i}\geq\min_{i=1}^{n}\|A_{i}\|_{2}^{2}n^{-4}$. $\blacksquare$

Proof: First we consider each entry of the matrix updates. For example, the variance of any entry in $\mathcal{H}((u^{T}y)^{4})=12(u^{T}y)^{2}yy^{T}$ can be bounded by $\|y\|_{2}^{8}$, which we can bound by $\mathop{\bf E\/}[\|y\|_{2}^{8}]\leq O(\mathop{\bf E\/}[\|Ax\|_{2}^{8}+\|\eta\|_{2}^{8}])$. This can be bounded by $O(n^{4}(\|A\|_{2}^{8}+\|\Sigma\|_{2}^{4}))$. This is also an upper bound for the variance (of any entry) of any of the other matrix updates when computing $\mathcal{H}(\widehat{P}(u_{0}))$.

Suppose that $(1-\epsilon)AD_{A}(u_{0})A^{T}\preceq\widehat{M}\preceq(1+\epsilon)AD_{A}(u_{0})A^{T}$, and let $\widehat{M}=BB^{T}$. Then there is a rotation matrix $R^{*}$ such that $\|B^{-1}AD_{A}(u_{0})^{1/2}-R^{*}\|_{F}\leq\sqrt{n}\epsilon$.

The intuition is: if any of the singular values of $B^{-1}AD_{A}(u_{0})^{1/2}$ are outside the range $[1-\epsilon,1+\epsilon]$, we can find a unit vector $x$ where the quadratic forms $x^{T}AD_{A}(u_{0})A^{T}x$ and $x^{T}\widehat{M}x$ are too far apart (which contradicts the condition of the lemma). Hence the singular values of $B^{-1}AD_{A}(u_{0})^{1/2}$ can all be set to one without changing the Froebenius norm of $B^{-1}AD_{A}(u_{0})^{1/2}$ too much, and this yields a rotation matrix.

Proof: Let $M=AD_{A}(u_{0})A^{T}$ and let $C=AD_{A}(u_{0})^{1/2}$, and so $M=CC^{T}$ and $\widehat{M}=BB^{T}$. The condition $(1-\epsilon)M\preceq\widehat{M}\preceq(1+\epsilon)M$ is well-known to be equivalent to the condition that for all vectors $x$, $(1-\epsilon)x^{T}Mx\leq x^{T}\widehat{M}x\leq(1+\epsilon)x^{T}Mx$.

Suppose for the sake of contradiction that $S=B^{-1}C$ has a singular value outside the range $[1-\epsilon,1+\epsilon]$. Assume (without loss of generality) that $S$ has a singular value strictly larger than $1+\epsilon$ (and the complementary case can be handled analogously). Hence there is a unit vector $y$ such that $y^{T}SS^{T}y>1+\epsilon$. But since $BSS^{T}B^{T}=CC^{T}$, if we set $x^{T}=y^{T}B^{-1}$ then we have $x^{T}\widehat{M}x=x^{T}BB^{T}x=y^{T}y=1$ but $x^{T}Mx=x^{T}CC^{T}x=x^{T}BSS^{T}B^{T}x=y^{T}SS^{T}y>1+\epsilon$. This is a contradiction and so we conclude that all of the singular values of $B^{-1}C$ are in the range $[1-\epsilon,1+\epsilon]$.

Let $U\Sigma V^{T}$ be the singular value decomposition of $B^{-1}C$. If we set all of the diagonal entries in $\Sigma$ to $1$ we obtain a rotation matrix $R^{*}=UV^{T}$. And since the singular values of $B^{-1}C$ are all in the range $[1-\epsilon,1+\epsilon]$, we can bound the Froebenius norm of $B^{-1}C-R^{*}$: $\|B^{-1}C-R^{*}\|_{F}\leq\sqrt{n}\epsilon$, as desired. $\blacksquare$

Our algorithm (and notation)

In this section we describe our overall algorithm. It uses as a blackbox the denoising and quasi-whitening already described above, as well as a routine for computing all local maxima of some “well-behaved” functions which is described later in Section 4.

Notation: Placing a hat over a function corresponds to an empirical approximation that we obtain from random samples. This approximation introduces error, which we will keep track of.

Step 2: Take $2N$ samples $y_{1},y_{2},...,y_{N},y_{1}^{\prime},y_{2}^{\prime},...,y_{N}^{\prime}$, and let

which is an empirical estimation of $P^{\prime}(u)$.

Step 3: Use the procedure AllOPT($\widehat{P}^{\prime}(u),\beta,\delta^{\prime},\beta^{\prime},\delta^{\prime}$) of Section 4 to compute all $n$ local maxima of the function $\widehat{P}^{\prime}(u)$.

Step 4: Let $R$ be the matrix whose rows are the $n$ local optima recovered in the previous step. Use procedure Recover of Section 5 to find $A$ and $\Sigma$.

Explanation: Step 1 uses the transformation $B^{-1}$ computed in the previous Section to quasi-whiten the data. Namely, we consider the sequence of samples $z=B^{-1}y$, which are therefore of the form $R^{\prime}Dx+\eta^{\prime}$ where $\eta=B^{-1}\eta$, $D=D_{A}(u_{0})$ and $R^{\prime}$ is close to a rotation matrix $R^{*}$ (by Lemma 2.11). In Step 2 we look at $\kappa_{4}((u^{T}z))$, which effectively denoises the new samples (see Lemma 2.2), and thus is the same as $\kappa_{4}(R^{\prime}D^{-1/2}x)$. Let $P^{\prime}(u)=\kappa_{4}(u^{T}z)=\kappa_{4}(u^{T}B^{-1}y)$ which is easily seen to be $E[(u^{T}R^{\prime}D^{-1/2}x)^{4}]$. Step 2 estimates this function, obtaining $\widehat{P}^{\prime}(u)$. Then Step 3 tries to find local optima via local search. Ideally we would have liked access to the functional $P^{*}(u)=(u^{T}R^{*}x)^{4}$ since the procedure for local optima works only for true rotations. But since $R^{\prime}$ and $R^{*}$ are close we can make it work approximately with $\widehat{P}^{\prime}(u)$, and then in Step 4 use these local optima to finally recover $A$.

Proof: In Step 1, by Lemma 2.11 we know once we use $z=B^{-1}y$, the whitened function $P^{\prime}(u)$ is inverse polynomially close to $P^{*}(u)$. Then by Lemma 5.3, the function $\widehat{P^{\prime}}(u)$ we get in Step 2 is inverse polynomially close to $P^{\prime}(u)$ and $P^{*}(u)$. Theorem 4.6 and Lemma 5.5 show that given $\widehat{P^{\prime}}(u)$ inverse polynomially close to $P^{*}(u)$, Algorithm 2: : AllOPT finds all local maxima with inverse polynomial precision. Finally by Theorem 5.6 we know $A$ and $W$ are recovered correctly up to additive $\epsilon$ error in Frobenius norm. The running time and sampling complexity of the algorithm is polynomial because all parameters in these Lemmas are polynomially related. $\blacksquare$

Note that here we recover $A$ up to a permutation of the columns and sign-flips. In general, this is all we can hope for since the distribution of $x$ is also invariant under these same operations. Also, the dependence of our algorithm on the various norms (of $A$ and $\Sigma$) seems inherent since our goal is to recover an additive approximation, and as we scale up $A$ and/or $\Sigma$, this goal becomes a stronger relative guarantee on the error.

Framework for iteratively finding all local maxima

Given that we can only ever hope for an additive approximation to a local maximum, one should be concerned about how the error accumulates when our goal is to find all local maxima. In fact, a naive strategy is to project onto the subspace orthogonal to the directions found so far, and continue in this subspace. However, such an approach seems to accumulate errors badly (the additive error of the last local maxima found is exponentially larger than the error of the first). Rather, the crux of our analysis is a novel method for bounding how much the error can accumulate (by refining old estimates).

Our strategy is to first find a local maximum in the orthogonal subspace, then run the local optimization algorithm again (in the original $n$-dimensional space) to “refine” the local maximum we have found. The intuition is that since we are already close to a particular local maxima, the local search algorithm cannot jump to some other local maxima (since this would entail going through a valley).

Throughout this section, we will assume that we are given oracle access to a function $f(u)$ and its gradient and Hessian. The procedure is also given a starting point $u_{s}$, a search range $\beta$, and a step size $\delta$. For simplicity in notation we define the following projection operator.

$\mbox{Proj}_{\perp u}(v)=v-(u^{T}v)u$, $\mbox{Proj}_{\perp u}(M)=M-(u^{T}Mu)uu^{T}$.

The basic step the algorithm is a modification of Newton’s method to find a local improvement that makes progress so long as the current point $u$ is far from a local maxima. Notice that if we add a small vector to $u$, we do not necessarily preserve the norm of $u$. In order to have control over how the norm of $u$ changes, during local optimization step the algorithm projects the gradient $\nabla f$ and Hessian $\mathcal{H}(f)$ to the space perpendicular to $u$. There is also an additional correction term $-\partial/\partial_{u}f(u)\cdot\|\xi\|^{2}/2$. This correction term is necessary because the new vector we obtain is $(u+\xi)/\left\lVert(u+\xi)\right\rVert_{2}$ which is close to $u-\|\xi\|_{2}^{2}/2\cdot u+\xi+O(\beta^{3})$. Step 2 of the algorithm is just maximizing a quadratic function and can be solved exactly using Lagrangian Multiplier method. To increase efficiency it is also acceptable to perform an approximate maximization step by taking $\xi$ to be either aligned with the gradient $\mbox{Proj}_{\perp u}\nabla f(u)$ or the largest eigenvector of $\mbox{Proj}_{\perp u}(\mathcal{H}(f(u)))$.

The algorithm is guaranteed to succeed in polynomial time when the function is Locally Improvable and Locally Approximable:

A function $f(u)$ is locally approximable, if its third order derivatives exist and for any $u$ and any direction $v$, the third order derivative of $f$ at point $u$ in the direction of $v$ is bounded by $0.01\delta/\beta^{3}$.

The analysis of the running time of the procedure comes from local Taylor expansion. When a function is Locally Approximable it is well approximated by the gradient and Hessian within a $\beta$ neighborhood. The following theorem from showed that the two properties above are enough to guarantee the success of a local search algorithm even when the function is only approximated.

If $|f(u)-f^{*}(u)|\leq\delta/8$, the function $f^{*}(u)$ is $(\gamma,\beta,\delta)$-Locally Improvable, $f(u)$ is $(\beta,\delta)$ Locally Approximable, then Algorithm 1 will find a vector $v$ that is $\gamma$ close to some local maximum. The running time is at most $O((n^{2}+T)\max f^{*}/\delta)$ where $T$ is the time to evaluate the function $f$ and its gradient and Hessian.

2 Finding all local maxima

Now we consider how to find all local maxima of a given function $f^{*}(u)$. The crucial condition that we need is that all local maxima are orthogonal (which is indeed true in our problem, and is morally true when using local search more generally in ICA). Note that this condition implies that there are at most $n$ local maxima.Technically, there are $2n$ local maxima since for each direction $u$ that is a local maxima, so too is $-u$ but this is an unimportant detail for our purposes. In fact we will assume that there are exactly $n$ local maxima. If we are given an exact oracle for $f^{*}$ and can compute exact local maxima then we can find all local maxima easily: find one local maximum, project the function into the orthogonal subspace, and continue to find more local maxima.

The following theorem gives sufficient conditions under which the above algorithm finds all local maxima, making precise the intuition given at the beginning of this section.

Orthogonal Local Maxima: The function has $n$ local maxima $v^{*}_{i},$ and they are orthogonal to each other.

Locally Improvable: $f^{*}$ is $(\gamma,\beta,\delta)$ Locally Improvable.

Improvable Projection: The projection of the function to any subspace spanned by a subset of local maxima is $(\gamma^{\prime},\beta^{\prime},\delta^{\prime})$ Locally Improvable. The step size $\delta^{\prime}\geq 10\delta$.

Lipschitz: If two points $\left\lVert u-u^{\prime}\right\rVert_{2}\leq 3\sqrt{n}\gamma$, then the function value $|f^{*}(u)-f^{*}(u^{\prime})|\leq\delta^{\prime}/20$.

Attraction Radius: Let $Rad\geq 3\sqrt{n}\gamma+\gamma^{\prime}$, for any local maximum $v^{*}_{i}$, let $T$ be $\min f^{*}(u)$ for $\left\lVert u-v^{*}_{i}\right\rVert_{2}\leq Rad$, then there exist a set $U$ containing $\left\lVert u-v^{*}_{i}\right\rVert_{2}\leq 3\sqrt{n}\gamma+\gamma^{\prime}$ and does not contain any other local maxima, such that for every $u$ that is not in $U$ but is $\beta$ close to $U$, $f^{*}(u)<T$.

If we are given function $f$ such that $|f(u)-f^{*}(u)|\leq\delta/8$ and $f$ is both $(\beta,\delta)$ and $(\beta^{\prime},\delta^{\prime})$ Locally Approximable, then Algorithm 2 can find all local maxima of $f^{*}$ within distance $\gamma$.

To prove this theorem, we first notice the projection of the function $f$ in Step 3 of the algorithm should be close to the projection of $f^{*}$ to the remaining local maxima. This is implied by Lipschitz condition and is formally shown in the following two lemmas. First we prove a “coupling” between the orthogonal complement of two close subspaces:

Proof: Let $S_{1}$ be $span\{v_{1},v_{2},...,v_{k}\}$, $S_{2}$ be $span\{v^{*}_{1},v^{*}_{2},...,v^{*}_{k}\}$ and $S_{1}^{\perp}$, $S_{2}^{\perp}$ be their orthogonal subspaces respectively. We first prove that for any unit vector $v\in S_{1}^{\perp}$, there is another unit vector $v^{\prime}\in S_{2}^{\perp}$ so that $v^{T}v^{\prime}\geq 1-4n\gamma^{2}$. In fact, we can take $v^{\prime}$ to be the unit vector along the projection of $v$ in $S_{2}^{\perp}$. To bound the length of the projection, we instead bound the length of projection to $S_{2}$. Since we know $v_{i}^{T}v^{\prime}=0$ for $i\leq k$ and $\left\lVert v_{i}-v^{*}_{i}\right\rVert\leq\gamma$, it must be that $(v^{*}_{i})^{T}v^{\prime}\leq 2\gamma$ when $\gamma<0.01$. So the projection of $v^{\prime}$ in $S_{2}$ has length at most $2\sqrt{n}\gamma$ and hence the projection of $v^{\prime}$ in $S_{2}^{\perp}$ has length at least $1-4n\gamma^{2}$.

The second equality cannot hold if these vectors are not orthogonal. And for any $w$,

So we conclude that the distance between these two vectors is at most $3\sqrt{n}\gamma$. $\blacksquare$

Using this lemma we see that the projected function is close to the projection of $f^{*}$ in the span of the rest of local maxima:

Let $g^{*}$ be the projection of $f^{*}$ into the space spanned by the rest of local maxima, then $|g^{*}(w)-g(w)|\leq\delta/8+\delta^{\prime}/20\leq\delta^{\prime}/8$.

Proof: The proof is straight forward because $|g^{*}(w)-g(w)|\leq|f^{*}(u)-f(u)|+|f^{*}(u)-f^{*}(u^{\prime})|$ for some $\left\lVert u-u^{\prime}\right\rVert_{2}\leq 3\sqrt{n}\gamma$, we know the first one is at most $\delta/8$ and the second one is at most $\delta^{\prime}/20$ by Lipschitz Condition. $\blacksquare$

Now we are ready to prove the main theorem.

Proof: [Theorem 4.6] By Theorem 4.4 the first column is indeed $\gamma$ close to a local maximum. We then prove by induction that if $v_{1}$, $v_{2}$, …, $v_{k}$ are $\gamma$ close to different local maxima, then $v_{k+1}$ must be close to a new local maximum.

By Lemma 4.8 we know $g_{k+1}$ is $(\gamma^{\prime},\beta^{\prime},\delta^{\prime})$ Locally Improvable, and because it is a projection of $f$ its derivatives are also bounded so it is $(\beta^{\prime},\delta^{\prime})$ Locally Approximable. By Theorem 4.4 $u^{\prime}$ must be $\gamma^{\prime}$ close to local maximum for the projected function. Then since the projected space is close to the space spanned by the rest of local maxima, $u^{\prime}$ is in fact $\gamma^{\prime}+3\sqrt{n}\gamma$ close to $v^{*}_{k+1}$ (here again we are reindexing the local maxima wlog.).

Now we use the Attraction Radius property, since $u$ is currently in $U$, $f^{*}(u)\geq T$, and each step we go to a point $u^{\prime}$ such that $\left\lVert u^{\prime}-u\right\rVert\leq\beta$ and $f^{*}(u^{\prime})>f^{*}(u)\geq T$. The local search in Algorithm 1 can never go outside $U$, therefore it must find the local maximum $v^{*}_{k+1}$. $\blacksquare$

Local search on the fourth order cumulant

Next, we prove that the fourth order cumulant $P^{*}(u)$ satisfies the properties above. Then the algorithm given in the previous section will find all of the local maxima, which is the missing step in our main goal: learning a noisy linear transformation $Ax+\eta$ with unknown Gaussian noise. We first use a theorem from to show that properties for finding one local maxima is satisfied.

Also, for notational convenience we set $d_{i}=2D_{A}(u_{0})_{i,i}^{-2}$ and let $d_{min}$ and $d_{max}$ denote the minimum and maximum values (bounds on these and their ratio follow from Claim 2.9). Using this notation $P^{*}(u)=\sum_{i=1}^{n}d_{i}(u^{T}R^{*}_{i})^{4}$.

When $\beta<d_{min}/10d_{max}n^{2}$, the function $P^{*}(u)$ is $(3\sqrt{n}\beta,\beta,P^{*}(u)\beta^{2}/100)$ Locally Improvable and $(\beta,d_{min}\beta^{2}/100n)$ Locally Approximable. Moreover, the local maxima of the function is exactly $\{\pm R^{*}_{i}\}$.

Proof: The proof appears in . Here for completeness we show the proof using our notations.

First we establish that $P^{*}(u)$ is Locally Improvable.Observe that this desirada is invariant under rotation, so we need only prove the theorem for $P^{*}(v)=\sum_{i=1}^{n}d_{i}v_{i}^{4}$. The gradient of the function is $\nabla P^{*}(v)=4(d_{1}v_{1}^{3},d_{2}v_{2}^{3},...,d_{n}v_{n}^{3})$. The inner product of $\nabla P^{*}(v)$ and $v$ is exactly $4\sum_{i=1}^{n}d_{i}v_{i}^{4}=4P^{*}(v)$. Therefore the projected gradient $\phi=\mbox{Proj}_{\perp v}\nabla P^{*}(v)$ has coordinate $\phi_{i}=4v_{i}(d_{i}v_{i}^{2}-P^{*}(v))$. Furthermore, the Hessian $H=\mathcal{H}(P^{*}(v))$ is a diagonal matrix whose $(i,i)^{th}$ entry is $12d_{i}v_{i}^{2}$.

Consider the case in which $\left\lVert\phi\right\rVert\geq P^{*}(v)\beta/4$. We can obtain an improvement to $P^{*}(v)\beta^{2}/100$ because we can take $\xi$ in the direction of $\phi$ and with $\left\lVert\xi\right\rVert_{2}=\beta/20$. The contribution of the Hessian term is nonnegative and the third term $-2P^{*}(u)\left\lVert\xi\right\rVert_{2}^{2}$ is small in comparison.

Hence, we can assume $\left\lVert\phi\right\rVert\leq P^{*}(v)\beta/4$. Now let us write out the expression of $\left\lVert\phi\right\rVert^{2}$

In particular every term $v_{i}^{2}(d_{i}v_{i}^{2}-P^{*}(v))^{2}$ must be at most $\beta^{2}(P^{*}(v))^{2}/16.$. Thus for any $i$, either $v_{i}^{2}\leq\beta^{2}$ or $(d_{i}v_{i}^{2}-P^{*}(v))^{2}\leq(P^{*}(v))^{2}/16$.

If there are at least 2 coordinates $k$ and $l$ such that $(d_{i}v_{i}^{2}-P^{*}(v))^{2}\leq(P^{*}(v))^{2}/16$, then we know for these two coordinates $v_{i}^{2}\in[0.75P^{*}(v)/d_{i},1.25P^{*}(v)/d_{i}]$. We choose the vector $\xi$ so that $\xi_{k}=\tau v_{l}$ and $\xi_{l}=-\tau v_{k}$. Wlog assume $\xi\cdot\phi\geq 0$ otherwise we use $-\xi$. Take $\tau$ so that $\tau^{2}(v_{l}^{2}+v_{k}^{2})=\beta^{2}$. Clearly $\left\lVert\xi\right\rVert=\beta$ and $\xi\cdot v=0$ so $\xi$ is a valid solution. Also $\tau^{2}$ is lower bounded by $\beta^{2}/(v_{l}^{2}+v_{k}^{2})\geq\frac{4}{5}\frac{\beta^{2}}{P^{*}(u)(1/d_{l}+1/d_{k})}$.

In the remaining case, all of the coordinates except for at most one satisfy $v_{i}^{2}\leq\beta^{2}$. Since we assumed $\beta^{2}<\frac{1}{n}$, there must be one of the coordinate $v_{k}$ that is large, and it is at least $1-n\beta^{2}$. Thus the distance of this vector to the local maxima $e_{k}$ is at most $3\sqrt{n}\beta$. $\blacksquare$

We then observe that given enough samples, the empirical mean $\widehat{P}^{\prime}(u)$ is close to $P^{*}(u)$. For concentration we require every degree four term $z_{i}z_{j}z_{k}z_{l}$ has variance at most $Z$.

$Z=O(d_{min}^{2}\lambda_{min}(A)^{8}\|\Sigma\|_{2}^{4}+d_{min}^{2})$.

Proof: We will start by bounding $\mathop{\bf E\/}[(z_{i}z_{j}z_{k}z_{l})^{2}]\leq\mathop{\bf E\/}[(z_{i}^{8}+z_{j}^{8}+z_{k}^{8}+z_{l}^{8})]$. Furthermore $\mathop{\bf E\/}[z_{i}^{8}]\leq O(\mathop{\bf E\/}[(B^{-1}Ax)_{i}^{8}+(B^{-1}\eta)_{i}^{8}])$. Next we bound $\mathop{\bf E\/}[(B^{-1}\eta)_{i}^{8}]$, which is just the eighth moment of a Gaussian with variance at most $\|B^{-1}\Sigma B^{-T}\|_{2}\leq\|B^{-1}\|_{2}^{2}\|\Sigma\|_{2}\leq d_{min}^{1/2}\lambda_{min}(A)^{-2}\|\Sigma\|_{2}$. Hence we can bound this term by$O(\|B^{-1}\Sigma B^{-T}\|_{2}^{4})=O(d_{min}^{2}\lambda_{min}(A)^{8}\|\Sigma\|_{2}^{4})$. Finally the remaining term $\mathop{\bf E\/}[(B^{-1}Ax)_{i}^{8}]$ can be bounded by $O(d_{min}^{2})$ because the variance of this random variable is only larger if we instead replace $x$ by an $n$-dimensional standard Gaussian. $\blacksquare$

Given $2N$ samples $y_{1},y_{2},...,y_{N},y_{1}^{\prime},y_{2}^{\prime},...,y_{N}^{\prime}$, suppose columns of $R^{\prime}=B^{-1}AD_{A}(u_{0})^{1/2}$ are $\epsilon$ close to the corresponding columns of $R^{*}$, with high probability the function $\widehat{P}^{\prime}(u)$ is $O(d_{max}n^{1/2}\epsilon+n^{2}(N/Z\log n)^{-1/2})$ close to the true function $P^{*}(u)$.

Proof: $\widehat{P}^{\prime}(u)$ is the empirical mean of $F(u,y,y^{\prime})=-(u^{T}B^{-1}y)^{4}+3(u^{T}B^{-1}y)^{2}(u^{T}B^{-1}y^{\prime})^{2}$. In Section 2 we proved that $P^{\prime}(u)=\mathop{\bf E\/}_{y,y^{\prime}}F(u,y,y^{\prime})=\sum_{i=1}^{n}2D^{-1/2}_{i,i}(u^{T}R_{i})^{4}=\sum_{i=1}^{n}\lambda_{i}(u^{T}R_{i})^{4}$. First, we demonstrate that $P^{\prime}(u)$ is close to $P^{*}(u)$, and then using concentration bounds we show that $\widehat{P}^{\prime}(u)$ is close to $P^{\prime}(u)$ (with high probability) over all $u$.

The first part is a simple application of Cauchy-Schwartz:

The first inequality uses the fact that $((u^{T}R^{\prime}_{i})^{2}+(u^{T}R^{*}_{i})^{2})\leq 3$, the second inequality uses the fact that when $\epsilon$ is small enough, $\left\lVert u^{T}R^{\prime}\right\rVert_{2}\leq 2$.

Next we prove that the empirical mean $\widehat{P}^{\prime}(u)$ is close to $P^{\prime}(u)$. The key point here is we need to prove this for all points $u$ since a priori we have no control over which directions local search will choose to explore. We accomplish this by considering $\widehat{P}^{\prime}(u)$ as a degree-4 polynomial over $u$ and prove that the coefficient of each monomial in $\widehat{P}^{\prime}(u)$ is close to the corresponding coefficient in $P^{\prime}(u)$. This is easy: the expectation of each coefficient of $F(u,y,y^{\prime})$ is equal to the correct coefficient, and the variance is bounded by $O(Z)$. The coefficients are also sub-Gaussian so by Bernstein’s inequality the probability that any coefficient of $\widehat{P}^{\prime}(u)$ deviates by more than $\epsilon^{\prime}$ (from its expectation) is at most $e^{-\Omega(\epsilon^{\prime 2}N/Z)}$. Hence when $N\geq O(Z\log n/\epsilon^{\prime 2})$ with high probability all the coefficients of $\widehat{P}^{\prime}(u)$ and $P^{\prime}(u)$ are $\epsilon^{\prime}$ close. So for any $u$:

Therefore $\widehat{P}^{\prime}(u)$ and $P^{*}(u)$ are $O(d_{max}n^{1/2}\epsilon+n^{2}(N/Z\log n)^{-1/2})$ close. $\blacksquare$

This proof can also be used to show that the derivatives of the function $\widehat{P}^{\prime}(u)$ is concentrated to the derivatives of the true function $P^{*}(u)$ because the derivatives are only related to coefficients, therefore $\widehat{P}^{\prime}(u)$ is also $(\beta,d_{min}\beta^{2}/100n)$ Locally Approximable.

The other properties required by Theorem 4.6 are also satisfied:

For any $\|u-u^{\prime}\|_{2}\leq r$, $|P^{*}(u)-P^{*}(u^{\prime})|\leq 5d_{max}n^{1/2}r$. All local maxima of $P^{*}$ has attraction radius $Rad\geq d_{min}/100d_{max}$.

Proof: The Lipschitz condition follows from the same Cauchy-Schwartz as appeared above. When two points $u$ and $u^{\prime}$ are of distance $r$, $|P^{*}(u)-P^{*}(u^{\prime})|\leq 5d_{max}n^{1/2}r$. Finally for the Attraction Radius, we know when $3\sqrt{n}\gamma+\gamma^{\prime}\leq d_{min}/100d_{max}$, we can just take the set $U$ to be $u^{T}R^{*}_{i}\geq 1-d_{min}/50d_{max}$. For all $u$ such that $u^{T}R^{*}_{i}\in[1-d_{min}/25d_{max},1-d_{min}/50d_{max}]$ (which contains the $\beta$ neighborhood of $U$), we know the value of $P^{*}(u)\leq T$. $\blacksquare$

Applying Theorem 4.6 we obtain the following Lemma (the parameters are chosen so that all properties required are satisfied):

Let $\beta^{\prime}=\Theta((d_{min}/d_{max})^{2})$, $\beta=\min\{\gamma n^{-1/2},\Omega((d_{min}/d_{max})^{4}n^{-3.5})\}$, then the procedure Recover($f,\beta,$ $d_{min}\beta^{2}/100n$ $,\beta^{\prime},d_{min}\beta^{\prime 2}/100n$) finds vectors $v_{1},v_{2},...,v_{n}$, so that there is a permutation matrix $\Pi$ and $k_{i}\in\{\pm 1\}$ and for all $i$: $\left\lVert v_{i}-(R\Pi\mbox{Diag}(k_{i}))^{*}_{i}\right\rVert_{2}\leq\gamma$.

After obtaining $\widehat{R}=\left[v_{1},v_{2},...,v_{n}\right]$ we can use Algorithm 3 to find $A$ and $\Sigma$:

Given a matrix $\widehat{R}$ such that there is permutation matrix $\Pi$ and $k_{i}\in\{\pm 1\}$ with $\|\widehat{R}_{i}-k_{i}(R^{*}\Pi)_{i}\|_{2}\leq\gamma$ for all $i$, Algorithm 3 returns matrix $\widehat{A}$ such that $\|\widehat{A}-A\Pi\mbox{Diag}(k_{i})\|_{F}\leq O(\gamma\left\lVert A\right\rVert_{2}^{2}n^{3/2}/\lambda_{min}(A))$. If $\gamma\leq O(\epsilon/\left\lVert A\right\rVert_{2}^{2}n^{3/2}\lambda_{min}(A))\times\min\{1/\left\lVert A\right\rVert_{2},1\}$, we also have $\|\widehat{\Sigma}-\Sigma\|_{F}\leq\epsilon$.

Recall that the diagonal matrix $D_{A}(u)$ is unknown (since it depends on $A$), but if we are given $R^{*}$ (or an approximation) and since $P^{*}(u)=\sum_{i=1}^{n}d_{i}(u^{T}R^{*}_{i})^{4}$, we can recover the matrix $D_{A}(u)$ approximately from computing $P^{*}(R^{*}_{i})$. Then given $D_{A}(u)$, we can recover $A$ and $\Sigma$ and this completes the analysis of our algorithm.

Proof: By Lemma 2.11 we know the columns of $R^{\prime}$ is close the the columns of $R$ (the parameters will be set so that the error is much smaller than $\gamma$), thus $\|\widehat{R}_{i}-k_{i}(R^{\prime}\Pi)_{i}\|_{2}\leq\gamma$. Applying Lemma 5.3 we obtain: $|\widehat{P}^{\prime}(\widehat{R}_{i})-P^{*}(\widehat{R}_{i})|\ll\gamma$. Furthermore, when $\|\widehat{R}_{i}-k_{i}R^{*}_{\Pi^{-1}(i)}\|_{2}\leq\gamma$ we know that $P^{*}(\widehat{R}_{i})/d_{\Pi^{-1}(i)}\in[1-3\gamma,1+3\gamma]$ (here we are abusing notation and use the permutation matrix as a permutation). Hence $\widehat{D}_{A}(u)_{i,i}/\left(D_{A}(u)\right)_{\Pi^{-1}(i),\Pi^{-1}(i)}\in[1-3\gamma,1+3\gamma]$. We have: