text:
\section{Analytic Formula for Population Risks}\label{sec:formula}

\subsection{Basics on Hermite Polynomials}\label{sec:basic-hermite}

In this section, we briefly review Hermite polynomials and Fourier analysis on Gaussian space.  Let $H_j$ be the probabilists' Hermite polynomial~\cite{wiki-hermite}, and let $h_j = \frac{1}{\sqrt{j!}}H_j$ be the normalized Hermite polynomials. The normalized Hermite polynomial forms a complete orthonormal basis in the function space  $L^2(\mathbb R, e^{-x^2/2})$ in the following sense\footnote{We denote by $L^2(\mathbb R, e^{-x^2/2})$ the weighted $L_2$ space, namely, $L^2(\mathbb R, e^{-x^2/2})\triangleq \left\{f : \int_{-\infty}^{\infty} f(x)^2e^{-x^2/2}dx < \infty \right\}$ }. For two functions $f,g$ that map $\mathbb R$ to $\mathbb R$, define the inner product  $\inner{f,g}$ with respect to the Gaussian measure as
\begin{align}
\inner{f,g} = \mathop{\mathbb E}\displaylimits_{x\sim \mathbb N(0,1)}\left[f(x)g(x)\right]\,.\nonumber
\end{align}
The polynomials $h_0,\dots, h_m,\dots$ are orthogonal to each other under this inner product:
\begin{align}
\inner{h_i,h_j} = \delta_{ij} \,.\nonumber
\end{align}
Here $\delta_{ij} = 1$ if $i = j$ and otherwise $\delta_{ij} = 0$. Given a function $\sigma\in L^2(\mathbb R, e^{-x^2/2})$ , let the $k$-th Hermite coefficient of $\sigma$ be defined as
\begin{align}
\hat{\sigma}_k = \inner{\sigma, h_k}\,.\nonumber\end{align}

Since $h_0,\dots, h_m,\dots,$ forms a complete orthonormal basis, we have the expansion that
\begin{align}
\sigma(x) = \sum_{k\in \mathbb{N}}\hat{\sigma}_k h_k(x)\,.\nonumber
\end{align}
We will leverage several other nice properties of the Hermite polynomials in our proofs. The following claim connects the Hermite polynomial to the coefficients of Taylor expansion of a certain exponential function. It can also serve as a  definition of Hermite polynomials.
\begin{claim}[{\cite[Equation 11.8]{o2014analysis}}]
	We have that for $t,z\in \mathbb R$,
	\begin{align}
	\exp(tz-\frac{1}{2}t^2) = \sum_{k=0}^{\infty} \frac{1}{k!}H_k(z) t^k \,. \nonumber
	\end{align}
	\label{claim:hermite-gf}
\end{claim}
The following Claims shows that the expectation $\mathop{\mathbb E}\displaylimits\left[h_n(x)h_m(y)\right]$ can be computed easily when $x,y$ are (correlated) Gaussian random variables.
\begin{claim}[{\cite[Section 11.2]{o2014analysis}}]~\label{prop:rhocorrelated}
Let $(x,y)$ be $\rho$-correlated standard normal variables (that is, both $x$,$y$ have marginal distribution $\mathbb N(0,1)$ and $\mathop{\mathbb E}\displaylimits[xy] = \rho$). Then,
	\begin{align}
	\mathop{\mathbb E}\displaylimits\left[h_m(x)h_n(y)\right] = \rho^n \delta_{mn}\,. \nonumber
	\end{align}
\end{claim}
\noindent As a direct corollary, we can compute $	\mathop{\mathbb E}\displaylimits_{x\sim \mathcal{N}(0,I_{d\times d} )}\left[\sigma(u^{\top}x)\gamma(v^{\top}x)\right]$ by expanding in the Hermite basis and applying the Claim above.
\begin{claim}\label{claim:1}
Let $\sigma,\gamma$ be two functions from $\mathbb R$ to $\mathbb R$ such that  $\sigma^2 , \gamma^2 \in L^2(\mathbb R, e^{-x^2/2})$.  Then, for any unit vectors $u,v\in \mathbb R^d$, we have that
	\begin{align}
	\mathop{\mathbb E}\displaylimits_{x\sim \mathcal{N}(0,I_{d\times d} )}\left[\sigma(u^{\top}x)\gamma(v^{\top}x)\right] & = \sum_{i\in \mathbb{N}} \hat{\sigma}_i\hat{\gamma}_i \inner{u,v}^i\,. \nonumber
	\end{align}
\end{claim}
\begin{proof}[Proof of Claim~\ref{claim:1}]

	\Authornotecolored{blue}{T}{TODO: make the Fourier expansion a bit more rigorous}
	Let $s = u^{\top}x$ and $t= v^{\top}x$. Then $s,t$ are two spherical standard normal random variables that are $\inner{u,v}$-correlated, and we have that
	\begin{align}
	\mathop{\mathbb E}\displaylimits_{x\sim \mathcal{N}(0,I_{d\times d} )}\left[\sigma(u^{\top}x)\gamma(v^{\top}x)\right] & = \mathop{\mathbb E}\displaylimits\left[\sigma(s)\gamma(t)\right] \,.\nonumber
	\end{align}
	We expand $\sigma(s)$ and $\gamma(t)$ in the Fourier basis and obtain that
	\begin{align}
	\mathop{\mathbb E}\displaylimits\left[\sigma(s)\gamma(t)\right] & = \mathop{\mathbb E}\displaylimits\left[\sum_{i \in \mathbb{N}}\hat{\sigma}_{i} h_i(s)\sum_{j\in \mathbb{N}} \hat{\gamma}_j h_j(t)\right]\nonumber \\
	& = \sum_{i,j}\hat{\sigma}_i\hat{\gamma}_j \mathop{\mathbb E}\displaylimits\left[h_i(s)h_j(t)\right] \nonumber\\
	& = \sum_{i} \hat{\sigma}_i\hat{\gamma}_i \inner{u,v}^i \tag{by Claim~\ref{prop:rhocorrelated}}
	\end{align}
\end{proof}

\subsection{Analytic Formula for population risk $f$ and $f'$}\label{sec:landscape:f}

In this section we prove Theorem~\ref{thm:population_risk} and Theorem~\ref{thm:fprime}, which both follow from the following more general Theorem.

\begin{theorem}\label{thm:population_general}

	Let $\gamma, \sigma \in L^2(\mathbb R, e^{-x^2/2})$, and $\hat{y} = a^{\top}\gamma(B x)$ with parameter $a\in \mathbb R^\ell$ and $B\in \mathbb R^{\ell\times d}$. Define the population risk $f_{\gamma}$ as
	\begin{align}
	f_{\gamma}(a,B) = \mathop{\mathbb E}\displaylimits\left[\norm{y-\hat{y}}^2\right]\,.\nonumber
	\end{align}
	Suppose $B = \begin{bmatrix}
	b_1^\top\\
	\vdots\\
	b_\ell^\top
	\end{bmatrix}$ and $B^{\star} = \begin{bmatrix}
	{b^{\star}_1}^\top\\
	\vdots\\
	{b^{\star}_m}^\top
	\end{bmatrix}$  and $b_i$'s and $b^{\star}_i$'s have unit $\ell_2$ norm. Then,

	\begin{align}
	f(a,B)		&=  \sum_{k\in\mathbb{N}}  \Norm{\hat{\sigma}_k\sum_{i\in [m]} a^{\star}_i {b^{\star}_i}^{\otimes k} - \hat{\gamma}_k\sum_{i\in [\ell]} a_i b_i^{\otimes k} }_F^2 + \textup{const}, \nonumber
	\end{align}
	where $\hat{\sigma}_k, \hat{\gamma}_k$ are the $k$-th Hermite coefficients of the function $\sigma$ and $\gamma$ respectively.
\end{theorem}

We can see that Theorem~\ref{thm:population_risk} follows from choosing $\gamma=\sigma$ and Theorem~\ref{thm:fprime} follows from choosing $\gamma = \hat{\sigma}_2h_2 + \hat{\sigma}_4 h_4$. The key intuition here is that we can decompose $\sigma$ into a weighted combination of Hermite polynomials, and each Hermite polynomial influence the population risk more or less independently (because they are orthogonal polynomials with respect to the Gaussian measure).
\begin{proof}[Proof of Theorem~\ref{thm:population_general}] We have
	\begin{align}
	f_{\gamma}&  = \mathop{\mathbb E}\displaylimits\left[\Norm{\hat{y}-y}^2\right]  = \mathop{\mathbb E}\displaylimits\left[\Norm{{a^{\star}}^{\top}\sigma(B^{\star} x)-a^{\top}\gamma(B x)}^2\right]\nonumber\\
	& = \mathop{\mathbb E}\displaylimits\left[\Norm{\sum_{i\in [m]} a^{\star}_i\sigma({b^{\star}_i}^{\top}x)- \sum_{i\in [\ell]}a_i\gamma(b_i^{\top}x)}^2\right]\nonumber\\
	& = \sum_{i\in [m],j\in [\ell]}\mathop{\mathbb E}\displaylimits\left[a^{\star}_ia^{\star}_j\sigma({b^{\star}_i}^{\top}x)\sigma({b^{\star}_j}^{\top}x)\right] + \sum_{i\in [m],j\in [\ell]}\mathop{\mathbb E}\displaylimits\left[a_ia_j\gamma(b_i^{\top}x)\gamma(b_j^{\top}x)\right] \nonumber\\
	& - 2\sum_{i\in [m],j\in [\ell]}\mathop{\mathbb E}\displaylimits\left[a^{\star}_ia_j\sigma({b^{\star}_i}^{\top}x)\gamma(b_j^{\top}x)\right] \nonumber	\\
	& = \sum_{i,j\in [m]} a^{\star}_ia^{\star}_j \sum_{k\in\mathbb{N}} \hat{\sigma}_k^2\inner{b^{\star}_i,b^{\star}_j}^k+\sum_{i,j\in [\ell]} a_ia_j \sum_{k\in\mathbb{N}} \hat{\gamma}_k^2\inner{b_i,b_j}^k \nonumber \\
	& - 2\sum_{i\in [m], j\in [\ell]} a^{\star}_ia_j \sum_{k\in\mathbb{N}} \hat{\sigma}_k\hat{\gamma}_k\inner{b^{\star}_i,b_j}^k \tag{by Claim~\ref{claim:1}}\\
	&=  \sum_{k\in\mathbb{N}} \Norm{\hat{\sigma}_k\sum_{i\in [m]} a^{\star}_i {b^{\star}_i}^{\otimes k} - \hat{\gamma}_k\sum_{i\in [\ell]} a_i b_i^{\otimes k} }_F^2 \,.\nonumber
	\end{align}
\end{proof}
\subsection{Analytic Formula for population risk $G$}

In this section we show that the population risk $G(\cdot)$ (defined as in equation~\eqref{eqn:GB}) has the following analytical formula:
\begin{align}
G(B) &= 2\sqrt{6}|\hat{\sigma_4}|\cdot \sum_{i\in [d]} a^{\star}_i \sum_{j,k\in [d], j\neq k}\inner{b^{\star}_i, b_j}^2 \inner{b^{\star}_i,b_k}^2 \nonumber\\
& -\frac{|\hat{\sigma}_4|\mu }{\sqrt{6}}\sum_{i,j\in [d]} a^{\star}_i\inner{b^{\star}_i,b_j}^4+ \lambda \sum_{i=1}^m(\norm{b_i}^2-1)^2 \,.\nonumber
\end{align}
The formula will be crucial for the analysis of the landscape of $G(\cdot)$ in Section~\ref{sec:landscape}. The formula follows straightforwardly from the following two theorems and the definition~\eqref{eqn:GB}.
\begin{theorem}\label{thm:formula}
	Let $\phi(\cdot,\cdot,\cdot)$ be defined as in equation~\eqref{eqn:def-phi}, we have that
	\begin{align}
	\mathop{\mathbb E}\displaylimits\left[y \cdot \sum_{j,k\in [d], j\neq k}\phi(b_j,b_k, x)\right]  & = 2\sqrt{6}\hat{\sigma}\cdot \sum_{i\in [d]} a^{\star}_i \sum_{j,k\in [d], j\neq k}\inner{b^{\star}_i, b_j}^2 \inner{b^{\star}_i,b_k}^2 \,.\nonumber
	\end{align}
\end{theorem}

\begin{theorem}\label{thm:4thpowerformula}

	Let $\varphi(\cdot,\cdot)$ be defined as in equation~\eqref{eqn:def-varphi}, then we have that
	\begin{align}
	\mathop{\mathbb E}\displaylimits\left[y\cdot \sum_{j\in [d]}\varphi(b_j, x)\right]  = \frac{\hat{\sigma}_4}{\sqrt{6}}\sum_{i,j\in [d]} a^{\star}_i\inner{b^{\star}_i,b_j}^4 \,.\nonumber
	\end{align}
\end{theorem}

\noindent In the rest of the section we prove Theorem~\ref{thm:formula} and ~\ref{thm:4thpowerformula}.

We start with a simple but fundamental lemma. Essentially all the result in this section follows from expanding the two sides of equation~\eqref{eqn:crucial} below.
\begin{lemma}\label{lem:expabst}
	Let $u, v\in \mathbb R^d$ be two fixed vectors and $x\sim \mathbb N(0,I_{d\times d})$. Then, for any $s,t\in \mathbb R$,
	\begin{align}
		\exp(\inner{u,v}st) =  \mathop{\mathbb E}\displaylimits\left[\exp(u^\top x t - \frac 1 2 \norm{u}^2 t^2)\exp(v^\top x s - \frac 1 2 \norm{v}^2 s^2)\right] \,. \label{eqn:crucial}
	\end{align}
\end{lemma}
\begin{proof}
Using the fact that $\mathop{\mathbb E}\displaylimits\left[\exp(v^\top x)\right] = \exp(\frac{1}{2}\norm{v}^2)$, we have that,
\begin{align}
&\mathop{\mathbb E}\displaylimits\left[\exp(u^\top x t - \frac 1 2 \norm{u}^2 t^2)\exp(v^\top x s - \frac 1 2 \norm{v}^2 s^2)\right] \nonumber\\
& = \mathop{\mathbb E}\displaylimits\left[\exp((tu+sv)^\top x)\right]\exp(-\frac{1}{2}\norm{u}^2t^2 -\frac{1}{2}\norm{v}^2s^2) \nonumber\\
& = \exp(\frac{1}{2}\norm{tu+sv}^2-\frac{1}{2}\norm{u}t^2-\frac{1}{2}\norm{v}^2s^2) \tag{by the formula $\mathop{\mathbb E}\displaylimits\left[\exp(v^\top x)\right] = \exp(\frac{1}{2}\norm{v}^2)$} \\
& = \exp(\inner{u,v}st)\,.\nonumber
\end{align}
\end{proof}
Next we extend some of the results in the previous section to the setting with different scaling (such as when $v$ in Claim~\ref{claim:1} is no longer a unit vector.)
\begin{lemma}\label{lem:ab4delta4k}
	Let $u$ be a fixed unit vector and $v$ be an arbitrary vector in $\mathbb R^d$.  Let $\varphi(v,x) = \frac{1}{8} \norm{v}^4 - \frac{1}{4}(v^{\top}x)^2 \norm{v}^2 +\frac{1}{24} (v^{\top}x)^4$.
\begin{align}
\inner{u,v}^4\delta_{4,k} = \mathop{\mathbb E}\displaylimits\left[H_k(u^{\top}x)\varphi(v, x)\right]
\end{align}
\end{lemma}
As a sanity check, we can verify that when $v$ is a unit vector, $\varphi(v,x) = \sqrt{24}h_4(v^\top x)$ and th Lemma reduces to a special case of Claim~\ref{prop:rhocorrelated}.
\begin{proof}
	Let $A, B$ be formal power series in variable $s,t$ defined as $A = \exp(\inner{u,v}st)$ and $B=\mathop{\mathbb E}\displaylimits\left[\exp(u^\top x t - \frac 1 2 \norm{u}^2 t^2)\exp(v^\top x s - \frac 1 2 \norm{v}^2 s^2)\right] $. We refer the readers to~\cite{wiki:power_series} for more backgrounds of power series. For casual readers, one can just think of $A$ as $B$ as two power series obtained by expanding the $exp(\cdot)$ via Taylor expansion. For a formal power series $A$ in variable $x$, let $[x^\alpha]A$ to denote coefficient in front of the monomial $x^\alpha$. By Lemma~\ref{lem:expabst}, we have that $A=B$, and thus
	\begin{align}
	\co{s^4t^k} A = \co{s^4t^k}B\,,
	\end{align}
	which implies that
	\begin{align}
		\frac{1}{24}\inner{u,v}^4\delta_{4,k} &= \mathop{\mathbb E}\displaylimits\left[\co{t^k}\left(\exp(u^\top x t- \frac{1}{2}t^2)\right) \cdot \co{s^4}\left(\exp(v^\top x s - \frac 1 2 \norm{v}^2 s^2)\right)\right]\nonumber\\
		&= \mathop{\mathbb E}\displaylimits\left[\frac{1}{k!}H_k(u^\top x)\varphi(v,x)\right]\,.
			\end{align}
	where the last line is by the fact that $\varphi(v,x) = \co{s^4}\exp(v^\top x s - \frac 1 2 \norm{v}^2 s^2)$. This can be verified by applying Claim \ref{claim:hermite-gf} with $t=s \norm{v}$ and $z=\frac{v^T x}{\norm{v}}$, and noting that $H_4 (x) = x^4-6x^2+3$.
\end{proof}

\noindent Now we are ready prove Theorem~\ref{thm:4thpowerformula} using Lemma~\ref{lem:ab4delta4k}.

\begin{proof}[Proof of Theorem~\ref{thm:4thpowerformula}]
	Using the fact that $\sigma(v^\top x) = \sum_{k=0}^{\infty} \hat{\sigma}_kh_k(v^\top x)$, we have that
	\begin{align}
	\mathop{\mathbb E}\displaylimits\left[y\cdot \sum_{j}\varphi(b_j, x)\right] & = \sum_{i,j\in [d]}a^{\star}_i \mathop{\mathbb E}\displaylimits\left[\sigma({b^{\star}_i}^\top x) \varphi(b_j,  x)\right] \nonumber\\
	& = \sum_{i,j\in [d]} a^{\star}_i \sum_{k}^{\infty} \mathop{\mathbb E}\displaylimits\left[\hat{\sigma}_k h_k({b^{\star}_i}^\top x) \varphi(b_j, x)\right] \nonumber \\
	& = \sum_{i,j\in [d]}a^{\star}_i \mathop{\mathbb E}\displaylimits\left[\hat{\sigma}_4 h_4({b^{\star}_i}^\top x) \varphi(b_j,x)\right] = \frac{\hat{\sigma}_4}{\sqrt{6}}\sum_{i,j\in [d]} a^{\star}_i  \inner{b^{\star}_i,b_j}^4 \tag{by Lemma~\ref{lem:ab4delta4k} and $h_j = \frac{1}{\sqrt{j}}H_j$}
	\end{align}
\end{proof}
\noindent Towards proving Theorem~\ref{thm:formula}, we start with the following Lemma. Inspired by the proofs above, we design a function $\phi(v,w,x)$ such that we can estimate $\inner{u,v}^2\inner{u,w}^2$ by taking expectation of $\mathop{\mathbb E}\displaylimits\left[\sigma(u^\top x)\phi(v,w,x)\right]$.
\begin{lemma}\label{lem:ab2ac2}
Let $a$ be a fixed unit vector in $\mathbb R^d$ and $v,w$ two fixed vectors in $\mathbb R^d$. Let $\varphi(\cdot,\cdot)$ be defined as in Lemma~\ref{lem:ab4delta4k}.  	Define $\phi(v,w,x)$ as
\begin{align}
\phi(v,w,x) & = \varphi(v+w,x)+\varphi(v-w,x) -2\varphi(v,x) - 2\varphi(w,x)\\
& = \frac{1}{2} \norm{v}^2\norm{w}^2 + \inner{v,w}^2 - \frac{1}{2}\norm{w}^2 (v^\top x)^2 - \frac{1}{2}\norm{v}^2(w^{\top}x)^2  \\
& - 2(v^\top x)(w^\top x) v^\top w + \frac{1}{2}(v^\top x)^2(w^\top x)^2\,. \nonumber\end{align}
Then, we have that
\begin{align}
 \mathop{\mathbb E}\displaylimits\left[\sigma(u^\top x) \phi(v,w,x)\right] = 2\sqrt{6}\hat{\sigma}_4 \inner{u,v}^2\inner{u,w}^2\nonumber\,.
\end{align}
	\end{lemma}
\begin{proof}
Using the fact that $\inner{u,v+w}^2 + \inner{u,v-w}^4 - 2\inner{u,v}^2 - 2\inner{u,w}^4 = 12\inner{u,v}^2\inner{u,w}^2$ and Lemma~\ref{lem:ab4delta4k}, we have that
\begin{align}
12\inner{u,v}^2\inner{u,w}^2\delta_{4,k} &= \mathop{\mathbb E}\displaylimits\left[H_k(u^\top x) (\varphi(v+w,x)+\varphi(v-w,x) -2\varphi(v,x) - 2\varphi(w,x))\right] \nonumber\\
& = \mathop{\mathbb E}\displaylimits\left[H_k(u^\top x) \phi(v,w,x)\right] \label{eqn:1}\,.
\end{align}

\noindent Using the fact that $\sigma(u^\top x) = \sum_{k=0}^{\infty} \hat{\sigma}_kh_k(u^\top x)$, we conclude that
\begin{align}
 \mathop{\mathbb E}\displaylimits\left[\sigma(u^\top x) \phi(v,w,x)\right] & =  \sum_{k=0}^{\infty} \hat{\sigma}_k\mathop{\mathbb E}\displaylimits\left[h_k(u^\top x)\phi(v,w,x)\right]  \nonumber\\
 & =  \frac{\hat{\sigma}_4}{\sqrt{6}}\mathop{\mathbb E}\displaylimits\left[H_4(u^\top x)\phi(v,w,x)\right] \tag{by Lemma~\ref{lem:ab4delta4k} and $h_j = \frac{1}{\sqrt{j}}H_j$} \\
 &=  2\sqrt{6}\hat{\sigma}_4 \inner{u,v}^2\inner{u,w}^2\tag{by Lemma~\ref{lem:ab4delta4k} again}
\end{align}
\end{proof}
\noindent Now we are ready to prove Theorem~\ref{thm:formula} by using Lemma \ref{lem:ab2ac2} for every summand.
\begin{proof}[Proof of Theorem~\ref{thm:formula}]
	We have that
	\begin{align}
	\mathop{\mathbb E}\displaylimits\left[y \cdot \sum_{j,k}\phi(b_j,b_k, x)\right] & = \sum_{i}a^{\star}_i \sum_{j,k}\mathop{\mathbb E}\displaylimits\left[\sigma({b^{\star}_i}^{\top}x)\phi(b_j,b_k,x)\right]\nonumber\\
	& =2 \sqrt{6}\hat{\sigma}_4\sum_{i}a^{\star}_i  \sum_{j,k}\inner{b^{\star}_i, b_j}^2 \inner{b^{\star}_i,b_k}^2	\,. \tag{by Lemma~\ref{lem:ab2ac2}}\end{align}
\end{proof}

\section{Recovering the Linear Layer}\label{sec:linear}

We will show in this section that if we have are given a $\delta$-approximation of $B^{\star}$, then it is easy to recover $a^{\star}$. The key observation here is that the correlation between the $\inner{b^{\star}_i,x}$ and the output $y$ is exactly proportional to $a^{\star}_i$. We also note that there could be multiple other ways to recover $a^{\star}$, e.g., using linear regression with the $\sigma(Bx)$ as input and the $y$ as output. We chose this algorithm mostly because of the ease of analysis.
\begin{algorithm}\caption{Recovering $a^{\star}$}\label{alg:recovera}
	{\bf Input: } A matrix $B$ with unit row norms that is row-wise $\delta$-close to $B^{\star}$ in Euclidean distance. \\
	{\bf Return: } Let $a'_i = 2\widehat{\E}[y\inner{x,b_i}]$ where $\widehat{\E}$ means the empirical average.  Set $a_i \leftarrow |a'_i|$ and $b_i \leftarrow b_i\mbox{sgn}(a'_i)$
\end{algorithm}

\begin{lemma*}[Restatement of Lemma~\ref{lem:recoveraintro}]
		Given a matrix $B$ whose rows are $\delta$-close to $B^\star$ in Euclidean distance up to permutation and sign flip with $\delta \le 1/(2\kappa^\star)$. Then, we can give estimates $a,B'$ (using e.g., Algorithm~\ref{alg:recovera}) such that there exists a permutation $P$ where  $\|a-P a^\star\|_\infty \le \delta a^{\star}_{\max}$ and $B'$ is row-wise $\delta$-close to $PB^{\star}$.
\end{lemma*}

\noindent To see why this simple algorithm works for recovering $a^{\star}$, we need the following simple claim.

\begin{claim}\label{clm:correlation}
For any vector $v$ we have
$$
\E[y\inner{x,v}] = \frac{1}{2} \sum_{i=1}^m a^{\star}_i\inner{b^{\star}_i, v}.
$$
\end{claim}

The proof of this claim follows immediately from the property of Hermite polynomials. Now we are ready to prove the corollary.

\begin{proof}
Without loss of generality we assume $B$ is close to a sign flip of $B^{\star}$. The unknown permutation does not change the proof.

Since $b_i$ is $\delta$ close to $B^{\star}_i$, let $u$ be the vector where $u_j = \inner{b^{\star}_j, b_i-b^{\star}_i}$, we have
$$a'_i = \sum_{i=1}^m a^{\star}_i\inner{b^{\star}_i, b_i} = \sum_{i=1}^m a^{\star}_i(\inner{b^{\star}_i, b^{\star}_i}+\inner{b^{\star}_i,b_i-b^{\star}_i}) = a^{\star}_i + \inner{a^{\star}_i,u} \in a^{\star}_i \pm a^{\star}_{max}\delta.$$
Therefore $a'_i$ is always positive, $a_i$ is in the desirable range and $\|B'_i - B^{\star}_i\|\le \delta$.

Similarly, if $-b_i$ is $\delta$ close to $B^{\star}_i$, we have $a'_i\in -a^{\star}_i\pm a^{\star}_{max}\delta$, and the conclusion still holds.

\end{proof}

For the settings considered in Section~\ref{sec:general}, the vectors $b^{\star}_i$ are not necessarily orthogonal. In this case we use the following algorithm:

\begin{algorithm}\caption{Recovering $a^{\star}$ for general case}\label{alg:recovergeneral}
	{\bf Input: } A matrix $B$ with unit row norms, and $B$ is $\delta$-close to $B^{\star}$ in spectral norm up to permutation and sign flip. \\
 Let $u_i = 2\widehat{\E}[y\inner{x,b_i}]$ where $\widehat{\E}$ means the empirical average.
 \\
 Let $a' = (BB^\top)^{-1} u$.\\
 {\bf Return:}
  Set $a_i \leftarrow |a'_i|$ and $b_i \leftarrow b_i\mbox{sgn}(a'_i)$
\end{algorithm}

\begin{lemma}\label{lem:recovergeneral}
		Given a matrix $B$ whose rows have unit norm, and $\|B-SPB^{\star}\|\le \delta$ for some permutation matrix $P$ and diagonal matrix $S$ with $\pm 1$ entries on diagonals.If $\frac{\sigma_{min}^2(B)}{4\sqrt{2}\kappa^\star\sqrt{m}}$, we can give estimates $a,B'$ (using e.g., Algorithm~\ref{alg:recovergeneral}) such that $\|a-Pa^\star\| \le \frac{2\sqrt{2}a^{\star}_{max}\sqrt{m}}{\sigma_{min}^{-2}(B)}\cdot \delta$ and $\|B' - PB^{\star}\|\le \delta$.
\end{lemma}

\begin{proof}
We again use Claim~\ref{clm:correlation}: in this case we know the vector $u$ satisfies $u = B(B^{\star})^\top a^{\star}$. As a result, for the vector $a'$, we have
$$
a' = (BB^\top)^{-1} (B(B^{\star})^\top)a^{\star} = (B^\dag)^\top (B^{\star})^\top a^{\star} = (B^{\star} B^\dag)^\top a^{\star}.
$$

By assumption we know $B = SPB^{\star} + E$ where $\|E\| \le \delta$. By the perturbation of matrix inverse (Theorem~\ref{thm:matpert}), we know if $\|E\|\le \delta \le\sigma_{min}(B)/2$, then $B^\dag = (B^{\star})^\dag P^{-1}S^{-1} + E'$ where $\|E'\| \le 2\sqrt{2}\sigma_{min}(B)^{-2} \delta$. Therefore
$$
a' = (P^{-1}S^{-1} +E')^\top a^{\star} = S^{-\top}P^{-\top} a^{\star} + (E')^\top a^{\star} = SPa +(E')^\top a^{\star}.
$$
(Here the last equality is because for both permutation matrix $P$ and sign flip matrix $S$, $P^{-\top} = P$ and $S^{-\top} = S$.) Therefore, coordinates of $a'$ are permutation and sign flips of $a^{\star}$, up to an error term $(E')^\top a^{\star}$.

