\section{Introduction}

\iffalse
It is known that most neural networks that perform well in classification tasks are overparameterized, i.e., the number of parameters in the neural network exceeds the number of data points in the training set used to train the network parameters.
In fact, it has even been shown experimentally that these neural networks can memorize random labels \cite{zhang2016understanding} but nevertheless they generalize well when trained with real data.
Therefore, it is of great interest to understand why overparameterized neural networks work well in practice despite the fact that conventional wisdom learned from a direct application of Rademacher complexity results in statistical learning theory \cite{mohri} tells us that overparameterization is typically undesirable.
\ruoyu{This paragraph is very similar to those papers
with implicit regularization; one naturally expect
to explain this phenomenon by implicit regularization.
But the current paper is NOT that ambitious, and does NOT intend to explain the phenomenon. If so, then the beginning tone is somewhat misleading?
I tend to rewrite this paragraph,
but leave it for now: before rewriting, want
to get your feedback. }
\fi

Machine learning practitioners are ultimately interested in test error. Nevertheless, due to the difficulty of directly analyzing test error, theoreticians often decompose the test error into two parts: training error and generalization gap. In the past few years, a lot of research has been devoted to a single aspect of the problem, either analyzing
 the training loss or analyzing the generalization gap.
Using optimization theory, a number of works
\cite{kawaguchi2016deep,freeman2016topology,choromanska2015loss,nguyen2017loss,nguyen2018loss,liang2018understanding, liang2018adding, liang2019revisiting}
have analyzed local minima in deep neural networks,
or prove converge to global minima of the training loss,
but do not provide bounds on generalization error.
Using statistical learning theory, \cite{bartlett2017spectrally,neyshabur2018towards}
analyzed the generalization gap, but did not
 show how to achieve small training error.
There are some recent works that intend to prove small test error for neural nets.
Most works explore the NTK-type analysis
 to prove a generalization bound (e.g.,
 \cite{arora2019fine,cao2019generalization}), but these works require a very large width, say, at least $\Omega(n^2)$. References
 \cite{ji2019, chen2019much} showed that the requirement for the width can be small if the data satisfies certain assumptions
 (i.e., the margin is large enough), but the width can still be large for other  data distributions.
One notable exception is
 \cite{allen2018learning}, which used a different initialization (scaling by $1/m$ instead of $1/\sqrt{m}$) and thus allow the trajectory to go beyond the
 kernel regime.
 However, their width is still at least $\Omega(n^r)$ (this is not explicitly stated; see Section \ref{subsec: related work} for detailed discussion).

   In this paper, we are interested in generic data distributions and mildly overparamterized neural nets with width in the order of $n$, where $n$ is the number of samples. The question we explore is the following: is there a polynomial time
  algorithm to find a point with small test error,
  for a mildly overparameterized neural nets?

We add an explicit $\ell_2$ regularization term to the loss function and first study the landscape of the resulting regularized empirical loss function.
We show that all local minima of this loss function, and other directional local minima (i.e., those points which are local minima only in certain specific directions), have good memorization and generalization performance.
In the course of identifying these properties, we also show that all of these points have a key common feature: at least one neuron is inactive at these local stationary points if the number of neurons is larger than the number of points in the dataset.
We exploit this property to design a
variant of gradient descent algorithm which finds a good set of parameters for the neural network. Specifically, we show that when the algorithm gets stuck at a directional local minimum which does not have good generalization performance, by perturbing the parameters of the inactive neuron, one can find a descent direction to reduce the value of the loss function.
We prove that for convolutional neural net,
the algorithm runs in polynomial time in the number of data points and the input dimension.
For a fully connected neural network, under an additional
assumption on the data distribution, we provide
another polynomial time algorithm that finds a point with
small test error.

\subsection{Relationship to prior work}\label{subsec: related work}

We now relate our work to prior work in a few different categories:
\begin{itemize}

 \item \emph{Neural Tangent Kernel and Related Results:}
 By exploring the connection between gradient descent on versions of the neural network linearized around the initialization and kernel regression on a kernel called the Neural Tangent Kernel (NTK),
 researchers proved the convergence of gradient descent to global minima
 for ultra-wide networks
 \cite{jacot2018neural,du18,arora2019fine}.
 The width we require in this work is in the order of $O(n)$,
 which is much smaller than the requirement of the NTK-type works.

 \item \emph{Optimization and Generalization}.
As mentioned earlier, a number of NTK-type results
provide small test error for neural nets.
Following the convention in NTK theory, most works (e.g.
\cite{arora2019fine,cao2019generalization})
require a very large width, say, at least $\Omega(n^2)$.
Recently,
 \cite{ji2019}\cite{chen2019much}
 showed that the width can be smaller than $O(n)$
 if the dataset is separable
 with a large enough margin.
More specifically, the number of parameters required is inversely proportional to the margin by which the binary labelled dataset can be separated,
thus for certain dataset the width can be as small
as $\text{poly}(\log (n))$.
However, if the dataset contains points that are very close to each other with different labels or if the dataset contains very small amount of mislabeling errors, then the number of parameters required could be very large or infinity.
One of our interests is in obtaining results for datasets which are not separable due to the fact, in practice, that a small amount of labeling errors can render a dataset non-separable.

A notable different work is
\cite{allen2018learning}, which did not utilize NTK-type analysis and proved small test error by running gradient descent.
However, their required width is still large (at least
$\Omega(n^4)$ or $\Omega(1/\epsilon^8)$ where $\epsilon$
is the desired test error).
Note that their theorems did not explicitly state the dependence of the width $m$ on $n$ or $1/\epsilon$,
and we derive the dependence as follows.
In their proof of Lemma B.4 in page 40,
the second and third equation after equation (B.4)
together enforce an requirement on $m$:
based on the discussions there, to let the first term
in the second equation be smaller than
$O(\epsilon)$, we need
$ \frac{ \sqrt{m} }{ \epsilon^2 } \frac{1}{\epsilon m} < \epsilon $,
 i.e., $m > \frac{1}{ \epsilon^8 }$.
The relation between $N$ and $\epsilon$
can be found is discussed
in their Remark B.7: they
are aware of a proof that can show $N$
 scales as $1/\epsilon^2$~\footnote{ Based on their Appendix B.5, technically their proof requires $ n > 1/ \epsilon^4$, thus the test error can only scale as $1/n^{1/4}$ which is worse than the typical test error $1/n^{1/2}$.}.
By using this desired scaling, their width $m$ would be at least $ \Omega( n^4)$.
Another difference is that their data are generated by a target network, while our data are generic.

\item \emph{Escaping saddle points}.
 There are some works which
 shed light on how easy or difficult it is to escape saddle points and reach a local or global minimum \cite{laurent2017multilinear,tian2017analytical,soltanolkotabi2019theoretical,mei2018mean,brutzkus2017globally,zhong2017recovery,li2017convergence,brutzkus2017sgd,wang2018learning,du2018power,oymak2019towards,janzamin2015beating, mondelli2018connection}.
Our work also handle certain points that
a variant of GD may get stuck at, and we provide a problem-dependent algorithm to escape such points.
Note that these points may not be the saddle point
of the original loss function, since we consider
a special variant of GD that only operates
on a subset of parameters.

   \iffalse
   \item \textit{Optmization in over-paramterized networks.}
   The dynamics of gradient descent performed on this linearized model in the limit when the number of neurons goes to infinity leads to the neural-tangent kernel, which is one of the dominant approximations to neural networks studied in the overparameterized regime in recent literature. While the early results were for squared loss functions, related results are also available for logistic loss which is more suitable for binary classification \cite{ziweipolylog,chen2019much}.
   The use of the logistic loss function leads to far fewer neurons in the neural network to achieve good generalization performance. In fact, in the case of logistic loss functions, theory seems to suggest that overparameterization is not required. However, it should be noted that the number of parameters required is inversely proportional to the margin by which the binary labelled dataset can be separated. Thus, if the dataset contains points that are very close to each other with different labels or if the dataset contains very small amount of mislabeling errors, then the number of parameters required could be very large or infinity.
   \fi

    \item \emph{Explicit Regularization}: In \cite{wei2019regularization}, the authors show that explicit regularization can decrease the number of samples required to learn a distribution. In \cite{wei2020implicit}, the dual role of dropout as both an implicit and explicit regularization is discussed. In \cite{liang2019revisiting},
    an regularizer is added to eliminate the decreasing path to infinity for ReQU (rectified quadratic unit) network. These works did not
    directly provide a bound on the test error.

    \item \emph{Computational Complexity Results:} There has been a long line of work pointing out the difficulty of training neural networks, we refer the readers to small sample of such works in \cite{goel2020superpolynomial,vu1998infeasibility,blum1989training}. We do not directly address this issue in our paper, as our polynomial-time algorithm works only for two cases: convolutional neural networks where each neuron only takes as input a subset of the input vector to the neural network,
    or fully connected neural networks under extra
     assumption on the data.

    \item \emph{Implicit Regularization}.
    One approach to understanding the performance of overparameterized networks is to show
    that the iterates generated by gradient descent exhibit some implicit bias.

   For instance, for linear regression, gradient descent
   converges to a minimum norm solution for overparamaterized  models; for matrix factorization,
   gradient descent converges to a low-rank solution
   for overparameterized models \cite{gunasekar2018implicit,LiEtAl2017Algorithmic}.
 Researchers have also studied implicit
 bias in deep neural networks.
 \cite{arora2019implicit,ji2018gradient}
 analyzed deep linear networks;
 for instance, \cite{ji2018gradient}
 showed that the weight matrices will converge
 to low-rank matrices for scalar output case.
    A few recent works
    \cite{ji2020directional,lyu2019gradient}
  analyzed deep non-linear homogeneous networks,
  but they assume the loss to be smaller than a certain threshold. They did not show that an algorithm can find a point below
  the required threshold, thus their results can be viewed
  as local analysis (near the  global minima).

\end{itemize}

\section{Preliminaries}\label{sec::prelim}