When $\delta \le \frac{\sigma_{min}^2(B)}{4\sqrt{2}\kappa^\star\sqrt{m}}$, we know $\|(E')^\top a^{\star}\| \le \|E'\| a^{\star}_{max}\sqrt{m} \le a^{\star}_{min}/2$, therefore the signs are all recovered correctly. After fixing the sign, we have $\|a - Pa\| \le \|(E')^\top a^{\star}\| \le \frac{2\sqrt{2}\deltaa^{\star}_{max}\sqrt{m}}{\sigma_{min}^{-2}(B)}$, and $\|B' - PB^{\star}\| \le \delta$.
\end{proof}

\section{Sample Complexity}
\label{sec:sample}

In this section we will show that our algorithm only requires polynomially many samples to find the desired solution. Note that we did not try to optimize the polynomial dependency.

\begin{theorem}[Theorem~\ref{thm:samplecomplexity} Restated]
In the setting of Theorem~\ref{thm:main}, suppose we use $N$ empirical samples to approximate $G$ and obtain function $\widehat{G}$. There exists a fixed polynomial such that if $N \ge \mbox{poly}(d, a^{\star}_{max}/a^{\star}_{min},1/\epsilon)$, with high probability for any point $B$ with

	 $\lambda_{min}(\nabla^2 \widehat{G}(B)) \ge -\tau_0/2$
	 and $\|\nabla \widehat{G}(B)\| \le \epsilon/2$, then  $B$ can be written as $B = DP +E $ where $P$ is a permutation matrix, $D$ is a diagonal matrix and $|E|_\infty \le O(\epsilon/(\hat{\sigma}_4a^\star_{\min}))$.
\end{theorem}
\newcommand{\widehat{G}}{\widehat{G}}
In order to bound the sample complexity, we will prove a uniform convergence result: we show that with polynomially many samples, the gradient and Hessian of $\widehat{G}$ are point-wise close to the gradient and Hessian of $G$, therefore any approximate local minimum of $\widehat{G}$ must also be an approximate local minimum of $G$.

However, there are two technical issues in showing the uniform convergence result. The first issue is that when the norm of $B$ is very large, both the gradient and Hessian of $G$ and $\widehat{G}$ are very large and we cannot hope for good concentration. We deal with this issue by showing when $B$ has a large norm, the empirical gradient  $\nabla \widehat{G}(B)$ must also have large norm, and therefore it can never be an approximate local minimum (we do this later in Lemma~\ref{lem:normbound}).
The second issue is that our objective function involves high-degree polynomials over Gaussian variables $x,y$, and is therefore not sub-Gaussian or sub-exponential. We use a standard truncation argument to show that the function does not change by too much if we restrict to the event that the Gaussian variables have bounded norm.

\newcommand{\mathcal{F}}{\mathcal{F}}
\begin{lemma}\label{lem:truncate}
Suppose $P'(B)+R(B) = \E_{(x,y)}[f(x,y,B)]$ where $f$ is a polynomial of degree at most $5$ in $x, y$ and at most $4$ in $B$. Also assume that the sum of absolute values of coefficients is bounded by $\Gamma$. For any $\epsilon \le \Gamma/2$, let $R = Cd\log (a^{\star}_{max}\Gamma/\epsilon)$ for a large enough constant $C$, let $\mathcal{F}$ be the event that $\|x\|^2 \le R$, and let $G_{trunc} = \E_{(x,y)}[f(x,y,B)1_\mathcal{F}]$. For any $B$ such that $\|b_i\|\le 2$ for all rows, we have
$$
\|\nabla G(B) - \nabla G_{trunc}(B)\| \le \epsilon,
$$
and
$$
\|\nabla^2 G(B) - \nabla^2 G_{trunc}(B)\| \le \epsilon,
$$
\end{lemma}

\begin{proof}
By standard $\chi^2$ concentration bounds, for large enough $C$ and any $z > R$, the probability that $\|x\|^2 \ge z$ is at most $\exp(-10z)$.

By simple calculation, it is easy to check that $\|\nabla_B f(x,y,B)\| \le 4\Gamma d^{1.5} a^{\star}_{max}\|x\|^5$, and $\|\nabla^2_B f(x,y,B)\| \le 12 \Gamma d^2 a^{\star}_{max}\|x\|^5$. We know $\|\nabla G(B) - \nabla G_{trunc}(B)\|  = \|\E[\nabla_B f(x,y,B)(1-1_\mathcal{F})\|$. The expectation between $\|x\|^2 \in [2^i R, 2^{i+1}R]$, for $i = 0,1,2,...$, is always bounded by $4\Gamma d^{1.5} a^{\star}_{max}\|2^{i+1}R\|^5 \exp(-2^iR) < \epsilon/2^{i+1}$. Therefore
$$
\|\nabla G(B) - \nabla G_{trunc}(B)\| \le \sum_{i=0}^\infty \epsilon/2^{i+1} \le \epsilon.
$$
The bound for the Hessian follows from the same argument.
\end{proof}

Finally, we combine this truncation with a result of \cite{mei2016landscape} that proves universal convergence of gradient and Hessian. For completeness here we state a version of their theorem with bounded gradient/Hessian:

\begin{theorem}[Theorem 1 in \cite{mei2016landscape}]\label{thm:concentration}
Let $f(\theta)$ be a function from $\mathbb R^p\to \mathbb R$ and $\hat{f}$ be its empirical version. If the norm of the gradient and Hessian of a function is always bounded by $\tau$, for variables in a ball of radius $r$ in $p$ dimensions, there exists a universal constant $C_0$ such that for $C = C_0\max\{\log r\tau/\delta, 1\}$, the following hold:
\begin{enumerate}
\item[(a)] The sample gradient converges to the population gradient. Namely if $N \ge Cp\log p$ we have
$$
\Pr[\sup_{\|\theta\|\le r}\|\nabla f(\theta) - \nabla \hat{f}_\theta\| \le \tau \sqrt{\frac{Cp\log n}{n}}]\ge 1-\delta.
$$
\item[(b)] The sample Hessian converges to the empirical Hessian. Namely if $N \ge Cp\log p$ we have
$$
\Pr[\sup_{\|\theta\|\le r}\|\nabla f(\theta) - \nabla \hat{f}_\theta\| \le \tau \sqrt{\frac{Cp\log n}{n}}]\ge 1-\delta.
$$
\end{enumerate}
\end{theorem}

As an immediate corollary of this theorem and Lemma~\ref{lem:truncate}, we have

\begin{corollary}\label{cor:uniformconverge}
In the setting of Theorem~\ref{thm:samplecomplexity}, for every $B$ whose rows have norm at most 2, we have with high probability,
$$
\|\nabla G(B) - \nabla \widehat{G}(B)\| \le \epsilon/2,
$$
and
$$
\|\nabla^2 G(B) - \nabla^2 \widehat{G}(B)\| \le \tau_0/2.
$$
\end{corollary}

\begin{proof}
On the other hand, for all such matrices $B$, by Lemma~\ref{lem:truncate} we know the gradient and Hessian of $G$ is close to the gradient and Hessian of $G_{trunc}$.
$$
\|\nabla G(B) - \nabla G_{trunc}(B)\| \le \epsilon/4,
$$
and
$$
\|\nabla^2 G(B) - \nabla^2 G_{trunc}(B)\| \le \tau_0/4.
$$

Now, the gradient and Hessian for individual samples for estimating $G_{trunc}$ are bounded by some $\mbox{poly}(d,1/\epsilon)$, therefore by Theorem~\ref{thm:concentration} we know the gradient and Hessian of $\widehat{G}$ are close to those of $G_{trunc}$. When $N \ge \mbox{poly}(d,1/\epsilon)$ for a large enough polynomial, we have with high probability, for all $B$ with all rows $\|b_i\|\le 2$,

$$
\|\nabla G_{trunc}(B) - \nabla \widehat{G}(B)\| \le \epsilon/4,
$$
and
$$
\|\nabla^2 G_{trunc}(B) - \nabla^2 \widehat{G}(B)\| \le \tau_0/4.
$$
The corollary then follows from triangle inequality.
\end{proof}

Finally we handle the case when $B$ has a row with large norm. We will show that in this case $\nabla \widehat{G}(B)$ must also be large, so $B$ cannot be an approximate local minimum.

\begin{lemma}\label{lem:normbound}
If $b_i$ is the row with largest norm and $\|b_i\| \ge 2$, then when $N\ge \mbox{poly}(d, a^{\star}_{max}/a^{\star}_{min})$ for some fixed polynomial, we have with high probability $\inner{\nabla \widehat{G}(B), b_i} \ge c\lambda\|b_i\|^4$ for some universal constant $c > 0$.
\end{lemma}

\begin{proof}
The proof of this Lemma is very similar to Claim~\ref{clm:ball}. Note that by equation~\eqref{eqn:GB} there are three terms in $\widehat{G}(B)$:\sloppy
(1) $\sign(\hat \sigma_4) \hat{\mathop{\mathbb E}\displaylimits}\left[y \cdot \sum_{j,k\in [d], j\neq k}\phi(b_j,b_k, x)\right]$,
(2) $- \mu \sign(\hat \sigma_4) \hat{\mathop{\mathbb E}\displaylimits}\left[y\cdot\sum_{j \in [d]} \varphi(b_j,x)\right]$, (3) $\lambda \sum_{i=1}^m(\norm{b_i}^2-1)^2$. Here $\hat{E}$ is the empirical average over the samples.

Note that the first two terms are homogeneous degree 4 polynomials over $B$, and the third term does not depend on the sample. By argument similar to Corollary~\ref{cor:uniformconverge}, we know for any $B$ where $b_i$ has the largest row norm, with the number of samples we choose the gradient of the first two terms is $ca^{\star}_{min}\|b_i\|^3$ close to the gradient of their expectations, where $c < 0.01$ is a small constant.

By Theorem~\ref{thm:main:formula}, we know the expectation of the first two terms are equal to $A1(B) = 2\sqrt{6}|\hat{\sigma}_4|\cdot \sum_{i\in [d]} a^{\star}_i \sum_{j,k\in [d], j\neq k}\inner{b^{\star}_i, b_j}^2 \inner{b^{\star}_i,b_k}^2$ and $A2(B)=  - \frac{|\hat{\sigma}_4|\mu }{\sqrt{6}}\sum_{i,j\in [d]} a^{\star}_i\inner{b^{\star}_i,b_j}^4$. Here the gradient of the first term always have positive correlation with $b_i$, so we can ignore it. For the second term, we know the gradient
$$
\frac{\partial}{\partial b_i}[A2(B)] = - \frac{|\hat{\sigma}_4|\mu}{\sqrt{6}} \sum_{j} a^{\star}_j\inner{b^{\star}_j,b_i}^3b^{\star}_j.
$$
Taking the inner-product with $b_i$, and use the fact that $b^{\star}_i$ form an orthonormal basis, we know
$$
\inner{\frac{\partial}{\partial b_i}[A2(B)],b_i} \ge - ca^{\star}_{min}\|b_i\|^4.
$$

On the other hand, when $\|b_i\|\ge 2$, we have for the third term
$$
\inner{\frac{\partial}{\partial b_i}[\lambda (\norm{b_i}^2-1)^2],b_i} \ge \lambda (\|b_i\|-1)^4 \ge \lambda\|b_i\|^4/16.
$$
Since $\lambda$ is larger than $a^{\star}_{max}$, we know the negative contribution from $A2$ and the difference between the empirical version and $G$ are both negligible. Therefore we have $\inner{\nabla \widehat{G}(B),b_i} \ge c\lambda \|b_i\|^4$ as desired.
\end{proof}

Now we are ready to prove Theorem~\ref{thm:samplecomplexity}:

\begin{proof}
By Lemma~\ref{lem:normbound}, any point $B$ with $\nabla \widehat{G}(B) \le \epsilon$ must have $\|b_i\|\le 2$ for all $i$.
Now by Corollary~\ref{cor:uniformconverge}, we know the point $B$ we have must satisfy
$$
\|\nabla G(B)\| \le \epsilon; \nabla^2 G(B) \succeq -\tau_0I.
$$
By point 3 in Theorem~\ref{thm:main}, this implies the guarantee on $B$.
\end{proof}
\section{Spurious Local minimum for function $P'$}
\label{sec:example}

In this section we give an example where the function $P'$ does have spurious local minimum.

In this example, $d = 4$, and the true vectors are the standard basis vectors $b^{\star}_i = e_i$. We will set $a^{\star}_1 = 1$, and $a^{\star}_2 = a^{\star}_3 = a^{\star}_4 = 2+\delta$ (where $\delta>0$ is an arbitrary positive constant).

The spurious local minimum that we consider is $b_1 = b_2 = e_1 = b^{\star}_1$, $b_3 = e_2 = b^{\star}_2$, $b_4 = \frac{\sqrt{2}}{2} e_3+\frac{\sqrt{2}}{2} e_4$. That is,
$$
B = \left(\begin{array}{cccc}1 & 0 & 0 & 0\\
1 & 0 & 0 & 0\\
0 & 1 & 0 & 0\\
0 & 0 & \frac{\sqrt{2}}{2} & \frac{\sqrt{2}}{2}\end{array}\right).
$$

The objective $P'(B) = 1$ and the only non-zero term is $a^{\star}_1 \inner{b^{\star}_1,b_1}^2 \inner{b^{\star}_1,b_2}^2$. In order to improve the objective locally, we need to change either $b_1$ or $b_2$, otherwise the term $a^{\star}_1 \inner{b^{\star}_1,b_1}^2 \inner{b^{\star}_1,b_2}^2$ is still 1, and all other terms ($a^{\star}_i\inner{b^{\star}_i,b_j}^2\inner{b^{\star}_i,b_k}^2$) are non-negative.

Assume we have a local perturbation $B'$, where $b'_1 = \sqrt{1-\epsilon_1^2} e_1 + \epsilon_1 u_1$, $b'_2 = \sqrt{1-\epsilon_2^2} e_1 + \epsilon_2 u_2$. Here $u_1, u_2$ are unit vectors that are orthogonal to $e_1$. Also, since this is a local perturbation, we make sure $\epsilon_1,\epsilon_2 \le \epsilon$, and $b_3(2) \ge 1-\epsilon$, $[b_4(3)]^2, [b_4(4)]^2 \ge 0.5 - \epsilon$. We will show that when $\epsilon$ is small enough, the objective function $P'(B') \ge 1$.

To see this, notice that the term $a^{\star}_1 \inner{b^{\star}_1,b_1}^2 \inner{b^{\star}_1,b_2}^2$ is now equal to $(1-\epsilon_1^2)(1-\epsilon_2^2)$. On the other hand, for $b_1$, we have
\begin{align*}
\sum_{i=2}^4 a^{\star}_i \sum_{k=3}^4 \inner{b^{\star}_i,b_1}^2 \inner{b^{\star}_i,b_k}^2 & = \epsilon_1^2 (2+\delta)\sum_{i=2}^4 \sum_{k=3}^4 \inner{b^{\star}_i,u_1}^2 \inner{b^{\star}_i,b_k}^2 \\
& = \epsilon_1^2 (2+\delta)\sum_{i=2}^4  \inner{b^{\star}_i,u_1}^2 (\sum_{k=3}^4\inner{b^{\star}_i,b_k}^2) \\
& \ge \epsilon_1^2 (2+\delta)\sum_{i=2}^4  \inner{b^{\star}_i,u_1}^2 \cdot \min_{i=2}^4\{\sum_{k=3}^4\inner{b^{\star}_i,b_k}^2\}\\
& \ge \epsilon_1^2 (2+\delta)(0.5-\epsilon).
\end{align*}

Similarly we have the same equation for $b_2$. Note that all the terms we analyzed are disjoint, therefore
$$
P'(B') \ge (1-\epsilon_1^2)(1-\epsilon_2^2) + \epsilon_1^2 (2+\delta)(0.5-\epsilon) + \epsilon_2^2 (2+\delta)(0.5-\epsilon).
$$
By removing higher order terms of $\epsilon$, it is easy to see that $P'(B') \ge 1$ when $\epsilon$ is small enough. Therefore $B$ is a local minima of $P'$.
\section{More backgrounds on Hermite Polynomials} \label{sec:backgrounds}

\section{Introduction}
Scalable optimization has been playing crucial roles in the success of deep learning, which has immense applications in artificial intelligence.  Remarkably, optimization issues are often addressed through designing new models that make the resulting training objective functions easier to be optimized. For example, over-parameterization~\cite{livni2014computational}, batch-normalization~\cite{ioffe2015batch}, and residual networks~\cite{he2016deep,he2016identity} are often considered as ways to improve the optimization landscape of the resulting objective functions.

How do we design models and objective functions that allow efficient optimization with guarantees?  Towards understanding this question in a principled way, this paper studies learning neural networks with one hidden layer. Roughly speaking, we will show that when the input is from Gaussian distribution and under certain simplifying assumptions on the weights, we can design an objective function $G(\cdot)$, such that\\
\vspace{-0.1in}

\noindent [a] all local minima of $G(\cdot)$ are global minima \\

\vspace{-0.1in}

\noindent [b] all the global minima are the desired solutions, namely, the ground-truth parameters (up to permutation and some fixed transformation).
\vspace{0.1in}

We note that designing such objective functions is challenging because 1) the natural $\ell_2$ loss objective does have bad local minimum, and 2) due to the permutation invariance\footnote{Permuting the rows of $B^{\star}$ and the coordinates of $a^{\star}$ correspondingly preserves the functionality of the network.}, the objective function inherently has to contain an exponential number of isolated local minima.

\subsection{Setup and known issues with proper learning}

We aim to learn a neural network with a one-hidden-layer using a non-convex objective function. We assume input $x$ comes from Gaussian distribution and the label $y$ comes from the model
\begin{align}
y =  {a^{\star}}^{\top}\sigma({B^{\star}}x) + \xi \label{eqn:model}
\end{align}
where $a^{\star} \in \mathbb R^{d}, B^{\star}\sim \mathbb R^{m \times d}$ are the ground-truth parameters, $\sigma(\cdot)$ is a element-wise non-linear function, and $\xi$ is a noise vector with zero mean. Here we can without loss of generality assume $x$ comes from spherical Gaussian distribution $\mathbb N(0,I_{d\times d})$. \footnote{This is because if $x\sim N(0,\Sigma)$, then we can whiten the data by taking $x' = \Sigma^{-1/2}x$ and define ${B^{\star}}' = B\Sigma^{1/2}$. We note that ${B^{\star}}'x' = Bx$ and therefore we main the functionality of the model. }

For technical reasons, we will further assume $m \le d$ and that $a^{\star}$ has non-negative entries.

The most natural learning objective is perhaps the $\ell_2$ loss function, given the additive noise.
Concretely, we can parameterize with training parameters $a\in \mathbb R^d, B\sim \mathbb R^{m \times d}$ of the same dimension as $a^{\star}$ and $B^{\star}$ correspondingly,
\begin{align}
\hat{y} = a^{\top}\sigma(Bx)\,,
\end{align}
and then use stochastic gradient descent to optimize the $\ell_2$ loss function. When we have enough training examples, we are effectively minimizing  the following population risk with stochastic updates,
\begin{align}
f(a,B) = \mathop{\mathbb E}\displaylimits\left[\norm{\hat{y} - y}^2\right]\,.
\end{align}

However,  {\bf empirically} stochastic gradient descent {\bf cannot} converge to the ground-truth parameters in the synthetic setting above when $\sigma(x) = \textup{ReLU}(x) = \max\{x,0\}$,  even if we have access to an infinite number of samples, and $B^{\star}$ is a orthogonal matrix. Such empirical results have been reported in~\cite{livni2014computational} previously, and we also provide our version in Figure~\ref{fig:relu} of Section~\ref{sec:experiments}.  This is consistent with observations and theory that over-parameterization is crucial for training neural networks successfully~\cite{livni2014computational, hardt2016gradient,soudry2016no}.

These empirical findings suggest that the population risk $f(a,B)$ has spurious local minima with inferior error compared to that of the global minimum. This phenomenon occurs even if we assume we know $a^{\star}$ or $a^{\star} = \mathbf{1}$ is merely just the all one's vector.  Empirically, such landscape issues seem to be alleviated by over-parameterization.
By contrast, our method described in the next section does not require over-parameterization and might be suitable for applications that demand the recovery of the true parameters.

\subsection{Our contributions}

Towards learning with the same number of training parameters as the ground-truth model, we first study the landscape of the population risk $f(\cdot)$ and give an analytic formula for it  --- as an explicit function of the ground-truth parameter and training parameter with the randomness of the data being marginalized out.  The formula in equation~\eqref{eqn:popluation_risk_tensor} shows that $f(\cdot)$ is implicitly attempting to solve simultaneously a finite number of low-rank tensor decomposition problems with commonly shared components.

Inspired by the formula, we design a new training model whose associated loss function --- named $f'$ and formally defined in equation~\eqref{eqn:fprime} --- corresponds to the loss function for decomposing a matrix (2-nd order tensor) and a 4-th order tensor (Theorem~\ref{thm:fprime}). Empirically,
stochastic gradient descent on $f'$
learns the network as shown in experiment section (Section~\ref{sec:experiments}).

Despite the empirical success of $f'$, we still lack a provable guarantee on the landscape of $f'$.
The second contribution of the paper is to design a more sophisticated objective function $G(\cdot)$ whose landscape is provably nice --- all the local minima of $G(\cdot)$ are proven to be global, and they correspond to the permutation of the true parameters. See Theorem~\ref{thm:main}.

Moreover, the value and the gradient of $G$ can be estimated using samples, and there are no constraints in the optimization. These allow us to use straightforward stochastic gradient descent (see guarantees in \cite{ge2015escaping,jin2017escape}) to optimize $G(\cdot)$ and converge to a local minimum, which is also a global minimum (Corollary~\ref{cor:algorithm}).

Finally, we also prove a finite-sample complexity result. We will show that with a polynomial number of samples, the empirical version of $G$ share almost the same landscape properties as $G$ itself (Theorem~\ref{thm:samplecomplexity}). Therefore, we can also use an empirical version of $G$ as a surrogate in the optimization.

\subsection{Related work}

The work of Arora et al.~\cite{arora2014provable} is one of the early results on provable algorithms for learning deep neural networks, where the authors give an algorithm for learning deep generative models with sparse weights. Livni et al.~\cite{livni2014computational}, Zhang et al.~\cite{zhang2016l1,zhang2017learnability}, and Daniely et al.~\cite{daniely2016toward} study the learnability of special cases of neural networks using ideas from kernel methods. Janzamin et al.~\cite{janzamin2015beating} give a polynomial-time algorithm for learning one-hidden-layer neural networks with twice-differential activation function and known input distributions, using the ideas from tensor decompositions.

A series of recent papers study the theoretical properties of non-convex optimization algorithms for one-hidden-layer neural networks.  Brutzkus and Globerson~\cite{brutzkus2017globally} and Tian~\cite{tian2017analytical} analyze the landscape of the population risk for one-hidden-layer neural networks with Gaussian inputs under the assumption that the weights vector associated to each hidden variable (that is, the filters) have disjoint supports. Li and Yuan~\cite{li2017convergence} prove that stochastic gradient descent recovers the ground-truth parameters when the parameters are known to be close to the identity matrix. Zhang et al.~\cite{DBLP:journals/corr/ZhangPS17} studies the optimization landscape of learning one-hidden-layer neural networks with a specific activation function, and they design a specific objective function that can recover a single column of the weight matrix. Zhong et al.~\cite{zhong2017recovery} studies the convergence of non-convex optimization from a good initializer that is produced by tensor methods. Our algorithm works for a large family of activation functions (including ReLU) and any full-rank weight matrix.
To our best knowledge, we give the first global convergence result for gradient-based methods for our general setting
\footnote{The work of ~\cite{janzamin2015beating,zhong2017recovery} are closely related, but they require tensor decomposition as the algorithm/initialization.

The optimization landscape properties have also been investigated on simplified neural networks models. Kawaguchi~\cite{kawaguchi2016deep} shows that the landscape of deep neural nets does not have bad local minima but has degenerate saddle points. Hardt and Ma~\cite{hardt17identity} show that re-parametrization using identity connection as in residual networks~\cite{he2016deep} can remove the degenerate saddle points in the optimization landscape of deep linear residual networks. Soudry and Carmon~\cite{soudry2016no} showed that an over-parameterized neural network does not have bad differentiable local minimum. Hardt et al.~\cite{hardt2016gradient} analyze the power of over-parameterization in a linear recurrent network (which is equivalent to a linear dynamical system.)

The optimization landscape has also been analyzed for other machine learning problems, including SVD/PCA phase retrieval/synchronization, orthogonal tensor decomposition, dictionary learning, matrix completion, matrix sensing \cite{baldi1989neural,srebro2003weighted, ge2015escaping,sun2015nonconvex, bandeira2016low,ge2016matrix,bhojanapalli2016global,ge2017no}. Our analysis techniques build upon that for tensor decomposition in~\cite{ge2015escaping} --- we add two additional regularization terms to deal with spurious local minimum caused by the weights $a^{\star}$ and to remove the constraints.

\subsection{Notations: }

We use $\mathbb{N},\mathbb{R}$ to denote the set of natural numbers and real numbers respectively. We use $\norm{\cdot}$ to denote the Euclidean norm of a vector and spectral norm of a matrix. We use $\norm{\cdot}_F$ to denote the Frobenius/Euclidean norm of a matrix or high-order tensor. For a vector $x$, let $\norm{x}_0$ denotes its infinity norm and for a matrix $A$, let $|A|_0$ be a shorthand for $\norm{\textup{vec}(A)}_0$ where $\textup{vec}(A)$ is the vectorization of $A$. For a vector $x$, let $\sl{x}$ denotes the second largest absolute values of the entries for $x$. We note that $\sl{\cdot}$ is not a norm.

We use $A\otimes B$ to denote the Kronecker product of $A$ and $B$, and $A^{\otimes k}$ is a shorthand for $A\otimes \cdots \otimes A$ where $A$ appears $k$ times. For vectors $a\otimes b$ and $a^{\otimes k}$ denote the tensor product. We use $\lambda_{\max}(\cdot), \lambda_{\min}(\cdot)$ to denote the largest and smallest eigenvalues of a square matrix. Similarly, $\sigma_{\max}(\cdot)$ and $\sigma_{\min}(\cdot)$ are used to denote the largest and smallest singular values. We denote the identity matrix in dimension $d\times d$ by $I_{d\times d}$, or $I$ when the dimension is clear from the context.

In the analysis, we rely on many properties of Hermite polynomials. We use $h_j$ to denote the $j$-th normalized Hermite polynomial. These polynomials form an orthonormal basis. See Section~\ref{sec:basic-hermite} for an introduction of Hermite polynomials.

We will define other notations when we first use them.
\section{Simulation}\label{sec:experiments}
\begin{figure}[!t]
	\centering
	\begin{subfigure}[b]{0.4\textwidth}
		\includegraphics[width=\textwidth]{plots/sgd_relu_d=50_its=50k_obj.eps}
	\end{subfigure}
	\begin{subfigure}[b]{0.4\textwidth}
		\includegraphics[width=\textwidth]{plots/sgd_relu_d=50_its=50k_paramErr.eps}
	\end{subfigure}
		\caption{Data are generated by a network with ReLU activation without noise. The training model uses the same architecture. Left: the estimated population risk doesn't converge to zero. Right: the parameter error using the surrogate in equation~\eqref{eqn:surrogate}.}\label{fig:relu}
\end{figure}
In this section, we provide simple simulation results that verify that minimizing $G(B)$ with SGD recovers a permutation of $B^{\star}$; however, minimizing Equation \eqref{eqn:population_risk} with SGD results in finding spurious local minima. Based on the formula for the population risk in Equation~\eqref{eqn:popluation_risk_tensor}, we also verified empirically the conjecture that SGD would successfully recover $B^{\star}$ using the activation functions $\gamma(z)= \hat \sigma_2 h_2 (z) + \hat \sigma_4 h_4 (z)$,\footnote{We also observed that using $\gamma(z) =\frac12 |z|$ also works but due to the space limitation we don't report the experimental results here. }  even if the data were generated via a model with ReLU activation. (See Section~\ref{sec:main-formula} for the rationale behind such conjectures.)

For all of our experiments, we chose $B^{\star} =I_{d \times d}$ with dimension $d=50$ and $a^{\star} = \mathbf{1}$ for simplicity, and the data is generated from a one-hidden-layer network with ReLU activation without noise. We use stochastic gradient descent with fresh samples at each iteration, and we plot the (expected) population error (that is, the error on a fresh batch of examples).

\begin{figure}
	\centering
	\begin{subfigure}[b]{0.4\textwidth}
		\includegraphics[width=\textwidth]{plots/sgd_h2_h4_obj_shifted.eps}

	\end{subfigure}
	~

	\begin{subfigure}[b]{0.4\textwidth}
		\includegraphics[width=\textwidth]{plots/sgd_h2_h4_paramErr.eps}

	\end{subfigure}
	\caption{The labels are generated from a network with ReLU activation. We learn with $\hat \sigma_2 h_2 + \hat\sigma_4 h_4$ activation. Left: the test loss subtracted by the theoretical global minimum value. Right: the error in parameter space measured by equation~\eqref{eqn:surrogate}}\label{fig:h2h4}

\end{figure}

\begin{figure}[!h]
	\centering
	\begin{subfigure}[b]{0.4\textwidth}
		\includegraphics[width=\textwidth]{plots/our_obj_objectiveValue.eps}
	\end{subfigure}
	~

	\begin{subfigure}[b]{0.4\textwidth}
		\includegraphics[width=\textwidth]{plots/our_obj_paramErr.eps}

	\end{subfigure}
	\caption{Learning with objective function $G(\cdot)$. Left: the test loss. Right: the error in parameter space measured by equation~\eqref{eqn:surrogate}. }\label{fig:our}

\end{figure}

To test whether SGD converges to a matrix $B$ which is equivalent to $B^{\star}$ up to permutation of rows, we use a  surrogate error metric to evaluate whether ${B^{\star}}^{-1}B$ is close to a permutation matrix. Given a matrix $Q$ with row norm 1, let
\begin{align}e(Q) = \min\{1-\min_i \max_j |Q_{ij}|, 1-\min_j \max_i |Q_{ij}|\}.\label{eqn:surrogate}\end{align}
Then we have that if $e(Q)\le \epsilon$ for some $\epsilon < 1/3$, then it implies that $Q$ is $\sqrt{2\epsilon}$-close to a permutation matrix in infinity norm. On the other direction, we know that if $e(Q) > \epsilon$, then $Q$ is not $\epsilon$-close to any permutation matrix in infinity norm. The latter statement also holds when $Q$ doesn't have row norm $1$.

Figure~\ref{fig:relu} shows that without over-parameterization, using ReLU as an activation function, SGD doesn't converge to zero test error and the ground-truth parameters. We decreased step-size by a factor of $4$ every $5000$ number of iterations after the error plateaus at $10000$ iterations. For the final $5000$ iterations, the step-size is less than $10^{-9}$, so we can be confident that the non-zero objective value is not due to the variance of SGD.
We see that none of the five runs of SGD converged to a global minimum.

Figure~\ref{fig:h2h4} shows that using $\hat{\sigma}_2h_2 + \hat{\sigma}_4h_4$ as the activation function, SGD with projection to the set of matrices $B$ with row norm 1 converges to the ground-truth parameters. We also plot the loss function which  converges the value of a global minimum. (We subtracted the constant term in equation~\eqref{eqn:fh2h4} so that the global minimum has loss 0.)

Figure~\ref{fig:our} shows that using our objective function $G(B)$, the iterate converges to the ground truth matrix $B^{\star}$. The fact that the parameter error goes up and down is not surprising,  because the algorithm first gets close to a saddle point and then breaks ties and converges to a one of the global minima.

Finally we note that using the loss function $G(\cdot)$ seems to require significantly larger batch (and sample complexity) to reduce the variance in the gradients estimation. We used batch size 262144 in the experiment for $G(\cdot)$. However, in contrast, for the $\hat{\sigma}_2h_2+\hat{\sigma_4}h_4$ we used batch size 8192 and for relu we used batch size 256.

\section{Conclusion}

In this paper we first give an analytic formula for the population risk of the standard $\ell_2$ loss, which empirically may converge to a spurious local minimum. We then design a novel population loss that is guaranteed to have no spurious local minimum.

Designing  objective functions with well-behaved landscape is an intriguing and fruitful direction.
We hope that our techniques can be useful for characterizing and designing the optimization landscape for other settings.

We conjecture that the objective $\alpha f_2 +\beta f_4$ has no spurious local minimum when $\alpha,\beta$ are reasonable constants and the ground-truth parameters are in general position\footnote{See equation~\eqref{eqn:f4} for the definition of $f_k$ and Theorem~\ref{thm:fprime} for how to access $\alpha f_2 + \beta f_4$ in the setting of one-hidden-layer neural nets.}.  We provided empirical evidence to support the conjecture.

Our results assume that the input distribution is Gaussian. Extending them to other input distributions is a very interesting open problem.

\section{Handling Non-Orthogonal Weights}

\label{sec:general}

\newcommand{F}{F}

In this section, we first show that when the weight vectors $\{b^{\star}_i\}'s$ are not orthonormal, the local optimum of a slight variant of $G(B)$ still allow us to recover $B^{\star}$. The main observation is that the set of local minima are preserved (in a certain sense) by linear transformation of the variables. We design an objective function $F(B)$ that is equivalent to $G(B)$ up to a linear transformation. This allows us to use Theorem~\ref{thm:main} as a black box to characterize all the local minima of $F$.

\subsection{Local Minimum after a Linear Transformation}

Given a function $f(y)$, we say function $g(\cdot)$ is a linear transformation of $f(\cdot)$ if there is a matrix $W$ such that $g(x) = f(Wx)$. If $W$ has full rank, the local minima of $f$ are closely related to the local minima of $g$.

We recall some standard notation in calculus first. We use $\nabla f(t)$ to denote the gradient of $f$ evaluated at $t$. For example, $\nabla f(Wx)$ is a shorthand for $\frac{\partial f(y)}{\partial y} \vert_{y=Wx}$, and similarly $\nabla^2 f(Wx)$ is  $\frac{\partial^2 f(y)}{(\partial y)^2} \vert_{y=Wx}$.

The following theorem then connects the gradients and Hessians of $f(Wx)$ and $g(x)$. Essentially, it shows that the set of local minima and saddle points have a 1-1 mapping between $f$ and $g$, and the corresponding norms/eigenvalues only differ multiplicatively by quantities related to the spectrum of $W$.

\begin{theorem}\label{thm:lineartransform} Let $W\in \mathbb R^{d\times m} (d\ge m)$ be a full rank matrix. Suppose $g:\mathbb R^m\rightarrow \mathbb R$ and $f:\mathbb R^d\rightarrow \mathbb R$ are twice-differentiable functions such that $g(x) = f(Wx)$ for any $x\in \mathbb R^m$. Then, for all $x\in \mathbb R^m$, the following three properties hold:
\begin{enumerate}
\item $\sigma_{min}(W) \|\nabla f(Wx)\| \le \|\nabla g(x)\| \le \sigma_{\max}(W)  \|\nabla f(Wx)\|$.
\item If $\lambda_{min}(\nabla^2 g(x)) < 0$, then
$$\sigma_{\max}(W) ^2 \lambda_{min}(\nabla^2 f(Wx))\le \lambda_{min}(\nabla^2 g(x))\le \sigma_{min}(W)^2 \lambda_{min}(\nabla^2 f(Wx)).$$
\item The point $x$ satisfies the first and second order optimality condition for $g$ iff $y=Wx$ also satisfy the first and second order optimality condition for $f$.
\end{enumerate}
\end{theorem}

\begin{proof}
The proof follows from the relationship between the gradients of $g$ and the gradients of $f$. By basic calculus, we have
$$
\nabla g(x) = \frac{\partial f(Wx)}{\partial x} = W^\top \frac{\partial f(y)}{\partial y}\Big \vert_{y = Wx} = W^\top \nabla f(Wx)
$$
which immediately implies bullet 1. Similarly, we can compute the second order derivative:
$$
\nabla^2 g(x) = W^\top [\nabla^2 f(Wx)] W.
$$
\sloppy To simplify notation, let $A = \nabla^2 f(Wx)$. Let $x=\arg\min_{\norm{x}=1} x^\top W^\top AWx$, and $y = (Wx)/\|W x\|$.  Therefore
 $$\lambda_{min}(A) \le y^\top A y \le \lambda_{min}(W^\top AW)/\norm{Wx}^2\le \lambda_{\min} (W^\top A W) / \norm{W}^2 .$$

On the other hand, let $y$ be the unit vector that minimizes $y^\top A y$, we know $y$ is in column span of $W$ because $f$ is only defined on the row span, so there must exist a unit vector $x$ such that $W x = \lambda y$ where $\lambda \ge \sigma_{min}(W)$. For this $x$ we have $\lambda_{min}(W^\top AW) \le x^\top W^\top AW x = \lambda^2 \lambda_{min}(A) \le \sigma_{min}^2(W)\lambda_{min}(A)$. This finishes the proof for 2.

Finally, notice that $W$ is full rank, so $\nabla g(x) = W^\top \nabla f(Wx) = 0$ iff $\nabla f(Wx) = 0$. Also, $\nabla^2 g(x) = W^\top [\nabla^2 f(Wx)] W \succeq 0$ iff $\nabla^2 f(Wx) \succeq 0$.
\end{proof}

\subsection{Objective for Non-Orthogonal Weights}

Now we will design a new objective function that can be linearly transformed to the orthonormal case. The main idea is to view the rows of $B^{\star}$ as the new basis that we work on (which is not necessarily orthogonal). Note that this is already the case for the first two terms of the objective function $G(B)$, we change the objective function as follows:
More concretely, we define

\begin{align}
F_{\alpha,\mu,\lambda}(B) &= 2\sqrt{6}\hat{\sigma}\cdot \sum_{i\in [d]} \alpha_i \sum_{j,k\in [d]}\inner{b^{\star}_i, b_j}^2 \inner{b^{\star}_i,b_k}^2 \nonumber\\
& -\frac{\hat{\sigma}_4\mu }{\sqrt{6}}\sum_{i,j\in [d]} \alpha_i\inner{b^{\star}_i,b_j}^4+ \lambda \sum_{j=1}^m((\sum_{i=1}^m \alpha_i\inner{b_j,b^{\star}_i}^2 - 1)^2-1)^2 \,.\nonumber
\end{align}

Note that the only change in the objective is the regularizer for the norm of $b_j$. It is now replaced by $((\sum_{i=1}^m \alpha_i\inner{b_j,b^{\star}_i}^2 - 1)^2-1)^2$, which tries to ensure the ``norm'' of $b_j$ in the basis defined by row of $B^{\star}$ to be 1. The objective function that we will optimize corresponds to choosing $\alpha_i = a^{\star}_i$.

Similar as before, this function can be computed as expectations

\begin{align}
F_{a^{\star},\mu,\lambda}(B) &= \mathop{\mathbb E}\displaylimits\left[y \cdot \sum_{j,k\in [d], j\neq k}\phi(b_j,b_k, x)\right] - \mu \mathop{\mathbb E}\displaylimits\left[y\cdot\sum_{j \in [d]} \varphi(b_j,x)\right] \nonumber \\
&  + \lambda \E_{(x,y),(x',y')}[\sum_{i=1}^m y\cdot \phi_2(b_i,x)\cdot y'\cdot \phi_2(b_i,x')], \label{eqn:ob}
\end{align}
where $(x',y')$ is an independent sample, and $\phi_2(v,x) = (v^\top x)^2 - \|v\|^2$.

Intuitively, if we can find a linear transformation that makes $\{b^{\star}_i\}$'s orthonormal, that will reduce the problem to the orthonormal case. This is in fact the whitening matrix:

\newcommand{o}{o}
\newcommand{\mathcal{B}}{\mathcal{B}}
\newcommand{\bar{b}}{\bar{b}}

Let $M = \sum_{i=1}^m a^{\star}_i b^{\star}_i (b^{\star}_i)^\top$  be the weighted covariance matrix of $b^{\star}_i$'s. Suppose the SVD of $M$ is $UDU^\top$ and let $W = UD^{-1/2}$. We apply the transformation $W^\top$ to the vectors $\sqrt{a^{\star}_i}b_i$'s and obtain that $o_i =  W^\top \sqrt{a^{\star}_i}b^{\star}_i$. We can verify that $o_i$'s are orthogonal vectors because

\begin{align}
\sum_{i\in [m]}o_io_i^\top = W^\top M W = I
\end{align}

For notational convenience, let's extend the definition of the $G(\mathcal{B})$ in equation by using the putting the relevant information in the subscript
	\begin{align}
G_{\alpha,\beta,\lambda, o}(\mathcal{B}) &= \sqrt{6}\hat{\sigma}\cdot \sum_{i\in [d]} a^{\star}_i \sum_{j,k\in [d], j\neq k}\inner{o_i, \bar{b}_j}^2 \inner{o_i,\bar{b}_k}^2  - \frac{\hat{\sigma}_4\mu }{\sqrt{6}}\sum_{i,j\in [d]} a^{\star}_i\inner{o_i,\bar{b}_j}^4 \,.\nonumber\nonumber\\
&+ \lambda \sum_{i=1}^m(\norm{\bar{b}_i}^2-1)^2 \nonumber
\end{align}
(That is, the index $o$ denotes the ground-truth solution with respect to which $G$ is defined.)

The next Theorem shows that we can rotate the objective function $F$ properly so that it matches the objective $G$ with a ground-truth vector $o_i$'s.

\begin{theorem}\label{lem:fglinear}
 Let $W$ be defined as above, and let $1/a^{\star}$ be the vector whose $i$-th entry is $1/a^{\star}_i$. Then, we have that
$$
G_{1/a^{\star},\mu,\lambda,o_i}(\mathcal{B}) = F_{a^{\star},\mu,\lambda}( \mathcal{B} W^\top) ).
$$
Note this can be interpreted as a linear transformation as in vector format $\mathcal{B} W^\top$ is equal to $\mathcal{B}\cdot (W^\top\otimes I_{d\times d})$.
\end{theorem}

\begin{proof}
The equality can be obtained by straightforward calculation. We note that since $\mathcal{B} = \begin{bmatrix}
\bar{b}_1^\top\\
\vdots\\
\bar{b}_m^\top
\end{bmatrix}$, the rows of $ \mathcal{B}\cdot (W^\top\otimes I_{d\times d})$ are $W\bar{b}_1,\dots, W\bar{b}_m$.

Therefore, we have that
\begin{align}
&\quad F_{a^{\star},\mu,\lambda}(\mathcal{B}\cdot (W^\top\otimes I_{d\times d}) )\\
 &= 2\sqrt{6}\hat{\sigma}\cdot \sum_{i\in [d]} a^{\star}_i \sum_{j\neq k\in [d]}\inner{b^{\star}_i, W\bar{b}_j}^2 \inner{b^{\star}_i,W\bar{b}_k}^2 \nonumber\\
&\quad -\frac{\hat{\sigma}_4\mu }{\sqrt{6}}\sum_{i,j\in [d]} a^{\star}_i\inner{b^{\star}_i,W\bar{b}_j}^4+ \lambda \sum_{j=1}^m(\sum_{i=1}^m a^{\star}_i\inner{W\bar{b}_j,b^{\star}_i}^2 - 1)^2 \,.\nonumber \\
& = 2\sqrt{6}\hat{\sigma}\cdot \sum_{i\in [d]} \frac{1}{a^{\star}_i} \sum_{j,k\in [d]}\inner{\sqrt{a^{\star}_i}W^\top b^{\star}_i, \bar{b}_j}^2 \inner{\sqrt{a^{\star}_i}W^\topb^{\star}_i,\bar{b}_k}^2 \nonumber\\
&\quad -\frac{\hat{\sigma}_4\mu }{\sqrt{6}}\sum_{i,j\in [d]} \frac{1}{a^{\star}_i}\inner{\sqrt{a^{\star}_i}W^\topb^{\star}_i,\bar{b}_j}^4+ \lambda \sum_{j=1}^m(\sum_{i=1}^m \inner{\bar{b}_j,\sqrt{a^{\star}_i}W^\topb^{\star}_i}^2 - 1)^2 \,.\nonumber \\
& = 2\sqrt{6}\hat{\sigma}\cdot \sum_{i\in [d]} \frac{1}{a^{\star}_i} \sum_{j,k\in [d]}\inner{o_i, \bar{b}_j}^2 \inner{o_i,\bar{b}_k}^2 \nonumber\\
&\quad -\frac{\hat{\sigma}_4\mu }{\sqrt{6}}\sum_{i,j\in [d]} \frac{1}{a^{\star}_i}\inner{o_i,\bar{b}_j}^4+ \lambda \sum_{j=1}^m(\sum_{i=1}^m \inner{\bar{b}_j,o_i}^2 - 1)^2 \,.\tag{by the definition of $o_i$'s}
\end{align}
\end{proof}

From Theorem~\ref{thm:main} we can immediately get the following Corollary (note that the only difference is that the coefficients now are $1/a^{\star}_i$ instead of $a^{\star}_i$). Recall $a^{\star}_{max} = \max_i a^{\star}_i$ and
$a^{\star}_{min} = \min_i a^{\star}_{min}$, we have

\begin{corollary}\label{cor:G}
Let $\kappa_a = a^{\star}_{max}/a^{\star}_{min}$.
Let $c$ be a sufficiently small universal constant (e.g. $c=0.01$ suffices). Assume $\mu \le c/\kappa_a$ and $\lambda \ge (ca^{\star}_{min})^{-1}$. The function $G_{1/a^{\star},\mu,\lambda,o_i}(\cdot)$ defined as in Theorem~\ref{lem:fglinear} satisfies that
	\begin{enumerate}
		\item[1.] A matrix $\mathcal{B}$ is a local minimum of $G$ if and only if $\mathcal{B}$ can be written as $\mathcal{B} = PDO$ where $O$ is a matrix whose rows are $o_i$'s, $P$ is a permutation matrix and $D$ is a diagonal matrix with  $D_{ii} \in \left\{\pm 1 \pm O(\mu/\lambdaa^{\star}_{min})\right\}$.
		 \item[2.] Any saddle point $\mathcal{B}$ has a strictly negative curvature in the sense that $\lambda_{\min}(\nabla^2 G(\mathcal{B}))\ge -\tau_0$ where $\tau_0 = c\min\{\mu /(\kappa_a a^{\star}_{max}d), \lambda\}$
		\item [3.]  Suppose $\mathcal{B}$ is an approximate local minimum in the sense that $\mathcal{B}$ satisfies $$\norm{\nabla g(\mathcal{B})}\le \epsilon \textup{  and  }\lambda_{\min}(\nabla^2 g(\mathcal{B})) \ge - \tau_0$$
		Then $\mathcal{B}$ can be written as $\mathcal{B} = PDO +E $ where $P$ is a permutation matrix, $D$ is a diagonal matrix and $|E|_\infty \le O(\epsilona^{\star}_{max}/\hat{\sigma}_4)$.
	\end{enumerate}
\end{corollary}

Finally, we can combine the theorem above and Theorem~\ref{thm:main-general} to give a guarantee for optimizing $F$.
Let $\Gamma$ be a diagonal matrix with $\Gamma_{ii} = \sqrt{a^\star_i}$. Let $M = {B^{\star}}^\top {\Gamma}^2 B^{\star}$ and $\kappa(M) = \|M\|/\sigma_{min}(M)$.

\begin{theorem}\label{thm:main-nonortho}
Let $c$ be a sufficiently small universal constant (e.g. $c=0.01$ suffices). Let $\kappa_a = a^{\star}_{max}/a^{\star}_{min}$. Assume $\mu \le c/\kappa_a$ and $\lambda \ge 1/(c\cdot a^\star_{\min})$. The function $F(\cdot)$ defined as in Theorem~\ref{lem:fglinear} satisfies that
	\begin{enumerate}
		\item[1.] A matrix $B$ is a local minimum of $F$ if and only if $B$ satisfy $B^{-\top} = P  D\GammaB^{\star}$ where $P$ is a permutation matrix, $\Gamma$ is a diagonal matrix with $\Gamma_{ii} = \sqrt{a^\star_i}$, and $D$ is a diagonal matrix with  $D_{ii} \in \left\{\pm 1 \pm O(\mu /\lambdaa^{\star}_{min})\right\}$. Furthermore, this means that all local minima of $F$ are also global.
		 \item[2.] Any saddle point $B$ has a strictly negative curvature in the sense that $\lambda_{\min}(\nabla^2 F(B))\ge -\tau_0$ where $\tau_0 = c\min\{\mu /(\kappa_a d a^\star_{\max}), \lambda\}\sigma_{min}(M)$.
		\item [3.]  Suppose $B$ is an approximate local minimum in the sense that $B$ satisfies $$\norm{\nabla F(B)}\le \epsilon \textup{  and  }\lambda_{\min}(\nabla^2 F(B)) \ge - \tau_0$$
		Then $B$ can be written as $B^{-\top} = P D\Gamma B^{\star} +E $ where $\Gamma, D, P$ are as in 1, the error term $\|E\| \le O(\epsilona^{\star}_{\max}\sqrt{md}\cdot \kappa(M)^{1/2}/\hat{\sigma}_4)$ (when $\epsilona^{\star}_{\max}\sqrt{md}\cdot \kappa(M)^{1/2}/\hat{\sigma}_4 < c$).
	\end{enumerate}
\end{theorem}

\begin{proof}
Note that we can immediately apply Theorem~\ref{thm:main} to $G_{1/a^{\star},\mu,\lambda,o_i}(B)$ to characterize all its local minima. See Corollary~\ref{cor:G}.

Next we will transform the properties for local minima of $G$ (stated in Corollary~\ref{cor:G}) to $F$ using Theorem~\ref{thm:lineartransform}. First we note that the transformation matrix $W$ and $M$ are closely related:
\begin{align}WW^\top = M, \sigma_{min}(W)^2 = 1/\|M\|, \|W\|^2 = 1/\sigma_{min}(M).\end{align}
This is because according to the definition of $W$, the SVD of $M$ is $M = UDU^\top$ and $W = UD^{-1/2}$, so $WW^\top = UD^{-1}U^\top = M^{-1}$. The claims of the singular values follow immediately from the SVD of $M$ and $W$.

As a result, all local minimum of $F$ are of the form $\mathcal{B} W^\top$ where $\mathcal{B}$ is a local minimum of $G$. For $B = \mathcal{B} W^\top$, the gradient and Hessian of $F(B)$ and $G(\mathcal{B})$ are also related by Theorem~\ref{thm:lineartransform}.

Let us first prove 1. By Corollary~\ref{cor:G}, we know every local minimum of $G$ is of the form $\mathcal{B} = PDO$. According to the definition of $O$ in Theorem~\ref{lem:fglinear}, we know each row vector $o_i$ is equal to $W^\top (a^{\star}_i)^{1/2} b^{\star}_i$, therefore $O = \Gamma B^{\star} W$. As a result, all local minima of $G$ are of the form $\mathcal{B} = P D \GammaB^{\star} W$. By Theorem~\ref{thm:lineartransform} and Theorem~\ref{lem:fglinear}, we know all local minima of $F$ must be of the form $B =\mathcal{B} W^\top =   P  D\Gamma B^{\star} WW^\top =  P  D \GammaB^{\star} M^{-1}$.

Now we try to compute $B^{-\top}$. To do that observe that  $ [\Gamma B^{\star}] M^{-1} [\Gamma B^{\star}]^\top  =  I$. Therefore $[\Gamma B^{\star} M^{-1}]^{-\top} = \Gamma B^{\star}$, and for any local minimum $B$, we have
 \begin{align}B^{-\top} & = (P  D\GammaB^{\star} M^{-1})^{-\top} = P^{-\top}D^{-\top}(\Gamma B^{\star} M^{-1})^{-\top} \nonumber
 \\
& = P^{-\top}D^{-\top}\Gamma B^{\star}.
\nonumber\end{align}
 Note that $P^\top$ is still a permutation matrix, and $D^{-\top}$ is still a matrix whose diagonal entries are $\{\pm 1\pm O(\mu/\lambdaa^{\star}_{\min})\}$, so this is exactly the form we stated in 1. More concretely, the rows of $B^{-\top}$ are permutations of $\sqrt{a^{\star}_i}b^{\star}_i$.

For bullet 2, it follows immediately from Property 2 in Theorem~\ref{thm:lineartransform}. Note that by property 2,
$$\lambda_{min}(\nabla^2 F(\mathcal{B} W^\top)) \le \frac{\lambda_{min}(\nabla^2 G(\mathcal{B}))}{\|W\|^2} = \lambda_{min}(\nabla^2 G(\mathcal{B}))\sigma_{min}(M).$$

Finally we will prove 3. Let $\mathcal{B} = BW^{-\top}$, so that $G(\mathcal{B}) = F(B)$. We will prove properties of $B$ using the properties of $\mathcal{B}$ from Corollary~\ref{cor:G}.

First we observe that by Theorem~\ref{thm:lineartransform},
$$
\lambda_{min}(\nabla^2 G(\mathcal{B})) \ge \|W\|^2 \lambda_{min}(\nabla^2 F(B)) \ge -c\min\{\mu/(\kappa_a da^{\star}_{\max},\lambda\}.$$
\sloppy
Therefore the second order condition for Claim 3 in Corollary~\ref{cor:G} is satisfied.
Now when $\|\nabla F(B)\| \le \epsilon$, we have $\|\nabla G(\mathcal{B})\| \le \epsilon \|W\| = \epsilon/\sigma_{min}(M)^{1/2}$. By Corollary~\ref{cor:G}, we know $\mathcal{B}$ can be expressed as $PDO + E'$ where $D$ is the diagonal matrix, $P$ is a permutation matrix and $|E'|_\infty \le O(\epsilona^{\star}_{\max}/(\hat{\sigma}_4\sigma_{min}(M)^{1/2}))$. We will apply perturbation Theorem~\ref{thm:matpert} for matrix inversion. Since $\sigma_{min}(PDO) \ge 1/2$, we know when $\|E'\| \le 1/4$,
$$\|(PDO+E')^{-1} - (PDO)^{-1}\| \le 8\sqrt{2}\|E'\|.$$
Here $\|E\|$ is bounded by $\|E\|_F \le \sqrt{md} |E'|_\infty \le O(\epsilona^{\star}_{\max}\sqrt{md}/(\hat{\sigma}_4\sigma_{min}(M)^{1/2}))$, which is smaller than $1/4$ when $\epsilon$ is small enough.

The corresponding point in $F$ is $B = \mathcal{B} W^\top$, and in 1 we have already proved $(PDOW^\top)^{-\top}$ is of the form we want, therefore we can define $E = B^{-\top} - (PDOW^\top)^{-\top} = (\mathcal{B} - PDO)^{-\top}W^{-1}$, and
$$
\|E\| = \|W^{-1}\|\|(PDO+E')^{-1} - (PDO)^{-1}\| = O(\epsilona^{\star}_{\max}\sqrt{md}\cdot \kappa(M)^{1/2}/\hat{\sigma}_4).
$$
This finishes the proof.
\end{proof}

\subsection{Handle Undercomplete Case}\label{sec:undercomplete}

The objective function $F$ can handle the case when the weights $b^{\star}_i$'s are not orthogonal, but still requires the number of components $m$ to be equal to the number of dimensions $d$. In this section we show how to use similar ideas for the case when the number of components is smaller than the dimension ($m < d$).

\newcommand{\mathcal S}{\mathcal{S}}
\newcommand{\mathcal{F}}{\mathcal{F}}

Note that all the terms in $F(B)$ only depends on the inner-products $\inner{b_j, b^{\star}_i}$. Let $\mathcal S$ be the span of $\{b^{\star}_i\}$'s and $P_\mathcal S$ be the projection matrix to this subspace, it is easy to see that $F(B)$ satisfies
$$
F(B) = F(BP_\mathcal S ).
$$
That is, the previous objective function only depends on the projection of $B$ in the space $\mathcal S$. Using similar argument as Theorem~\ref{thm:main-nonortho}, it is not hard to show the only local optimum in $\mathcal S$ satisfies the same conditions, and allow us to recover $B^{\star}$. However, without modifying the objective, the local optimum of $F(B)$ can have arbitrary components in the orthogonal subspace $\mathcal S^\perp$.

In order to prevent the components from $\mathcal S^\perp$, we add an additional $\ell_2$ regularizer: define $\mathcal{F}_{\alpha,\mu,\lambda,\delta}$ as follows:

\begin{align}
\mathcal{F}_{\alpha,\mu,\lambda,\delta}(B) = F_{\alpha,\mu,\lambda}(B) + \frac{\delta}{2} \|B\|_F^2 \label{eqn:def:f}
\end{align}

Intuitively, since the first term $F_{\alpha,\mu,\lambda}(B)$ only cares about the projection $ BP_\mathcal S$, minimizing $\|B\|_F^2$ will remove the components in the orthogonal subspace of $\mathcal S$. We will choose $\delta$ carefully to make sure that the additional term does not change the local optima of $F_{\alpha,\mu,\lambda}(B)$ by too much, while still ensuring a small projection on $\mathcal S^\perp$.

In this case we will consider pseudo-inverse instead of inverse. In particular, for a $m\times d$ matrix $B$, define its pseudo-inverse $B^\dag$ to be the matrix such that $BB^\dag = I_{m\times m}$ and $B^\dag B$ is the projection to the row span of $B$.

Let $M = \sum_{i=1}^m a^{\star}_i b^{\star}_i (b^{\star}_i)^\top$, $\kappa(M) = \|M\|/\sigma_{m}(M)$.
\begin{theorem}\label{thm:main-under}
For any desired accuracy $\epsilon_0$, we can choose parameters $\epsilon,\delta,\tau_0,\mu,\lambda$, such that for the objective function $\mathcal{F}_{a^{\star},\mu,\lambda,\delta}(B)$, for any $B$ such that
$$
\|\nabla \mathcal{F}(B)\| \le \epsilon, \quad \nabla^2 \mathcal{F}(B) \ge -\tau_0/2,
$$
we have $[B^{\dag}]^\top = B^{\star} D\Gamma P+E$ where $\Gamma$ is a diagonal matrix with entries $\sqrt{a^{\star}_i}$, $D$ is a diagonal matrix with entries close to 1, $P$ is a permutation matrix and $\|E\| \le \epsilon_0$.

To choose the parameters, let $c$ be a sufficiently small universal constant (e.g. $c=0.01$ suffices).  Assume $\mu \le c/\kappa^\star$ and $\lambda \ge 1/(c\cdot a^\star_{\min})$. Let $\tau_0 = c\min\{\mu /(\kappa d a^\star_{\max}), \lambda\}\sigma_{min}(M)$. Let $\delta \le \min\{\frac{c\hat{\sigma}_4\epsilon_0}{a^{\star}_{max}\cdot m\sqrt{d} \kappa^{1/2}(M)}, \tau_0/2\}$, and $\epsilon = \min\{\lambda\sigma_{min}(M)^{1/2},c\delta/\sqrt{\|M\|}, c\epsilon_0\delta\sigma_{min}(M)\}$.
\end{theorem}

We first show that if the gradient is small, then the point cannot have a large component in $\mathcal S^\perp$.

\begin{lemma}\label{lem:outside}
If $\|\nabla \mathcal{F}_{a^{\star},\mu,\lambda,\delta}(B)\| \le \epsilon$, then $\|P_{\mathcal S^\perp} B\|_F \le \epsilon/\delta$.
\end{lemma}

\begin{proof}
Since $F_{a^{\star},\mu,\lambda}(B)$ only depends $BP_\mathcal S$, we know $ \nabla F_{a^{\star},\mu,\lambda}(B)P_{\mathcal S^\perp} = 0$. Therefore $\epsilon \ge \|\nabla \mathcal{F}_{a^{\star},\mu,\lambda,\delta}(B)P_{\mathcal S^\perp}\|_F = \|(\delta B)P_{\mathcal S^\perp}\|_F = \delta \|P_{\mathcal S^\perp} B\|_F$, and we have $\| BP_{\mathcal S^\perp}\|_F \le \epsilon/\delta$ as desired.
\end{proof}

Next we show that if the gradient of $\mathcal{F}_{a^{\star},\mu,\lambda,\delta}(B)$ is small, and $\delta$ is also small, then the gradient of $F_{a^{\star},\mu,\lambda}(B)$ can be bounded.

\begin{lemma} \label{lem:inside}
In the setting of Theorem~\ref{thm:main-under}, if $\|\nabla \mathcal{F}_{a^{\star},\mu,\lambda,\delta}(B)\| \le \epsilon \le \lambda\sigma_{min}(M)^{1/2}$, then we have $$\|\nabla F_{a^{\star},\mu,\lambda}(B)\| \le \epsilon + \delta \sqrt{2m/\sigma_{min}(M)}.$$
\end{lemma}

Towards proving Lemma~\ref{lem:inside}, we first bound the norm of $B$ by the following claim:

\begin{claim}\label{clm:ball} If $\|\nabla \mathcal{F}_{a^{\star},\mu,\lambda,\delta}(B)\|\le \lambda\sigma_{min}(M)^{1/2}$, then each row $b_i$ must satisfy $b_i^\top M b_i \le 2$.
\end{claim}
\begin{proof} We prove by contradiction. Assume towards contradiction that there is a column $b_i$ such that $b_i^\top M b_i \ge 2$. We consider the quantity,
$$
\langle \frac{\partial  \mathcal{F}_{a^{\star},\mu,\lambda,\delta}}{\partial b_i}(B), b_i \rangle.
$$

Note that $\mathcal{F}_{a^{\star},\mu,\lambda,\delta}$ has 4 terms:
(1) $2\sqrt{6}\hat{\sigma}\cdot \sum_{i\in [d]} \alpha_i \sum_{j,k\in [d]}\inner{b^{\star}_i, b_j}^2 \inner{b^{\star}_i,b_k}^2 $, (2)
$-\frac{\hat{\sigma}_4\mu }{\sqrt{6}}\sum_{i,j\in [d]} \alpha_i\inner{b^{\star}_i,b_j}^4$, (3) $\lambda \sum_{j=1}^m((\sum_{i=1}^m \alpha_i\inner{b_j,b^{\star}_i}^2 - 1)^2-1)^2$, (4) $\frac{\delta}{2}\|B\|_F^2$.

Among these 4 terms, the first, third and forth terms all contribute positively to this inner-product (because when $b_i$ is moved to $(1-\epsilon) b_i$ all those terms clearly decrease). Term 2 $-\frac{\hat{\sigma}_4\mu}{\sqrt{6}} \sum_{i,j\in[m]}a^{\star}_i\inner{b^{\star}_i,b_j}^4$ contribute negatively. Therefore we can ignore terms 1 and 4:
$$
\langle \frac{\partial  \mathcal{F}_{a^{\star},\mu,\lambda,\delta}}{\partial b_i}(B), b_i \rangle \ge \langle \frac{\partial}{\partial b_i}[-\frac{\hat{\sigma}_4\mu}{\sqrt{6}} \sum_{i\in[m]}a^{\star}_i\inner{b^{\star}_i,b_j}^4 + \lambda (\sum_{i\in[m]}a^{\star}_i\inner{b^{\star}_i,b_j}^2 - 1)^2], b_i \rangle.
$$

Let $b_i^\top M b_i = C \ge 2$, we know $\sum_{i\in[m]}a^{\star}_i\inner{b^{\star}_i,b_j}^4 \le \frac{1}{a^{\star}_{min}}\sum_{i\in[m]}(a^{\star}_i)^2\inner{b^{\star}_i,b_j}^4 \le C^2/a^{\star}_{min}$. Therefore,
$$
\langle \frac{\partial}{\partial b_i}[-\frac{\hat{\sigma}_4\mu}{\sqrt{6}} \sum_{i\in[m]}a^{\star}_i\inner{b^{\star}_i,b_j}^4], b_i\rangle = -\frac{4\hat{\sigma}_4\mu}{\sqrt{6}} \sum_{i\in[m]}a^{\star}_i\inner{b^{\star}_i,b_j}^4 \ge -\frac{4\hat{\sigma}_4\mu}{\sqrt{6}} \cdot \frac{C^2}{a^{\star}_{min}}
$$

On the other hand,
$$
\langle \frac{\partial}{\partial b_i}[\lambda (\sum_{i\in[m]}a^{\star}_i\inner{b^{\star}_i,b_j}^2 - 1)^2, b_i\rangle = 4\lambda(b_i^\top M b_i - 1)(b_i^\top M b_i) = 4\lambda C(C-1).
$$
By the choice of $\lambda, \mu$, we can see that the negative term is negligible, and we know
$$
\langle \frac{\partial  \mathcal{F}_{a^{\star},\mu,\lambda,\delta}}{\partial b_i}(B), b_i \rangle \ge 2\lambda C(C-1)
$$
Since $b_i^\top M b_i = C$, we have $\|b_i\|\le \sqrt{C/\sigma_{min}(M)}$. Therefore the norm of the gradient is at least $2\lambda C(C-1)/\|b_i| \ge 2\sqrt{2}\lambda \sigma_{min}(M)$, this contradicts with the assumption. The norm of the rows must all be bounded.
\end{proof}

\begin{proof}[Proof of Lemma~\ref{lem:inside}]
We have that $b_i^\top M b_i \le 2$ implies $\|b_i\|^2 \le 2/\sigma_{min}(M)$. The norm of the whole matrix is bounded by$\|B\|_F \le \sqrt{\sum_{i=1}^m \|b_i\|^2} \le \sqrt{2m/\sigma_{min}(M)}$, so by triangle inequality we have
$$
\|\nabla F_{a^{\star},\mu,\lambda}(B)\| \le \|\nabla \mathcal{F}_{a^{\star},\mu,\lambda,\delta}(B)\| + \delta\|B\|_F \le \epsilon + \delta \sqrt{2m/\sigma_{min}(M)}.
$$
\end{proof}
\noindent Finally we are ready to prove Theorem~\ref{thm:main-under}.

\begin{proof} [Proof of Theorem~\ref{thm:main-under}]
We will separate $B$ into two components $B_\mathcal S = BP_{\mathcal S} $ and $B_\perp =  BP_{\mathcal S^\perp}$.

We will first show that $B_{\mathcal S}$ is close to the desirable solution. To do that we will use Theorem~\ref{thm:main-nonortho} \footnote{If we restrict all the vectors to the subspace $\mathcal S$, we can still apply Theorem~\ref{thm:main-nonortho} as long as we replace all inverses with pseudo-inverses.}. By the choice of $\epsilon, \delta$, we know from Lemma~\ref{lem:inside} that $\|\nabla F_{a^{\star},\mu,\lambda}(B_\mathcal S)\| \le 2\delta \sqrt{2m/\sigma_{min}(M)}$. Also, $\nabla^2 F_{a^{\star},\mu,\lambda}(B_\mathcal S) \ge \nabla^2 \mathcal{F}_{a^{\star},\mu,\lambda}(B) - \delta \ge -\tau_0$. Therefore we know $B_{\mathcal S}$ must be of the form
$$
[B_{\mathcal S}^{\dag}]^\top = P D\Gamma B^{\star}+ E_1,
$$
where $\|E_1\| < \epsilon_0/2$. Also at the same time from the proof of Theorem~\ref{thm:main-nonortho} we know $B_{\mathcal S} = (PDO+E')W^\top$ where $PDO+E'$ has singular values close to 1. Therefore $\sigma_{min}(B_{\mathcal S}) \ge \sigma_{min}(W)/2 = 1/2\sqrt{\|M\|}$.

By Lemma~\ref{lem:outside} we know $\|B_\perp\|_F \le \epsilon/\delta$. We apply inverse matrix perturbation (Theorem~\ref{thm:matpert}) again, using $B = B_{\mathcal S}+B_\perp$, therefore we know
$$
B^{\dag} = B_{\mathcal S}^{\dag} + E_2,
$$
where $\|E_2\| \le O(\|B_\perp\|_F/\sigma_{min}^2(B_{\mathcal S})) \le \epsilon_0/2$.

Combining these two perturbations we know
$$
[B^{\dag}]^\top = PD\Gamma B^{\star}+ E_1+E_2^\top,
$$
and the error term $E_1+E_2^\top$ has spectral norm at most $\epsilon_0$.
\end{proof}

\subsection{Toolbox: Matrix Perturbation}

In the proof we used the following theorem for the perturbation of matrices.

\begin{theorem}[Stewart and Sun~\cite{stewart1990matrix}]\label{thm:matpert}
Consider the perturbation of a matrix $A$: if $B = A+E$,then we have
$$
\|B^{\dag} - A^{\dag}\| \le \sqrt{2} \|A^{\dag}\| \|B^{\dag}\| \|E\|.
$$

As a corollary, if $\|E\| \le \sigma_{min}(A)/2$, then we have
$$
\|B^{\dag} - A^{\dag}\| \le 2\sqrt{2} \sigma_{min}(A)^{-2} \|E\|.
$$
\end{theorem}

\subsection{Local Optimality of All the Variables} \label{sec:landscape:multi-columns}

In this section we prove Theorem~\ref{thm:main-general}. Results in Subsection~\ref{sec:landscape:single-column} have established that if $B$ is a local minimum, then each row $b_s$ of $B$ has to be a scaling of a basis vector. In this section we show that these basis vectors need to be distinct from each other. The following proposition summaries such a claim (with a weak error analysis).
\begin{proposition}

	In the setting of Theorem~\ref{thm:main-general}, suppose
 $B$ satisfies $$\norm{\nabla g(B)}\le \epsilon \textup{  and  }\lambda_{\min}(\nabla^2 g(B)) \ge - \tau $$ for parameters \sloppy $\tau,\epsilon$ satisfying $0\le \tau \le c\min\{\mu \beta_{\min}/(\kappa_\alpha d),\lambda\}$ and $\epsilon\le c\min\{\alpha_{\min}, \sqrt{\tau^3/\beta_{\min}}\}$. Then, the matrix $B$ can be written as $$B =  DP + E\,,$$
 where $D$ is diagonal such that $\forall i, |D_{ii}|\in [1/4,2]$, and $P$ is a permutation matrix, and $|E|_\infty \le \delta$ with $\delta = \left(\frac{\tau}{\mu \beta_{\min}}\right)^{1/2} $.
 \label{prop:main-prop}
\end{proposition}

As alluded before, in the proof we will first apply the results in Section~\ref{sec:landscape:single-column} to show that when $B$ is a local minimum, each row $b_s$ has a unique large entry. Then we will show that the largest entries of each row sit on different columns. The key intuition behind the proof is that if two rows, say row $s,t$, have their large entries on the same column, then it means that there exists a column--- say column $k$ --- that doesn't contain largest entry of any row. Then either row $s$ or $t$ will violate Lemma~\ref{lem:subproblem-ub-alpha}. Or in other words, either row $s$ or $t$ can move their mass into the column $k$ to decrease the function value. This contradicts the assumption that $B$ is a local minimum.
\begin{proof}
As pointed in the paragraph below equation~\eqref{eqn:def:h}, when we restrict our attention to a particular row of $B$ and fix the rest of the rows the function $G_{\alpha,\beta,\mu}$ reduces to the function $h(\cdot)$ in equation~\eqref{eqn:def:h} so that we can apply lemmas in Section~\ref{sec:landscape:single-column}.

Concretely, fix an index $s\in [d]$ and let $x=b_s$. For all $i\in [d]$, let $\bar \alpha_i = \alpha_i \sum_{j \neq s} (b_j ^\top e_i)^2$, and $\bar \beta_i = \mu \beta_i $. Then we have that
\begin{align}
G_{\alpha,\beta,\mu}(B)  = \sum_{i=1}^d \bar{\alpha}_i x_i ^2 - \sum_{i} \bar{\beta}_i x_i ^4 + \lambda \big( \norm{x}^2-1\big)^2  \label{eqn:12}
\end{align}

We view the function above as $h(x)$. Now we apply Lemma~\ref{lem:second-largest-and-norm} (by replacing $\alpha,\beta$ in Lemma~\ref{lem:second-largest-and-norm} by $\bar{\alpha},\bar{\beta}$). The assumption that $\lambda_{\min}(\nabla^2 g_{\alpha,\beta,\mu}(B)) \ge -\tau$ implies that $\lambda_{\min}(\nabla^2 h(x)) \ge -\tau $ since $\nabla^2 h(x)$ is a submatrix of $\nabla ^2 g(B)$.\sloppy~Moreover, $\norm{\nabla h(x)} \le \norm{\nabla G(B)}\le \epsilon \le \sqrt{\tau^3/(\mu \beta_{\min})}$

Hence by Lemma~\ref{lem:second-largest-and-norm}, we have that the second largest entry of $|b_s|$ satisfies
\begin{align}\forall s, \sl{b_s} \le \delta .
\label{eq:second-largest-delta}
\end{align}
where $\delta \triangleq  \left(\frac{\tau}{\mu \beta_{\min}}\right)^{1/2}$ for the ease of notation. We can check that $\delta \le \frac{1}{4\sqrt{\kappa_\alpha d}}$ by the assumption. Therefore, we have essentially shown that each row of $B$ has only one single large entry, since the second largest entry is at most $\delta$.

Next we show that each row of $B$ has largest entries on distinct columns. For each row $j\in [d]$, let $i_j = \arg\max_{i} | e_i ^\top b_j|$ be the index of the largest entry of $b_j$. We will show that $i_1, \ldots, i_d$ are distinct.

For the sake of contradiction, suppose they are not distinct, that is, there are two distinct rows $s,t$ that have the same largest entries on column $l$, that is, we assume that $i_s =i_t =l$.
This implies that $\{i_1,\ldots, i_d\}\neq [d]$ and let $k\in [d]$ be the index such that $k \notin \{i_1,\ldots, i_d\}$.
We note that by the assumption $\delta = \left(\frac{\tau}{\mu \beta_{\min}}\right)^{1/2} \le \frac{1}{4 \sqrt{\kappa_\alpha d}}\le \frac{1}{4\sqrt{d}}$. We first bound from above $\bar \alpha_k$
\begin{align}
\bar \alpha_k &= \alpha_k \sum_{j \neq s} (b_j^\top e_k)^2 \le \alpha_k d \delta^2\le \frac{1}{16}\alpha_{\min}. \tag{by $\delta \le \frac{1}{4\sqrt{\kappa_\alpha d}}$ }
\end{align}

Assume in addition without loss of generality that $|b_s^\top e_l|\le |b_t^\top e_l|$. Let
\begin{align}
z_l \triangleq \sum_{j \neq s} (b_j^\top e_{l})^2
\label{eqn:13}
\end{align}
be the sum of squares of the entries on the column $l$ without entry $b_j ^\top e_l$, and that  $\bar \alpha_l =\alpha_l z_l$ . We first prove that $z_l \ge 1/3$.

For the sake of contradiction, assume $z_l < 1/3.$  Then we have that $
\bar \alpha_l = \alpha_l z \le \frac{1}{3} \alpha_l \,.$
This implies that $\lambda \ge 4\max\{\bar{\alpha}_{l}, \tau\}$, and since $l$ is the index of the largest column of $b_s$  we can invoke Lemma~\ref{lem:norm} and conclude that $\norm{b_s}^2\ge 1/2$. This further implies that
\begin{align}
(b_s^\top e_l)^2 \ge \norm{b_s}^2 -d\sl{b_s}\ge 1/2 - d\delta^2\ge 1/3 \tag{by $\delta \le 1/(4\sqrt{d})$}
\end{align}
Since we have assumed that $|b_s^\top e_l|\le |b_t^\top e_l|$. Then we obtain that
\begin{align}
z_l \ge |b_t^\top e_l|^2 \ge |b_s^\top e_l|^2 \ge 1/3\,,\nonumber
\end{align}
which contradicts the assumption. Therefore, we conclude that $z_l\ge 1/3$.  Then we are ready to bound $\bar{\alpha}_l$ from below:
\begin{align}
\bar \alpha_l &= \alpha_l z_l \ge \frac{1}{3}\alpha_{l} \,.\nonumber
\end{align}
The display above and Equation \eqref{eqn:13}  implies that
\begin{align}
\bar{\alpha}_l - \bar{\alpha}_k \ge \frac{1}{4}\alpha_{\min}. \label{eqn:10}
\end{align}

Note that $l$ is the largest entry in absolute value in the vector $b_s$. We will apply Lemma \ref{lem:subproblem-ub-alpha}. We fix every row of $B$ except $b_s$ and consider the objective as a function of $b_s$ only.  Again let $\bar \alpha_i = \alpha_i \sum_{j \neq s} (b_j ^\top e_i)^2$, and $\bar \beta_i = \mu \beta_i $ and we have the equation~\eqref{eqn:12}.  (Note that now $\bar{\alpha}$ depends on the choice of $s$ which we fixed.)  Lemma~\ref{lem:subproblem-ub-alpha} gives us that
\begin{align*}
\bar \alpha_{l}& \le \bar \alpha_k +2\epsilon+2\tau+4\bar{\beta_\ell}.
\end{align*}
Since $\epsilon \le \frac{1}{50}\alpha_{\min}$, $\tau \le \frac{1}{50}\alpha_{\min}$ and $\bar{\beta}_l = \mu \beta_l \le \frac{1}{50}\alpha_{\min}$, we obtain that
\begin{align}
\bar \alpha_{l}& \le \bar \alpha_k + \frac{1}{5}\alpha_{\min}
\end{align}
which contradicts equation~\eqref{eqn:10}.  Thus we have established that $i_1,\dots, i_d$ are distinct.

Finally, let $Q$ be the matrix that only contain the largest entries (in absolute value) of each columns of $B$. Since $i_1,\dots, i_d$ are distinct, we have that $Q$ contains exactly one entry per row and per column. Therefore $Q$ can be written as $DP$ where $P$ is a permutation matrix and $D$ is a diagonal matrix. Moreover, we have that $\norm{b_s}_{\infty}^2 \ge \norm{b_s}^2 - d\sl{b_s}^2 \ge  1/4 $ and $\norm{b_s}^2 \le 2$. Therefore, the largest entry of each row has absolute value between $1/4$ and $2$. Therefore $|D|_{ii} \in [1/4,2]$. Let $E = B-PD$. Then we have that $|E|_\infty \le \max_s \sl{b_s}\le \delta$,which completes the proof.

\end{proof}

\noindent Applying Lemma~\ref{lem:exact-recovery}, we can further strengthen Proposition~\ref{prop:main-prop} with better error bounds and better control of the largest entries of each column.
\begin{proposition}[Strengthen of Proposition~\ref{prop:main-prop}]\label{prop:strengthen}
	In the setting of Proposition~\ref{prop:main-prop}. Suppose in addition that $\tau$ satisfies
	$\tau \le c\mu \beta_{\min}^2/\beta_{\max}$.  Then, the matrix $B$ can be written as $$B =  DP + E\,,$$
	where $P$ is a permutation matrix, $D$ is diagonal such that
	$$\forall i\in [d],~~ \frac{1}{1-\frac{\mu\beta_i}{\lambda}} \left(1-  \frac{18d \epsilon^2}{\beta_{\min}^2}-\frac{\epsilon}{2\lambda}\right)\le |D_{ii}|^2\le \frac{1}{1-\frac{\mu\beta_i}{\lambda}}\left(1+\frac{\epsilon}{2\lambda}\right)$$  and $$|E|_\infty \le\frac{3\epsilon}{\beta_{\min}}.$$

\end{proposition}

\begin{proof}
	By Proposition~\ref{prop:main-prop}, we know that $|E|_\infty\le \delta = \left(\frac{\tau}{\mu \beta_{\min}}\right)^{1/2}$. Now we use Lemma~\ref{lem:exact-recovery} to strength the error bound.

As we have done in the proof of Proposition~\ref{prop:main-prop}, we again fix an arbitrary $s\in [d]$ and all the rows except $b_s$ and view $G_{\alpha,\beta,\mu}$ as a function of $b_s$. For all $i\in [d]$, let $\bar \alpha_i = \alpha_i \sum_{j \neq s} (b_j ^\top e_i)^2$, and $\bar \beta_i = \mu \beta_i $ and view $G_{\alpha,\beta,\mu}$ as a function of the form $h(x)$ with $\alpha,\beta$ replaced by $\bar{\alpha},\bar{\beta}$, namely,
\begin{align}
h(x) = \sum_{k} \bar \alpha_k x_k ^2- \sum_{k} \bar{\beta}_k x_k ^4 + \lambda \big( \norm{x}^2-1\big)^2 + \textup{const} \nonumber
\end{align}

We will verify the condition of Lemma~\ref{lem:exact-recovery}. Let $i$ be the index of the largest entry in absolute value of the vector $b_s$. Since we have shown that the largest entry in each row sits on different columns, and the second largest entry is always less than $\delta$, we have that,
\begin{align}
\bar \alpha_i &= \alpha_i \sum_{j \neq s} (b_j^\top e_i)^2 \le \alpha_i d \delta^2\le \frac{1}{16}\alpha_{\min}. \tag{by $\delta \le \frac{1}{4\sqrt{\kappa_\alpha d}}$ }
\end{align}
For any $k\neq i$, we know that the column $k$ contains some entry $(k,j_k)$ which is the largest entry of some row, and we also have that $j_k\neq s$ since the largest entry of row $s$ is on column $i$. Therefore, we have that
\begin{align}
\bar \alpha_k &= \alpha_k \sum_{j \neq s} (b_j ^\top e_k)^2\ge \alpha_l (b_{j_k} ^\top e_k)^2 \ge \alpha_l ( \norm{b_{k}}^2 - d \delta^2)\nonumber\\
&\ge \frac{1}{3}\alpha_{l} \tag{by $\delta \le 1/(4\sqrt{d})$}
\end{align}
Therefore, $\bar{\alpha_k} \ge \bar{\alpha_i}$ for any $k\neq i$ and thus $i = \textup{argmin}_k |\bar{\alpha}_k|$.
By the fact that $|E|_\infty\le \delta$, we have that $\norm{x_{-i}}_\infty\le \delta \le 0.1\min\{ 1/\sqrt{d}, \sqrt{\beta_{\min}/(\beta_{\max})}\}$. Now we are ready to apply Lemma~\ref{lem:exact-recovery} and obtain that $\sl{b_s}\le \frac{3\epsilon}{\beta_{\min}}$. Applying the argument for every row $s$ gives  $|E|_\infty\le \frac{3\epsilon}{\beta_{\min}}$.

Finally, we give the bound for the entires in $D$.
 Let $v$ be a short hand for $\nabla h(b_s)$ which is equal to the $s$-th column of $\nabla G(B)$. Since $B$ is an $\epsilon$-approximate stationary point, then we have that $\norm{v}\le \epsilon$ and by straightforward calculation of the gradient, we have
\begin{align*}
v_i &=2 \bar \alpha_i x_i -4 \mu \beta_i x_i ^3 + 4\lambda( \sum_{j=1}^d x_j ^2 -1 ) x_i \,.
\end{align*}
Since $x_i\neq 0$, dividing by $x_i$ gives,
\begin{align*}
0 &=2\bar \alpha_i -4 \mu \beta_i x_i ^2 + 4\lambda( \sum_{j=1}^d x_j ^2 -1) - \frac{v_i }{x_i}\\
&= (4\lambda - 4 \mu \beta_i ) x_i ^2+ 4\lambda \sum_{j \neq i} x_j ^2 +2 \bar \alpha_i - 4\lambda - \frac{v_i}{x_i}
\end{align*}
Rearranging the equation above gives,
\begin{align*}
x_i ^2 &= \frac{1}{4 \lambda - 4 \mu \beta_i} \left(4\lambda-2 \bar \alpha_i  -4\lambda\sum_{j \neq i} x_j ^2 - \frac{v_i}{x_i}\right)\\
&=\frac{1}{1- \frac{\mu \beta_i}{\lambda}} \left( 1- \frac{\bar \alpha_i}{2\lambda}- \sum_{j \neq i} x_j ^2 -\frac{v_i}{4\lambda x_i}\right)
\end{align*}
To upper bound $x_i ^2 $ , we note that $|v_i| < \epsilon$, $\bar \alpha_i >0$, and $\sum_{j\neq i} x_j ^2 >0$, so
\begin{align}
x_i ^2  &\le \frac{1}{\left(1-\frac{\mu\beta_i}{\lambda}\right)}\left(1+\frac{\epsilon}{2\lambda}\right) \le 1 + \frac{2\mu \beta_i +\epsilon}{\lambda} \tag{since $\lambda \ge 4\mu \beta_i$}
\end{align}
For the lower bound of $x_i ^2$, we note that $|E|_\infty  \le \delta = \frac{3 \epsilon}{\beta_{\min}}$ implies $\sum_{j \neq i} x_j ^2 \le d \delta^2$. Moreover, we have proved that each rows has largest entry at different columns. Also note that the largest entry of row $b_s$ is on column $i$. Therefore, we have $\bar \alpha_i = \alpha_i \sum_{j \neq s} (b_j ^T e_i)^2\le \alpha_{\max} d \delta^2$. Using these two estimates and $\delta = \frac{3\epsilon}{\beta_{\min}}$, we have
\begin{align*}
x_i ^2 &\ge \frac{1}{1-\frac{\mu\beta_i}{\lambda}}( 1- (\frac{\alpha_{\max}}{2\lambda} +1) d \delta^2 - \frac{\epsilon}{2\lambda})\\
&= \frac{1}{1-\frac{\mu\beta_i}{\lambda}}\left(1-  \frac{18d \epsilon^2}{\beta_{\min}^2}-\frac{\epsilon}{2\lambda}\right)
\end{align*}
\end{proof}

\noindent Finally we are ready to prove Theorem~\ref{thm:main-general} by applying Proposition~\ref{prop:main-prop}.

\begin{proof}[Proof of Theorem \ref{thm:main-general}]
By setting  $\epsilon=0, \tau = 0$ in Proposition \ref{prop:main-prop}, we have that any local minimum $B$ satisfies that $B = DP$ where $P$ is a permutation matrix and $D$ is a diagonal and the precise diagonal entries of $D$. It can be verified  that all these points have the same function value, so that they are all global minimizers.

Towards proving the second bullet, we note that a saddle point $B$ satisfies that $\nabla G(B) =0$. We will prove that $\lambda_{\min}(\nabla^2 G(B))\le  -\tau_0$. For the sake of contradiction, suppose $\lambda_{\min}(\nabla^2 G(B)) \ge - \tau_0$. Then setting $\epsilon = 0$ and $\tau = \tau_0$ in  Propostion~\ref{prop:strengthen}, we have that $B = DP$ and $D_{ii} =  \left\{\pm \sqrt{\frac{1}{1-\mu \beta_i/\lambda}}\right\}$, which by bullet 1 implies that $B$ is a local minimum. This contradicts the assumption that $B$ is a saddle point.

The 3rd bullet is a just a rephrasing of Proposition~\ref{prop:strengthen}.
\end{proof}

\section{Overview: Landscape Design and Analysis}\label{sec:overview}

In this section, we present a general overview of ideas behind the design of objective function $G(\cdot)$. Inspired by the formula~\eqref{eqn:popluation_risk_tensor}, in Section~\ref{sec:overview:samples}, we envision a family of possible objective functions for which we have unbiased estimators via samples. In Section~\ref{sec:overview:localmin}, we pick a specific function that feeds our needs: a) it has no spurious local minimum;  b) the global minimum corresponds to the ground-truth parameters.

\subsection{Which objective can be estimated by samples?} \label{sec:overview:samples}

Recall that in equation~\eqref{eqn:population_risk} of Theorem~\ref{thm:population_risk} we give an analytic formula for the straightforward population risk $f$. Although the population risk $f$ doesn't perform well empirically, the lesson that we learn from it help us design better objective functions. One of the key fact that leads to the proof of Theorem~\ref{thm:population_risk} is that for any continuous and bounded function $\gamma$, we have that
\begin{align}
\mathop{\mathbb E}\displaylimits\left[y\cdot \gamma(b_i^\top x)\right] = \sum_{k \in\mathbb{N}}\hat{\gamma}_k\hat{\sigma}_k\left(\sum_{j\in [d]} a^{\star}_j \inner{b^{\star}_j,b_i}^k\right) \,.\nonumber
\end{align}
Here $\hat{\sigma}_k$ and $\hat{\gamma}_k$ are the $k$-th Hermite coefficient of the function $\sigma$ and $\gamma$. That is, letting $h_k$ the $k$-th normalized probabilists' Hermite polynomials~\cite{wiki-hermite} and $\inner{\cdot,\cdot}$ be the standard inner product between functions, we have $\hat{\sigma}_k = \inner{h_k, \sigma}$.

Note that $\gamma$ can be chosen arbitrarily to extract different terms. For example, by choosing $\gamma = h_k$, we obtain that
\begin{align}
\mathop{\mathbb E}\displaylimits\left[y\cdot h_k(b_i^\top x)\right] = \hat{\sigma}_k \sum_{j\in [d]} a^{\star}_j \inner{b^{\star}_j,b_i}^k \,.\label{eqn:6}
\end{align}
That is, we can always access functions forms that involves weighted sum of the powers of $\inner{b^{\star}_i,b_i}$, as in RHS of equation~\eqref{eqn:6}.

Using a bit more technical tools in Fourier analysis (see details in Section~\ref{sec:formula}), we claim that most of the symmetric polynomials over variables $\inner{b^{\star}_i,b_j}$ can be estimated by samples:

\begin{claim}[informal]\label{claim:informal} For an arbitrary polynomial $p()$ over a single variable, there exits a corresponding function $\phi^p$ such that
	\begin{align}
	\mathop{\mathbb E}\displaylimits\left[y\cdot \phi^p(B,x)\right] = \sum_j a^{\star}_j \sum_i p(\inner{b^{\star}_j,b_i})
	\label{eqn:14}
	\end{align}
	Moreover, for an any polynomial $q(\cdot,\cdot)$ over two variables, there exists corresponding $\phi^q$ such that
	\begin{align}
	\mathop{\mathbb E}\displaylimits\left[y \cdot \phi^q(B,x)\right] = \sum_j a^{\star}_j \sum_{i,k} q(\inner{b^{\star}_j,b_i}, \inner{b^{\star}_k,b_i})\label{eqn:15}
	\end{align}
\end{claim}
\noindent We will not prove these two general claims.
Instead, we only focus on the formulas in Theorem~\ref{thm:formula} and Theorem~\ref{thm:4thpowerformula}, which are two special cases of the claims above.

Motivated by Claim~\ref{claim:1}, in the next subsection, we will pick an objective function which has no spurious local minimum among those functional forms on the right-hand sides of equation~\eqref{eqn:14} and~\eqref{eqn:15}.

\subsection{Which objective has no spurious local minima?}\label{sec:overview:localmin}

As discussed briefly in the introduction, one of the technical difficulties to design and analyze objective functions for neural networks comes from the permutation invariance --- if a matrix $B$ is a good solution, then any permutation of the rows of $B$ still gives an equally good solution (if we also permute the coefficients in $a$ accordingly). We only know of a very limited number of objective functions that guarantee to enjoy permutation invariance and have no spurious local minima~\cite{ge2015escaping}.

We start by considering the objective function used in~\cite{ge2015escaping},
\begin{align}
\min ~& P(B)= \sum_i \sum_{j\neq k} \inner{b^{\star}_i, b_j}^2\inner{b^{\star}_i,b_k}^2 \nonumber \\
s.t. ~& \forall i\in [d], \norm{b_i} = 1
\end{align}
Note that here we overload the notation by using $b^{\star}_i$'s to denote a set of fixed vectors that we wanted to recover and using $b_i$'s to denote the variables. Careful readers may notice that $P(B)$ doesn't fall into the family of functions that we described in the previous section (that is, RHS equation of~\eqref{eqn:14} and~\eqref{eqn:15}), because it lacks the weighting $a^{\star}_i$'s. We will fix this issue later in the subsection.  Before that we first summarize the nice properties of the landscape of $P(B)$.

For the simplicity of the discussion, let's assume $B^{\star}=\begin{bmatrix}(b^{\star}_1)^\top\\ \vdots \\ (b^{\star}_d)^\top\end{bmatrix}$ forms an orthonormal matrix in the rest of the subsection. Then, any permutation and sign-flip of the rows of $B^{\star}$ leads to a global minimum of $P(\cdot)$ --- when $B = SQB^{\star}$ with a permutation matrix $Q$ and a sign matrix $S$ (diagonal with $\pm 1$), we have that $P(B) =0$ because one of $\inner{b^{\star}_i, b_j}^2$ and $\inner{b^{\star}_i,b_k}^2$ has to be zero for all $i,j,k$\footnote{Note that $B^{\star}$ is orthogonal, and $j\neq k$}).

It turns out that these permutations/sign-flips of $B^{\star}$ are also the only local minima\footnote{We note that since there are constraints here, by local minimum we mean the local minimum on the manifold defined by the constraints. } of function $P(\cdot)$. To see this, notice that $P(B)$ is a degree-2 polynomial of $B$. Thus if we pick an index $s$ and fix every row except for $b_s$, then $P(B)$ is a quadratic function over unit vector $b_s$ \--- reduces to an smallest eigenvector problem. Eigenvector problems are known to have no spurious local minimum. Thus the corresponding function (w.r.t $b_s$) has no spurious local minimum. It turns out the same property still holds when we treat all the rows as variables and add the row-wise norm constraints (see proof in \cite{ge2015escaping}).

However, there are two issues with using objective function $P(B)$. The obvious one is that it doesn't involve the coefficients $a^{\star}_i$'s and thus doesn't fall into the forms of equation~\eqref{eqn:15}. Optimistically, we would hope that for nonnegative $a^{\star}_i$'s the weighted version of $P$ below would also enjoy the similar landscape property
\begin{align}
P'(B) = \sum_i a^{\star}_i\sum_{j\neq k} \inner{b^{\star}_i, b_j}^2\inner{b^{\star}_i,b_k}^2 \nonumber
\end{align}
When $a^{\star}_i$'s are positive, indeed the global minimum of $P'$ are still just all the permutations of the $B^{\star}$.\footnote{This is the main reason why we require $a^{\star}\ge 0$.} However, when $\max a^{\star}_i > 2\min a^{\star}_i$, we found that $P'$ starts to have spurious local minima . It seems that spurious local minimum often occurs when a row of $B$ is a linear combination of a smaller number of rows of $B^{\star}$. See Section~\ref{sec:example} for a concrete example.

To remove such spurious local minima, we add a regularization term below that pushes each row of $B$ to be close to one of the rows of $B^{\star}$,
\begin{align}
R(B) = -\mu \sum_i a^{\star}_i \sum_{j} \inner{b^{\star}_i,b_j}^4
\end{align}
We see that for each fixed $j$, the part in $R(B)$ that involves $b_j$ has the form
\begin{align}
-\mu \sum_{i} a^{\star}_i\inner{b^{\star}_i,b_j}^4  = -\mu \inner{\sum_i a^{\star}_i{b^{\star}_i}^{\otimes 4}, b_j^{\otimes 4}}
\end{align}
This is commonly used objective function for decomposing tensor $\sum_i a^{\star}_i{b^{\star}_i}^{\otimes 4}$. It's known that for orthogonal $b^{\star}_i$'s,  the only local minima are $\pmb^{\star}_1,\dots, \pmb^{\star}_d$~\cite{ge2015escaping}. Therefore, intuitively $R(B)$ pushes each of the $b_i$'s towards one of the $b^{\star}_i$'s. \footnote{However, note that $R(B)$ by itself doesn't work because it does not prevent the solutions where all the $b_i$'s are equal to the same $b^{\star}_j$.} Choosing $\mu$ to be small enough, it turns out that $P'(B) + R(B)$ doesn't have any spurious local minimum as we will show in Section~\ref{sec:landscape}.

Another issue with the choice of $P'(B) +R(B)$ is that we are still having a constraint minimization problem. Such row-wise norm constraints only make sense when the ground-truth $B^{\star}$ is orthogonal and thus has unit row norm. A straightforward generalization of $P(B)$ to non-orthogonal case requires some special constraints that also depend on the covariance matrix $B^{\star}{B^{\star}}^\top$, which in turn requires a specialized procedure to estimate. Instead, we move the constraints into the objective function by considering adding another regularization term that approximately enforces the constraints.

It turns out the following regularizer suffices for the orthogonal case,
\begin{align}
S(B) = \lambda \sum_i (\norm{b_i}^2-1)^2\,.
\end{align}
Moreover, we can extend this easily to the non-orthogonal case (see Section~\ref{sec:general}) without estimating any statistics of $B^{\star}$ in advance.
We note that $S(B)$ is not the Lagrangian multiplier and it does change the global minima slightly. We will take $\lambda$ to be large enough so that $\norm{b_i}$ has to be close to 1. As a summary, we finally use the unconstrained objective
\begin{align}
\min G(B) \triangleq P'(B) + R(B) + S(B)\nonumber
\end{align}
Since $R(B)$ and $S(B)$ are degree-4 polynomials of $B$, the analysis of $G(B)$ is much more delicate, and we cannot use much linear algebra as we could for $P'(B)$. See Section~\ref{sec:landscape} for details.

Finally we note that a feature of this objective $G(\cdot)$ is that it only takes $B$ as variables. We will estimate the value of $a^{\star}$ after we recover the value of $B$. (see Section~\ref{sec:linear}).
\section{Main Results}

\subsection{Connecting $\ell_2$ Population Risk with Tensor Decomposition} \label{sec:main-formula}

We first show that a natural $\ell_2$ loss for the one-hidden-layer neural network can be interpreted as simultaneously decomposing tensors of different orders.

A straightforward approach of learning the model~\eqref{eqn:model} is to parameterize the prediction by
\begin{align}
\hat{y} = a^{\top}\sigma(Bx)\,,
\end{align}
where $a\in \mathbb R^d, B\sim \mathbb R^{m \times d}$ are the training parameters. Naturally, we can use $\ell_2$ as the empirical loss, which means the population risk is
\begin{align}
f(a,B) = \mathop{\mathbb E}\displaylimits\left[\norm{\hat{y} - y}^2\right]\,. \label{eqn:population_risk}
\end{align}

\noindent Throughout the paper, we use ${b^{\star}_1}^\top,\dots, {b^{\star}_m}^\top$ to denote the row vectors of $B^{\star}$ and similarly for $B$. That is, we have $B = \begin{bmatrix}
b_1^\top\\
\vdots\\
b_m^\top
\end{bmatrix}$ and $B^{\star} = \begin{bmatrix}
{b^{\star}_1}^\top\\
\vdots\\
{b^{\star}_m}^\top
\end{bmatrix}$. Let $a_i$ and $a^{\star}_i$'s be the coordinates of $a$ and $a^{\star}$ respectively.

\noindent We give the following analytic formula for the population risk defined above.

\begin{theorem}\label{thm:population_risk} Assume vectors $b_i,b^{\star}_i$'s are unit vectors. Then, the population risk $f$ defined in equation~\eqref{eqn:population_risk} satisfies that
	\begin{align}
	f(a,B)        &=  \sum_{k\in\mathbb{N}} \hat{\sigma}_k^2 \Norm{\sum_{i\in [m]} a^{\star}_i {b^{\star}_i}^{\otimes k} - \sum_{i\in [m]} a_i b_i^{\otimes k} }_F^2 + \textup{const}\,. \label{eqn:popluation_risk_tensor}
	\end{align}
	where $\hat{\sigma}_k$ is the $k$-th Hermite coefficient of the function $\sigma$. See section~\ref{sec:basic-hermite} for a short introduction of Hermite polynomial basis. \footnote{  When $\sigma = ReLU$, we have that $\hat \sigma_0=\frac{1}{\sqrt{2\pi}}$, $\hat \sigma_1=\frac{1}{2}$. For $n\ge 2$ and even, $\hat \sigma_n = \frac{((n-3)!!)^2 }{\sqrt{2 \pi n!}}$. For $n \ge 2$ and odd, $\hat \sigma_n =0$.}

\end{theorem}

\noindent \textit{Connection to tensor decomposition: }We see from equation~\eqref{eqn:popluation_risk_tensor} that the population risk of $f$ is essentially an average of infinite number of loss functions for tensor decomposition. For a fixed $k\in \mathbb{N}$, we have that the $k$-th summand in equation~\eqref{eqn:popluation_risk_tensor} is equal to (up to the scaling factor $\hat{\sigma}_k^2$)
\begin{align}f_k \triangleq \Norm{T_k- \sum_{i\in [m]} a_i b_i^{\otimes k} }_F^2\,. \label{eqn:f4}\end{align}
where $T_k = \sum_{i\in [m]} a^{\star}_i {b^{\star}_i}^{\otimes k}$ is a $k$-th order tensor in $(\mathbb R^{d})^{\otimes k}$.  We note that the objective $f_k$ naturally attempts to decompose the $k$-order rank-$m$ tensor $T_k$ into $m$ rank-1 components $a_1b_i^{\otimes k},\dots, a_mb_m^{\otimes k}$.

The proof of Theorem~\ref{thm:population_risk} follows from using techniques in Hermite Fourier analysis, which is deferred to Section~\ref{sec:landscape:f}.

\paragraph{Issues with optimizing $f$:} It turns out that optimizing the population risk using stochastic gradient descent is empirically difficult. Figure~\ref{fig:relu} shows that in a synthetic setting where the noise is zero, the test error empirically doesn't converge to zero for sufficiently long time with various learning rate schemes, even if we are using fresh samples in iteration. This suggests that the landscape of the population risk has some spurious local minimum that is not a global minimum. See Section~\ref{sec:experiments} for more details on the experiment setup.

\paragraph{An empirical fix:} Inspired by the connection to tensor decomposition objective described earlier in the subsection, we can design a new objective function that takes exactly the same form as the tensor decomposition objective function $f_2+f_4$. Concretely, let's define
\begin{align}
\hat{y}' =  a^\top \gamma(Bx)
\end{align}
where $\gamma = \hat{\sigma_2}h_2 + \hat{\sigma}_4 h_4$  and $h_2(t) = \frac{1}{\sqrt{2}}(t^2-1)$ and $h_4(t) = \frac{1}{\sqrt{24}}(t^4-6t^2+3)$ are the 2nd and 4th normalized probabilists' Hermite polynomials~\cite{wiki-hermite}. We abuse the notation slightly by using the same notation to denote the its element-wise application on a vector. Now for each example we use $\norm{\hat{y}'-y}^2$ as loss function. The corresponding population risk is
\begin{align}
f'(a,B) = \mathop{\mathbb E}\displaylimits\left[\norm{\hat{y}'-y}^2\right]\,.\label{eqn:fprime}
\end{align}

Now by an extension of Theorem~\ref{thm:population_risk}, we have that the new population risk is equal to the $\hat{\sigma}_2^2f_2+\hat{\sigma}_4^2f_4$.

\begin{theorem}\label{thm:fprime}
	Let $f'$ be defined as in equation~\eqref{eqn:fprime} and $f_2$ and $f_4$ be defined in equation~\eqref{eqn:f4}. Assume $b_i,b^{\star}_i$'s are unit vectors. Then, we have \begin{align}f' =\hat{\sigma}_2^2 f_2 + \hat{\sigma}_4^2 f_4 +\textup{const}\label{eqn:fh2h4}\end{align}
\end{theorem}

\noindent It turns out stochastic gradient descent on the objective $f'(a,B)$ (with projection to the set of matrices $B$ with row norm 1) converges empirically to the ground truth $(a^{\star},B^{\star})$ or one of its equivalent permutations. (See Figure ~\ref{fig:h2h4}.) However, we don't know of any existing work for analyzing the landscape of the objective $f'$ (or $f_k$ for any $k \ge 3$). We conjecture that the landscape of $f'$ doesn't have any spurious local minimum under certain mild assumptions on $(a^{\star},B^{\star})$. Despite recent attempts on other loss functions for tensor decomposition~\cite{ge2017on}, we believe that analyzing $f'$ is technically challenging and its resolution will be potentially enlightening for the understanding landscape of loss function with permutation invariance.  See Section~\ref{sec:experiments} for more experimental results.

\subsection{Landscape design for orthogonal $B^{\star}$}

The population risk defined in equation~\eqref{eqn:fprime} --- though works empirically for randomly generated ground-truth $(a^{\star},B^{\star})$ --- doesn't have any theoretical guarantees. It's also possible that when $(a^{\star},B^{\star})$ are chosen adversarially or from a different distribution, SGD no longer converges to the ground-truth.

To solve this problem,  we design another objective function $G(\cdot)$, such that the optimizer of $G(\cdot)$ still corresponds to the ground-truth, and $G()$ has provably nice landscape --- all local minima of $G()$ are global minima.

In this subsection, for simplicity,  we work with the case when $B^{\star}$ is an orthogonal matrix and state our main result. The discussion of the general case is deferred to the end of this Section and Section~\ref{sec:general}.

We define our objective function $G(B)$ as
\begin{align}
G(B) &\triangleq \sign(\hat \sigma_4) \mathop{\mathbb E}\displaylimits\left[y \cdot \sum_{j,k\in [d], j\neq k}\phi(b_j,b_k, x)\right] - \mu \sign(\hat \sigma_4) \mathop{\mathbb E}\displaylimits\left[y\cdot\sum_{j \in [d]} \varphi(b_j,x)\right] \nonumber\\
&+ \lambda \sum_{i=1}^m(\norm{b_i}^2-1)^2 \label{eqn:GB}
\end{align}
where $\varphi(\cdot,\cdot)$ is defined as
\begin{align}
\varphi(v,x) = \frac{1}{8} \norm{v}^4 - \frac{1}{4}(v^{\top}x)^2 \norm{v}^2 +\frac{1}{24} (v^{\top}x)^4\,.\label{eqn:def-varphi}
\end{align}
and $\phi(\cdot, \cdot, \cdot)$ is defined as
\begin{align}
\phi(v,w,x)
& = \frac{1}{2} \norm{v}^2\norm{w}^2 + \inner{v,w}^2 - \frac{1}{2}\norm{w}^2 (v^\top x)^2 - \frac{1}{2}\norm{v}^2(w^{\top}x)^2 \nonumber \\
&+ 2(v^\top x)(w^\top x) v^\top w + \frac{1}{2}(v^\top x)^2(w^\top x)^2\,. \label{eqn:def-phi}\end{align}

The rationale behind of the choices of $\phi$ and $\varphi$ will only be clearer and relevant in later sections. For now, the only relevant property of them is that both are smooth functions whose derivatives are easily computable.

We remark that we can sample $G(\cdot)$ using the samples straightforwardly --- it's defined as an average of functions of examples and the parameters. We also note that only parameter $B$ appears in the loss function. We will infer the value of $a^{\star}$ using straightforward linear regression after we get the (approximately) accurate value of $B^{\star}$.

Due to technical reasons, our method only works for the case when $a^{\star}_i > 0$ for every $i$. We will assume this throughout the rest of the paper. The general case is left for future work. Let $a^\star_{\max} = \max a^{\star}_i$, $a^\star_{\min} = \min a^{\star}_i$, and $\kappa^\star = \max a^{\star}_i/\min a^{\star}_i$. Our result will depend on the value of $\kappa^\star.$ Essentially we treat $\kappa^\star$ as an absolute constant that doesn't scale in dimension. The following theorem characterizes the properties of the landscape of $G(\cdot)$.

\begin{theorem}\label{thm:main}
	Let $c$ be a sufficiently small universal constant (e.g. $c=0.01$ suffices) and suppose the activation function $\sigma$ satisfies $\hat{\sigma}_4 \neq 0$. Assume $\mu \le c/\kappa^\star$, $\lambda \ge c^{-1}a^\star_{\max}$, and $B^{\star}$ is an orthogonal matrix. The function $G(\cdot)$ defined as in equation~\eqref{eqn:GB} satisfies that
	\begin{enumerate}
		\item[1.] A matrix $B$ is a local minimum of $G$ if and only if $B$ can be written as $B = DPB^{\star}$ where $P$ is a permutation matrix and $D$ is a diagonal matrix with  $D_{ii} \in \left\{\pm 1 \pm O(\mu a^\star_{\max}/\lambda)\right\}$.\footnote{More precisely, $|D_{ii}| = \sqrt{\frac{1}{1-\mu|\hat{\sigma}_4|a^{\star}_i/(\sqrt{6}\lambda)}}$} Furthermore, this means that all local minima of $G$ are also global.
		\item[2.] Any saddle point $B$ has a strictly negative curvature in the sense that $\lambda_{\min}(\nabla^2 G(B))\ge -\tau_0$ where $\tau_0 = c\min\{\mu a^\star_{\min}/(\kappa^\star d), \lambda\}$

		\item [3.]  Suppose $B$ is an approximate local minimum in the sense that $B$ satisfies $$\norm{\nabla G(B)}\le \epsilon \textup{  and  }\lambda_{\min}(\nabla^2 G(B)) \ge - \tau_0$$
		Then $B$ can be written as $B = PDB^{\star} +EB^{\star} $ where $P$ is a permutation matrix, $D$ is a diagonal matrix satisfying the same bound as in bullet 1, and $|E|_\infty \le O(\epsilon/(\hat{\sigma}_4a^\star_{\min}))$.

		As a direct consequence, $B$ is $O_d(\epsilon)$-close to a global minimum in Euclidean distance, where $O_d(\cdot)$ hides polynomial dependency on $d$ and other parameters.
	\end{enumerate}
\end{theorem}

The theorem above implies that we can learn $B^{\star}$ (up to permutation of rows and sign-flip) if we take $\lambda$ to be sufficiently large and optimize $G(\cdot)$ using stochastic gradient descent. In this case, the diagonal matrix $D$ in bullet 1 is sufficiently close to identity (up to sign flip) and therefore a local minimum $B$ is close to $B^{\star}$ up to permutation of rows and sign flip. The sign of each $b^{\star}_i$ can be recovered easily after we recover $a$ (see Lemma~\ref{lem:recoveraintro} below.)

\sloppy Stochastic gradient descent converges to a local minimum \cite{ge2015escaping} (under the additional property as established in bullet 2 above), which is also a global minimum for the function $G(\cdot)$. We will prove the theorem in Section~\ref{sec:landscape} as a direct corollary of Theorem~\ref{thm:main-general}.
The technical bullet 2 and 3 of the theorem is to ensure that we can use stochastic gradient descent to converge to a local minimum as stated below.\footnote{In the most general setting, converging to a local minimum of a non-convex function is NP-hard.}
\begin{corollary}\label{cor:algorithm}
	In the setting of Theorem~\ref{thm:main}, we can use stochastic gradient descent to optimize function $G(\cdot)$ (with fresh samples at each iteration) and converge to an approximate global minimum $B$ that is $\epsilon$-close to a global minimum  in time $\poly(d,1/\epsilon)$.
\end{corollary}

\noindent After approximately recovering the matrix $B^{\star}$, we can also recover the coefficient $a^{\star}$ easily. Note that fixing $B$, we can fit $a$ using simply linear regression. For the ease of analysis, we analyze a slightly different algorithm. The lemma below is proved in Section~\ref{sec:linear}.

\begin{lemma}\label{lem:recoveraintro}
	Given a matrix $B$ whose rows have unit norm, and are $\delta$-close to $B^\star$ in Euclidean distance up to permutation and sign flip with $\delta \le 1/(2\kappa^\star)$. Then, we can give estimates $a,B'$ (using e.g., Algorithm~\ref{alg:recovera}) such that there exists a permutation $P$ where  $\|a-P a^\star\|_\infty \le \delta a^{\star}_{\max}$ and $B'$ is row-wise $\delta$-close to $PB^{\star}$.
\end{lemma}

The key step towards analyzing objective function $G(B)$ is the following theorem that gives an analytic formula for $G(\cdot)$.

\begin{theorem}\label{thm:main:formula} The function $G(\cdot)$ satisfies
	\begin{align}
	G(B) &= 2\sqrt{6}|\hat{\sigma}_4|\cdot \sum_{i\in [d]} a^{\star}_i \sum_{j,k\in [d], j\neq k}\inner{b^{\star}_i, b_j}^2 \inner{b^{\star}_i,b_k}^2  - \frac{|\hat{\sigma}_4|\mu }{\sqrt{6}}\sum_{i,j\in [d]} a^{\star}_i\inner{b^{\star}_i,b_j}^4\nonumber\nonumber\\
	&+ \lambda \sum_{i=1}^m(\norm{b_i}^2-1)^2 \label{eqn:GB:formula}
	\end{align}
\end{theorem}

\noindent Theorem~\ref{thm:main:formula} is proved in Section~\ref{sec:formula}. We will motivate our design choices with a brief overview in Section~\ref{sec:overview} and formally analyze the landscape of $G$ in Section~\ref{sec:landscape} (see Theorem~\ref{thm:main-general}).
\paragraph{Finite sample complexity bounds}: Extending Theorem~\ref{thm:main}, we can characterize the landscape of the empirical risk $\widehat{G}$, which implies that stochastic gradient on $\widehat{G}$ also converges approximately to the ground-truth parameters with polynomial number of samples.
\begin{theorem}\label{thm:samplecomplexity}
	In the setting of Theorem~\ref{thm:main}, suppose we use $N$ empirical samples to approximate $G$ and obtain empirical risk $\widehat{G}$. There exists a fixed polynomial $\mbox{poly}(d, 1/\epsilon)$ such that if $N \ge \mbox{poly}(d, 1/\epsilon)$, then with high probability the landscape of $\widehat{G}$ very similar properties to that of $G$.

	\sloppy Precisely, if $B$ is an approximate local minimum in the sense that
	$\lambda_{min}(\nabla^2 \widehat{G}(B)) \ge -\tau_0/2$ and $\|\nabla \widehat{G}(B)\| \le \epsilon/2$, then  $B$ can be written as $B = DPB^{\star} +EB^{\star} $ where $P$ is a permutation matrix, $D$ is a diagonal matrix and $|E|_\infty \le O(\epsilon/(\hat{\sigma}_4a^\star_{\min}))$.
\end{theorem}

All of the results above assume that $B^{\star}$ is orthogonal. Since the local minimum are preserved by linear transformation of the input space, these results can be extended to the general case when $B^{\star}$ is not orthogonal but full rank (with some additional technicality) or the case when the dimension is larger than the number of neurons ($m < d$). See Section~\ref{sec:general} for details.

\section{Landscape of Population Risk $G(\cdot)$}\label{sec:landscape}

In this section we prove Theorem~\ref{thm:main}. Since the landscape property is invariant with respect to rotations of parameters, without loss of generality we assume $B^{\star}$ is the identity matrix $I$ throughout this section. (See Section~\ref{sec:general} for a precise statement for the invariance.)  Recall that by Theorem~\ref{thm:main:formula}, the population risk $G(\cdot)$ in the case of $B^{\star} = I$ is equal to
\begin{align}
G(B) & = 2\sqrt{6} |\hat\sigma_4|\sum_{i} a^{\star} _i \sum_{j\neq k} (b_j ^\top e_i)^2 (b_k ^\top e_i)^2 \nonumber\\
& -\frac{|\hat{\sigma}_4|\mu }{\sqrt{6}} \sum_{i=1}^d a^{\star}_i \sum_{j=1}^d (b_j ^\top e_i)^4 + \lambda \sum_{j=1}^d \big ( \norm{b_j}^2 -1\big)^2\,.
\label{eq:obj}
\end{align}
In the rest of section we work with the formula above for $G(\cdot)$ instead of the original definition. In fact, for future reference, we study a more general version of the function $G$. For nonnegative vectors $\alpha,\beta$ and nonnegative number $\mu$, let $G_{\alpha,\beta,\mu}$ be defined as
\begin{align}G_{\alpha,\beta,\mu} (B) = \sum_{i=1}^d \alpha_i \sum_{j \neq k} (b_j ^\top e_i)^2 ( b_k ^\top e_i)^2- \mu \sum_{i=1}^d \beta_i \sum_{j=1}^d (b_j ^\top e_i)^4+\lambda \sum_{j=1}^d (\norm{b_j}^2-1)^2\label{eqn:def:galphabetamu}
\end{align}
Here $e_i$ denotes the $i$-th natural basis vector. We see that $G$ is sub-case of $G_{\alpha,\beta,\mu}$ and we prove the following extension of Theorem~\ref{thm:main}. Let $\alpha_{\max} = \max_i \alpha_i$ and $\alpha_{\min} = \min_i \alpha_i$.

\begin{theorem}\label{thm:main-general}
	Let $\kappa_\alpha = \alpha_{\max}/\alpha_{\min}$ and $c$ be a sufficiently small universal constant (e.g. $c=10^{-2}$ suffices).  Suppose $\mu \le c\alpha_{\min}/\beta_{\max}$ and $\lambda\ge  4\max(\mu \beta_{\max}, \alpha_{\max})$.
	Then,	the function $G_{\alpha,\beta,\mu}(B)$ defined as in equation~\eqref{eqn:def:galphabetamu} satisfies that
	\begin{enumerate}
		\item[1.] A matrix $B$ is a local minimum of $G_{\alpha,\beta,\mu}$ if and only if $B$ can be written as $B =DP$ where $P$ is a permutation matrix and $D$ is a diagonal matrix with  $D_{ii} \in \left\{\pm \sqrt{\frac{1}{1-\mu \beta_i/\lambda}}\right\}$.
		\item[2.] Any saddle point $B$ has strictly negative curvature in the sense that $\lambda_{\min}(\nabla^2 G_{\alpha,\beta,\mu}(B))\le -\tau_0$ where $\tau_0 = c\min\{\mu \beta_{\min}/(\kappa_\alpha d),\mu \beta_{\min}^2/\beta_{\max}, \lambda\}$
		\item [3.]  Suppose $B$ is an approximate local minimum in the sense that $B$ satisfies $$\norm{\nabla g(B)}\le \epsilon \textup{  and  }\lambda_{\min}(\nabla^2 g(B)) \ge - \tau_0$$
		Then $B$ can be written as $B = DP +E $ where $P$ is a permutation matrix, $D$ is a diagonal matrix with the entries satisfying 	$$ \frac{1}{1-\frac{\mu\beta_i}{\lambda}} \left(1-  \frac{18d \epsilon^2}{\beta_{\min}^2}-\frac{\epsilon}{2\lambda}\right)\le D_{ii}^2\le \frac{1}{1-\frac{\mu\beta_i}{\lambda}}\left(1+\frac{\epsilon}{2\lambda}\right)$$ and $E$ is an error matrix satisfying $$|E|_\infty \le 3\epsilon/\beta_{\min}.$$
		As a direct consequence, $B$ is $O_d(\epsilon)$-close to a global minimum in Euclidean distance, where $O_d(\cdot)$ hides polynomial dependency on $d$ and other parameters.
	\end{enumerate}
\end{theorem}

\noindent Here we recall that $|E|_\infty$ denotes the largest entries in the matrix $E$. Theorem~\ref{thm:main} follows straightforwardly from Theorem~\ref{thm:main-general} by setting $\alpha=2\sqrt{6}|\hat{\sigma}_4| a^{\star}$ and $\beta = |\hat{\sigma}_4| a^{\star}/\sqrt{6}$. In the rest of the section we prove Theorem~\ref{thm:main-general}.

Note that our variable $B$ is a matrix of dimension $d\times d$ and we use $b_i$ to denote the rows of $B$, that is, $B = \begin{bmatrix}b_1^\top \\ \vdots \\ b_d^\top \end{bmatrix}$.  Naturally, towards analyzing the properties of a local minimum $B$, the first step is that we pick a row $b_s$  of $B$ and treat only $b_s$ as variables and others rows as fixed. We will show that local optimality of $b_s$ will imply that $b_s$ is equal to one of the basis vector $e_j$ up to some scaling factor. This step is done in Section~\ref{sec:landscape:single-column}. Then in Section~\ref{sec:landscape:multi-columns} we show that the local optimality of all the variables in $B$ implies that each of the rows of $B$ corresponds to different basis vector, which implies that $B$ is a permutation matrix (up to scaling of the rows).
\subsection{Step 1: Analysis of Local Optimality of a Single Row} \label{sec:landscape:single-column}
Suppose we fix $b_1,\cdots, b_{s-1}, b_{s+1}, \cdots, b_d$, and optimize only over $b_s$, we obtain the objective $h$ of the following form:
\begin{align}
h_{\alpha,\beta,\lambda}(x) =\sum_{i=1}^d \alpha_i x_i ^2 - \sum_{i=1}^d \beta_i x_i ^4 + \lambda \big( \norm{x}^2-1\big)^2 \label{eqn:def:h}
\end{align}
We can see that setting $\alpha_i = a^{\star}_i \sum_{k \neq s} (b_k ^\top e_i)^2,\, \beta_i = a^{\star}_i, \text{ and } x= b_s$ gives us the original objective $G(B)$. In this subsection, we will work with $h(\cdot)$ and analyze the properties of the local minima of $h(\cdot)$.

The following lemma shows that a local minimum $x$ of the objective $h(\cdot)$ must be a scaling of a basis vector. Recall that $\sl{x}$ denotes the second largest absolute value of the entries of $x$. The lemma deals generally an approximate local minimum, though we suggest casual readers simply think of $\epsilon, \tau =0$ in the lemma.
\begin{lemma}\label{lem:second-largest-and-norm}
Let $h(\cdot)$ be defined in equation~\eqref{eqn:def:h} with non-negative vectors $\alpha$ and $\beta$ in $\mathbb R^d$.  Suppose parameters $\epsilon,\tau\ge 0$ satisfy that $\epsilon\le \sqrt{\tau^3/\beta_{\min}}$. If some point $x$ satisfies $\norm{\nabla h(x)} \le \epsilon$ and $\lambda_{\min} ( \nabla^2 h(x)) \ge - \tau $, then we have
$$
\sl{x} \le \sqrt{\frac{\tau}{\beta_{\min}}}.
$$

\end{lemma}

\begin{proof}Without loss of generality, we can take $\epsilon = \sqrt{\tau^3/\beta_{\min}}$ which means $\tau = \epsilon^{2/3}\beta_{\min}^{1/3}$.
The gradient and Hessian of function $h(\cdot)$ are
\begin{align}
\nabla h(x) &= 2 \textup{diag}(\alpha) x - 4 \textup{diag}(\beta) x^{\odot 3} +\gamma x\nonumber\\
\nabla ^2 h(x) &= 2 \textup{diag}(\alpha) - 12 \textup{diag}( \beta \odot x^{\odot 2}) + \gamma I + 8 \lambda xx^\top.\label{eqn:hessian}
\end{align}
where $\gamma \triangleq 4\lambda (\norm{x}^2-1)$.

Let $S=\{i: |x_i | \ge \delta\}$ be the indices of the coordinates that are significantly away from zero, where $\delta = \left( \frac{\epsilon}{\beta_{\min}} \right)^{1/3}$. Since $\norm{\nabla h(x)}\le \epsilon$, we have that $|\nabla h(x)_i| \le \epsilon$ for every $i\in [d]$, which implies that
\begin{align}
\forall i \in [d], \left|2 \alpha_i x_i + \gamma x_i - 4\beta_i x_i ^3\right|\le \epsilon\label{eqn:4}
\end{align}
which further implies that
\begin{align}
\forall i\in S, \left|2 \alpha_i +\gamma  -4 \beta_i x_i ^2\right| \le  \frac{\epsilon}{\delta}\label{eqn:3}
\end{align}
If $|S| = 1$, then we are done because $\sl{x} \le \delta$. Next we prove that $|S|\ge 2$. For the sake of contradiction, we
assume that $|S| \ge 2$. Moreover, WLOG, we assume that $|x|_1\ge |x|_2 $ are the two largest entries of $|x|$ in absolute values.

We take $v\in \mathbb R^d$ such that $v_{1} = -x_2/\sqrt{x_1 ^2 +x_2 ^2}$, and $v_2 = x_1/\sqrt{x_1 ^2 +x_2 ^2}$, and $v_j = 0$ for $j \ge 2$. Then we have that $v^\top x=0$ and $\norm{v} = 1$. We evaluate the quadratic form and have that
\begin{align}
v^\top \nabla^2 h(x) v &= v^\top (2 \textup{diag}(\alpha) +\gamma I_{d} )v -12 v^\top \textup{diag}(\beta \odot x^{\odot 2}) v \tag{since $v^\top x = 0$}\\
&= (2\alpha_1 + \gamma) v_1 ^2+ (2\alpha_2 + \gamma) v_2 ^2 - 12\beta_1 v_1 ^2 x_1 ^2 -12 \beta_2 v_2 ^2 x_2 ^2\nonumber\\
&\le- 8\beta_1 v_1 ^2 x_1 ^2 -8 \beta_2 v_2 ^2 x_2 ^2+ \frac{\epsilon}{\delta} \tag{by equation~\eqref{eqn:3} and $\norm{v}=1$}\\
&\le -8 (\beta_1+\beta_2)  \frac{x_2 ^2x_1 ^2 }{x_1 ^2 +x_2^2}+ \frac{\epsilon}{\delta}\nonumber\\
&\le - 8 \beta_{\min} x_2 ^2 + \frac{\epsilon}{\delta}\,.\nonumber
\end{align}
Recall that $\delta = \left( \frac{\epsilon}{\beta_{\min}} \right)^{1/3}$. Then we conclude that $$v^\top \nabla^2 h(x) v \le - 6 \beta_{\min}^{1/3} \epsilon^{2/3}= -6\tau.$$ This contradicts with the assumption that $\lambda_{\min} (\nabla^2 h(x)) \ge  -\beta_{\min}^{1/3} \epsilon^{2/3} = \tau$ and that $\norm{v} = 1$. Therefore we have $|S|=1$ and
\begin{align*}
\sl{x} \le \delta = \left( \frac{\epsilon}{\beta_{\min}}\right)^{1/3}  \le \sqrt{\frac{\tau}{\beta_{\min}}} \tag{ using $\epsilon \le \sqrt{\tau^3/\beta_{\min}}$}
\end{align*}
\end{proof}

\noindent For future reference, we can also show that for a sufficiently strong regularization term (sufficiently large $\lambda$),  the norm of a local minimum $x$ should be bounded from below and above by $1/2$ and $2$. This are rather coarse bounds that suffice for our purpose in this subsection. In Section~\ref{sec:landscape:multi-columns} we will show that all the rows of a local minimum $B$ of $G$ have norm close to 1.
\begin{lemma}\label{lem:norm}
	In the setting of Lemma~\ref{lem:second-largest-and-norm},
\begin{enumerate}
	\item Suppose in addition that  $\lambda\ge 4\max(\beta_{\max}, \tau)$ and  $\epsilon \le 0.1\beta_{\min} d^{-3/2}$, then   $$\norm{x}^2 \le 2\,.$$

\item  Let $i^\star = \arg\max_{i} |x_i|$. In addition to the previous conditions in bullet 1, assume that $\lambda\ge 4\alpha_{i^\star }$. Then, $$\norm{x}^2 \ge \frac12\,.$$
\end{enumerate}
\end{lemma}

We remark that we have to state the conditions for the upperbounds and lowerbounds separately since they will be used with these different conditions.

\begin{proof}
Let $S=\{i: |x_i | \ge \delta\}$ be the indices of the coordinates that are significantly away from zero, where $\delta = \left( \frac{\epsilon}{\beta_{\min}} \right)^{1/3}$.	We first show that $\norm{x}^2\le 2$. We divide into two cases:

\begin{enumerate}
	\item[1] $S$ is empty. Since $\epsilon \le 0.1\beta_{\min} d^{-3/2}$, then  $\delta \le \frac{\sqrt{2}}{\sqrt{d}}$. We conclude that $\norm{x}^2\le 2$.
	\item[2] $S$ is non-empty. For $i \in S$, recall equation~\eqref{eqn:3} which implies that
	\begin{align*}
	4\lambda( \norm{x}^2-1)&\le \frac{\epsilon}{\delta}+4\beta_i x_i ^2 \\
	&\le  \frac{\epsilon}{\delta} +4 \beta_{\max} \norm{x}^2 \nonumber\\
	\norm{x}^2& \le  \frac{\epsilon}{4\lambda \delta} +\frac{ \beta_{\max}}{ \lambda} \norm{x}^2+1
	\end{align*}

	Since $\lambda\ge 4\beta_{\max}$, so $\lambda \ge \beta_{\min}^{1/3} \epsilon^{2/3}\ge \frac{3\epsilon}{4\delta}$, and thus from the display above  we have  that $\norm{x}^2\le 2$.
\end{enumerate}

Next we show that $\norm{x}^2 \ge \frac12$. Again we divide into two cases:
\begin{enumerate}
\item $S$ is empty.
For the sake of contradiction, assume that $\norm{x}^2 \le \frac12$, then $\gamma \le -2 \lambda$. We show that there is sufficient negative curvature. Recall that
\begin{align*}
\nabla^2 h(x) &= 2 \textup{diag} (\alpha) -12 \textup{diag}(\beta\odot x^{\odot 2}) +\gamma I + 8 \lambda xx^\top\\
&\preceq  2 \textup{diag} (\alpha) -12 \textup{diag}(\beta\odot x^{\odot 2})-2\lambda I + 8 \lambda xx^\top
\end{align*}
Choose index $j^\star$ so that $\alpha_{j ^\star} =\alpha_{\min}$, then
\begin{align*}
e_{j^\star} ^\top \nabla^2  h(x) e_{j^\star} &= 2 \alpha_{\min} - 12 \beta_{j^\star} x_{j^\star} ^2 - 2\lambda +8 \lambda x_{j^\star}^2\\
&\le 2 \alpha_{\min} + 8\lambda \delta^2 - 2\l\\
&\le 2\alpha_{\min} - \lambda ( 2- 8 \delta^2)\\&
 \le 2\alpha_{\min} -\frac43 \lambda  \tag{by  $\delta^2 \le \frac{1}{12}$}\\
 &\le -\frac{5}{6}\lambda \le -3\tau \tag{ by $ \lambda \ge 4\max\{\alpha_{\min},\tau\}$}
\end{align*}
This contradicts with the fact that $\lambda_{\min} (\nabla^2 h(x)) \ge -\tau$. Thus when $S$ is empty, $\norm{x}^2 \ge \frac12$.

\item $S$ is non-empty. Recall that $i^\star =\arg\max_i |x_i|$, and by definition $i^\star \in S$. Using Equation \eqref{eqn:3}
\begin{align*}
\gamma &\ge - 2 \alpha_{i^\star} - \frac{\epsilon}{\delta}\\
\end{align*}
which implies that
\begin{align*}
\norm{x}^2& \ge 1 - \frac{\alpha_{i^\star}}{\lambda}- \frac{\epsilon}{4 \lambda \delta}.
\end{align*}
Since $\lambda\ge 4 \alpha_{i^\star}$, and $ \lambda \ge \beta_{\min}^{1/3} \epsilon^{2/3}\ge \frac{\epsilon}{\delta}$, we conclude that $\norm{x}^2 \ge 1/2$.
\end{enumerate}
\end{proof}
We have shown that a local minimum $x$ of $h$ should be a scaling of the basis vector $e_{i^\star}$. The following lemma strengthens the result by demonstrating that not all basis vector can be a local minimum --- the corresponding coefficient $\alpha_{i^\star }$ has to be reasonably small for $e_{i^\star}$ being a local minimum. The key intuition here is that if $\alpha_{i^\star}$ is very large compared to other entries of $\alpha$, then if we move locally the mass  of $e_{i^\star}$ from entry $i^{\star}$ to some other index $j$, the objective function will be likely to decrease because $\alpha_jx_j^2$ is likely to be smaller than $\alpha_{i^\star}x_{i^\star}^2$. (Indeed, we will show that such movement will cause a second-order decrease of the objective function in the proof.)
\begin{lemma}
	In the setting of Lemma~\ref{lem:second-largest-and-norm}, let $i^\star = \arg\max_{i} |x_i|$. If $\norm{\nabla h(x)} \le \epsilon$, and $\lambda_{\min}(\nabla^2 h(x)) > -\tau$ for $0\le \tau\le 0.1\beta_{\min}/d$ and $\epsilon \le \sqrt{\tau^3/\beta_{\min}}$, then
\begin{align}
\alpha_{i^\star} \le   \alpha_{\min} +2\epsilon+2\tau+4\beta_{i^\star}\,. \nonumber\end{align}
\label{lem:subproblem-ub-alpha}
\end{lemma}
\begin{proof}
For the ease of notation, assume WLOG that $i^\star = 1$. Let $\delta = (\tau/\beta_{\min})^{1/2}$. By the assumptions, we have that $\delta \le \frac{1}{\sqrt{6d}}$. By Lemma \ref{lem:second-largest-and-norm}, we have $\norm{x}^2 \ge \frac12$, which implies that
\begin{align}
x_1 ^2 \ge \norm{x}^2 - (d-1) \sl{x}^2\ge \frac12 -d \sl{x}^2 \ge 1-d\delta^2 \ge \frac{1}{3}\,.\label{eqn:x1}
\end{align}

Define $v=-\left(\frac{x_k}{x_1}\right) e_1 + e_k$.  Since $x_1$ is the largest entry of $x$, we can verify that
$
1\le \norm{v} ^2 = 1+ \frac{x_k^2}{x_1 ^2}\le 2
$.  By the assumption, we have that
\begin{align}v^\top \nabla ^2 h(x) v \ge - \tau \norm{v}^2 \ge -4 \tau\,.\label{eqn:5}\end{align}

On the other hand, recall the form of Hessian (equation~\eqref{eqn:hessian}), by straightforward algebraic manipulation, we have that
\begin{align}
v^\top \nabla^2 h(x) v &= v^\top \left(2 \textup{diag}(\alpha) - 12 \textup{diag}( \beta \odot x^{\odot 2}) + \gamma I + 8 \lambda xx^\top.\right) v\nonumber\\
& =2\alpha_1 (\frac{x_k }{x_1})^2 + 2\alpha_k-12 \beta_1 x_k ^2 - 12 \beta_k x_k ^2  + \gamma (\frac{x_k}{x_1})^2  +\gamma\tag{by $v^\top x = 0$}\\
&\le (2\alpha_1+\gamma) (\frac{x_k }{x_1})^2-12 (\beta_1+\beta_k) x_k ^2  + 2\alpha_k +(4\beta_1 x_1 ^2 -2 \alpha_1 +\frac{\epsilon}{|x_1|})\tag{by equation~\eqref{eqn:3}}\\
&\le(4 \beta_1 x_1^2 +\frac{\epsilon}{|x_1|} ) (\frac{x_k }{x_1})^2-12 (\beta_1+\beta_k) x_k ^2  + 2\alpha_k  +(4\beta_1 x_1 ^2 -2 \alpha_1 +\frac{\epsilon}{|x_1|})
\tag{by equation~\eqref{eqn:3}}\\
& = -8 \beta_1 x_k ^2   - 12 \beta_k x_k ^2 +4\beta_1 x_1 ^2 + 2\alpha_k-2 \alpha_1+ 4\epsilon \tag{by $|x_k|\le |x_1|$ and $|x_1|^2\ge 1/3$}\\
& \le   2\alpha_k-2 \alpha_1+ 4\epsilon + 8\beta_1\,. \tag{by $\norm{x}\le 2$ using Lemma~\ref{lem:second-largest-and-norm}}
\end{align}

Combining equation~\eqref{eqn:5} and the equation above gives
\begin{align*}
\alpha_1
&\le \alpha_k + 2 \epsilon + 2\tau +4\beta_1 \,.
\end{align*}

Since $k$ is arbitrary we complete the proof.
\end{proof}

The previous lemma implies that it's very likely that the local minimum $x$ can be written as $x =  x_{i^\star}e_{i^\star}$ and the index $i^\star$ is also likely to be the argmin of $\alpha$. The following technical lemma shows that when this indeed happens, then we can strengthen Lemma~\ref{lem:second-largest-and-norm} in terms of the error bound's dependency on $\epsilon$ and $\tau$. In Lemma~\ref{lem:second-largest-and-norm}, we have that $\sl{x}$ is bounded by a function of $\tau$. Here we strengthen the bound to be a function that only depends on $\epsilon$. Thus as long as $\tau$ be small enough so that we can apply Lemma~\ref{lem:second-largest-and-norm} and Lemma~\ref{lem:subproblem-ub-alpha} to  meet the condition of the lemma below, then we get an error bound that goes to zero as $\epsilon$ goes to zero. This translates to the error bound in bullet 3 of Theorem~\ref{thm:main-general} where the bound on $E$ only depends on $\epsilon$. For casual readers we suggest to skip this Lemma since its precise functionality will only be clearer in the proof of Theorem~\ref{thm:main-general}.
\begin{lemma}
In the setting of Lemma~\ref{lem:second-largest-and-norm}, in addition we assume that  $i = \textup{argmin}_k |\alpha_k|$ and that $x$ can be written as $x= x_{i} e_{i}+x_{-i}$ satisfying
\begin{align}
\norm{x_{-i}}_{\infty} \le 0.1\min\{ 1/\sqrt{d}, \sqrt{\beta_{\min}/(\beta_{\max})}\}\nonumber\,.
\end{align}
Then, we can strengthen the bound to
\begin{align}
\norm{x_{-i}}_\infty \le \frac{3\epsilon}{\beta_{\min}}.\nonumber
\end{align}
\label{lem:exact-recovery}
\end{lemma}
\begin{proof}
WLOG, let $i=1$. Let $x_j $ be the second largest entry of $x$ in absolute value. Define $v_1= 4\beta_1 x_1 ^2 -2 \alpha_1 -\gamma $, and similarly $v_j = 4\beta_j x_j ^2 - 2\alpha_j -\gamma$. Since $\norm{\nabla h(x) } \le \epsilon$, by equation~\eqref{eqn:3}, we have that $|v_1| \le  \frac{\epsilon}{|x_1|}$ and $|v_2| = \frac{\epsilon}{|x_j|}$. Subtracting $ 4 \beta_1 x_1 ^2 = 2\alpha_1 +\gamma +v_1$ and $ 4 \beta_j x_j ^2 = 2\alpha_j +\gamma+v_j $, we obtain,
\begin{align}
4 \beta_1 x_1 ^2 &= 4 \beta_j x_j ^2 -2 (\alpha_j -\alpha_1) +v_1-v_j \nonumber\\
&\le 4 \beta_j x_j ^2 +(v_1-v_j)\tag{since $\alpha_j -\alpha_1\ge0$}
\end{align}

Since $ \norm{x}^2 \ge \frac12$, then $x_1 ^2 \ge \frac12 - d\delta^2 \ge \frac13$. Since $|x_j| \le \delta$,
\begin{align*}
4 \beta_j x_j ^2 \le 4 \beta_{\max}\delta^2
\end{align*}
Combining the above two displays,
\begin{align}
(v_1-v_j) &\ge 4 \beta_1 x_1 ^2 - 4 \beta_{j} \delta ^2
\nonumber\\
&\ge 4 \beta_1 x_1 ^2 - 4 \beta_{\max} \delta^2\nonumber\\
&\ge \frac43 \beta_1 - \frac23 \beta_{\min} \tag{ using $x_1 ^2 \ge \frac13$ and $ \delta \le \sqrt{\beta_{\min}/(6\beta_{\max})}$}\nonumber\\
&\ge \frac23 \beta_{\min}
\end{align}
Since $|v_1| \le \frac{\epsilon}{|x_1|} $ and $ |v_2| = \frac{\epsilon}{|x_2|}$,
\begin{align}
2 \frac{\epsilon}{|x_j|} \ge \frac23 \beta_{\min},
\end{align}
and re-arranging gives $|x_j| \le 3\frac{\epsilon}{\beta_{\min}} $.
\end{proof}

redpajama_set_name: RedPajamaArXiv