\textbf{Notation for single layered network with ReLU activation.} In this paper, we consider both fully-connected neural network (FNN) and convolutional neural network (CNN) with rectified linear units (ReLUs). For simplicity of notation, we use the following way the express the outputs of both FNNs and CNNs. Given an input vector of dimension $d$, let $\{\phi_k\}_{k\in \mathbb{N}}$ be a fixed vector series where each vector $\phi_k$ is a binary vector of dimension $d$, i.e., $\phi_k\in\{0,1\}^d$.
We use $\odot$ to denote the Hadamard product of two vectors. Specifically, for a vector $\phi\in\{0, 1\}^d$ and an vector $x\in\mathbb{R}^d$, the Hadamard product  $x\odot \phi $ is a vector of dimension $d$ where the $k$-th coordinates of the vector $x\odot \phi$ is the product of the $k$-th coordicates of vectors $x$ and $\phi$, i.e.,
$$(x\odot \phi )(k)=x(k)\phi(k)= \left\{\begin{matrix}x(k) &\text{if }\phi(k)=1,\\ 0&\text{if }\phi(k)=0.\end{matrix}\right.$$
Given a series $\{\phi_k\}_{k\in \mathbb{N}}$, we define the output of neural network $f$ as
$$f(x;\bm{\theta})=\sum_{j=1}^ma_j(\bm{w}_j^\top (x\odot \phi_j))_+,$$
where $m$ denotes the number of neurons in the neural network and the vector $\bm{\theta}$ consists of all parameters (i.e., $a_j$s and $\bm{w}_j$s) in the neural network. Next, we will show how to choose the series $\{\phi_k\}$ to express the outputs of FNN and CNN respectively.

\textbf{FNN.} To represent the output of a fully-connected neural network, we can choose $\phi_k=\bm{1}_d$ for all $k\in\mathbb{N}$. This indicates $x\odot \phi_k=x$ for all $x\in\mathbb{R}^d$ and $k\in\mathbb{N}$ and thus indicates that the output of $f$ becomes
$$f(x;\bm{\theta})=\sum_{j=1}^ma_j(\bm{w}_j^\top (x\odot \phi_j))_+=\sum_{j=1}^ma_j(\bm{w}_j^\top x)_+,$$
which is exactly the same as the notation of the output of a fully-connected neural network. We also note here that, in this case, the series $\{\phi_k\}_k$ is a periodic series with a period of one.

\textbf{CNN.}
To represent the output of a convolutional neural network with filters of size $r\in[1, d]$, we can choose a periodic series $\{\phi_k\}_{k\in\mathbb{N}}$ satisfying that for each integer $k\in[1, d-r+1]$,
$$\phi_k(j)=\left\{\begin{matrix}1 &\text{if } k\le j\le k+r-1,\\ 0&\text{otherwise,}\end{matrix}\right.$$
 and that the periodic vector series $\{\phi_k\}_{k\in\mathbb{N}}$ has a period of $d-r+1$, i.e., $\phi_k=\phi_{k+d-r+1}$, for any $k\in\mathbb{N}$. Now it is straight-forward to see that the output of $f$ denotes a convolutional neural network with filters of size $r$.

\textbf{Loss and error.} We use $\mathcal{D}=\{(x_{i},y_{i})\}_{i=1}^{n}$ to denote a dataset containing $n$ samples, where $x_{i}\in\mathbb{R}^{d}$ and $y_{i}\in\{-1,1\}$ denote the feature vector and the label of the $i$-th sample, respectively. Given a neural network $f(x;\bm{\theta})$ parameterized by $\bm{\theta}$ and a univariate loss function $\ell:\mathbb{R}\rightarrow\mathbb{R}$, in binary classification tasks, we define the regularized  empirical loss $L_{n}(\bm{\theta};\bm{\lambda})$ as a linear combination of a regularizer $V(\bm{\theta};\bm{\lambda})$ parameterized by a vector $\bm{\lambda}$ and the average loss of the network $f$ on a sample in the dataset. We define the training error (also called the misclassification error) $R_{n}(\bm{\theta};f)$ as the misclassification rate of the network $f$ on the dataset $D$, i.e.,
\begin{equation}
L_{n}(\bm{\theta};\bm{\lambda})=\sum_{i=1}^{n}\ell(-y_{i}f(x_{i};\bm{\theta}))+V(\bm{\theta};\bm{\lambda})
\end{equation}
and
\begin{equation} R_{n}(\bm{\theta};f)=\frac{1}{n}\sum_{i=1}^{n}\mathbb{I}\{y_{i}\neq \text{sgn}(f(x_{i};\bm{\theta}))\},
\end{equation}
where $\mathbb{I}$ is the indicator function. Given an underlying distribution $\mathbb{P}_{\bm{X}\times Y}$, we define the test error $R(\bm{\theta}; f)$ of a neural network as the misclassification error of the neural network on the underlying distribution, i.e.,
$$R(\bm{\theta}; f)=\mathbb{P}_{\bm{X}\times Y}(Y\neq \text{sgn}(f(X;\bm{\theta}))).$$

\section{Assumptions}\label{sec::main-results}

In this section, we introduce several assumptions on the univariate loss function and dataset.

\begin{assumption}[Loss function]\label{assumption::loss-1}
 Assume that the univariate loss function $\ell$ is convex, non-decreasing and twice differentiable. Assume that both function $\ell$ and its derivative $\ell$ are $1$-Lipschitiz and that there exists a positive real number $a\in\mathbb{R}^+$ such that $\ell'(z)\le e^{az}$ holds for any $z\in\mathbb{R}$.
\end{assumption}
\textbf{Remark: } The Lipschitz constant need not be $1,$ we assume it for the simplicity of notation. In fact, if a general univariate loss $\ell$ and its derivative $\ell'$ is $L$-Lipschitz, then we can normalize the loss function by setting $\ell_{new}=\ell/L$ to satisfy the assumption. Similarly, if $\ell'(z)\le be^{az}$ for some positive $a, b\in\mathbb{R}$ and all $z\in\mathbb{R}$, then we can also normalize the loss by setting $\ell_{new}=\ell/b$ to satisfy the assumption.

Given a periodic series of vectors $\bm{\phi}=(\phi_k)_{k\ge 1}$, we define the class of functions $\mathcal{H}_{\bm{\phi}}$  as follows,
\begin{equation}\label{H def}
\mathcal{H}_{\bm{\phi}}=\bigcup_{m=1}^{\infty}\left\{h:x\mapsto \sum_{j=1}^{m}a_j(\bm{w}_j^\top (x\odot \phi_j))_+\Bigg|\sum_{j=1}^m|a_j|=1, \|\bm{w}_j\|_2=1,j=1,...,m\right\}.
\end{equation}
Now we introduce the assumption on the dataset. For the problem where we train the neural network to memorize all points in the dataset, we assume that the most of samples in the dataset can be separated by a neural network in the function class $\mathcal{H}_{\bm{\phi}}$ with a positive margin.

\begin{assumption}[Dataset]\label{assumption::dataset-1}
Assume that $\|x_i\|_2\le 1$ holds for all $i\in[n]$. For a given series of vectors $\bm{\phi}=(\phi_k)_{k\ge 1}$, assume that there exists a  number $E\in[0, n]$, a margin $\gamma\in(0,1]$ and a ReLU network $h\in\mathcal{H}_{\bm{\phi}}$ such that
$$\sum_{i=1}^n\mathbb{I}\{y_i h(x_i)\ge \gamma\}\ge n-E.$$
\end{assumption}

However, for the problem where we train the neural network to achieve good performance on the underlying distribution, we assume that with high probability, data samples drawn from the underlying distribution can be separated by a neural network in the function class $\mathcal{H}_{\bm{\phi}}$ with  a positive margin.

\begin{assumption}[Data Distribution]\label{assumption::data-distribution}
Assume $\mathbb{P}_{X\times Y}(\|X\|_2\le 1)=1$.
For a given series of vectors $\bm{\phi}=(\phi_k)_{k\ge 1}$, assume that there exists a real number $\varepsilon\in [0,0.5)$, a margin $\gamma>0$ and a ReLU network $h\in\mathcal{H}_{\bm{\phi}}$ such that $\mathbb{P}_{X\times Y}(Yh(X)\ge \gamma)=1-\varepsilon $. Assume that the samples in the dataset are independently drawn from the underlying distribution $\mathbb{P}_{X\times Y}$.

\end{assumption}

\textbf{Remark:} In fact, the above assumption covers the several cases shown below
\begin{itemize}[leftmargin=*]
    \item \textbf{Clearly separable cases.} If all data points can be separated by a neural network $h\in\mathcal{H}_{\bm{\phi}}$ with a large margin $\gamma>0$, then $\sum_{i=1}^n\mathbb{I}\{y_i h(x_i)\ge \gamma\}\ge n$ or $\mathbb{P}_{X\times Y}(Yh(X)\ge \gamma)=1$.
    \item \textbf{Almost clearly separable cases.} If most of samples drawn from the underlying distribution can be separated by a neural network $h\in\mathcal{H}$ with a large margin $\gamma>0$.
\end{itemize}

\textbf{Remark:}
Assumption~\ref{assumption::dataset-1}
is different from Assumption \ref{assumption::data-distribution}
as the latter requires independent samples
while the former does not.
We distinguish Assumption~\ref{assumption::dataset-1} and
Assumption~\ref{assumption::data-distribution} because
some of our results require
Assumption ~\ref{assumption::dataset-1}
while some of our results require
Assumption \ref{assumption::data-distribution}.
Further, it is straightforward to show that under Assumption~\ref{assumption::data-distribution}, when the number of samples is sufficiently large, Assumption~\ref{assumption::dataset-1} holds with high probability.

\section{All Local Minima Memorize and Generalize Well}\label{sec::single}
In this section, we will show that all local minima of the empirical loss have good memorization and generalization performance. Recall that for a single-layered ReLU network consisting of $m$ neurons, the output of the neural network is defined as $f(x;\bm{\theta})=\sum_{j=1}^{m}a_{j}\left(\bm{w}_{j}^{\top}(x\odot \phi_j)\right)_+$. Now we define the empirical loss as

\begin{equation}\label{eq::loss-single}
L_{n}(\bm{\theta};\bm{\lambda})=\sum_{i=1}^{n}\ell(-y_{i}f(x_{i};\bm{\theta}))+\frac{1}{2}\sum_{j=1}^{m}\lambda_{j}\left[a_{j}^{2}+\|\bm{w}_{j}\|^{2}_{2}\right],
\end{equation}
where all regularizer coefficients $\lambda_{j}$'s are positive numbers and the vector $\bm{\lambda}=(\lambda_{1},...,\lambda_{m})$ consists of all regularizer coefficients. We note that after adding the regularizer, the empirical loss $L_n$ is coercive (i.e., $L_n(\bm{\theta})\rightarrow\infty$ as $\|\bm{\theta}\|_2\rightarrow\infty$) and always has a  global minimum. Now we present the following theorem to show that if the network size is larger than the dataset size ($m\ge n+1$ for FNN; for CNN, the threshold is $(n+1)(d-r+1)$)
and the regularizer coefficient vector $\bm{\lambda}$ is carefully chosen, then every local minimum of the empirical loss ${L}_{n}(\bm{\bm{\theta}}^{*};\bm{\lambda})$ achieves zero training error on the dataset $\mathcal{D}$.

\begin{theorem}\label{thm::single}
Let $m\ge (n+1)(d-r+1)$ and $\lambda_0\in(0, n)$. Under Assumption~\ref{assumption::loss-1} and \ref{assumption::dataset-1}, there exists  a zero measure set $\mathcal{C}\subset\mathbb{R}^{m}$ such that for any $\bm{\lambda}\in(\lambda_0/2,\lambda_{0})^{m}\setminus \mathcal{C}$, both of the following statements are true:
\begin{itemize}[leftmargin=*]
\item[(1)] the empirical loss $L_{n}(\bm{\theta};\bm{\lambda})$ is coercive.
\item[(2)] every local minimum $\bm{\theta}^{*}$ of the loss ${L}_{n}(\bm{\bm{\theta}}^{*};\bm{\lambda})$ achieves a training error at most $\frac{\lambda_0+2E}{ \ell'(0)\gamma n}$, i.e., $R_{n}(\bm{\theta^{*}};f)\le\frac{\lambda_0+2E}{ \ell'(0)\gamma n}$.
\end{itemize}
\end{theorem}
\textbf{Remarks:} The proof is provided in Appendix~\ref{appendix::proof::memo}. (1) The first part shows that the empirical loss is coercive and thus eliminate the possibility that a descent algorithm can diverge to infinity. (2) If we are choosing $\lambda_0=\sqrt{n}$ and if the portion of samples that cannot be separated with a large margin is small (i.e., $E/n$ is small), then the training error at every local minimum is of order $\mathcal{O}\left(\frac{1}{\gamma \sqrt{n}}+\frac{E}{\gamma n}\right)$ and thus very small when the number of samples $n$ is very large. (3) For the fully connected neural network where $r=d$, the amount of neurons we need is just $n+1$.

One may wonder whether over-parameterization
leads to overfitting.
The next theorem states that the test
error of every local minimum is bounded above.

\begin{theorem}\label{thm::single-gen}
Let $m\ge (n+1)(d-r+1)$ and $\lambda_0\in(0, n)$. Assume that the data samples in the dataset are independently drawn from an underlying distribution $\mathbb{P}_{X\times Y}$. Under Assumption~\ref{assumption::loss-1} and \ref{assumption::dataset-1}, with probability at least $1-\delta$, there exists  a zero measure set $\mathcal{C}\subset\mathbb{R}^{m}$ such that for any $\bm{\lambda}\in(\lambda_0/2,\lambda_{0})^{m}\setminus \mathcal{C}$, both of the following statements are true:
\begin{itemize}[leftmargin=*]
\item[(1)] the empirical loss $L_{n}(\bm{\theta};\bm{\lambda})$ is coercive.
\item[(2)] every local minimum $\bm{\theta}^{*}$ of the loss ${L}_{n}(\bm{\bm{\theta}}^{*};\bm{\lambda})$ achieves a test error at most
\begin{align*}
\mathbb{P}(Yf(X;\bm{\theta}^*)<0)=\mathcal{O}\left(\frac{\lambda_0}{\gamma n}+\frac{E\ln n}{\gamma\lambda_0\sqrt{n}}+\frac{E}{\gamma n}+\frac{\ln n}{\gamma\sqrt{n}}+\sqrt{\frac{\log(1/\delta)}{n}}\right).
\end{align*}
\end{itemize}
\end{theorem}
The proof is provided in Appendix~\ref{appendix::proof::gen}. (1) It is easy to see from the first two terms in the upper bound that the regularizer coefficient $\lambda_0$ cannot be too small or too large since the upper bound goes to infinity when  $\lambda_0\rightarrow 0$ or  $\lambda_0\rightarrow \infty$. In fact,  if we choose $\lambda_0=\sqrt{n}\ln n$, then the upper bound becomes $\mathcal{O}\left(\frac{\ln n}{\gamma \sqrt{n}}+\frac{E}{\gamma {n}}+\sqrt{\frac{\log(1/\delta)}{n}}\right)$ and the upper bound goes to $\mathcal{O}(E/(\gamma n))$ as $n$ goes to infinity. When two classes of samples can be well-separated, i.e., $E=0$, the upper bound goes to zero as $n$ goes to infinity. (2) For the fully connected neural network where $r=d$, the amount of neurons we need is just $n+1$.

Remark: A major difference with
earlier works on neural nets is that we directly
bound the generalization error
of the local minimum,
while those works provided bounds on the generalization gap (the gap between training error and test error)
and did not provide bounds on the training error.

\section{Can Other Points that Memorize and Generalize Well?}
\label{sec::other}
In the previous section, we show that every local minimum of the empirical loss memorizes and generalizes well. However, some existing work has shown that finding a local minimum can be hard. Therefore, in this section, we will show that a set of points including but not limited to local minima can also have good memorization and generalization performance. In the next section, we will show that there exists an algorithm finding a point in this set with polynomial number of computations.

Recall that for a single-layered ReLU network consisting of $m$ neurons, the output of the neural network is defined as $f(x;\bm{\theta})=\sum_{j=1}^{m}a_{j}\left(\bm{w}_{j}^{\top}(x\odot \phi_j)\right)_+$ and we define the empirical loss as
\begin{equation}
L_{n}(\bm{\theta};\bm{\lambda})=\sum_{i=1}^{n}\ell(-y_{i}f(x_{i};\bm{\theta}))+\frac{1}{2}\sum_{j=1}^{m}\lambda_{j}\left[a_{j}^{2}+\|\bm{w}_{j}\|^{2}_{2}\right].
\end{equation}

\begin{theorem}[Memorization]\label{thm::other-memo}
Let $m\ge (n+1)(d-r+1)$. Under Assumption~\ref{assumption::loss-1} and \ref{assumption::dataset-1},  if there exists a number $C>0$ such that for any $k=1,...,d-r+1$,
$$\max_{\bm{u}\in\mathbb{B}^d}\left|\sum_{i=1}^n\ell'(-y_i f(x_i;\bm{\theta}^*))y_i(\bm{u}^\top (x_i\odot \phi_k))_+\right|\le C,$$
then the neural network with parameters $\bm{\theta}^*$ achieves a training error at most $\frac{C+2E}{ \ell'(0)\gamma n}$.
\end{theorem}
\textbf{Remark:} The proof is provided in Appendix~\ref{appendix::proof::thm::other-memo}. (1) For every local minimum, we can show that  the inequality holds with $C=\lambda_0$. Thus Theorem \ref{thm::other-memo}
is a generalization of Theorem \ref{thm::single}.
(2) This Theorem shows that for any set of parameters satisfying the above inequalities with a small constant $C$, then this set of parameters can achieve small training error.

Next, we will present the theorem showing that some points other than local minima can also achieve good generalization performance.

\begin{theorem}[Generalization]\label{thm::other-gen}
Let $m\ge (n+1)(d-r+1)$ and choose $\bm{\lambda}\in[\lambda_0/2, \lambda_0]^m$ for some $\lambda_0>0$. Under Assumptions~\ref{assumption::loss-1} and \ref{assumption::dataset-1}, with probablity at least $1-\delta$, for any set of parameters $\bm{\theta}^*$ satisfying all of the following conditions,
\begin{itemize}
\item[(1)] $\|\bm{w}_j^*\|_2= |a_j^*| $ holds for any $j=1,...,m$,
\item[(2)] $\|\nabla_{\bm{a}}L_n(\bm{\theta}^*;\bm{\lambda})\|_2\le \varepsilon$,
\item[(3)] there exists a positive number $C>0$ such that  for any $k=1,...,d-r+1$
$$\max_{\bm{u}\in\mathbb{B}^d}\left|\sum_{i=1}^n\ell'(-y_i f(x_i;\bm{\theta}^*))y_i(\bm{u}^\top (x_i\odot \phi_k))_+\right|\le C,$$
\end{itemize}
the testing error of the neural network $f(\cdot;\bm{\theta}^*)$ has an upper bound of
$$ \mathbb{P}(Yf(X;\bm{\theta}^*)<0)=\mathcal{O}\left(\frac{C+E}{\gamma {n}}+\frac{(C+E)\ln n}{\gamma \lambda_0\sqrt{n}}+\frac{\varepsilon \|\bm{a}^*\|_2}{\lambda_0\sqrt{n}}+\sqrt{\frac{\log(1/\delta)}{n}}\right)
.$$
\end{theorem}
\textbf{Remarks: }The proof is provided in Appendix~\ref{appendix::proof::thm::other-gen}. (1) At every local minimum, $\varepsilon=0$ by the first condition of the critical point and $C=\lambda_0$ by the definition of the local minimum. (2) This theorem shows that if a set parameter is very closed to the critical point (i.e., $\varepsilon$ is small) and the inequalities in the third conditions holds with $C>0$, then by setting $\lambda_0=\sqrt{n}\ln n$, we can still achieve an
upper bound of an order of $\mathcal{O}\left(\frac{C+E}{\gamma {n}}+\frac{\ln n}{\gamma \sqrt{n}}+\sqrt{\frac{\log(1/\delta)}{n}}\right)$.

An algorithm to find a good set of parameters which generalizes well is presented in Algorithm~\ref{algo::main}, which further calls two other algorithms Algorithm \ref{algo::coeff} and \ref{algo::perturb}. In the next section, we present the intuition behind the algorithms and prove that they converge in polynomial time to a good solution.

\begin{algorithm}[t]\label{algo::main}
\SetAlgoLined
\KwResult{ Parameter vector $\bm{\theta}^*=(a_1^*,...,a_m^*,a^*_1\bm{u}^*_1,...,a_m^*\bm{u}_m^*)$ }
\tcp{Initialization}
Choose $|a_j(0)|\le 1$ and $\|\bm{u}_j(0)\|_2=1$ for each $j\in[m]$ and $\bm{\lambda}(0)=\lambda_0\bm{1}_m$ and $\lambda_0\ge \sqrt{n}$ \;
Choose $C\ge \lambda_0$ and $K=\max\{\lceil L_n(\bm{\theta}(0);\bm{\lambda}(0)), 2n\rceil\}$\;
 \tcp{Setting Up The Loop Index }
 $k = 0$\;
 \tcp{Loop}
 \While{ $k\le L_n(\bm{\theta}(0);\bm{\lambda}(0))$ }{
 \iffalse

$$\min_{p_1,...,p_n\in\mathbb{R}}\sum_{j:\phi_j=\phi_s}\left(\lambda_j(k) -\sum_{i=1}^np_iy_i\text{sgn}(a_j(k))(\text{sgn}(a_j(k))\bm{u}_j^\top(k) (x_i\odot \phi_s))_+\right)^2\ge \left(\frac{\lambda_0}{8K}\right)^2,$$
$\bm{\lambda}(k)\preccurlyeq\bm{\lambda}(k-1)$ and $\|\bm{\lambda}(k)-\bm{\lambda}(k-1)\|_\infty\le \frac{\lambda_0}{2K}$\;
\fi
 Call Algorithm~\ref{algo::coeff} to update regularizer coefficient $\bm{\lambda}(k)$\;
  \tcp{ Running The Gradient Descent Algorithm}
  Initialize $\bm{\Theta}(0)=(\alpha_j(0), \alpha_j(0)\bm{u}_j(k))_j$ such that $\alpha_j( 0)=a_j(k)$ for each $j\in[m]$\;
  Setting $t=0$; $\eta_k=\frac{1}{72\max[L_n(\bm{\Theta}(0)), 2n]}$\;
  \While{ $\|\nabla_{\bm{\alpha}}L_{m}(\bm{\Theta}(t);\bm{\lambda}(k))\|_{2} >\frac{\lambda_0}{16K\sqrt{n}}$}
  {$\bm{\alpha}(t+1)=\bm{\alpha}(t)-\eta_k\nabla_{\bm{\alpha}}L_n(\bm{\Theta}( t);\bm{\lambda}(k))$, $t\leftarrow t+1$
  }
  \tcp{Gradient descent algorithm terminates at the $T_k$-th step}
  Call Algorithm~\ref{algo::perturb} to perturb the inactive neuron from $\bm{\Theta}(T_k)$ to $\bm{\theta}(k+1)$ and check the termination condition\;
  $k\leftarrow k+1$

 }
 \caption{Gradient Descent Algorithm}
\end{algorithm}

\begin{algorithm}[t]\label{algo::coeff}
\SetAlgoLined
{
\tcp{Update The Regularizer Coefficients}
  Choose $\bm{\lambda}(k)$ such that   for each $s\in[d-r+1]$
$$\min_{p_1,...,p_n\in\mathbb{R}}\sum_{j:\phi_j=\phi_s}\left(\lambda_j(k) -\sum_{i=1}^np_iy_i\text{sgn}(a_j(k))(\text{sgn}(a_j(k))\bm{u}_j^\top(k) (x_i\odot \phi_s))_+\right)^2\ge \left(\frac{\lambda_0}{8K}\right)^2,$$
$\bm{\lambda}(k)\preccurlyeq\bm{\lambda}(k-1)$ and $\|\bm{\lambda}(k)-\bm{\lambda}(k-1)\|_\infty\le \frac{\lambda_0}{2K}$\;

 }
 \caption{Updating Regularizer Coefficients}
\end{algorithm}

\begin{algorithm}[t]\label{algo::perturb}
\SetAlgoLined
{\KwResult{ Parameter vector $\bm{\theta}(k+1)$. }
\tcp{Perturb the inactive neuron}
   For each $s\in[d-r+1]$, finding the unit vector $\bm{v}_s^*$ such that
  $$\left|\sum_{i=1}^n \ell_i'(T_k)y_i ({\bm{v}_s^*}^\top (x_i\odot \phi_s))_+\right|\ge \max_{\bm{u}\in\mathbb{B}^d}\left|\sum_{i=1}^n \ell_i'(T_k)y_i (\bm{u}^\top (x_i\odot \phi_s))_+\right|-C,$$
  where $\ell_i'(T_k)\triangleq\ell'(-y_if(x_i;\bm{\Theta}(T_k)))$\;
  \tcp{Check termination conditions}
  \eIf{ $\left|\sum_{i=1}^n \ell_i'(T_k)y_i ({\bm{v}_s^*}^\top (x_i\odot \phi_s))_+\right|\le 5\lambda_0,\forall s\in[d-r+1]$}{
  Break the loop and output $\bm{\theta}^*=\bm{\Theta}(T_{k})$\;
   }{
   Finding an index $s\in[d-r+1]$ such that $\left|\sum_{i=1}^n \ell_i'(T_k)y_i ({\bm{v}_s^*}^\top (x_i\odot \phi_s))_+\right|> 5\lambda_0$\;
   Finding an index $j\in[m]$ such that $\phi_j=\phi_s$ and $|a_j|\le \frac{1}{\sqrt{4n}}$\;
   $a_j(k+1)= \text{sgn}\left(\sum_{i=1}^n \ell_i'(T_k)y_i ({\bm{v}_s^*}^\top (x_i\odot \phi_s))_+\right)\sqrt{\frac{1}{\lambda_0}}$\;
   $\bm{u}_j(k+1)=\text{sgn}(a_j(k+1))\bm{v}_s^*$\;
   $a_r(k+1)= a_r(k), \bm{u}_r(k+1)=\bm{u}_r(k)$ for $r\neq s$\;
  }
 }
 \caption{Perturb Inactive Neurons And Check Termination Condition}
\end{algorithm}

\section{A Polynomial Time Training Algorithm}

\subsection{Algorithm}\label{sec::intuition}
In this subsection, we provide intuition as to why Algorithm~\ref{algo::main} finds parameters satisfying all conditions in Theorem~\ref{thm::other-gen}.  .

To satisfy the first condition in Theorem~\ref{thm::other-gen} (i.e., $\|\bm{w}_j^*\|_2= |a_j^*| $ for any $j=1,...,m$), we just need to set $\bm{w}_j= \alpha_j\bm{u}_j, a_j=\alpha_j$ for some unit vector $\bm{u}_j$ and update the scalar $\alpha_j$ and the unit vector $\bm{u}_j$, respectively. Now the parameter vector  becomes $\bm{\Theta}=(\alpha_j, \alpha_j\bm{u}_j)_{j=1,\dots, m}$. Considering the update rule for each $\alpha_j$, i.e.,
\begin{align*}
\alpha_j(t+1)&=\alpha_j(t)-\eta\nabla_{\alpha_j}L_n(\bm{\Theta}(t))\\
&=\left[1-2\eta\lambda_j+2\eta\sum_{i=1}^n\ell'_i(-y_i f(x_i;\bm{\Theta}(t)))y_iy_i\text{sgn}(\alpha( t))(\text{sgn}(\alpha(t))\bm{u}_j^\top (x_i\odot \phi_j))_+\right]\alpha_j(t)
\end{align*}
and choosing a sufficiently small step size $\eta\le \frac{1}{8n}$, we observe that $\text{sgn}(\alpha_j(t+1))=\text{sgn}(\alpha_j(t))=\text{sgn}(\alpha_j(0))$ holds for all $t=0,1,2...$. This indicates that (1) we can run the gradient descent algorithm  on the vector $\bm{\alpha}$ for a sufficiently long time to ensure that the gradient norm $\|\nabla_{\bm{\alpha}}L_n(\bm{\Theta}(t);\lambda)\|_2$ is sufficiently small and that  (2) the sign of each $\alpha_j$ does not change with respect to the time. Furthermore, using the fact that $\text{sgn}(\alpha_j(t))=\text{sgn}(\alpha_j(0))$ for all $t$, we obtain that for each unique template feature vector $\phi_k$, $1\le k\le d-r+1$ (the series $\{\phi_k\}_{k\ge 1}$ is a periodic series of a period equal to $d-r+1$), we have
\begin{align*}
\|\nabla_{\bm{\alpha}}L_n(\bm{\Theta}(t))\|^2_2
&=\sum_{j=1}^m\left(\lambda_j -\sum_{i=1}^n\ell'(1-y_if(x_i;\bm{\Theta}(t)))y_i\text{sgn}(\alpha_j(0))(\text{sgn}(\alpha_j(0))\bm{u}_j^\top (x_i\odot \phi_j))_+\right)^2\alpha_j^2(t)\\
&\ge \sum_{j:\phi_j=\phi_k}\left(\lambda_j -\sum_{i=1}^n\ell'(1-y_if(x_i;\bm{\Theta}(t)))y_i\text{sgn}(\alpha_j(0))(\text{sgn}(\alpha_j(0))\bm{u}_j^\top (x_i\odot \phi_k))_+\right)^2\alpha_j^2(t).
\end{align*}
Thus, if for each $k\in[d-r+1]$, we can choose regularizer coefficients $\lambda_j$s satisfying
\begin{equation}\label{eq::intuition::1}
\min_{p_1,...,p_n\in\mathbb{R}}\sum_{j:\phi_j=\phi_k}\left(\lambda_j -\sum_{i=1}^np_iy_i\text{sgn}(\alpha_j(0))(\text{sgn}(\alpha_j(0))\bm{u}_j^\top (x_i\odot \phi_k))_+\right)^2\ge \text{const}^2
\end{equation}
then when the gradient descent norm is small, for each unique template feature vector $\phi_k$, $k\in[d-r+1]$, the parameters of one of the neurons is also small, i.e.,
$$\|\nabla_{\bm{\alpha}}L_n(\bm{\Theta}(t))\|_2\le \delta\times\text{const}\implies \max_{k\in[d-r+1]}\min_{j:\phi_j=\phi_k} |\alpha_j|=\max_{k\in[d-r+1]}\min_{j:\phi_j=\phi_k}\|\alpha_j\bm{u}_j\|_2\le\delta.$$

Therefore, we only need to handle the third condition in Theorem~\ref{thm::other-gen}. From the prior analysis, we know that after running the gradient descent algorithm for a long time, for each template feature vector $\phi_k$, there exists a neuron with small parameters. This indicates that  we can perturb the inactive neuron in the following way to make the empirical loss decrease further. Since
 the loss function $\ell$ is gradient Lipschitz, we obtain the following inequality
$$L_n(\tilde{\bm{\theta}};\bm{\lambda})-L_n(\bm{\Theta};\bm{\lambda})\le -\sum_{i=1}^n \ell'(-y_i f(x_i;\bm{\theta}))y_i \tilde{\alpha}_j(\tilde{\alpha}_j\tilde{\bm{u}}_j^\top (x_i\odot \phi_j))_++\lambda_j\tilde{\alpha}_j^2+2n\tilde{\alpha}_j^4+n\alpha_j^2+2n\alpha_j^4,$$
when we only perturb the $j$-th neuron from $(\alpha_j, \alpha_j\bm{u}_j)$ to $(\tilde{\alpha}_j,\tilde{\alpha}_j\tilde{\bm{u}}_j)$ for an arbitrary unit vector $\tilde{\bm{u}}_j$. Since we already know that after running the gradient descent algorithm for a sufficiently long time, for each unique feature vector $\phi_k$, there exists a neuron of index $k_j$ and with $|\alpha_{k_j}|\le \delta= 1/\sqrt{4n}$, then using the fact that $\lambda_0\ge \sqrt{n}$,  we have
$$L_n(\tilde{\bm{\theta}};\bm{\lambda})-L_n(\bm{\theta};\bm{\lambda})\le -\frac{1}{\lambda_0}\sum_{i=1}^n \ell'(-y_i f(x_i;\bm{\theta}))y_i \text{sgn}(\tilde{\alpha}_{k_j})(\text{sgn}(\tilde{\alpha}_{k_j})\tilde{\bm{u}}_{k_j}^\top (x_i\odot \phi_k))_++4,$$
when we only perturb the $k_j$-th neuron from $(\alpha_{k_j}, \alpha_{k_j}\bm{u}_{k_j})$ to $(\tilde{\alpha}_{k_j},\tilde{\alpha}_{k_j}\tilde{\bm{u}}_{k_j})$ for an arbitrary unit vector $\tilde{\bm{u}}_{k_j}$ and $|\tilde{\alpha}_{k_j}|=1/\sqrt{\lambda_0}$. Therefore, for any arbitrary unit vector $\bm{v}_k$, setting $\text{sgn}(\tilde{\alpha}_{k_j})=\text{sgn}\left(\sum_{i=1}^n \ell'(-y_i f(x_i;\bm{\theta}))y_i (\bm{v}_k^\top (x_i\odot \phi_k))_+\right)$ and $\tilde{\bm{u}}_{k_j}=\text{sgn}(\tilde{\alpha}_{k_j})\bm{v}_k$, we have
$$L_n(\tilde{\bm{\theta}};\bm{\lambda})-L_n(\bm{\theta};\bm{\lambda})\le -\frac{1}{\lambda_0}\left|\sum_{i=1}^n \ell'(-y_i f(x_i;\bm{\theta}))y_i (\bm{v}_k^\top (x_i\odot \phi_k))_+\right|+4.$$
Assume that we have an algorithm which always finds a vector $\bm{v}_k^*$ for each $k\in[d-r+1] $satisfying $$\left|\sum_{i=1}^n \ell'(-y_i f(x_i;\bm{\theta}))y_i ({\bm{v}_k^*}^\top (x_i\odot \phi_k))_+\right|\ge \max_{\bm{u}\in\mathbb{B}^d}\left|\sum_{i=1}^n \ell'(-y_i f(x_i;\bm{\theta}))y_i (\bm{u}^\top (x_i\odot \phi_k))_+\right|-C$$
for some positive constant $C>0$.
If for each $k\in[d-r+1]$, $\bm{v}_k^*$ satisfies
$$\left|\sum_{i=1}^n \ell'(-y_i f(x_i;\bm{\theta}))y_i ({\bm{v}_k^*}^\top (x_i\odot \phi_k))_+\right|\le 5\lambda_0$$
then we terminate the algorithm; otherwise, we set the new $\alpha_{k_j}$ and $\bm{u}_{k_j}$ to $\tilde{\alpha}_{k_j}$ and $\tilde{\bm{u}}_{k_j}$, respectively, and rerun the gradient descent.
Therefore, this means that before the algorithm terminates, each time when we perturb the inactive neuron, the empirical loss  decreases by a constant
$$L_n(\tilde{\bm{\theta}};\bm{\lambda})-L_n(\bm{\theta};\bm{\lambda})\le -1.$$
It further indicates that when the algorithm terminates, we should have
$$ \max_{\bm{u}\in\mathbb{B}^d}\left|\sum_{i=1}^n \ell'(-y_i f(x_i;\bm{\theta}))y_i (\bm{u}^\top (x_i\odot \phi_1))_+\right|\le 5{\lambda_0}+C$$
which makes the third condition in Theorem~\ref{thm::other-gen} hold. We note here that if the algorithm continues, then in the next iteration, we  have to update the regularizer coefficient vector $\bm{\lambda}$ to make the inequality~\eqref{eq::intuition::1} hold since we update the unit vector $\bm{u}_k$s. Furthermore, when we update the coefficient vector $\bm{\lambda}$, we also need to ensure the empirical loss does not increase and the coefficient vector $\bm{\lambda}$ always locates in the region $[\lambda_0/2, \lambda_0]^m$.

\subsection{Performance and Complexity Analysis}
In this subsection, we present the result on the performance and computational complexity of Algorithm~\ref{algo::main}. We first show that Algorithm~\ref{algo::main} finds the output $\bm{\theta}^*$ with small testing error if the number of samples is sufficiently large.

\begin{theorem}\label{thm::algo}
Let $m\ge (n+1)(d-r+1)$ and $\lambda_0=\sqrt{n}\ln n$. Assume that the data samples in the dataset are independently drawn from an underlying distribution $\mathbb{P}_{X\times Y}$. Under Assumption~\ref{assumption::loss-1} and \ref{assumption::dataset-1}, with probability at least $1-\delta$, the output $\bm{\theta}^*$ of the Algorithm~\ref{algo::main} satisfies
\begin{align*}
\mathbb{P}(Yf(X;\bm{\theta}^*)<0)=\mathcal{O}\left(\frac{C+E}{\gamma n}+\frac{\ln n}{\gamma \sqrt{n}}+\sqrt{\frac{\log(1/\delta)}{n}}\right).
\end{align*}
\end{theorem}
\textbf{Remark:} The proof is provided in Appendix~\ref{appendix::proof::thm-algo}. This upper bound on the testing error depends on a constant $C$, which denotes how well we find a unit vector $\bm{v}_s^*$  for each $s\in[d-r+1]$ such that
  $$\left|\sum_{i=1}^n \ell'(-y_if(x_i;\bm{\theta}))y_i ({\bm{v}_s^*}^\top (x_i\odot \phi_s))_+\right|\ge \max_{\bm{u}\in\mathbb{B}^d}\left|\sum_{i=1}^n \ell'(-y_if(x_i;\bm{\theta}))y_i (\bm{u}^\top (x_i\odot \phi_s))_+\right|-C.$$
If we can perfectly solve the optimization problem on the right hand side of the inequality, then $C=0$. Therefore, it is straightforward to see that the complexity and testing performance of Algorithm~\ref{algo::main} depends on the method we use to solve the above optimization problem.

For each $s\in[d-r+1]$, an easy way to solve the optimization problem
$$\max_{\bm{u}\in\mathbb{B}^d}G(\bm{u};\phi_s,\bm{\theta})\triangleq\left|\sum_{i=1}^n \ell'(-y_if(x_i;\bm{\theta}))y_i (\bm{u}^\top (x_i\odot \phi_s))_+\right|$$
is  exhaustive search. Specifically, for each $s\in[d-r+1]$, we can search every possible vector in the following discrete set,
$$\mathcal{A}_s=\left\{\bm{u}\in\mathbb{B}^d:\bm{u}(k)=0\text{ for }k\notin I_s\text{ and }\bm{u}(k)\in\bigcup_{h=0}^{\lceil \lambda_0\rceil}\left\{\frac{h\lambda_0}{n}\right\}\text{ for } k\in I_s\right\},$$
where the index set $I_s\triangleq\{j\in[d]: \phi_s(d)=1\}$, and set the output $\bm{v}_s^*$ be the vector $\bm{u}$ maximizing $G$. Since by assumption, $|I_s|\le r$ for  each $s\in[d-r+1]$, then this exhaustive search algorithm has a total computation complexity of $\mathcal{O}((d-r+1)(n/\lambda_0)^r)=\mathcal{O}\left({dn^{r/2}}\right)$ since $\lambda_0\ge \sqrt{n}$. Furthermore, the approximation error of the output of this exhaustive search algorithm is
$$\max_{j\in[d-r+1]}\left|G(\bm{v}^*_j;\phi_j,\bm{\theta})-\max_{\bm{u}\in\mathbb{B}^{d}}G(\bm{u};\phi_j,\bm{\theta})\right|\le \frac{\lambda_0}{n}\sum_{i=1}^n\ell'(-y_i f(x_i;\bm{\theta}))\le \lambda_0.$$
Combining this inequality with Theorem~\ref{thm::algo}, we have the following Corollary.
\begin{corollary}\label{cor::1}
Let $m\ge (n+1)(d-r+1)$ and $\lambda_0=\sqrt{n}\ln n$. Assume that the data samples in the dataset are independently drawn from an underlying distribution $\mathbb{P}_{X\times Y}$. Under Assumption~\ref{assumption::loss-1} and \ref{assumption::dataset-1}, with probability at least $1-\delta$, the output $\bm{\theta}^*$ of  Algorithm~\ref{algo::main} satisfies
\begin{align*}
\mathbb{P}(Yf(X;\bm{\theta}^*)<0)=\mathcal{O}\left(\frac{E}{\gamma n}+\frac{\ln n}{\gamma \sqrt{n}}+\sqrt{\frac{\log(1/\delta)}{n}}\right).
\end{align*}
Furthermore, Algorithm~\ref{algo::main} has a computation complexity of $\mathcal{O}\left(dKm^2+dKn^{r/2}+nK^5/\lambda_0\right)$ where $K=\max\{\lceil L_n(\bm{\theta}(0);\bm{\lambda}(0)), 2n\rceil\}$.
\end{corollary}
\textbf{Remark:} The proof is provided in Appendix~\ref{appendix::proof::corollary-1}. We note here that for
CNN with a given template size $r$, the complexity of Algorithm~\ref{algo::main} is of polynomial order of dimension $d$, samples size $n$ and the network size $m$.
For instance, when the template size is $r = 3$
and $ m = (n+1) (d - 2 )$,
the complexity is $ \mathcal{O}( K n^2 d^3
 + d K n^{3/2} + \sqrt{n} K^5/\ln n ) .$
 For FNN where $r = d$,
 this complexity is not polynomial,
and we will discuss how to reduce the complexity
for FNN next.

\subsection{Discussion: Improved Complexity for FNN}
From Corollary~\ref{cor::1}, we can see that when $r=d$ (i.e., the neural network is a fully connected neural network), the algorithm has a computation complexity of $\mathcal{O}(n^{d/2})$. This is due to the fact that solving the following optimization problem requires a computational complexity of $\mathcal{O}(n^{d/2})$,
\begin{align}\label{eq::opt}
\max_{\bm{u}\in\mathbb{B}^d}\left|\sum_{i=1}^n \ell'(-y_if(x_i;\bm{\theta}))y_i (\bm{u}^\top x_i)_+\right|\end{align}
where the vector $\bm{\theta}$ is given.
In this subsection, we will show that when the underlying data distribution has some additional properties, the above optimization problem can be solved by an algorithm of much smaller complexity with high probability.

\begin{assumption}\label{assump::telgarsky}
Let $\gamma \ge 8\varepsilon>0$. Assume that there exists a function $\bar{\bm{v}}(\omega):\mathbb{R}^d\rightarrow \mathbb{R}^d$ with $\|\bar{\bm{v}}(\omega)\|_2\le 1$  such that
$$\mathbb{P}\left(Y\int\bar{\bm{v}}^\top({\omega})X\mathbb{I}\{X^\top \omega\ge 0\}\mu(d\omega)\ge \gamma\right)\ge 1-\varepsilon,$$
where $\mu$ is the uniform measure defined on the sphere $\mathbb{S}^{d-1}$.
\end{assumption}

\begin{algorithm}[t]\label{algo::opt}
\SetAlgoLined
{
\tcp{Fast Solver of Problem~\eqref{eq::opt}}
\tcp{At the initialization of Algorithm~\ref{algo::main}}
  Choose $\bm{w}_1,...,\bm{w}_{2M}$ from the sphere $\mathbb{S}^{d-1}$ uniformly at random and choose $r>0$\;
 \tcp{At each time of running Algorithm~\ref{algo::perturb}, perturb the inactive neuron in the following direction.}
For each $j\in[2M]$, calculate $$\bm{v}_j=\frac{\sum_{i=1}^n y_i\ell'(-y_i f(x_i;\bm{\theta})) x_i \mathbb{I}\{\bm{\omega}_j^\top x_i\ge 0\}}{\left\|\sum_{i=1}^n y_i\ell'(-y_i f(x_i;\bm{\theta})) x_i \mathbb{I}\{\bm{\omega}_j^\top x_i\ge 0\}\right\|_2}.$$
 }
 Then $\bm{v}^*=\bm{\omega}_{j^*} + ra_j \bm{v}_{j^*}^*$ where
 $$j^*=\arg\max_j \left|\sum_{i=1}^n y_i\ell'(-y_i f(x_i;\bm{\theta}))\left[(\bm{\omega}_j+ra_j \bm{v}^*_j)^\top x_i\right]_+\right|$$
 \caption{Finding The Optimal Descent Direction $\bm{v}^*$}
\end{algorithm}
Now we will show that under Assumption~\ref{assump::telgarsky}, the optimization problem defined by \eqref{eq::opt} can be solved by  Algorithm~\ref{algo::opt}. The performance of the above algorithm is guaranteed by the following lemma.
\begin{lemma}\label{lemma::algo-performance}
Let $\varepsilon_0=\max\{\varepsilon, n^{-1/3}\}$ and  $\bm{\omega}_1,...,\bm{\omega}_{2M}$ be i.i.d. random vectors uniformly distributed on the sphere $\mathbb{S}^{d-1}$. Under Assumption~\ref{assump::telgarsky}, if $$r<\frac{\gamma-4\varepsilon_0}{16d},\quad n\ge \frac{\ln (6/\delta)}{2\varepsilon_0^2}\quad\text{and}\quad M\ge \max\left\{\frac{\ln(4n/\delta)}{\gamma^2},\frac{4\ln(6n/\delta)}{r^2(\gamma-4\varepsilon_0)^2}\right\},$$ then
\begin{align*}
\mathbb{P}\left( \forall \bm{\beta}\in[0,1]^n:\max_{j\in[2M]}\left|\sum_{i=1}^n y_i\beta_i\left[(\bm{\omega}_j+ra_j \bm{v}_j)^\top x_i\right]_+\right|\ge \frac{r(\gamma-4\varepsilon_0)}{8}\max_{\bm{u}\in\mathbb{B}^d}\left|\sum_{i=1}^ny_i\beta_i(\bm{u}^\top x_i)\right|\right)\ge 1-\delta
\end{align*}
where
$$\bm{v}_j=\frac{\sum_{i=1}^n y_i\beta_i x_i \mathbb{I}\{\bm{\omega}_j^\top x_i\ge 0\}}{\left\|\sum_{i=1}^n y_i\beta_i x_i \mathbb{I}\{\bm{\omega}_j^\top x_i\ge 0\}\right\|_2}.$$
\end{lemma}
\textbf{Remark: }Lemma~\ref{lemma::algo-performance} ensures that if the radius $r$ is sufficiently small, the sample size $n$ is sufficiently large and the number of random vector $M$
is sufficiently large, then with high probability, the $v^*$ found by Algorithm~\ref{algo::opt} is an approximate solution of the problem~\eqref{eq::opt}.
\begin{corollary}\label{cor::2}
Let $m\ge n+1$ and $\lambda_0=\sqrt{n}\ln n$. Assume that the data samples in the dataset are independently drawn from an underlying distribution $\mathbb{P}_{X\times Y}$. Under Assumption~\ref{assumption::loss-1} and \ref{assump::telgarsky}, with probability at least $1-\delta$, the output $\bm{\theta}^*$ of  Algorithm~\ref{algo::main} satisfies
\begin{align*}
\mathbb{P}(Yf(X;\bm{\theta}^*)<0)=\mathcal{O}\left(\frac{\varepsilon}{\gamma }+\frac{\ln n}{r\gamma^2 \sqrt{n}}+\sqrt{\frac{\log(1/\delta)}{n}}\right).
\end{align*}
Furthermore, Algorithm~\ref{algo::main} has a computation complexity of $\mathcal{O}(dKm^2+ MK + nK^5/\lambda_0).$ where $K=\max\{\lceil L_n(\bm{\theta}(0);\bm{\lambda}(0)), 2n\rceil\}$.
\end{corollary}

\begin{proof}
The upper bound on $\mathbb{P}(Yf(X;\bm{\theta}^*)<0)$ follows directly from Theorem~\ref{thm::other-gen}, with $C=\frac{40\lambda_0}{r\gamma}$, $\varepsilon = \frac{\lambda_0}{16K\sqrt{n}}$ and $\|\bm{a}^*\|_2\le 2\sqrt{\frac{K}{\lambda_0}}$, where the upper bound on $\varepsilon$ and $\|\bm{a}^*\|_2$ follows from Claim~\ref{claim::algo::2} in Appendix~\ref{appendix::proof::thm::other-gen}. Furthermore, following the same analysis on the computational complexity in the proof of Corollary~\ref{cor::1}, we know that the complexity of choosing the regularization coefficient is $\mathcal{O}(dm^2)$ and the complexity of running gradient descent is $\mathcal{O}(nK^4/\lambda_0)$. Furthermore, it is easy to see that running Algorithm~\ref{algo::opt} requires time complexity of $\mathcal{O}(M)$. Since there are $K$ iterations in total, then the time complexity of the whole algorithm is $\mathcal{O}(dKm^2+ MK + nK^5/\lambda_0).$
\end{proof}

\iffalse
\textbf{Initialization}. Choose any $\bm{\theta}(0)$ such that $L_n(\bm{\theta}(0))<\infty$.

\textbf{Running}. Finding a series of parameters $\{\bm{\theta}(k)\}_{k\ge 1}$ and regularizer coefficient vectors $\{\bm{\lambda}(k)\}_{k\ge 0}\subset[\lambda_0/2, \lambda_0]$ such that there exists two positive constant $\phi_1$, $\phi_2$, $\phi_3$ and for each $k=1,2,...,$
\begin{itemize}
    \item[(1)] $\|\bm{w}_j(k)\|_2=|a_j(k)|$ for $j=1,...,m$,
    \item[(2)] $\|\nabla_{\bm{a}}L_n(\bm{\theta}(k);\bm{\lambda}(k))\|_2\le \sqrt{1/n}$
    \item[(3)]
    $$\max_{1\le j\le m}\left|\sum_{i=1}^n\ell'(-y_i f(x_i;\bm{\theta}(k)))y_i \left(\frac{\bm{w}_j^\top(k) x_i}{\|\bm{w}_j(k)\|_2}\right)_+\right|\ge \phi_1 \max_{\bm{u}\in\mathbb{B}^d}\left|\sum_{i=1}^n\ell'(-y_i f(x_i;\bm{\theta}(k)))y_i(\bm{u}^\top (x_i\odot \phi_j))_+\right|$$
    \item[(4)] $L_n(\bm{\theta}(k);\bm{\lambda}(k))\le L_n(\bm{\theta}(k-1);\bm{\lambda}(k-1))-\phi_2$.
\end{itemize}

\textbf{Termination}. We terminate the process if
$$\max_{1\le j\le m}\left|\sum_{i=1}^n\ell'(-y_i f(x_i;\bm{\theta}(k)))y_i \left(\frac{\bm{w}_j^\top(k) x_i}{\|\bm{w}_j(k)\|_2}\right)_+\right|\le \phi_3\lambda_0.$$
\begin{lemma}
We assume that the above process terminates on the $M$-th iteration. Then $M\le L_n(\bm{\theta}(0))/c_2$ and the parameter vector $\bm{\theta}(M)$ satisfies the following three conditions:
\begin{itemize}
    \item[(1)] $\|\bm{w}_j(M)\|_2=|a_j(M)|$ for $j=1,...,m$
    \item[(2)] $\|\nabla_{\bm{a}}L_n(\bm{\theta}(M);\bm{\lambda}(M))\|_2\le \sqrt{1/n}$
    \item[(3)] $$\max_{\bm{u}\in\mathbb{B}^d}\left|\sum_{i=1}^n\ell'(-y_i f(x_i;\bm{\theta}(M)))y_i(\bm{u}^\top (x_i\odot \phi_j))_+\right|\le \frac{\phi_3\lambda_0}{\phi_1}.$$
\end{itemize}

\end{lemma}

\textbf{Remark:} Using Theorem~\ref{thm::other-gen}, we are able to bound the testing error of $\bm{\theta}(M)$.

\begin{proof}
We assume that the process terminates on the $M$-th iteration. Clearly, $M\le L_n(\bm{\theta}(0))/\phi_2<\infty$ since the empirical loss $L_n$ is non-negative. Further, at the $M$-th iteration, we have
(1) $\|\bm{w}_j(M)\|_2=|a_j(M)|$ for $j=1,...,m$;  (2) $\|\nabla_{\bm{a}}L_n(\bm{\theta}(M);\bm{\lambda}(M))\|_2\le \sqrt{1/n}$.
(3) $$\max_{\bm{u}\in\mathbb{B}^d}\left|\sum_{i=1}^n\ell'(-y_i f(x_i;\bm{\theta}(M)))y_i(\bm{u}^\top (x_i\odot \phi_j))_+\right|\le \frac{\phi_3\lambda_0}{\phi_1},$$
where the part (3) follows from the condition that
$$\max_{1\le j\le m}\left|\sum_{i=1}^n\ell'(-y_i f(x_i;\bm{\theta}(M)))y_i \left(\frac{\bm{w}_j^\top(M) x_i}{\|\bm{w}_j(M)\|_2}\right)_+\right|\ge \phi_1 \max_{\bm{u}\in\mathbb{B}^d}\left|\sum_{i=1}^n\ell'(-y_i f(x_i;\bm{\theta}(M)))y_i(\bm{u}^\top (x_i\odot \phi_j))_+\right|$$
and termination condition
$$\max_{1\le j\le m}\left|\sum_{i=1}^n\ell'(-y_i f(x_i;\bm{\theta}(M)))y_i \left(\frac{\bm{w}_j^\top(M) x_i}{\|\bm{w}_j(M)\|_2}\right)_+\right|\le \phi_3\lambda_0.$$

\end{proof}

\subsection{Warm-up of Constructing Algorithms}
Now we will explain more in detail how we use the intuition presented in Section~\ref{sec::intuition} to construct the algorithm.

\textbf{Initialization}. We are setting $a_j(0)=\|\bm{w}_j(0)\|_2$ and $L_n(\bm{\theta}_0)\le 2n$.

\textbf{Balanced parameters and small gradient norm}. We first explain how we make the parameters balanced and gradient norm small. The main fact we exploit is that at the $k$-th outer iteration, if we are writing $\bm{w}(k,t)_j=c_j(k, t)\bm{u}(k)_j$ for some unit vectors $\bm{u}(k)_j$, initially setting $a_j(k,0)=c_j(k,0)$, fixing all unit vectors $\bm{u}_j$ and only running gradient descent on $\bm{a}$ and $\bm{c}$, i.e.,
$$\bm{a}(k, t+1)=\bm{a}(k, t)-\eta \nabla_{\bm{a}}L_n(\bm{\theta}(k, t)),\quad \bm{c}(k, t+1)=\bm{c}(k, t)-\eta \nabla_{\bm{c}}L_n(\bm{\theta}(k, t))$$
then in the whole gradient descent process, $a_j(k, t)=c_j(k, t)$, $t=1,2,...$. The other useful fact used in our analysis is that during the whole gradient process, the sign of $a_j(k, t)$ and $c_j(k, t)$ does not change with respect to $t$.  Furthermore, since we are running the gradient descent on $\bm{a}=(a_1,...,a_m)$ and the empirical loss function is locally gradient Lipschitz, then for sufficiently small constant step size, the gradient norm $\min_t\|\nabla_{\bm{a}}L_n(\bm{\theta}(k, t))\|_2$ converges to zero at a rate of $\mathcal{O}(n/\sqrt{t})$.

\textbf{Inactive neuron.} By simple calculation, we know that
\begin{align*}
\|\nabla_{(\bm{a}, \bm{c})}L_n(\bm{\theta})\|^2_2&= \sum_{j=1}^m\left(\lambda_j a_j-\sum_{i=1}^n\ell'(1-y_if(x_i;\bm{\theta}))y_i(a_j\bm{u}_j^\top (x_i\odot \phi_j))_+\right)^2\\
&= \sum_{j=1}^m\left(\lambda_j -\sum_{i=1}^n\ell'(1-y_if(x_i;\bm{\theta}))y_i\text{sgn}(a_j)(\text{sgn}(a_j)\bm{u}_j^\top (x_i\odot \phi_j))_+\right)^2a_j^2
\end{align*}
and that in each outer iteration, during the gradient descent process, $\text{sgn}(a_j(k, t))$ and $\bm{u}_j(k)$ does not change with respect to the time $t$. This means that if before running the gradient descent, we carefully choose $\bm{\lambda}(k)\preccurlyeq\bm{\lambda}(k-1)$  so that
\begin{equation*}
\min_{p_1,...,p_n\in[0,1]}\sum_{j=1}^m\left(\lambda_j(k)-\text{sgn}(a_j(k, 0))\sum_{i=1}^np_iy_i(\text{sgn}(a_j(k, 0))\bm{u}_j(k)^\top (x_i\odot \phi_j))_+\right)^2\ge \text{const}(n)
\end{equation*}
then
$$\min_j a^2_j(k, t)\le \frac{\|\nabla_{(\bm{a}, \bm{c})}L_n(\bm{\theta}(k, t))\|^2_2}{\text{const}(n)}$$
where const$(n)$ denotes some constant depending only on $n$. This indicates that when the gradient norm is suffciently small, one of the neuron must have a suffciently small $|a_j|=\|\bm{w}_j\|_2=|c_j|$.

\textbf{Perturbing the inactive neuron}. If one of the neuron $j$ has a sufficiently small $|a_j|=\|\bm{w}_j\|_2=|c_j|\le \varepsilon$, now stop the gradient descent and perturb it. Without loss of generality, we assume that  $|a_1|=\|\bm{w}_1\|_2=|c_1|\le \varepsilon$ for some sufficiently small $\varepsilon$, then for a given unit vector $\bm{v}$, if we set
\begin{align*}
\tilde{a}_{1}&=\tilde{c}_{1}=\delta\text{sgn}\left(\sum_{i=1}^n\ell'(-y_i f(x_i;\bm{\theta}))y_i (\bm{v}^\top (x_i\odot \phi_j))_+\right),\\
\tilde{\bm{u}}_{1}&=\delta\bm{v}\text{sgn}\left(\sum_{i=1}^n\ell'(-y_i f(x_i;\bm{\theta}))y_i (\bm{v}^\top (x_i\odot \phi_j))_+\right)
\end{align*}
and by Taylor expansion, we should have
\begin{align}\label{eq::warmup-1}
L_n(\tilde{\bm{\theta}})-L_n(\bm{\theta})\le -\delta^2 \left|\sum_{i=1}^n\ell'(-y_i f(x_i;\bm{\theta}))y_i (\bm{v}^\top (x_i\odot \phi_j))_+\right|+\lambda_0\delta^2 +n\varepsilon^2 + \frac{n(\varepsilon^2+\delta^2)^2}{8}.
\end{align}
Therefore, to find the descent perturbation direction, we have to solve the following optimization problem
$$\max_{\bm{v}\in\mathbb{B}^d}G(\bm{u}_j;\alpha,\phi_j)\triangleq\left|\sum_{i=1}^n\ell'(-y_i f(x_i;\bm{\theta}))y_i (\bm{v}^\top (x_i\odot \phi_j))_+\right|.$$
We note here that by definition, $\|\phi_j\|_1\le r$. Then it is easy to see that for each $j\in[d-r+1]$, this maximization problem can be solved by exhaustively searching in the space of parameters $\bm{v}(k)$'s where $\phi_j(k)\neq 0$. Specifically, for each given vector $\phi_j$, we can exhaustively search the parameters $\{\bm{v}(k):\phi_j(k)=1\}$ in the set $\{1/\lceil \sqrt{n}\rceil, 2/\lceil \sqrt{n}\rceil,...,1\}$ and output the answer $\bm{v}_j$ to be vector achieves the maximum value on $\left|\sum_{i=1}^n\ell'(-y_i f(x_i;\bm{\theta}))y_i (\bm{v}^\top (x_i\odot \phi_j))_+\right|$. Then for each $\phi_j$, we can find an approximate solution $\bm{v}_j$ with an approximation error of $\sqrt{n}$ within $(\lceil \sqrt{n}\rceil)^{r}$ steps.

\textbf{Termination of the algorithm.} We terminate the algorithm when $\max_{j\in[d-r+1]}G(\bm{u}_j;\alpha,\phi_j)\le 2\lambda_{\max}$, this means that before the termination,
$$\max_{j\in[d-r+1]}G(\bm{u}_j;\alpha,\phi_j)> 2\lambda_{\max}.$$
This provides the constant decrements on the empirical loss when we perturb the inactive neuron. At the same time, when the algorithm terminates, we have
\begin{align*}
2\lambda_{\max}\ge \max_{j\in[d-r+1]}G(\bm{u}_j;\alpha,\phi_j)\ge \max_{j\in[d-r+1]}\max_{\bm{u}\in\mathbb{S}^{d-1}}G(\bm{u};\alpha,\phi_j)- \sqrt{n}.
\end{align*}
This means that before the algorithm terminates, the inequality
$$\max_{j\in[d-r+1]}G(\bm{u}_j;\alpha,\phi_j)> 2\lambda_{\max}$$
always holds and when the algorithm terminates, we should have
$$\max_{j\in[d-r+1]}\max_{\bm{u}\in\mathbb{S}^{d-1}}G(\bm{u};\alpha,\phi_j)\le 2\lambda_{\max}+\sqrt{n}\le 3\lambda_0,$$
where the last inequality follows from the fact that $\lambda_{\max}\ge \sqrt{n}$.
Recall the inequality \eqref{eq::warmup-1} that when are perturbing some inactive neuron with the mask vector $\phi_j$, we  have
\begin{align*}
L_n(\tilde{\bm{\theta}})-L_n(\bm{\theta})\le -\delta^2 \left|\sum_{i=1}^n\ell'(-y_i f(x_i;\bm{\theta}))y_i (\bm{v}^\top (x_i\odot \phi_j))_+\right|+\lambda_{\max}\delta^2 +n\varepsilon^2 + \frac{n(\varepsilon^2+\delta^2)^2}{8}.
\end{align*}
Setting $\delta=\sqrt{8\lambda_0/(9n)}$, $\varepsilon \le \delta\sqrt{\lambda_{\max}/(2n)}$ and $\lambda_0\ge \sqrt{n}$, then
\begin{align*}
L_n(\tilde{\bm{\theta}})-L_n(\bm{\theta})\le -1/18.
\end{align*}
This indicates that each time we perturb the inactive neuron, the empirical loss will decrease by $1/18$. Since the initial value of the empirical loss is bounded and the empirical loss is always non-negative, then the algorithm will terminate within $18L_n(\bm{\theta}_0)$ iterations.

Based on the above analysis, we can check that all of three conditions in Theorem~\ref{thm::other-gen} are satisfied and thus we can use Theorem~\ref{thm::other-gen} to show that the testing error of the obtained neural network is small when the number of samples $n$ is large enough.

\subsection{Algorithm}
Now we introduce the algorithm we construct to find a point which generalizes well. We first introduce several notations. In this subsection, we also rewrite each $\bm{w}_j$ as $c_j\bm{u}_j$ where $c_j=\|w_j\|_2$ and $\bm{u}_j$ is a unit vector. We use $a_j(k)$, $c_j(k)$ and $\bm{u}_j(k)$ to denote the value of $a_j$, $c_j$ and $\bm{u}_j$ at the beginning of the $k$-th outer iteration. We also define vector $\bm{q}_i(k)$ at the $k$-th outer iteration and for the $i$-th sample in the dataset as
\begin{equation}\label{eq::algo-1}
\bm{q}_i(k)=\left(\begin{matrix}y_i\text{sgn}(a_1(k))(\text{sgn}(a_1(k))\bm{u}_1^\top (x_i\odot \phi_1))_+ \\.\\.\\.\\y_i\text{sgn}(a_m(k))(\text{sgn}(a_m(k))\bm{u}_m^\top (x_i\odot \phi_m))_+\end{matrix}\right).
\end{equation}

\textbf{Initialization.} Setting $\lambda_{\max}>\sqrt{n}$ and for each neuron $j=1,...,m$, setting $a_j(0)=c_j(0)=\|\bm{w}_j(0)\|_2\le \sqrt{\frac{n}{4m\lambda_{\max}}}$, setting $\lambda_j(0)=\lambda_{\max}$ and arbitrarily choosing $\bm{u}_j(0)\in\mathbb{S}^{d-1}$.

\textbf{Running.} At each outer iteration $k=1,2,3....$, run:
\begin{itemize}[leftmargin=*]
\item[(1)] \textbf{Updating the regularization coefficients.} Given the values of $\{a_j(k)\}_{j}$ and $\{\bm{u}_j(k)\}_{j}$, computing vectors $\bm{q}_1(k),...,\bm{q}_n(k)$ using the equation~\eqref{eq::algo-1}.  Next, we find a unit vector $\bm{v}(k)$ such that $\bm{v}(k)^\top \bm{q}_i(k)=0$ holds for all $i\in[n]$. We note here that such a vector always exists since $q_i$'s are $m$-dimensional vectors and $m>n+1$. Now we update the coefficient vector $\bm{\lambda}$ in the following way.
\begin{itemize}
    \item If $|\bm{v}^\top(k)\bm{\lambda}(k-1)|\ge \frac{\lambda_{\max}}{144L_n(\bm{\theta}(0))}$, then we set $\bm{\lambda}(k)=\bm{\lambda}(k-1)$.
    \item If $|\bm{v}^\top(k)\bm{\lambda}(k-1)|< \frac{{\lambda_{\max}}}{144{L_n(\bm{\theta}(0))}}$ then
    \begin{itemize}
        \item If $\sum_{j=1}^mv_j^2(k)\mathbb{I}\{v_j(k)>0\}\ge 1/2$, then we define $$\bm{\lambda}(k)=\bm{\lambda}(k-1)-\frac{\lambda_{\max}}{36L_n(0, 0)}\left(\begin{matrix}v_1\mathbb{I}\{v_1(k)>0\}\\.\\.\\.\\ v_m\mathbb{I}\{v_m(k)>0\}\end{matrix}\right).$$
        \item Else, then we define
$$\bm{\lambda}(k)=\bm{\lambda}(k-1)+ \frac{\lambda_{\max}}{36L_n(0, 0)}\left(\begin{matrix}v_1\mathbb{I}\{v_1(k)<0\}\\.\\.\\.\\ v_m\mathbb{I}\{v_m(k)<0\}\end{matrix}\right).$$
    \end{itemize}
\end{itemize}
\item[(2)] \textbf{Finding critical points.} We use the gradient descent algorithm on parameters $\{a_j(k)\}_{j}$ and $\{c_j(k)\}_{j}$ to find points where $a_j=c_j$ and $\|\nabla_{\bm{a}}L_n(\bm{\theta};\bm{\lambda})\|_2$ is small. We note here that during the gradient descent, we fix the unit vectors $\{\bm{u}_j(k)\}_{j}$.  We initialize the gradient  descent algorithm in the following way: setting $\bm{a}(k, 0)=\bm{c}(k, 0)=\bm{a}(k)=\bm{c}(k)$. At the $(t+1)$-th step of the gradient descent algorithm, updating $\bm{a}$ and $\bm{c}$ using the following equations:
$$\bm{a}(k, t+1)=\bm{a}(k, t)-\eta \nabla_{\bm{a}}L_n(\bm{\theta}(k, t)),\quad \bm{c}(k, t+1)=\bm{c}(k, t)-\eta \nabla_{\bm{c}}L_n(\bm{\theta}(k, t))$$
where we choose the step size $\eta=1/(18n)$. We terminate at time $T$ when $$\|\nabla_{(\bm{a},\bm{c})}L_n(\bm{\theta}(k, T))\|_2\le 1/(288n).$$
\item[(3)] \textbf{Perturbing the inactive neurons.} For each index $j\in[m]$ such that $|a_j(k, T)|=|c_j(k, T)|\le \frac{1}{3\sqrt{n}}$,  we are
\begin{itemize}
    \item defining $I_j=\{k\in[d]:\phi_j(k)=1\}$, checking all possible combination of $\{\bm{v}(k)\}_{k\in I_j}\subset \{1/\lceil\sqrt{n}\rceil, 2/\lceil\sqrt{n}\rceil, ...,1\}$ and outputting the vector $\bm{v}_j$ to be the vector $\bm{v}$ maximizing
    $$\left|\sum_{i=1}^n\ell'(-y_i f(x_i;\bm{\theta}(k, T)))y_i (\bm{v}^\top (x_i\odot \phi_j))_+\right|,$$
    \item setting \begin{align*}
{a}_{j}(k+1)&={c}_{j}(k+1)=\text{sgn}\left(\sum_{i=1}^n\ell'(-y_i f(x_i;\bm{\theta}(k, T)))y_i (\bm{v}_j^\top (x_i\odot \phi_j))_+\right)\sqrt{\frac{8\lambda_\max}{9n}},\\
{\bm{u}}_{j}(k+1)&=\bm{v}\text{sgn}\left(\sum_{i=1}^n\ell'(-y_i f(x_i;\bm{\theta}(k, T)))y_i (\bm{v}_j^\top (x_i\odot \phi_j))_+\right)
\end{align*}
For the other neuron with $|a_j(k, T)|=|c_j(k, T)|> \frac{1}{3\sqrt{n}}$, we set $a_j(k+1)=a_j(k, T)$, $c_j(k+1)=c_j(k, T)$ and $\bm{u}_j(k+1)=\bm{u}(k)$.
\end{itemize}

\textbf{Termination.} We terminate the algorithm at the step $K$ if
 $$\max_{j\in m}\left|\sum_{i=1}^n\ell'(-y_i f(x_i;\bm{\theta}(K)))y_i (\bm{u}_j(K)^\top (x_i\odot \phi_j))_+\right|\le 2\lambda_max.$$
\end{itemize}
\fi

\section{Conclusions}

We consider a single hidden layer overparameterized neural network trained for binary classification. Our main results characterize local minima in certain directions which have good generalization performance and present a specialized version of gradient descent to find such points. Unlike much of the prior work on this topic, we are able to handle datasets which are not separable through the use of explicit $\ell_2$ regularization. Interestingly, the amount of overparameterization we need is just one more than the number of points in the dataset. While we do not consider the robustness of the output of the neural network to small changes in the input, it is interesting that the number of neurons we require for the landscape properties to hold in fully connected neural networks is just one more than the number conjectured to be required in \cite{bubeck2020law} for such type of robustness.

\bibliographystyle{plain}


\section{Keywords:} turbulence, outer planets, magnetosphere, weak turbulence, wave-particle interaction, aurora}
\end{abstract}

\vspace*{1cm}

\section{Introduction and origin of turbulence}

All four outer planets of the solar system, Jupiter, Saturn, Uranus and Neptune, possess strong internal dynamo magnetic fields with polar field strength of $\sim$ 2,000,000 nT,  40,000 nT, 100,000 nT and 90,000 nT, respectively
\citep{conn07,ness10}.
They are also fast rotators with periods of 9.9 hrs, 10.7 hrs, 17.2 hrs and 16.1 hrs, respectively \citep{seid02}. These two effects lead to large and rotationally dominated magnetospheres with average magnetopause standoff distances on the sub-solar side of  65, 20, 20, and 26 planetary radii, respectively \citep{vasy09,bage09,bage13}.

One of the most important difference between the magnetospheres of the outer planets and the magnetosphere of Earth is the presence of large mass sources well inside the outer planets' magnetospheres and the resultant effects generated by the associated radial mass transport.
In case of Jupiter, its inner Galilean moon Io, located deeply within the magnetosphere at 6 R$_J$, produces about 10$^3$ kg s$^{-1}$ of SO$_2$, which eventually gets ionized \citep{thom04}. An additional, but smaller source of mass is its moon Europa with about 50 kg s$^{-1}$ of H$_2$O \citep{bage20}. In case of Saturn, the mass source is the moon Enceladus, with a time-variable injection rate of 200-1000 kg s$^{-1}$ of H$_2$O \citep{doug09,saur11,fles10a,hedm13}. In case of Uranus and Neptune, also major moons orbit within their magnetosphere. However, their mass loss rates are not well constrained. They are considered smaller compared to those in the magnetospheres of Jupiter and Saturn based on the significantly smaller mass densities in the magnetospheres of Uranus and Neptune \citep[e.g.,][]{bage09,bage13}.

 \begin{figure}[h!]
\begin{center}
\includegraphics[width=20cm]{figure_turbulence_article_4_0
\end{center}
\caption{Sketch of turbulent wave packages in Jupiter's magnetosphere. The interaction of counter-propagating Alfv\'en wave packages generates a turbulent cascade.
Waves packages are reflected at the ionosphere of Jupiter or other density gradients in the magnetosphere.
}
\label{fig:sketch}
\end{figure}
 The centrifugal force in the fast rotating magnetospheres causes the plasma to be concentrated near the equatorial regions (called plasma sheets, see Figure \ref{fig:sketch}) and drives transport of plasma radially outward. Due to conservation of angular momentum the outward moving plasma does not to fully corotate with the planet. This generates two types of magnetospheric stresses: Stretched magnetic field lines in the radial direction  from the radial transport (see Figure \ref{fig:sketch})  and  bend backed magnetic field lines in azimuthal direction due to the subcorotating plasma. These magnetic stresses couple to the planets' ionospheres and cause transport of angular momentum between the planets' ionospheres and magnetospheres \citep{hill01,gold07}. The stress balance is however not in steady state, but constantly disturbed due to non-continuous radial transport of plasma, which is observed to occur through the interchange of mass loaded magnetic flux tubes with adjacent less loaded flux tubes further out \citep{kive97c}. The imbalance of magnetic stresses cause Alfv\'en waves to propagate along the magnetic field to achieve stress balance. The Alfv\'en wave packages are partially reflected at the ionospheres or other boundaries, such as large density gradients of the plasma sheet. The resulting counter-propagating Alfv\'en wave packages interact nonlinearly and thus generate a turbulent cascade \citep{saur02,saur18a}.
On the MHD scales the magnetospheric plasma sheet is thought to be a key region where the cascade is being driven \citep{saur02}.

The non-linear interaction of counter-propagating wave packages is the fundamental building block of MHD and plasma turbulence. It is established  observationally, theoretically and in numerical modeling \citep[e.g.,][]{gold95,howe13}.  There are several key differences between turbulence in the solar wind \citep[e.g.,][]{brun13}  and the turbulence in the magnetospheres of the outer planets. The magnetospheric fluctuations are small amplitude fluctuations $\delta B \ll B_0$ compared to the dominant planetary background magnetic field $B_0$. The magnetospheres are also  finite in size. This implies that turbulent scales are bound and a maximum length scale for the inertial range of turbulence exists.
It also causes reflection of the wave packages at the ionospheric boundaries. Additionally the magnetospheres are highly inhomogeneous with large densities and small magnetic field strength in the magnetospheric plasma sheets and vice versa at high latitudes (see Figure \ref{fig:sketch}). Therefore the generator region, i.e. where the turbulent cascade is driven, and the dissipation regions of the turbulent fluctuations do not need to be spatially collocated \citep[]{saur18a}.
Another interesting difference between the solar wind and the magnetospheres of the outer planets is that the composition of ions has a much larger spread in mass. The solar wind is mostly composed of H$^+$ and He$^{++}$ compared to H$^+$, O$^+$, S$^+$, H$_2$O$^+$ ... in Jupiter's or Saturn's magnetospheres. This causes a larger spread of gyro frequencies and inertial length scales and will influence the pathways of turbulent dissipations.

In the reminder of this review we will present observations of turbulent fluctuations in Jupiter's and Saturn's magnetospheres. Additionally, we discuss their impact on overall properties of these magnetospheres such as aurora or temperature structure.

\vspace*{1cm}

\section{Magnetospheric Turbulence}
Turbulence and its implications for the magnetospheres of Jupiter and Saturn have been studied by a number of authors.
There are no published studies about turbulence in the magnetospheres of Uranus and Neptune due to a lack of appropriate data.

\vspace*{1cm}

\subsection{Jupiter magnetosphere}
Turbulence in Jupiter's magnetosphere has been investigated on various levels of detail.
Before the Galileo epoch, \citet{nish76} considered whistler mode turbulence
for pitch angle diffusion of energetic electrons and ions. \citet{barb81a} and \citet{barb84}  investigated ion heating due to high-frequency
kilohertz plasma wave turbulence and low frequency Alfv\'enic turbulence, respectively.
\citet{glas95} suggests a Kolmogorov-type of turbulence
in the very low-frequency milli-Hertz range based on Voyager magnetic field measurements.

Magnetic field measurements by the Galileo orbiter led to a dedicated turbulence analysis by \citet{saur02}. The analysis covered the magnetic field fluctuations in Jupiter's middle magnetosphere (9-24 R$_J$) within the magnetohydrodynamic (MHD) range between  2.8 $\times$ 10$^{-4}$ -- 3.6 $\times$ 10$^{-2}$ Hz, i.e., below the ion cyclotron frequency. Even though the relative velocity between the magnetospheric plasma and the spacecraft is much smaller compared to spacecraft in the solar wind, the Taylor theorem can still be applied \citep{saur13}. The properties of the magnetic field fluctuations in Jupiter's middle magnetosphere differ from those in the solar wind. In the Jovian system, the ratio of the magnetic fluctuations $\delta B$ compared to the background magnetic field $B_0$ is less than 10$^{-1}$, while it is of order unity in the solar wind.
A weak turbulence cascade is expected in Jupiter's middle magnetosphere based on the ratio $\epsilon =\tau_A / \tau_{nl}$ between the Alfv\'en time $\tau_A$ and the non-linear interaction time $\tau_{nl}$. This ratio $\epsilon$ lies in the range of 0.06 to 0.34 and is indicative of weak turbulence \citep{ng97,galt00}.  Magnetic field fluctuations parallel to the background magnetic field exhibit a power law with
spectral index of -2 \citep{saur02} also indicative of weak turbulence \citep{ng97,galt00}. The other components often show deviations from clear power laws possibly due to other perturbations in the Jupiter's highly dynamic magnetosphere or locally non-fully developed turbulent cascades.

Turbulence in Jupiter's magnetosphere has subsequently been analyzed by \citet{tao15} who combined low and high-time resolution magnetic field data of the Galileo spacecraft within MHD and kinetic scales, i.e., within 1 $\times$ 10$^{-4}$ Hz to 1 Hz. The resultant magnetic field spectra show spectral breaks well organized by the ion scales (i.e., ion cyclotron frequency, the Doppler-shifted ion inertial length scales and the ion gyroradius, which assume similar values in Jupiter's middle magnetosphere). For frequencies smaller than those associated with the ion scales, the spectral indices lie in the range of -0.6 and -1.9 and for higher frequencies within the range of -1.7 and -2.5 \citep{tao15}. The authors also show that the turbulence is intermittent, in particular in the equatorial region of Jupiter's magnetospheric plasma.

\vspace*{1cm}

\subsection{Saturn magnetosphere}

Saturn's magnetosphere also hosts small scale turbulent magnetic field fluctuations as demonstrated by \citet{vonp14}. The average amplitudes of the fluctuations compared to the magnetospheric background field assumes values of $\delta B / B$ = 0.07. An exemplary spectrum is shown in Figure \ref{fig:spectrum}.  On MHD scales the spectral slope of the fluctuations perpendicular and parallel to the background magnetic field varies between -0.8
 and -1.7.
\begin{figure}[h!]
\begin{center}
\includegraphics[width=15cm]{Figure_von_Papen.pdf}
\end{center}
\caption{
 Power spectral density (PSD) of turbulent magnetic field fluctuations in Saturn's magnetosphere (from \citet{vonp14}). $P_\perp$ and $P_\parallel$ refer to the power spectral densities of fluctuations perpendicular and parallel to the background magnetic field, respectively. $\lambda_{W^+}$ and
$\rho_{W^+}$indicate the Doppler shifted  ion inertial length scale and gyroradius of water ions, respectively, and
$f_{c, W^+}$  and $f_{c,H^+}$ are the cyclotron frequencies of water ions and protons, respectively.}
\label{fig:spectrum}
\end{figure}
 %
At ion scales, i.e., near the ion inertial length scale or the Doppler-shifted ion gyroradius, a spectral break occurs, with a sub-ion, i.e., high-frequency, slope of -2.3 inside of radial distances of 9 $R_S$ (see Figure \ref{fig:spectrum}). Further outside in the magnetosphere the sub-ion spectrum steepens to -2.6. These slopes could be consistent with turbulence of kinetic Alfv\'en waves. The probability density function of the fluctuations have non-Gaussian tails and a power law increase of the flatness indicates that the turbulence is intermittent \citep{vonp14}.  The turbulence on the kinetic scales is estimated to be strong in Saturn's magnetosphere with $\epsilon > 1$, i.e., the non-linear time is shorter than the wave propagation time, which could result from Saturn's weaker planetary magnetic field compared to Jupiter's field.

In a follow-up study,  \citet{vonp16} investigated the spatial and temporal structure of the turbulence in Saturn's magnetosphere. In local time coordinates, enhanced fluctuations are seen at noon, possibly resulting from flux tube interchanges \citep[e.g.,][]{kive97c}, which are considered a prime source of the turbulence. In a frame rotating with Saturn, increased fluctuations are seen at 65$^{\circ}$ southern and 250$^{\circ}$ northern magnetic phase. The later correlation is enigmatic and is related to the still unresolved planetary period oscillations of Saturn \citep{espi00,prov09,gurn09,vonp16}. The oscillations might result from a coupling between Saturn's atmosphere, ionosphere and magnetosphere. The variability in the turbulence can thus be considered as an effect of the Alfv\'enic energy fluxes underlying these coupling processes. The spatial variability of the turbulent fluctuations in Saturn's magnetosphere has been subsequently studied by  \citet{kami17}. The authors also found enhanced fluctuations between 10 hr and 20 hr local time and the turbulence to be more quiet between 3 hr and 9 hr local time.

\vspace*{1cm}

\subsection{Large scale implications of turbulence}

In addition to serving as a laboratory for turbulence studies, the fluctuations in the outer planets' magnetospheres have major implications on overall properties of their planetary systems.

\vspace*{1cm}

\subsubsection{Magnetospheric heating}

The magnetospheres of Jupiter and Saturn possess ion temperatures which strongly increase
from the plasma source regions radially outwards by approximately two orders of magnitude
to values up to nearly 10$^8$ K \citep{bage11}.
These observations are considered a long standing puzzle since major energy input rates are required to explain these temperature increases. Otherwise the radial transport of magnetospheric plasma would lead to adiabatic cooling with distance.
Dissipation of turbulent magnetic field fluctuations is a very powerful energy source which can explain the temperature increase.

First suggestions of the importance of turbulent heating in Jupiter's magnetosphere in addition to pickup energization and subsequent radial transport date back to \citet{barb84}. Based on the observed spatial distribution of the magnetic field fluctuations and a model for a weak turbulent cascade, \citet{saur04a} calculated the dissipation rate in Jupiter's middle magnetosphere to a total amount of 5 $\times $ 10$^{12}$ W
and demonstrated that this rate and its spatial distribution can explain the observed radial temperature profile
\citep{fran02a}. To a similar conclusion came  \citet{ng18} who used a radial transport model based on advection  and not on flux tube diffusion as in \citet{saur04a}. The turbulent dissipation rates are consistent with the required heating rate independently extracted from Galileo spacecraft  observations by \citet{bage11}.  A possible dissipation mode could be ion cyclotron damping  in particular of the heavy sulfur and oxygen ions of Jupiter's magnetosphere \citep{saur18a}.
In Saturn's magnetosphere the dissipation of the magnetic field fluctuations based on a strong turbulence model could provide power
 on the order of 10$^{11}$ W \citep{vonp16,kami17,neup21}, which is roughly consistent with the energy needed to heat the magnetosphere to its observed temperatures \citep{bage11}.

\vspace*{1cm}

\subsubsection{Aurora}

Jupiter's main auroral oval lies on magnetic field lines that map to an equatorial distance of about 20 - 30 R$_J$ \citep[e.g.,][]{clar98,clar02,hill01}
The turbulent power of the magnetic fluctuations in Jupiter's magnetosphere maximizes at the same radial distance \citep{saur03}. This  distance is also the region where the corotation of Jupiter's magnetosphere breaks down and magnetospheric-ionospheric  coupling currents maximize
as derived by
\cite{hill79,hill01,cowl01}. The turbulent fluctuations result from the non-steady radial transport and the resultant imbalanced stresses between the magnetosphere and Jupiter's ionosphere. The stress balance is achieved by Alfv\'en waves propagating between both regions. When the Alfv\'en waves reach the polar regions of the magnetosphere, the perpendicular length scales of the Alfv\'en waves grow small due to Jupiter's increasing magnetic field. Just above the ionosphere, the electron inertial length scale assumes its large values of $\sim$50 km due to the low electron density. When the turbulent cascade reaches this scale, the fluctuations are subject to electron Landau damping resulting in stochastically accelerated electrons \citep{saur18a}. This process is consistent with the bi-directional and energetically broad band electron distributions  recently measured by the Juno spacecraft \citep[e.g.,][]{mauk17,clar18}. These electrons are one of the key contributors to Jupiter's main auroral oval. Magnetic field fluctuations measured by the Juno spacecraft near the auroral particle acceleration region showed a spectrum with a slope of -2.29 $\pm$0.09 within 0.5 -- 5 $\times$ 10$^{-1}$ Hz \citep{gers19}. They could be due to kinetic/inertial Alfv\'en waves, however, the underlying turbulent processes for this spectrum are not fully clear yet as the Alfv\'en velocity is close to the speed of light in this region. The aurora of Saturn is less well understood, but also within Saturn's magnetosphere bi-directional energetically broad band electron beams have been observed on field lines which are connected to the aurora \citep{saur06,mich09} and which appear to be driven by stochastic, turbulent processes \citep{vonp16}.

\vspace*{1cm}

\subsection{Turbulence in Alfv\'en wings of the Moons}
A new and evolving area of turbulence research in the outer planets' magnetosphere are the Alfv\'en wings generated by the moons. In Figure \ref{fig:sketch}, we exemplarily show Jupiter's moon Io. The moon is an obstacle to the rotating magnetospheric plasma and  generates Alfv\'en waves, which propagate along the rotating background magnetic field towards Jupiter. In the rest frame of each moon two standing waves develop (north and south, respectively).
They are referred to as Alfv\'en wings \citep{neub80,neub98,goer80,saur13}. The moons can thus be considered as gigantic Alfv\'en wave antenna. The Alfv\'en wings  are reflected at gradients of the background density and generate filamented structures \citep{chus05,hess11}. The magnetic field fluctuations in Io's Alfv\'en wings at high latitudes have been recently measured by the Juno spacecraft. Their spectrum is consistent with a power law spectral index of -2.35 $\pm$0.09 within 0.1 -- 800 Hz \citep{sula20}.
When the filamented waves approach kinetic scales, wave-particle interaction sets in and can energize electrons through electron Landau damping and ions through cyclotron damping \citep{hess10,szal18,szal20,saur18a}. Similar particle acceleration processes have also been observed in the Alfv\'en wings of Europa \citep{alle20} and Ganymede \citep{szal20a}.

\vspace*{1cm}

\subsection{Magnetosheaths' of Jupiter and Saturn}
Turbulence in the magnetosheaths of the outer planets has received very little attention to date. Turbulence in the magnetosheaths is distinct from turbulence within the outer planets' magnetospheres. The magnetosheaths are characterized by large amplitude magnetic fluctuations $\delta B/ B_0 \sim 1$ \citep{hadi15}. The ultimate source of free energy stems from the shocked solar wind plasma with Alfv\'en Mach numbers as high as $\sim$100 \citep{mast13},  which causes enhanced plasma temperatures and large plasma beta turbulence.

\cite{alex08} analyzed turbulent magnetic field fluctuations in Saturn's magnetosheath obtained the Cassini spacecraft and identified  Alfv\'en vortices similar to those found in Earth's magnetosheath.
Alfv\'en vortices are non-linear magnetic structures associated with sheared velocity perturbations propagating obliquely to the external field.  \citep{hadi15} also studied Saturn's magnetosheath and found power-spectral energy densities of the magnetic field fluctuations with a $\sim$$f^{-1}$ scaling at MHD scales and an $\sim$$f^{-2.6}$ scaling at sub-ion scales.  The authors also found based on the compressibility of magnetic and density fluctuations that compressible magnetosonic  slow modes at MHD scales dominate rather than the Alfv\'en mode. Magnetic field turbulence in Jupiter's magnetosheath was studied by \cite{bolz14} based on Ulysses measurements. \cite{bolz14} considered a similar frequency range as \cite{hadi15}, but fitted the power spectral density to only one power law instead of two and thus found an average slope of around $-2$.

\vspace*{1cm}

\section{Discussion: Outstanding issues and outlook}

For a more comprehensive understanding of turbulence in the outer planets' magnetospheres further measurements and studies are required. The Galileo spacecraft recorded only a small fraction of its mission in a high-frequency modus. The Juno spacecraft is currently exploring the high latitude regions of Jupiter's magnetosphere including the moons' Alfv\'en wings. With its instrumentation and its polar orbit detailed studies of turbulence and the associated wave-particle processes of turbulent dissipation are on their way.

New missions with appropriate instrumentation to Uranus and Neptune are being discussed \citep[e.g.,][]{arri14}, but no mission has been decided upon.
The orientations of the planetary magnetic moments at Uranus and Neptune are highly tilted compared to the planets' spin axis. Also the plasma densities and the magnetic field fluctuations in these magnetospheres are much smaller compared to those of Jupiter and Saturn. It will be interesting to see how these properties influences the possible turbulent nature of the field and plasma fluctuations.

Further studies are warranted on where and what types of turbulent cascades are driven in the magnetospheres of the outer planets. For example, are there cascades on kinetic scales outside the plasma sheets? Also a better understanding of the nature of the dissipation mechanisms of the turbulent fluctuations and the resultant ion or electron distribution functions is an important matter. Further observational studies might aim to better resolve the effects of the various dissipation and turbulent break scales expected to occur at the inertial lengths scales and the gyroradii of the electrons and the various ion species, respectively. Because the magnetospheres are highly inhomogeneous, the spatial variability of these parameters and their effects can be more easily studied compared to the solar wind.
Suggestions for the roles of Landau and ion cyclotron damping of Alfv\'enic turbulence in the magnetosphere of Jupiter have been made \citep[e.g.,][]{saur18a}. However, the role of other wave modes or non-resonant wave-particle interaction mechanisms warrant further analysis.

Numerical simulations with hybrid or particle-in-cell models can be a useful tool to better understand the dissipation mechanisms and the associated transport of particles and energy \citep[e.g.,][]{dela21}. Future models could focus on more
realistic magnetic field geometries with current sheets and inhomogeneous plasma densities along field lines.

Turbulence in the magnetosheaths of the outer planets have not received much attention. Dedicated  study might investigate the structure of turbulence downstream of the planets' bow shocks to test the universality of turbulent evolution in Earth's and other planets' magnetosheaths. In the much larger magnetosheath of the outer planets' turbulence might evolve toward a more fully developed turbulent states compared to Earth.
While turbulence in the middle magnetospheres of Jupiter and Saturn is unique due to the plasma of the moons and the fast rotating planets, the existence of turbulence in the magnetospheric tails and the magnetosheaths could be universal across all magnetized planets.

Turbulent fluctuations might  be used as a diagnostic tool to probe the interior of the icy moons in the outer solar system. Icy moons, such as Europa or Ganymede, possess subsurface water oceans \citep{kive00,kive02,saur10,saur15}. These saline and thus electrically conductive oceans have been detected through electromagnetic induction caused by time-variable  magnetospheric fields. For existing observations only the $\sim$10 hr rotation period of Jupiter could be applied. Turbulent fluctuations provide a broad range of frequencies which can be used in future induction studies to additionally probe the interior of the moons  (e.g., with JUICE or Europa Clipper data).

In summary, compared to the solar wind or the interstellar medium, the magnetospheres of the outer planets are complementary laboratories to study plasma turbulence.
These magnetospheres are distinct due the planets' large background magnetic fields, the bound and inhomogeneous nature of the turbulent systems, where the location of the turbulent cascade and the turbulent dissipation are not necessarily collocated.

\vspace*{1cm}

\section*{Conflict of Interest Statement}

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

\section*{Author Contributions}

J.S. selected the material for this review and wrote the manuscript.

\section*{Funding}
This project has received funding from the European Research Council
(ERC) under the European Union’s Horizon 2020 research and innovation
programme (grant agreement No. 884711).

\section*{Acknowledgments}
J.S. thanks Sascha Janser for helpful comments on the manuscript  and Filip Elekes for the graphics support.

\bibliographystyle{frontiersinSCNS_ENG_HUMS}