




\section{Tsaknakis-Spirakis Algorithm in Detail}\label{app:TSAlgo}
As is stated in the main body, the TS algorithm consists of two parts, a descent algorithm and an adjustment step. We state the TS algorithm in three parts: Calculating the steepest direction, descent algorithm and adjustment step. For interested readers, see \cite{DBLP:journals/corr/abs-2107-01471} and \cite{DBLP:journals/im/TsaknakisS08} for a thorough description and analysis. Again, in this section, we always abbreviate Dini stationary points as stationary points.
\subsection{Calculating the Steepest Direction}
We first present formulas of Dini directional derivatives. Let $S_R(y):=\suppmax(Ry)$ and $S_C(x):=\suppmax(\T{C}x)$. It has been shown by \cite{DBLP:journals/im/TsaknakisS08} that 
\begin{align*}
    Df_R(x,y,x',y')&=\max_{S_R(y)}(Ry')-\T{x'}Ry-\T{x}Ry'+\T{x}Ry-f_R(x,y),\\
    Df_C(x,y,x',y')&=\max_{S_C(x)}(\T{C}x')-\T{x'}Cy-\T{x}Cy'+\T{x}Cy-f_C(x,y).
\end{align*}
If $f_R(x,y)\neq f_C(x,y)$, we have
\[Df(x,y,x',y')=\begin{cases}
Df_R(x,y,x',y'),& f_R(x,y)>f_C(x,y),\\
Df_C(x,y,x',y'),& f_R(x,y)<f_C(x,y).
\end{cases}\]
If $f_R(x,y)=f_C(x,y)$, we have
\begin{align*}
    Df(x,y,x',y')&=\max\{Df_R(x,y,x',y'),Df_C(x,y,x',y')\}\\
    &=\max\{T_1(x,y,x',y'),T_2(x,y,x',y')\}-f(x,y),
\end{align*}
where
\begin{align*}
T_1(x,y,x',y')&=\max_{S_R(y)}(Ry')-\T{x'}Ry-\T{x}Ry'+\T{x}Ry,\\ T_2(x,y,x',y')&=\max_{S_C(x)}(\T{C}x')-\T{x'}Cy-\T{x}Cy'+\T{x}Cy.
\end{align*}
We introduce linear convex combinations via $\rho, w$ and $z$ to smoothen the maximum operators:
\begin{align*}
T(x,y,x',y',\rho,w,z):=&\rho(\T{w}Ry'-\T{x}Ry'-\T{x'}Ry+\T{x}Ry)\\
         &\quad+(1-\rho)(\T{x'}Cz-\T{x}Cy'-\T{x'}Cy+\T{x}Cy),
\end{align*}
where $\rho\in[0, 1]$, $w\in\Delta_m$, $\supp(w)\subseteq S_R(y)$, $z\in\Delta_n$, $\supp(z)\subseteq S_C(x)$.

We assume that $f_R(x,y)=f_C(x,y)$. Then calculating the steepest direction becomes a minimax problem:
\begin{gather*}
    \min_{x',y'}\max_{\rho,w,z} T(x,y,x',y',\rho,w,z),\\
    \begin{split}
        \text{s.t. } (x',y'),(w,z)&\in\Delta_m\times\Delta_n,\\
    \rho&\in[0,1],\\
    \supp(w)&\subseteq S_R(y),\\
    \supp(z)&\subseteq S_C(x).
    \end{split}
\end{gather*}
\cite{DBLP:journals/corr/abs-2107-01471} and \cite{DBLP:journals/im/TsaknakisS08} show that this minimax problem can be solved by dual linear programmings. Suppose $(\rho_0,w_0,z_0)$ is the optimal solution of the dual LP. This tuple is called a \emph{dual solution}.
\subsection{Descent Algorithm}
The descent algorithm is divided into 2 steps: seeking the descent direction and line search.

Seek the descent direction, i.e., to find the steepest direction of $f$. To do so, first fix one of $x,y$ and adjust the other to make $f_R=f_C$; then the problem goes back to the first part.

Having obtained the steepest direction $(x' - x,y' - y)$, the algorithm searches for a proper step size in that direction. This step is called line search in optimization theory.

The pseudo-code of the descent algorithm is presented in \Cref{algo:dp}.
\vspace{0.5em}
\begin{breakablealgorithm}
\renewcommand{\algorithmicrequire}{\textbf{Input}}
\newcommand{\OUTPUT}[1]{\textbf{Output} #1}
\newcommand{\GOTO}[1]{\textbf{goto} #1}
\renewcommand{\algorithmiccomment}[1]{// \textit{#1}}
\caption{Searching for a $\delta-$stationary point (\cite{DBLP:journals/corr/abs-2107-01471}, Appendix B)}
\label{algo:dp}
\begin{algorithmic}[1]
\REQUIRE Precision $\delta$, payoff matrices $R_{m\times n},C_{m\times n}$ and initial strategy $(x,y)\in\Delta_m\times\Delta_n$.
\bigskip
\STATE \COMMENT{seeking descent direction.}
\IF{$f_R(x,y)\neq f_C(x,y)$}\label{algline:DP1}
\IF{$f_R(x,y)>f_C(x,y)$}
\STATE Fix $y$ and solve the following LP with respect to $x$:
\begin{align*}
    &\min\limits_{x}\{\max(Ry)-\T{x}Ry\}, \\
    \text{s.t. }&f_R(x,y)\geq f_C(x,y), \\
    & x\in\Delta_m.
\end{align*} 
\ENDIF
\medskip
\IF{$f_R(x,y)<f_C(x,y)$}
\STATE Fix $x$ and solve the following LP with respect to $y$:
\begin{align*}
    &\min\limits_{y}\{\max(\T{C}x)-\T{x}Cy\}, \\
    \text{s.t. }&f_C(x,y)\geq f_R(x,y), \\
    & y\in\Delta_n.
\end{align*}

\ENDIF
\ENDIF
\medskip
\STATE Solve the following minimax problem with respect to $(x',y')$:
\begin{align*}
    &\min_{x',y'}\max_{\rho,w,z} T(x,y,x',y',\rho,w,z), \\
    \text{s.t. }& (x',y')\in\Delta_m\times\Delta_n, \\
    & \rho\in[0,1],\ \supp(w)\subseteq S_R(y)\text{ and }\supp(z)\subseteq S_C(x).
\end{align*}
\STATE Let the optimal value of LP be $V$ and the corresponding solution be $(x',y')$ and $(\rho, w, z)$.
\medskip
\IF{$V-f(x,y)\geq-\delta$}
\STATE \OUTPUT{$(x,y)$}\RETURN
\ENDIF
\bigskip
\STATE \COMMENT{line search.}
\STATE $\bar{S}_R(y)\leftarrow[m]\setminus S_R(y)$, $\bar{S}_C(x)\leftarrow[n]\setminus S_C(x)$.
\STATE 
\begin{align*}
    &\epsilon_1^*\leftarrow\min_{i\in\bar{S}_R(y)}\left(\frac{\max(Ry)-(Ry)_i}{\max(Ry)-(Ry)_i+(Ry')_i-\max_{S_R(y)}(Ry')}\right), \\
    &\epsilon_2^*\leftarrow\min_{j\in\bar{S}_C(x)}\left(\frac{\max(\T{C}x)-(\T{C}x)_j}{\max(\T{C}x)-(\T{C}x)_j+(\T{C}x')_j-\max_{S_C(x)}(\T{C}x')}\right), \\
    &\epsilon^*\leftarrow\min\{\epsilon_1^*,\epsilon_2^*,1\}.
\end{align*}
\STATE $H\leftarrow\min\{\T{(x'-x)}R(y'-y),\T{(x'-x)}C(y'-y)\}$.
\medskip
\IF{$H<0$}
\STATE $\epsilon^*\leftarrow\min\{\epsilon^*,|V-f(x,y)|/(2|H|)\}$.
\ENDIF
\medskip
\STATE $x\leftarrow x+\epsilon^*(x'-x)$, $y\leftarrow y+\epsilon^*(y'-y)$.
\STATE \GOTO{\Cref{algline:DP1}}
\end{algorithmic}
\end{breakablealgorithm}
The descent algorithm will find a $\delta$-stationary point in $O(\delta^{-2})$ steps and every step takes a time of $\mathrm{poly}(m,n)$. Thus the descent algorithm is a polynomial-time algorithm.

\subsection{Adjustment Step}
After finding a ($\delta$-)stationary point $(x^*,y^*)$ and the corresponding dual solution $(\rho^*,w^*,z^*)$, the algorithm makes a further adjustment. The algorithm then outputs the better solution between the stationary point and the adjusted strategy. The original TS algorithm in \cite{DBLP:journals/im/TsaknakisS08} presented one, and \cite{DBLP:journals/corr/abs-2107-01471} gave another two better adjustment methods. We state them below:
\begin{enumerate}[fullwidth, label=\textbf{Method \arabic*}. ]
    \item \begin{interposition}[(\cite{DBLP:journals/im/TsaknakisS08}, Section 4)]
Define
\begin{align*}
    \lambda&:=\min_{y':\supp(y')\subseteq S_C(x^*)}\{\T{(w^*-x^*)}Ry'\},\\
    \mu&:=\min_{x':\supp(x')\subseteq S_R(y^*)}\{\T{x'}C(z^*-y^*)\}.
\end{align*}
The adjusted strategy pair is
\begin{align*}
\AdjFi:=\begin{cases}\left(\frac{1}{1+\lambda-\mu}w^*+\frac{\lambda-\mu}{1+\lambda-\mu}x^*,z^*\right),&\lambda\geq\mu,\\
\left(w^*,\frac{1}{1+\mu-\lambda}z^*+\frac{\mu-\lambda}{1+\mu-\lambda}y^*\right),&\lambda<\mu.\end{cases}
\end{align*}
\end{interposition}
\item \begin{interposition}[(\cite{DBLP:journals/corr/abs-2107-01471}, Section 3, Method 2)]
Define
\begin{align*}
    \alpha^*&:=\argmin_{\alpha\in[0,1]} f(\alpha w^*+(1-\alpha)x^*,z^*),\\
    \beta^*&:=\argmin_{\beta\in[0,1]} f(w^*,\beta z^*+(1-\beta)y^*).
\end{align*}
The adjusted strategy pair is 
\begin{align*}
  \AdjSe:=\begin{cases}\left(\alpha^*w^*+(1-\alpha^*)x^*,z^*\right),& f_C(w^*,z^*)\geq f_R(w^*,z^*),\\
\left(w^*,\beta^*z^*+(1-\beta^*)y^*\right),& f_C(w^*,z^*)<f_R(w^*,z^*).
\end{cases}  
\end{align*}
\end{interposition}
\item \begin{interposition}[(\cite{DBLP:journals/corr/abs-2107-01471}, Section 3, Method 3)]
Define\footnote{The denominator of $p^*$ or $q^*$ may be zero. In this case, simply define $p^*$ or $q^*$ to be $0$.}
\begin{align*}
p^*&:=\frac{f_R(x^*,z^*)}{f_R(x^*,z^*)+f_C(w^*,z^*)-f_R(w^*,z^*)},\\
q^*&:=\frac{f_C(w^*,y^*)}{f_C(w^*,y^*)+f_R(w^*,z^*)-f_C(w^*,z^*)}.
\end{align*}
The adjusted strategy pair is defined as
\begin{align*}
  \AdjTh:=\begin{cases}\left(p^*w^*+(1-p^*)x^*,z^*\right),& f_C(w^*,z^*)\geq f_R(w^*,z^*),\\
\left(w^*,q^*z^*+(1-q^*)y^*\right),& f_C(w^*,z^*)<f_R(w^*,z^*).
\end{cases}  
\end{align*}
\end{interposition}
\end{enumerate}
All three methods take a polynomial time. \cite{DBLP:journals/corr/abs-2107-01471} shows that all three methods can attain the lower bound $b\approx 0.3393$. However, Method 2 performs the best in non-worst cases.
\section{Introduction}
Plenty of problems in different fields of computer science and mathematical programming are formalized into optimization problems, often nonconvex and nonsmooth ones. As R.T.Rockafellar \cite{RockafellarQuote} stated, ``In fact the great watershed in optimization isn't between linearity and nonlinearity, but convexity and nonconvexity.'' Convex optimization can be solved by systematical approaches. However, different nonconvex problems have different kinds of nonconvexity. Challenges of stability in nonconvex and nonsmooth optimizations arise widely in applications:
\begin{itemize}[fullwidth, listparindent=\parindent]
    \item \emph{Deep neural networks: geometric characterizations of the loss landscape and stability of nearly well-trained neural networks}.
    
    The triumph of AlexNet by Krizhevsky, Sutskever and Hinton \cite{DBLP:conf/nips/KrizhevskySH12} on the classification problem of ImageNet \cite{DBLP:conf/cvpr/DengDSLL009} stimulates a burst of empirical studies for deep neural networks. Deep neural networks show a powerful prediction capability just like some black magic. However, many difficulties arise when we train a neural network. Neural networks are classic nonconvex (and probably nonsmooth) functions. To resolve these difficulties, understanding nonconvex and nonsmooth optimization is a key step. The first question is the geometric characterization of loss landscape, especially the characterization of stationary points. The answer on smooth\footnote{A smooth function is the one differentiable of any order at every point.} models has been given by Cooper's works \cite{DBLP:journals/corr/abs-2005-04210,DBLP:journals/corr/abs-1804-10200}. However, there is no such result for nonsmooth models. Technical barriers exist to generalizing Cooper's works to nonsmooth cases. The second question arises from understanding the stability of nearly well-trained neural networks. Such understanding is crucial to fine-tuning and initialization. For example, can we guarantee that the gradient descent method finds a global minimum point if we initialize a neural network near the same global minimum point? Studies on such stability are scarce even for smooth models.
    
    \item \emph{Parametric stability and non-escapable solutions in continuous optimization}.
    
    Analysis for nonsmooth optimizations has been developed by Clarke's works \cite{clarke1973necessary, clarke1990optimization} and Rockafellar's work \cite{Rockafellar+2015} in the 1970s. Many optimization algorithms for nonsmooth objective functions with some regularity conditions have been devised: standard subgradient methods \cite{Shor1985}, bundle methods \cite{KiwielBundle, LemarechalBundle} and  gradient sampling methods \cite{GradientSamplingBurke2002, BurkeRobustGradientSampling2005}. Convergence of many methods in nonconvex situations has been proved in the past decades. Numerous works focus on parametric stability for specific problems \cite{Ikeda1991,KOJIMA198093, zlobec2001stable}. Well-studied topics of parametric stability include condition numbers \cite{trefethen1997numerical} and sensitivity analysis for linear programmings \cite{jansen1997sensitivity} in the last century. On the other hand, a non-escapable solution means that once the optimization procedure goes nearby, it will get close to that solution no matter the trajectory of the procedure. Such solutions are tolerant to precision parameter and algorithmic numerical errors. However, for general nonconvex and nonsmooth optimization problems, parametric stability and non-escapable solutions are barely studied. 
    
    \item \emph{Stability analysis of the Tsaknakis-Spirakis algorithm in approximate Nash equilibria}.
    
    Ever since the fundamental existence proof by Nash \cite{nash1951non}, Nash equilibria have become the key solution concept in game theory. Computer scientists have made great efforts on Nash equilibrium computations. Finding an approximate solution in general was proved to be PPAD-complete for 3NASH by Daskalakis, Goldberg and Papadimitriou \cite{DBLP:journals/im/TsaknakisS08} and for 2NASH by Chen, Deng and Teng \cite{DBLP:journals/jacm/ChenDT09}. It is conjectured that PPAD-hard problem cannot be solved efficiently. These results lead to many algorithms finding an $\epsilon$-approximate Nash equilibrium in polynomial time for some small constant $\epsilon>0$. For bi-matrix games,  a direct approach is to view approximation ratio $\epsilon$ as an objective function with players' strategies as the input. Then the state-of-the-art approximation algorithm proposed by Tsaknakis and Spirakis \cite{DBLP:journals/im/TsaknakisS08}\footnote{We call it TS algorithm for short.} follows. The TS algorithm reaches an approximation ratio of roughly $0.3393$, proven tight recently by \cite{DBLP:journals/corr/abs-2107-01471}. The empirical studies \cite{DBLP:journals/corr/abs-2107-01471,DBLP:conf/wea/FearnleyIS15,DBLP:conf/wine/TsaknakisSK08} on the TS algorithm show a far better approximation performance than the theoretical ratio. Whether there is a more robust and practical approximation ratio remains blank.
\end{itemize}

Facing these challenges, we consider optimization problems of minimizing a continuous (probably nonsmooth and nonconvex) function $f: \Omega\to\mathbb{R}$, where $\Omega$ is an area on $\mathbb{R}^n$. We demand that at every point in $\Omega$, all possible unit direction vectors form a closed set. Function $f$ is required to be directionally differentiable at every point $x\in \Omega$ in every possible direction. We also require that directional derivatives are continuous in unit direction vectors. Such functions, forming a function class $\mathscr{F}$, cover a wide range of nonconvex and nonsmooth problems arising from various areas. A \emph{stationary point} is a point with all possible directional derivatives nonnegative. We focus on two notions of stability on stationary points of $f\in\mathscr{F}$: \emph{parametric stability} and \emph{convergence stability}. They both characterize the local stability of a stationary point but from different perspectives. The word ``local'' means that the stability we study is restricted to small perturbations in a small neighborhood of a stationary point.

Parametric stability is, literally, a stability concept about parameters. Suppose $f\in\mathscr{F}$ can be described by a parameter $r$. A stationary point is \emph{parametric-stable} if a small perturbation on $r$ always leads to a small change in both the location of the new stationary point and its $f$ value. One widely studied model on parameter perturbations is smoothed analysis, proposed and developed by Spielman and Teng \cite{DBLP:journals/jacm/SpielmanT04}. By considering computational complexities, smoothed analysis becomes a beautiful method to explain why the simplex algorithm for linear programming usually takes polynomial time \cite{DBLP:journals/jacm/SpielmanT04}. In addition, Chen, Deng and Teng \cite{DBLP:conf/focs/ChenDT06} shows that no algorithm computing a Nash equilibrium in a two-player game algorithm runs in polynomial time in a smoothed meaning unless $\mathrm{PPAD}\subseteq \mathrm{RP}$. Parametric stability extends such an idea to general optimization problems and we here consider stability concepts rather than complexity. Similar concepts of stability considered in the literature are numerical stability \cite{atkinson2005theoretical,higham2002accuracy}, condition numbers \cite{trefethen1997numerical} and sensitivity analysis for linear programmings \cite{jansen1997sensitivity}. However, parametric stability studied in this work applies to more general continuous optimization problems.

Convergence stability is a stability concept about approximate stationary points. Most optimization algorithms can only find an approximate stationary point with a small absolute gradient value of $\delta$ (called $\delta$-stationary point) due to limited time and precision. Convergence stability is motivated from two natural stability problems:
\begin{enumerate}
    \item Do $\delta$-stationary points get close to an exact stationary point as $\delta\to 0$?
    \item Suppose the optimization algorithm $\mathcal{A}$ finds a $\delta$-stationary point as its solution. If we perturb the solution a little and run $\mathcal{A}$ starting from the perturbed point, will $\mathcal{A}$ still find a nearby approximate stationary point?
\end{enumerate}
These two problems elicit a more formal description of convergence stability. We say a stationary point $x^*$ is \emph{convergence-stable} for optimization algorithm $\mathcal{A}$ if $\mathcal{A}$ almost surely finds a $\delta$-stationary point arbitrarily close to $x^*$ for sufficiently small $\delta$ when we run $\mathcal{A}$ with any initial point near $x^*$.  

To get nontrivial results on convergence stability, we focus on the optimization algorithm $\mathcal{A}$ whose iteratively generated point sequence $\{x_k\}_{k\geq 0}$ starting from arbitrary initial point $x_0$ almost surely satisfies
\begin{itemize}
    \item \emph{ultimately decreasing}, that is, if $\bar{x}$ is a limit point\footnote{We here means there are infinitely many $k$ such that $x_k$ belongs to any neighborhood of $\bar{x}$. If $\{x_k\}$ is finite, stopping at $k_0$, then we assume $x_k=x_{k_0}$ for all $k>k_0$.} of $\{x_k\}$, then $f(\bar{x})\leq f(x_0)$,
    \item \emph{$\delta$-result-stationary}, that is, $\mathcal{A}$ finds a $\delta$-stationary point as a limit point of $\{x_k\}$,
    \item \emph{$\lambda$-path-bounded}, that is, there exists $\lambda>0$ such that for any $ k=0,1,\cdots$ and any $\theta\in[0,1]$, $f(\theta x_{k+1}+(1-\theta)x_{k})\leq f(x_0)+\lambda$.
\end{itemize}
These conditions cover a wide range of optimization algorithms, including gradient descent, bundle methods, subgradient methods and gradient sampling methods.\footnote{For lots of descent optimization algorithms, ultimately decreasing and $\delta$-result-stationary is proved. $\lambda$-path-bounded can be guaranteed by setting a sufficiently small step size.}

For typical functions in $\mathscr{F}$, we give full characterizations of parametric stability and convergence stability qualitatively. We summarize them as follows. We consider function $f\in\mathscr{F}$ with a domain being a compact\footnote{A set $S\subset\mathbb{R}^n$ is compact if $S$ is bounded and closed.} area and its stationary point $x^*$.
\begin{enumerate}
\it
    \item Suppose $f=\max_{i=1}^s\{f_i\}$, where $f_i$ has continuous differentials. Then $x^*$ is convergence-stable if and only if $x^*$ is the unique minimum point as well as the unique stationary point in some neighborhood of $x^*$.
    \item Suppose $f$ is finite compositions of analytic mappings\footnote{A mapping $f:\Omega\to\mathbb{R}^m$ is called analytic if for every $x\in \Omega$, there exists a neighborhood of $x$ in which $f$ is equals to its infinite Taylor expansion. Analytic property is closed under compositions and arithmetic operations.} and the maximum function on some cube $[-B,B]^n$. (Such functions form a class $\mathscr{A}$) Then $x^*$ is convergence-stable if and only if $x^*$ is the unique minimum point in some neighborhood of $x^*$.
    \item Suppose $f$ has a compact parameter space  $\mathcal{R}$ and a compact domain $\Omega$, and $f$ is continuous on $\Omega\times\mathcal{R}$. If $x^*$ is the unique minimum point in some neighborhood of $x^*$, then $x^*$ is parametric-stable. Inversely, if $x^*$ is a saddle point or local minimum points forming a simple smooth curve, there exists some $f\in\mathscr{F}$ so that $x^*$ is not parametric-stable for $f$.\footnote{The situation is indeed more complicated. The exact meaning of such existence will be discussed in \Cref{sec:prop-of-stability}.}
\end{enumerate}

These results are natural and easy to prove in convex and smooth settings. Things become subtle when $f$ is \emph{nonconvex} and \emph{nonsmooth}. Specifically, there exists a counterexample (\Cref{ex:nonex-for-converge-stable}) of $f$ in the form of finite compositions of smooth mappings\footnote{It means the mapping is differentiable of all orders. A smooth mapping is not necessarily analytic.} and the maximum function. Even though $x^*$ is the unique minimum point as well as the unique stationary point of $f$ in some neighborhood of $x^*$, $x^*$ is \emph{not} convergence-stable for $f$. Thus mere smoothness is not enough to eliminate pathological functions. Analyticity is rather necessary.

To eliminate the difference between our primitive intuitions and such subtlety, we need careful arguments. Our proofs of these three results exhibit an intrinsic connection between stability and geometry theory. Properties of stability are majorly determined by the geometric structure of the objective function regardless of what kinds of optimization algorithms are applied. Parametric stability is indeed some kind of \emph{deformation stability}: When the parameter slightly changes, the graph of $f$ near stationary point $x^*$ should not change drastically so $x^*$ will not slide away nor vanish. Convergence stability is linked to \emph{area partitioning}: Partition the domain into \emph{finitely many} areas so that a sequence of $\delta$-stationary points lie in some particular well-conditioned area. As $\delta\to 0$, the limit behavior of $\delta$-stationary points is thus confined. Area partitioning is natural for $f=\max_i\{f_i\}$: On area $\Omega_i$, $f=f_i$. For $f\in\mathscr{A}$, however, the partitioning involves deep geometry theory. We resort to \emph{semi-analytic sets} and \emph{analytic manifolds}. We carefully partition the domain into areas being both semi-analytic sets and analytic manifolds. A delicate area partitioning guarantees the isolated local minimum point to be an isolated stationary point. The technique we use is \emph{stratification theorems}, proposed by Whitney \cite{WhitneyTangents} for complex analytic sets and extended to real semi-analytic sets by {\L}ojasiewicz \cite{lojasiewicz1965ensembles, lojasiewicz1995semi}. Note that such stratification would not exist if a set is merely a smooth manifold.

Inspired by Bolte, Daniilidis and Lewis's work \cite{DBLP:journals/siamjo/BolteDL07}, we extend our proof of the second statement and characterize the geometric properties of all stationary points of a function $f\in\mathscr{A}$. Specifically, stationary points of $f$ form a finite union of connected subsets $A_i$ satisfying all the following conditions:
\begin{itemize}
    \item Each $A_i$ is a sub-analytic set as well as an analytic manifold (possibly a single point).
    \item Each point in $A_i$ is a stationary point.
    \item $f$ is a constant on the closure of $A_i$.
    \item If $x^*$ is a non-isolated local minimum point, then there exists a neighborhood $U$ of $x^*$, such that for every $A_i$ and every $x\in A_i\cap U$, $x$ is also a non-isolated local minimum point.
\end{itemize}

Our general theory has several applications. We try to resolve (or at least understand) the challenges mentioned at the beginning. The response to the challenges can be summarized as a conclusion of the stability theory in classical optimization and an opening of stability theory in deep learning. Hence we present the results in a reversed order. 
\begin{itemize}[fullwidth, listparindent=\parindent]
    \item \emph{There exists a stable (both parametric-stable and convergence-stable) tight instance for the TS algorithm. And hence $0.3393$ is a robust tight approximation ratio. That is, on such an instance, with respect to the random choice of initial points, the algorithm can find such a stationary point with positive probability in experiments.}
    
    This application is an example of applying our theory and doing quantitative analysis for a particular function. It shows the possibility to do calculations on a function with a complicated form. Such analysis allows us to algorithmically verify whether a solution found by the TS algorithm is both parametric-stable and convergence-stable. By verifying the existence of such a solution in a tight instance, we are able to prove the robust tight approximation ratio for the TS algorithm.

    \item \emph{For functions in $\mathscr{A}$ and most optimization algorithms with small enough step sizes and precision parameters, parametric stability and non-escapability properties are mainly determined by the local geometric properties of the stationary point found by the algorithms.}
    
    Most optimization algorithms need to be descent and convergent \cite{bagirov2014introduction}. Thus with small enough step sizes and precision parameters, they satisfy our requirements on optimization algorithms. Our general results can directly apply to them. Non-escapability is one interpretation of convergence stability. Thus such results are natural corollaries.
    
    \item \emph{We characterize geometric properties of all stationary points in any (possibly nonsmooth) neural network as presented above. Using these properties, we show that for most overparameterized smooth neural networks\footnote{Overparameterization means that the number of parameters is far greater than the number of training samples. Smoothness means that functions (operators) used to construct a neural network is smooth. Thus no maximum function is used.} trained with quadratic loss, all points with zero training loss fail to be convergence-stable for the gradient descent method. This indicates that slightly different initializations (even sufficiently close to a global minimum point) could result in qualitatively different trained parameters.}
    
    It is worth noting that function class $\mathscr{A}$ covers most deep neural networks and loss functions. Cooper \cite{DBLP:journals/corr/abs-1804-10200} shows that with a slight perturbation on the training input, the zero-training-loss points (if exist) of an overparameterized \emph{smooth} neural network trained with quadratic loss form a high dimensional smooth manifold (thus is non-isolated). Cooper also describes the locus of stationary points in a \emph{smooth} feedforward neural network and shows that most stationary points are highly degenerated and not isolated \cite{DBLP:journals/corr/abs-2005-04210}. Our work makes a further step to show the geometric properties of \emph{nonsmooth} neural networks as presented above. By combining our results on convergence stability and the characterization of zero-training-loss points, stability claim for global minimum points follows. Our geometric characterization of local minimum points allows us to extend the claim to nonsmooth deep neural networks with non-isolated local minimum points. 
\end{itemize}

This paper is organized as follows. In \Cref{sec:preliminaries}, we present basic concepts and notations, propose function class $\mathscr{F}$, and show some examples and basic properties about $\mathscr{F}$. In \Cref{sec:notion-of-stability}, we present the motivations and definitions of parametric stability and convergence stability, and demonstrate some examples and nonexamples of these two concepts. In \Cref{sec:prop-of-stability}, we reveal the deep connection between two notions of stability and geometry theory, and use them to prove sufficient and necessary conditions for stability. Our results are presented from simple ones to complicated ones. In \Cref{sec:app}, we give three applications of our theory: The first is a quantitative stability analysis of the TS algorithm; the second is a general conclusion on non-escapable solutions in optimization; the last is the latter shows the connection between our theory and deep neural networks in a qualitative approach. In \Cref{sec:discussion}, we summarize our results of two notions of stability and discuss further implications of our methodology and present some possible future directions.
\section{Missing Proofs in the Main Body}
\subsection{Proof of \Cref{lemma:equicontinuity}}\label{app:proof-equicontinuity}
Since $f(x;R)$ is continuous on the compact set $\Omega\times\mathscr{R}$, by Cantor's uniform continuity theorem, $f(x;R)$ is uniformly continuous on $\Omega\times\mathscr{R}$, i.e., for any $\epsilon>0$, there are some $\delta>0$ such that if $x,y\in \Omega$ and $R_1,R_2\in\mathscr{R}$, then
\[\norm{(x-y,R_1-R_2)}<\delta\Longrightarrow |f(x;R_1)-f(y;R_2)|<\epsilon.\]
Take $y=x$ and $\norm{R_1-R_2}<\delta$ and we have
\[\norm{R_1-R_2}<\delta\Longrightarrow |f(x;R_1)-f(x;R_2)|<\epsilon,\]
which is the definition of the equicontinuity.

\subsection{Proof of \Cref{lemma:Df>0circle}}\label{app:proof-Df>0circle}

Since $f\in\mathscr{F}$, the set of all valid direction unit vectors at $x^*$, denoted by $V$, is compact.
For every $v\in V$, there exists $l_v>0$ such that every point on the segment $I_v:=\{x^*+\theta v:\theta\in [0,l_v]\}$ is in $\Omega$. We have the Taylor expansion of $f$ on $I_v$,
\[f(x^*+\theta v)=f(x^*)+\frac{\partial f(x^*)}{\partial v}\theta+o(\theta).\] 
Since $\partial f(x^*)/\partial v>0$, there exists $r_v>0$ such that if $0<\theta<r_v$, then
\[f(x^*+\theta v)>f(x^*)+\frac{\theta}2\frac{\partial f(x^*)}{\partial v}.\] 

By assumption, $f$ satisfies Lipschitz condition. More precisely, $|f(x')-f(x'')|\leq Q\lVert x'-x''\rVert$ for any $x'$ and $x''$ in $\Omega$. Given $v\in V$, for any $v'$ such that $\lVert v'-v\rVert<\frac{1}{4Q}\frac{\partial f(x^*)}{\partial v}$, for any $0<\theta<r_v$, whenever $x^*+\theta v'\in \Omega$,
\begin{align*}
    f(x^*+\theta v')&\geq f(x^*+\theta v)-Q\theta\lVert v-v'\rVert\\
        &>f(x^*)+\frac{\theta}2\frac{\partial f(x^*)}{\partial v}-Q\theta\frac{1}{4Q}\frac{\partial f(x^*)}{\partial v}\\
        &=f(x^*)+\frac{\theta}4\frac{\partial f(x^*)}{\partial v}.
\end{align*}
Note that sets $B\left(v,\frac{1}{4Q}\frac{\partial f(x^*)}{\partial v}\right)$ form an open cover of $V$. By Heine–Borel finite covering lemma, there is a finite subcover 
\[\left\{B\left(v_1,\frac{1}{4Q}\frac{\partial f(x^*)}{\partial v_1}\right),\ldots, B\left(v_k,\frac{1}{4Q}\frac{\partial f(x^*)}{\partial v_k}\right)\right\}.\]
The lemma is proven by taking $r=\min\{r_{v_1},\ldots,r_{v_k}\}$ and $c=\frac{1}{4Q}\min\left\{\frac{\partial f(x^*)}{\partial v_1},\ldots,\frac{\partial f(x^*)}{\partial v_k}\right\}$.
\begin{comment}
Recall that $\partial f(x)/\partial s>0$ for every valid direction $s$ is equivalent to $Df(x,x')>0$ for all $x'\in \Omega\setminus\{x^*\}$.
Let 
\begin{align*}
    W:=\{&x^*+\theta_0(x'-x^*)\in \Omega: \text{$\theta_0$ is the supremum of $\theta$ such that }\\
    &x^*+\theta'(x'-x^*)\in \Omega\text{ for all }\theta'\in(0,\theta]\}.
\end{align*}
In other words, $W$ is the set of points that reach the farthest possible distance on some ray from origin $x^*$. We prove that $W$ is compact. $W$ is a subset of $\Omega$, thus bounded. It suffices to prove that $W$ is closed. Take a sequence $\{x_k\}_{k\geq 1}\subseteq W$ that converges to $x_0$. We prove that $x_0\in W$.

Since $f\in\mathscr{F}$, the set of all valid direction unit vectors at $x^*$, denoted by $V$, is compact. Clearly $x^*\notin W$. Thus $x_k\neq x^*$ and $v_k:=(x_k-x^*)/\norm{x_k-x^*}$ is well-defined. Since $V$ is compact, we can assume $v_k$ converges to $v_0\in V$ by taking a subsequence of $\{x_k\}$. By the definition of $W$, there is a unique point in $W\cap\{x^*+tv_0:t\geq 0\}$, denoted by $\Tilde{x}_0$. We prove that $\Tilde{x}_0=x_0$, and then $x_0\in W$. Alternatively, we prove that there are infinitely many $x_k$ in $B(\Tilde{x}_0,\delta)$ for any $\delta>0$. 

Define the cone at $x^*$:
\[CV(v,u):=\{x^*+tv'\in\mathbb{R}^n:t\geq 0, v'\in V,\norm{v-v'}\leq u\}.\]

Suppose the diameter of $\Omega$ is $d<+\infty$. By the definition of $W$, there is a unique point in $W\cap CV(v,0)$ for every $v\in V$. Then $x\in W\cap CV(v_0,\delta/(2d))$ implies $x\in B(\Tilde{x}_0,\delta)$. There are infinitely many $v_k$ such that $\norm{v_k-v_0}<\delta/(2d)$, so there are infinitely many $x_k\in B(\Tilde{x}_0, \delta)$. Thus we proved that $W$ is compact.

Since $Df(x^*,x)>0$ for all $x\in W$, by Weierstrass' minimum theorem, $\min_{x\in W} Df(x^*,x)=N>0$.

Given $x\in W$, define the segment \[I_x=\{x':x'=\theta x^*+(1-\theta) x\text{ for some }\theta\in[0,1]\},\] 
which is the collection of points on the direction from $x^*$ to $x$. On every $I_x$, we can do Taylor expansion, i.e., 
\[f(x^*+\theta(x-x^*))=f(x^*)+Df(x^*,x)\theta+o(\theta).\] 
Thus, there exists $r_x>0$ such that if $\theta<r_x$, then
\[f(x^*+\theta(x-x^*))>f(x^*)+\frac{1}{2}Df(x^*,x)\theta.\] 

By assumption, $f$ satisfies Lipschitz condition. More precisely, $|f(x')-f(x'')|\leq Q\lVert x'-x''\rVert$ for valid $x'$ and $x''$. Given $x\in W$, if $x'$ satisfies $\lVert x'-x\rVert<\frac{N}{4Q}$, then
\begin{align*}
    f(x^*+\theta(x'-x^*))&\geq f(x^*+\theta(x-x^*))-Q\theta\lVert x-x'\rVert\\
        &>f(x^*)+\frac{1}{2}Df(x^*,x)\theta-\frac{N}{4}\theta\\
        &\geq f(x^*)+\frac{N}{2}\theta-\frac{N}{4}\theta\\
        &=f(x^*)+\frac{N}{4}\theta.
\end{align*}
Note that sets $U\left(x,\frac{N}{4Q}\right)$ form an open cover of $W$. By Heine–Borel finite covering lemma, there is a finite subcover 
\[\left\{U\left(x_1,\frac{N}{4Q}\right),\ldots, U\left(x_k,\frac{N}{4Q}\right)\right\}.\]
The lemma is proven by taking $r=\min\{r_{x_1},\ldots,r_{x_k}\}$ and $c=N/4$.


\end{comment}

\subsection{Proof of \Cref{thm:analytic-max-is-convergence-stable}}\label{app:proof-analytic-max-convergence-stable}

To prove \Cref{thm:analytic-max-is-convergence-stable}, it suffices to prove the following lemma, which is the function class $\mathscr{A}$'s version of \Cref{lemma:max-type-delta-SP-distance} (in a stronger form).

\begin{lemma}\label{lemma:analytic-abs-delta-SP-converge}
Suppose $f(x)$ is finite compositions of analytic mappings and the maximum function defined on a cube $\Omega=[-B,B]^n$, i.e., $f\in\mathscr{A}$. If a stationary point $x^*$ is the unique local minimum point in $B(x^*,r)\cap\Omega$ for some radius $r>0$, then there exists $r_1>0$ such that for any $\varepsilon>0$, there exists $\delta_0>0$, such that every $\delta$-stationary point in $B(x^*,r_1)\cap\Omega$ is in $B(x^*,\varepsilon)\cap\Omega$.
\end{lemma}

The proof of \Cref{lemma:analytic-abs-delta-SP-converge} can be divided into three steps. First, we rewrite $f\in\mathscr{A}$ into finite compositions of analytic mappings and the absolute-value function. Then, we partition the domain $\Omega$ into so-called \emph{semi-analytic sets} and show deep connections between our problem and properties of semi-analytic sets. Finally, we utilize these connections to prove the lemma by the contradiction. 

We will prove \Cref{lemma:analytic-abs-delta-SP-converge} for semi-analytic domain $\Omega$ after introducing  semi-analytic sets.

\subsubsection{From the Maximum Function to the Absolute-Value Function}
The maximum function $\max\{f_i\}$ is an index selection function: Different variable $x$ would lead to different indices chosen by the maximum operator. However, there are too many index candidates. It brings verbose statements for the proof. We thus use the absolute-value function $|x|$ instead, which only chooses two possible items. A famous trick to convert maximum operator and absolute-value operator is
\[\max\{a,b\}=\frac{|a-b|+a+b}{2}.\]
Thus a series of finite compositions of analytic mappings and the maximum function are equivalent to a series of finite compositions of analytic mappings and the absolute-value function. To make notations simple, we will consider the \emph{canonical form} below.

Consider $f:\Omega\to \mathbb{R}$, where $\Omega\subseteq\mathbb{R}^n$ is a compact area. $f$ is defined by the following evaluation process:
\begin{align*}
    z_1:&=f_1(x),\\
    z_2:&=f_2(x,|z_1|),\\
    &\cdots\\
    z_i:&=f_i(x,|z_1|,|z_2|,\cdots,|z_{i-1}|),\\
    &\cdots\\
    z_m:&=f_m(x,|z_1|,|z_2|,\cdots,|z_m|),\\
    f(x):&=z_m,
\end{align*}
where $f_i:\mathbb{R}^{n+i-1}\to\mathbb{R}$. For convenience, we use $z_1(x),z_2(x),\cdots,z_m(x)$ to denote the value of these intermediate variables as functions of $x$. We further require that for every $i=1,\cdots,m+1$, $f_i$ is analytic on the compact area
\[\Omega_{i}:=\{(x,|z_1(x)|,|z_2(x)|,\cdots,|z_{i-1}(x)|):x\in\Omega\}\subseteq\mathbb{R}^{n+i-1}.\] 
It is clear that every function in the form of finite compositions of analytic mappings and the absolute-value function can be expressed in such canonical form.

\subsubsection{Semi-Analytic Geometry}
Before proposing the approach to partitioning $\Omega$, we introduce \emph{semi-analytic sets} and their geometric properties. We will use them for later proofs. For further readings, see \cite{lojasiewicz1995semi} and \cite{Semianalytic-and-subanalytic-sets}.

We say a set $L\subseteq\mathbb{R}^n$ to be \emph{semi-analytic} if in the neighborhood of each point $x_0\in\mathbb{R}^n$, $L$ is the finite union of sets of the form $\{x\in\mathbb{R}^n:f_i(x)=0,g_j(x)> 0, i=1,\cdots,i_0,j=1,\cdots,j_0\}$, where $f_i$ and $g_j$ are analytic functions from $\mathbb{R}^n$ to $\mathbb{R}^m$ in the neighborhood of $x_0$. 

We now instead suppose the domain $\Omega$ is a compact semi-analytic area and prove the theorem for functions on this kind of domains. 
We have the following properties of semi-analytic sets:

\begin{proposition}[\cite{lojasiewicz1965ensembles}]
\label{prop:semi-analytic-set-properties}
Let $L\subseteq\mathbb{R}^n$ be a semi-analytic set. The following statements hold: 
\begin{itemize}
    \item The closure $\bar{L}$ is a semi-analytic set.
    \item The finite union and finite intersection of a collection of semi-analytic sets is semi-analytic.
    \item $L$ is the locally finite union of connected semi-analytic sets (its connected components).
    \item If $L$ is connected, then any two points, $a_1$ and $a_2$, in $L$ can be joined by a semi-analytic curve. This means that there exists some embedding $\varphi: [0, 1]\to L$ whose image is a semi-analytic set in $\mathbb{R}^n$, with $\varphi(0) = a_0$ and $\varphi(1) = a_2$.
\end{itemize}
\end{proposition}

It is worth noting that statement ``locally finite union'' can be improved to ``finite union on compact area $\Omega$'' by Heine–Borel finite covering lemma. Also, a semi-analytic set restricted on semi-analytic domain $\Omega$, i.e., intersecting with $\Omega$, is still semi-analytic. We will use these tricks to improve every local property of semi-analytic sets to a global property on $\Omega$.

Let $g$ be a smooth mapping between two Euclidean spaces and $L$ is a submanifold on $g$'s domain. We use $g\vert_L$ to denote the restriction of $g$ on $L$ and $\mathrm{rank}_x(g)$ to denote the rank of Jacobi matrix of $g$ at point $x$. We have the following property:
\begin{proposition}[\cite{Semianalytic-and-subanalytic-sets}, Proof of Lemma 3.4]
\label{prop:semi-analytic-set-rank}
Let $L\subseteq\mathbb{R}^n$ be a semi-analytic set. If $L$ is smooth and connected, and $\varphi:\mathbb{R}^n\to\mathbb{R}^m$ is the projection mapping $(x_1,x_2,\cdots,x_n)\mapsto (x_1,x_2,\cdots,x_m)$, where $m<n$, then for every $k\leq m$, $\{x\in L:\mathrm{rank}_{x}(\pi\vert_L)\leq k\}$ is a semi-analytic set.
\end{proposition}

To present further properties, we need \emph{stratification theorems} motivated by Whitney \cite{WhitneyTangents}. The stratification tries to divide a manifold into a finite set of submanifolds based on some given sets. In the semi-analytic settings, we state such stratification as follows. Given an analytic manifold $M$, a \emph{normal stratification} is a finite decomposition of $M$ into semi-analytic strata, i.e. into subsets being both analytic submanifolds and semi-analytic sets. {\L}ojasiewicz proved the following theorem:
\begin{theorem}[Theorem on Normal Stratification, \cite{lojasiewicz1965ensembles}]
For arbitrarily given semi-analytic sets $E_1,\cdots, E_k$ and a point $a\in M$, there exists a normal stratification of its arbitrarily small neighbourhood $U$ which is consistent with those sets. It means that each of the sets $E_i\cap U$ is a union of some strata of this stratification.
\end{theorem}
In fact, it can be required that the strata satisfy more conditions. Let $M$ be a manifold in $\mathbb{R}^n$. $x\in M$. Denote tangent space at $x$ by $ T_x M$. Note that $ T_x M$ is a Euclidean space thus closed. Let $N$ be another manifold in $\mathbb{R}^n$. Define the distance between two tangent spaces contained in $\mathbb{R}^n$ as follows:
\[\delta( T_x M,  T_y N):=\max_{v\in T_x M,\norm{v}=1}\min_{u\in T_y N,\norm{u}=1}\norm{u-v}.\]
Note that the minimum and the maximum can be attained since the set of unit vectors in a tangent space is compact and the distance function is continuous.

A \emph{semi-analytic stratification} of an analytic manifold $M$ is, by definition, a locally finite decomposition $\mathcal{I}$ of this manifold into semi-analytic strata, fulfilling the following condition: For every stratum $\Gamma\in\mathcal{I}$, its boundary $\partial \Gamma = \bar{\Gamma}\setminus\Gamma$ is a finite union of strata of family $\mathcal{I}$ having dimensions less than $\dim\Gamma$.

Let us consider triples $(\Lambda,\Gamma,a)$, where $\Lambda, \Gamma$ are semi-analytic strata of an analytic manifold $M$, such that $\Lambda\subseteq\partial\Gamma$ and $a\in\Lambda$. A well-studied condition proposed by Whitney \cite{WhitneyTangents} is 
\[\delta( T_a\Lambda, T_x\Gamma)\to 0\quad\text{when}\quad\Gamma\ni x\to a.\]
We call this condition \emph{Whitney's condition (a)}. We have the following stratification theorem:
\begin{theorem}[Theorem on Semi-Analytic Stratification, \cite{lojasiewicz1995semi,Stratifications-Whitney}]
\label{thm:semi-analytic-stratification}
For every locally finite family of semi-analytic subsets of a manifold $M$, there exists a semi-analytic stratification, consistent with this family and such that for every pair $\Lambda, \Gamma$ of strata of this stratification, such that $\Lambda\subseteq\Gamma$, the Whitney's condition (a) is fulfilled at every point of the stratum $\Lambda$. 
\end{theorem}
We say a semi-analytic set is a \emph{manifold semi-analytic set} if it is also an analytic manifold. We have a direct corollary from \Cref{thm:semi-analytic-stratification} and \Cref{prop:semi-analytic-set-properties}:
\begin{proposition}\label{prop:semi-analytic-Whitney}
Let $L\subseteq\mathbb{R}^n$ be a semi-analytic set. Suppose $L$ is relatively compact, i.e., the closure $\bar{L}$ is compact. Then $L$ is the disjoint union of a finite collection $\{A_k\}$ of connected manifold semi-analytic subsets. Moreover, $A_k\cap\partial A_l\neq\varnothing$ implies that $A_k\subseteq\partial A_l$, and that $A_k$ and $A_l$ satisfy the Whitney's (a) condition: For any $a\in A_k$,
\[\lim_{A_l \ni x\to a}\delta( T_aA_k,  T_xA_l)=0.\]
\end{proposition}

\subsubsection{Area Partitioning}
Now we are ready to give the partitioning of $\Omega$. Define
\[\mathrm{sign}(t):=\begin{cases}
1,&t>0,\\
0,&t=0,\\
-1,&t<0.
\end{cases}
\]
Then we define $\sigma_i(x):=\mathrm{sign}(z_i(x))$ for $x\in\Omega,i=1,\cdots,m-1$, and $\sigma(x):=(\sigma_1(x),\cdots,\sigma_{m-1}(x))$. 

For $s=(s_1,\cdots,s_{m-1})\in\{-1,0,1\}^{m-1}$, we define a function $f^s(x)$ by the following evaluation process:
\begin{align*}
    z_1^s&:=f_1(x),\\
    z_2^s&:=f_2(x,s_1z_1),\\
    &\cdots\\
    z_i^s&:=f_i(x,s_1z_1,s_2z_2,\cdots,s_{i-1}z_{i-1}),\\
    &\cdots\\
    z_m^s&:=f_m(x,s_1z_1,s_2z_2,\cdots,s_{m-1}z_{m-1}),\\
    f^s(x)&:=z_m^s.
\end{align*}
Then we define 
\[Q_{s}:=\{x\in\Omega:\text{For every } i=1,\cdots,m-1,\sigma_i(x)=s_i\}.\]
The closure of $Q_{s}$ is 
\[\bar{Q}_{s}=\{x\in\Omega:\text{For every } i=1,\cdots,m-1,\sigma_i(x)=s_i\text{ or }\sigma_i(x)=0\}.\]
For $x\in\bar{Q}_{s},f(x)=f^s(x)$, and $f^s(x)$ is analytic on $\bar{Q}_{s}$. Moreover, both $Q_s$ and $\bar{Q}_s$ are relatively compact semi-analytic sets. We complete our partitioning via $Q_s$ and $\bar{Q}_s$.

\subsubsection{Limit Behaviour of $\delta$-Stationary Points as $\delta\to 0$}
Next, we present the limit behaviour of a sequence of $\delta$-stationary points as $\delta\to 0$ from a location perspective. We say $x_0\in\Omega$ is a \emph{limit-stationary point}, if there exists a point sequence $\{x_k\}_{k\geq 1}$ in $\Omega$ and a corresponding nonnegative sequence $\{d_k\}_{k\geq 1}$, such that 
\begin{itemize}
    \item $\lim_{k\to\infty}x_k=x_0$, 
    \item $\limsup_{k\to\infty}d_k=0$, and
    \item for $k=1,2,\cdots$, $x_k$ is an $d_k$-stationary point.
\end{itemize}

Suppose $x_0\in\Omega$ is such a limit-stationary point. Since $\{x_k\}$ is an infinite sequence, there is at least one $s\in\{-1,0,1\}^{m-1}$ such that for infinitely many $x_k$, $\sigma(x_k)=s$. Without loss of generality we assume that $\sigma(x_k)=s$ for all $k\geq 1$. This implies that $x_k\in Q_{s}$ and $x_0\in\bar{Q}_s$.

We first state the relation between the location of the limit-stationary point and the location of the $\delta$-stationary point sequence. By \Cref{prop:semi-analytic-Whitney}, $\bar{Q}_{s}$ can be decomposed as a disjoint union of a finite collection $\mathcal{A}$ of connected manifold semi-analytic subsets. There is $A\in\mathcal{A}$ so that $x_0\in A$, and there exists $X\in\mathcal A$ such that $X$ contains infinitely many items of $\{x_k\}$. Again without loss of generality we assume that $x_k\in X$ for all $k$. As $\lim_{k\to\infty}x_k=x_0$, by \Cref{prop:semi-analytic-Whitney} we know that either $A=X$, or $A\subseteq \partial X$. In either case the Whitney's (a) condition holds for $A$ and $X$, so $\lim_{k\to\infty}\delta( T_{x_0}A; T_{x_k}X)=0$. 


Let the graph of $f^s$ on $A$ be
\[A':=\{(x,y)\in\mathbb{R}^{n+1}:x\in A,y-f^{s}(x)=0\}.\]
Since $A$ is a connected manifold semi-analytic set in $\mathbb{R}^n$ and $f^{s}(x)$ is analytic on $\bar{Q}_s\supseteq A$, $A'$ is a connected manifold semi-analytic set in $\mathbb{R}^{n+1}$.

Our next step is, roughly speaking, to construct a set $A^*\subseteq A$ consisting of points in $A$ whose gradient on manifold $A$ is $0$. Let $\pi:\mathbb{R}^{n+1}\to\mathbb{R}$ be the projection $\pi(x,y)=y$. 
\[A^*:=\{x\in A:\mathrm{rank}_{(x,f^s(x))}(\pi\vert_{A'})=0\}.\]
By \Cref{prop:semi-analytic-set-rank}, $A^*$ is semi-analytic since $A'$ is a manifold.

A key observation is that
\begin{lemma}\label{lemma:x0-in-A*}
$x_0\in A^*$.
\end{lemma}
The lemma says that in some way, limit point $x_0$ has zero gradient on manifold $A$.
\begin{proof}
We prove the lemma by contradiction. For mapping $g$, let $J_g(x)$ be the Jacobi matrix of $g$ at $x$. Suppose that $x_0\notin A^*$, which means $\mathrm{rank}_{(x_0,f^s(x_0))}(\pi\vert_{A'})>0$. Furthermore, $\pi\vert_{A'}$ maps $A'$ to $\mathbb{R}$, a one-dimensional space, so $\mathrm{rank}_{(x_0,f^s(x_0))}(\pi\vert_{A'})\leq 1$. Thus
\[\mathrm{rank}_{(x_0,f^s(x_0))}(\pi\vert_{A'})=1.\]

Let $K=\dim A$. Then $\dim A'=K$ since $A'$ is the graph of $f^s$ on $A$, thus homeomorphic to $A$. Suppose the mapping $\psi:A\to A'$ such that $\psi(x)=(x,f^s(x))$. Note that $\pi\circ\psi=f^s$. Take a local chart $(U,\xi)$ of $A'$ at $(x_0,f^s(x_0))$, where $U$ is the homeomorphic neighborhood and $\xi$ is the analytic coordinate transformation to $\mathbb{R}^K$. Define $\hat{\xi}:=\xi\circ\psi$. Let $z:=\xi(x_0,f^s(x_0))=\hat{\xi}(x_0)\in\mathbb{R}^K$. Since $\xi$ is the analytic coordinate transformation, we have 
\[\mathrm{rank}_z(\pi\circ\xi^{-1})=\mathrm{rank}_{(x_0,f^s(x_0))}(\pi\vert_{A'})=1.\]
By the chain rule of differentials, that is to say 
\[J_{\pi\circ\xi^{-1}}(z)=J_{\pi\circ\psi\circ\hat{\xi}^{-1}}(z)=\mathrm{grad} {f^s}(x_0)J_{\hat{\xi}^{-1}}(z)\neq 0_K^T.\] 
Therefore there is some $v\in\mathbb{R}^K$ such that $\mathrm{grad} {f^s}(x_0)J_{\hat{\xi}^{-1}}(z)v<0$. Let
\[u:=\frac{J_{\hat{\xi}^{-1}}(z)v}{\norm{J_{\hat{\xi}^{-1}}(z)v}}.\]
We have $u\in  T_{x_0}A$, $\norm{u}=1$, and $\mathrm{grad} {f^s}(x_0)u<0$. Note that $\norm{\mathrm{grad} {f^s}(x)}$ is bounded on compact set $\bar{Q}_s$, so we have 
\begin{align*}
    &\liminf_{k\to\infty}\min_{T_{x_k}\mathbb{R}^n\ni u'\text{ and is valid}}\frac{\partial f^s(x_k)}{\partial u'}\\
\leq&\liminf_{k\to\infty}\min_{u'\in  T_{x_k}X}\frac{\partial f^s(x_k)}{\partial u'}\\
=&\liminf_{k\to\infty}\min_{u'\in  T_{x_k}X,\norm{u'}=1}\mathrm{grad} {f^s}(x_k)u'\\
\leq&\liminf_{k\to\infty}\min_{u'\in  T_{x_k}X,\norm{u'}=1}\left(\mathrm{grad} {f^s}(x_k)u+\norm{\mathrm{grad} {f^s}(x_k)}\cdot\norm{u'-u}\right)\\
\leq &\liminf_{k\to\infty}\left(\mathrm{grad} {f^s}(x_k)u+\norm{\mathrm{grad} {f^s}(x_k)}\cdot \max_{u''\in  T_{x_0}A,\norm{u''}=1}\min_{u'\in  T_{x_k}X,\norm{u'}=1}\norm{u'-u''}\right)\\
=&\liminf_{k\to\infty}\left(\mathrm{grad} {f^s}(x_k)u+\norm{\mathrm{grad} {f^s}(x_k)}\cdot\delta\left( T_{x_0}A,  T_{x_k}X\right)\right)\\
=&\liminf_{k\to\infty}\mathrm{grad} {f^s}(x_k)u\\
=& \mathrm{grad} {f^s}(x_0)u<0.
\end{align*}
This contradicts the fact $\liminf_{k\to\infty}d_k\geq 0$. Therefore $x_0\in A^*$.
\end{proof}

By \Cref{prop:semi-analytic-set-properties}, there are finitely many connected components of $A^*$. Assume that $\hat{A}\subseteq A^*$ is the connected component containing $x_0$. We have the following property:
\begin{lemma}\label{lemma:constant-on-hatA}
$f(x)$ is a constant on $\bar{\hat{A}}$, i.e., the closure of $\hat{A}$.
\end{lemma}
\begin{proof}
Let $x''$ be any point in $\hat{A}$. We need to show that for every $x\in\bar{\hat{A}}$, $f(x'')=f(x)$. Note that $\bar{\hat{A}}$ is a semi-analytic set. By \Cref{prop:semi-analytic-set-properties}, there is a semi-analytic curve $\varphi(t):[0,1]\to\bar{\hat{A}}$ such that $\varphi(0)=x,\varphi(1)=x''$, and $\varphi(t)\in\bar{\hat{A}}$ for $t\in(0,1]$. $\varphi(t)$ is analytic everywhere except for finitely many points, so without loss of generality we assume that $\varphi(t)$ is analytic on $(0,1)$. It suffices to prove that $f^s(\varphi(0))=f^s(\varphi(1))$. 

For any mapping $g$, let $J_g(x)$ be the Jacobi matrix of $g$ at $x$. We have
\[f^s(\varphi(0))-f^s(\varphi(1))=\int_0^1\mathrm{grad} {f^s}(\varphi(t))J_{\varphi}(t)dt.\] 
The (Lebesgue) integral is well-defined since $\mathrm{grad} {f^s}(\varphi(t))J_{\varphi}(t)$ is bounded on $[0,1]$. By the property of curves, $J_{\varphi}(t)\in T_{\varphi(t)}A$. We claim that for every $x\in A^*$ and every $u\in T_x A$, $\mathrm{grad} {f^s}(x)u=0$. Then it follows that $f^s(\varphi(0))-f^s(\varphi(1))=0$, therefore $f^s(x)=f^s(x'')$, i.e., $f(x)=f(x'')$.

Let $\dim A=K$. Take a local chart $(U,\xi)$ at $x\in A^*$ at $(x,f^s(x))$, where $U$ is the homeomorphic neighborhood and $\xi$ is the analytic coordinate transformation to $\mathbb{R}^K$. Define $\hat{\xi}:=\xi\circ\psi$. Let $z:=\xi(x,f^s(x))=\hat{\xi}(x)\in\mathbb{R}^K$. Recall the definition of $A^*$,
\[\mathrm{rank}_z\left(\pi\vert_{A'}\right)=0.\]
Since $\xi$ is the analytic coordinate transformation, we have $\mathrm{rank}_z(\pi\circ\xi^{-1})=\mathrm{rank}_z(\pi\vert_{A'})=0$. By the chain rule of differentials, that is to say 
\[J_{\pi\circ\xi^{-1}}(z)=J_{\pi\circ\psi\circ\hat{\xi}^{-1}}(z)=\mathrm{grad} {f^s}(x)J_{\hat{\xi}^{-1}}(z)= 0_K^T.\] For every $u\in T_x A=\mathbb{R}^K$, it can be written in $J_{\hat\xi^{-1}}v$ for some $v\in\mathbb{R}^K$, Thus
\[\mathrm{grad} f^s(x)u=\mathrm{grad} {f^s}(x)J_{\hat{\xi}^{-1}}(z)v=0.\]
\end{proof}

\subsubsection{Isolation Property of local minimum point $x^*$}


We have shown that, for every $s\in\{-1,0,1\}^{m-1}$, there are finitely many $A$'s in the corresponding decomposition $\mathcal A$ of $\bar{Q}_s$, and we have constructed a semi-analytic subset $A^*$ from each $A$. Let $P$ denote the union of all these $A^*$'s. $P$ is still a semi-analytic set, and has finitely many connected components by \Cref{prop:semi-analytic-set-properties}. As $\Omega=\bigcup_{s\in\{-1,0,1\}^{m-1}}Q_s$, for any $x_0\in\Omega$ which is a limit-stationary point, $x_0\in P$.

Now suppose a stationary point $x^*\in\Omega$ is the unique local minimum point in $B(x^*,r)$ for some radius $r$. Observe that any stationary point is also a limit-stationary point, so we have $x^*\in P$. We have the following lemma:
\begin{lemma}\label{lemma:x*-is-isolated}
$x^*$ is an isolated point of $P$. In other words, there is some $r_1>0$ such that $B(x^*,2r_1)\cap P=\{x^*\}$.
\end{lemma}
\begin{proof}
Suppose for the contrary that $x^*$ is not an isolated point of $P$. Note that there are only finitely many $A^*$'s and finitely many $\hat{A}$'s. In addition, $P$ is the union of $\hat{A}$'s. Thus there exists some infinite set $\hat{A}$ such that $x^*$ is a cluster point of $\hat{A}$. However, by \Cref{lemma:constant-on-hatA} we know that for any $x\in \hat{A}$, $f(x)=f(x^*)$. Thus there exists an infinite sequence $x_k\to x^*$ in $\hat{A}$ such that $f(x_k)=f(x^*)$ for any $k$. That contradicts the fact that $x^*$ is the unique local minimum point in $B(x^*,r)$. Therefore $x^*$ must be an isolated point of $P$. 
\end{proof}

At last, we prove \Cref{lemma:analytic-abs-delta-SP-converge} by contradiction. Suppose that there is some $\varepsilon_0>0$ such that for every $\delta>0$, there exists $x'\in (B(x^*,r_1)\setminus B(x^*,\varepsilon_0))\cap\Omega$, such that $x'$ is a $\delta$-stationary point. Equivalently, there is an infinite point sequence $\{x_k\}$ in $(B(x^*,r_1)\setminus B(x^*,\varepsilon_0))\cap\Omega$, and a corresponding sequence $\{d_k\}$, such that $\limsup_{k\to\infty}d_k=0$, and for $k=1,2,\cdots$, $x_k$ is an $d_k$-stationary point. Take any convergent subsequence of $\{x_k\}$, and assume it converges to some point $x_0$. By \Cref{lemma:x0-in-A*}, $x_0\in P$ and $\norm{x_0-x^*}\in[\varepsilon_0,r_1]$, which contradicts the fact $B(x^*,2r_1)\cap P=\{x^*\}$.

Now we finally complete the proof of \Cref{thm:analytic-max-is-convergence-stable}.


\subsection{Proof of \Cref{prop:geometric-char-of-sp}}\label{app:proof-geometric-char-of-sp}
We have shown that all the limit-stationary points (also, stationary points) are contained in a semi-analytic set $P$, which has some good properties:
\begin{enumerate}
    \item $P$ is a semi-analytic set, and can be decomposed into the union of a finite collection of connected manifold semi-analytic subsets.
    \item $P$ has finitely many connected components, and $f$ is constant on each component.
\end{enumerate}


However, one may notice that, not every point in $P$ is a limit-stationary point. Inspired by \cite{DBLP:journals/siamjo/BolteDL07} and \cite{coste2000introduction}, we introduce sub-analytic sets and sub-analytic functions, and give a better characterization on the locus of all the limit-stationary points, and stationary points, respectively.

We say a set $A\subseteq\mathbb{R}^n$ to be \emph{sub-analytic} if each point $x\in \mathbb{R}^n$ admits a neighborhood $U$, such that 
$$A\cap U=\pi_{n+m,n}(B),$$
where $B$ is a bounded semi-analytic subset of $\mathbb{R}^{n+m}$ for some $m$, and $\pi_{n+m,n}(x_1,\cdots,x_{n+m})=(x_1,\cdots,x_n)$ is the canonical projection from $\mathbb{R}^{n+m}$ to $\mathbb{R}^n$.


\begin{proposition}[\cite{Semianalytic-and-subanalytic-sets}, Section 3]\label{prop:sub-analytic-set-properties}
We have the following properties of sub-analytic sets:
\begin{itemize}
    \item Every semi-analytic set is sub-analytic.
    \item Sub-analytic sets are closed under finite union and intersection.
    \item The closure, interior, and complement of a sub-analytic set is sub-analytic.
    \item The Cartesian product of two sub-analytic sets is sub-analytic.
    \item If $A\subseteq \mathbb{R}^n$ is a sub-analytic set, then $\pi_i(A)$ is sub-analytic, here \[\pi_i(x_1,\cdots,x_n)=(x_1,\cdots,x_{i-1},x_{i+1},\cdots,x_n).\]
\end{itemize}
\end{proposition}

For a mapping $f:\Omega\to \mathbb{R}^m$ defined on $\Omega\subseteq \mathbb{R}^n$, the \emph{graph} of $f$, denoted by $\mathrm{graph} f$, is defined by
$$\mathrm{graph} f:=\{(x,y):x\in \Omega,y=f(x)\}\subseteq\mathbb{R}^{n+m}.$$
We say $f$ is a sub-analytic mapping, if $\mathrm{graph} f$ is a sub-analytic subset of $\mathbb{R}^{n+m}$. Similarly we say a function $f:\Omega\to \mathbb{R}$ is sub-analytic if $\mathrm{graph} f$ is a sub-analytic subset of $\mathbb{R}^{n+1}$.
We will use the following facts:
\begin{proposition}[\cite{denkowska1979certaines}]\label{prop:sub-analytic-mapping-properties}
\begin{itemize}
    \item If $f:A\to B$ is a sub-analytic mapping, and $E$ is a bounded sub-analytic subset of $A$, then $f(E)$ is a sub-analytic set.
    \item If $f:A\to B,g:B\to C$ are both sub-analytic mappings, and $f$ maps bounded sets to bounded sets, then the composition $g\circ f$ is sub-analytic.
    \item Let $A$ be a sub-analytic set. A function $f:A\to\mathbb{R}$ is sub-analytic if and only if the epigraph of $f$, defined as $\mathrm{epi}f:=\{(x,y)\in A\times\mathbb{R}:y\geq f(x)\}$, is sub-analytic.
\end{itemize}
\end{proposition}

An important observation in \cite{coste2000introduction} is that ``all first-order formula is definable'', in other words, sub-analytic sets are able to express quantifiers.

\begin{lemma}\label{lemma:sub-analytic-first-order-formula}
If $A\subseteq \mathbb{R}^n,B\subseteq \mathbb{R}^m,C\subseteq \mathbb{R}^{m+n}$ are sub-analytic sets, then $$\begin{aligned}
\{x\in A:\exists y\in B,(x,y)\in C\}&\text{ and}\\
\{x\in A:\forall y\in B,(x,y)\in C\}&
\end{aligned}$$ are both sub-analytic sets.
\end{lemma}
\begin{proof}
Firstly, observe that $\{x\in A:\exists y\in B,(x,y)\in C\}=\pi_{n+m,n}((A\times B)\cap C)$. As sub-analytic sets are closed under complement operation, and therefore under set minus operation, we can similarly write $\{x\in A:\forall y\in B,(x,y)\in C\}$ as $A\setminus\pi_{n+m,n}(C\setminus(A\times B))$. By \Cref{prop:sub-analytic-set-properties}, we immediately get that $\{x\in A:\exists y\in B,(x,y)\in C\}$ and $\{x\in A:\forall y\in B,(x,y)\in C\}$ are both sub-analytic.
\end{proof}

Now we consider $f\in\mathscr{A}$ defined on a compact semi-analytic set $\Omega\subseteq \mathbb{R}^n$, as previously described. We show that:
\begin{lemma}\label{lemma:f-is-sub-analytic}
$f(x)$ and $\partial f(x)=\min_{s}\frac{\partial f(x)}{\partial s}=\min_{x'\neq x}\frac{Df(x,x')}{\norm{x'-x}}$ are both sub-analytic functions.
\end{lemma}
\begin{proof}
Recall that $f(x)$ is defined with $f_1,\cdots,f_m$, by the evaluation process that for $i=1,\cdots,m$, $z_i(x):=f_i(x,|z_1(x)|,|z_2(x)|,\cdots,|z_{i-1}(x)|)$, and $f(x):=z_m(x)$.

$z_1(x)=f_1(x)$ is analytic on $\Omega$, so $\mathrm{graph} z_1$ is an semi-analytic subset of $\mathbb{R}^{n+1}$, thus $z_1(x)$ is sub-analytic. Clearly the absolute-value function is also sub-analytic. Additionally, because of continuity, each $z_i(x)$ is bounded on compact set $\Omega$. Therefore if $z_1(x),\cdots,z_{k}(x)$ are sub-analytic, by \Cref{prop:sub-analytic-mapping-properties} it follows that $z_{k+1}(x)=f_{k+1}(x,|z_1(x)|,|z_2(x)|,\cdots,|z_{k}(x)|)$ is also sub-analytic. By induction we get that $f(x)=z_m(x)$ is a sub-analytic function.

Then we prove that $\partial f(x)$ is also a sub-analytic function with the help of \Cref{lemma:sub-analytic-first-order-formula}. For convenience, we first show that the Dini derivative $Df(\cdot,\cdot)$ is sub-analytic. For simplicity, we assume that $\Omega$ is a closed convex polyhedron. We write $\mathrm{graph}{Df(\cdot,\cdot)}$ as a first-order logic expression
\[\{(x,x',\delta)\in\Omega\times\Omega\times\mathbb{R}:\forall \epsilon>0,\exists r>0,\forall t>0, t>r\vee(\delta-\epsilon)t< f(x+t(x'-x))-f(x)<(\delta+\epsilon)t\},\]
which is sub-analytic. As each $f_i$ is analytic, $Df(x,x')$ is bounded on $\Omega\times\Omega$.


The epigraph of $\partial f(x)$, $\mathrm{epi} \partial f:=\left\{(x,\delta)\in\Omega\times\mathbb{R}:\delta\geq \partial f(x)=\min_{x'\neq x}\frac{Df(x,x')}{\norm{x'-x}}\right\}$ can be written as 
\[\mathrm{epi} \partial f=\{(x,\delta)\in\Omega\times\mathbb{R}:\exists x'\in \Omega, \norm{x'-x}>0\wedge Df(x,x')\leq \delta\norm{x'-x}\}.\]
Therefore $\partial f(x)$ is a sub-analytic function, by \Cref{prop:sub-analytic-mapping-properties}.
We note that $\partial f(x)$ is also bounded.
\end{proof}

We are ready to prove our characterization result on stationary points in \Cref{prop:geometric-char-of-sp}.
Let $S:=\{x\in\Omega:x\text{ is a stationary point of }f\}$, we have
\[S=\{x\in\Omega: \partial f(x)\geq 0\}=\{x\in\Omega:\forall y<0, \partial f(x)>y\}.\]
By \Cref{lemma:f-is-sub-analytic} and \Cref{lemma:sub-analytic-first-order-formula}, $S$ is a sub-analytic set.

By the theorem on sub-analytic stratification \cite{lojasiewicz1995semi, denkowska1979certaines} similar to \Cref{thm:semi-analytic-stratification}, we know that $S$ can be decomposed as the union of a finite collection $\{S_i\}$ of sub-analytic subsets, so that each $S_i$ is either a connected analytic manifold, or a single point.

Note that $S\subseteq P$, where $P$ is the semi-analytic set described at the beginning of the proof. So each $S_i$ must be a subset of some connected component $P'$ of $P$, and therefore $\bar{S_i}\subseteq \bar{P'}$. Recall that $f$ is constant on the closure of every connected component of $P$, therefore for any $x_1,x_2\in \bar{S_i}$, we have $f(x_1)=f(x_2)$.

Suppose $x^*\in\Omega$ is a non-isolated local minimum point of $f$, then there exists a neighborhood $U$ of $x^*$, so that $f(x)\geq f(x^*)$ for all $x\in U$. $x^*$ is clearly a stationary point, so $x^*\in S$. $S$ has finitely many connected components. Furthermore, for each connected component $S'$, $f$ is constant on $\bar{S'}$. If $x^*\in \bar{S'}$, we have $f(x)=f(x^*)$ for $x\in S'$; otherwise, we have $\inf_{x\in S'}\lVert x-x^*\rVert>0$. Therefore there is $r>0$ such that for all $x\in S\cap B(x^*,r),f(x)=f(x^*)$. Let $U^*=U\cap B(x^*,r)$, then for every $x\in S\cap U^*$, there is a neighborhood $U'$ of $x$, such that $U'\subseteq U$. It follows that $f(x)=f(x^*)\leq \inf_{x'\in U'}f(x')$, i.e., $x$ is also a non-isolated local minimum point. Now we complete the proof.

We remark that, the set of all limit-stationary points, denoted as $S^{\lim}$, can be written as
\[S^{\lim}=\{x\in\Omega:\forall\delta>0,\forall r>0,\exists x'\in\Omega,\lVert x-x'\rVert<r\wedge \partial f(x')>-\delta\}.\]
Thus all limit-stationary points also form a sub-analytic set, and have similar geometric properties.

\subsection{Proof of \Cref{lemma:descent-monotone}}\label{app:proof-descent-monotone}
For simplicity, we use $S_R(y):=\suppmax(Ry)$ and $S_C(x):=\suppmax(\T{C}x)$. The following lemma shows that in a neighborhood of a point with $f_R=f_C$, function $f$ behaves like the maximum of two quadratic functions in every direction.
\begin{lemma}[\cite{DBLP:journals/im/TsaknakisS08}, Appendix 8.1]\label{lemma:exact-form-Df}
Suppose $f_R(x,y)=f_C(x,y)$. Define
\[I(y,y',\epsilon):=\epsilon\left(Ry'-\max_{S_R(y)}(Ry')1_m\right)-(1-\epsilon)\left(\max(Ry)1_m-Ry\right).\]
Similarly, define
\[I(x,x',\epsilon):=\epsilon\left(\T{C}x'-\max_{S_C(x)}(\T{C}x')1_n\right)-(1-\epsilon)\left(\max(\T{C}x)1_n-\T{C}x\right).\]
Define $Hf_R(x,y,x',y'):=\T{(x'-x)}R(y'-y)$, and $Hf_C:=\T{(x'-x)}C(y'-y)$.
Denote the complement sets by $\bar{S_R}(y):=[m]\setminus S_R(y)$, $\bar{S_C}(x):=[n]\setminus S_C(x)$. Define
\[r(x',y'):=\min\{\epsilon\in[0,1]:\max_{\bar{S_R}(y)}(I(y,y',\epsilon))\leq 0\text{ and }\max_{\bar{S_C}(x)}(I(x,x',\epsilon))\leq 0\}.\]
If we use $Df_R$ for short as $Df_R(x,y,x',y')$ and other abbreviations are similar, then for all $\epsilon\in[0,r(x',y')]$, \[f(x+\epsilon(x'-x),y+\epsilon(y'-y))=f(x,y)+\max\left\{\epsilon Df_R-\epsilon^2Hf_R,\epsilon Df_C-\epsilon^2Hf_C\right\}.\]
\end{lemma}
Now we prove \Cref{lemma:descent-monotone}.

If $f_R(x_k,y_k)=f_C(x_k,y_k)$, this property is guaranteed by the line search algorithm: It computes $\epsilon^*\leq r(x',y')$ such that $f(x+\epsilon(x'-x),y+\epsilon(y'-y))\, (\epsilon\in[0,\epsilon^*])$ is the maximum of two quadratic functions. Furthermore,  
\[\epsilon^*\leq\frac{\max\{Df_R(x,y,x',y'),Df_C(x,y,x',y')\}}{2\min\{Hf_R(x,y,x',y'),Hf_C(x,y,x',y')\}}.\]
Since $(x',y')$ is a descent direction, by the property of quadratic functions, $f(x+\epsilon(x'-x),y+\epsilon(y'-y))$ is a decreasing convex function when $\epsilon\in[0,\epsilon^*]$.

If $f_R(x,_k,y_k)\neq f_C(x_k,y_k)$, without loss of generality we assume $f_R(x,_k,y_k)>f_C(x_k,y_k)$. In this case $y_{k+1}=y_k$, and $x_{k+1}$ is the optimal solution of the LP problem 
\begin{gather*}
    \min_{x\in\Delta_m} f_R(x,y_k),\\
    \text{ s.t. }f_R(x,y_k)\geq f_C(x,y_k).
\end{gather*}
It is clear that both $x_{k}$ and $x_{k+1}$ are feasible solutions of the LP problem, so $$f_R(x_{k},y_k)\geq f_R(x_{k+1},y_k)=f_R(x_{k+1},y_{k+1}).$$
Moreover by the linearity of target function in LP, for $\theta\in(0,1)$,
\begin{align*}
    &f((1-\theta)x_k+\theta x_{k+1},(1-\theta)y_k+\theta y_{k+1})\\
    =&f((1-\theta)x_k+\theta x_{k+1},y_k)\\
    =&(1-\theta)f_R(x_{k},y_{k})+\theta f_R(x_{k+1},y_{k})\\
    \leq &f_R(x_{k},y_k).
\end{align*}

\subsection{Proof of \Cref{prop:strict-delta-SP-distance}}\label{app:proof-strict-delta-SP-distance}
Let boundary set 
\begin{align*}
    W:=\{(x',y')\in\Delta_m\times\Delta_n:(x+\theta(x'-x),y+\theta(y'-y))\notin\Delta_m\times\Delta_n\text{ for }\theta>1\}.
\end{align*}
In other words, $W$ is the set of points that reach the farthest possible distance on some ray from origin $x^*$. By the property of simplexes (convex polyhedrons), $W$ is compact. Recall that we use $S_R(y)=\suppmax(Ry)$ and $S_C(x)=\suppmax(\T{C}x)$. 

We enhance \Cref{lemma:exact-form-Df} so that there exists a uniform $r>0$, which does not depend on the direction $(x',y')$.
\begin{lemma}\label{lemma:exact-form-circle}
Use the notation in \Cref{lemma:exact-form-Df}. Then there exists $\gamma>0$ such that if $(x+\epsilon(x'-x),y+\epsilon(y'-y))\in\Delta_m\times\Delta_n$ satisfies $0<\epsilon<\gamma$, then
\[f(x+\epsilon(x'-x),y+\epsilon(y'-y))=f(x,y)+\max\left\{\epsilon Df_R-\epsilon^2Hf_R,\epsilon Df_C-\epsilon^2Hf_C\right\}.\]
\end{lemma}
\begin{proof}
It suffices to prove the case when $(x',y')\in W$. For every index $i\in\bar{S_R}(y)$,
\[I(y,y',\epsilon)_i=\epsilon\left((Ry')_i-\max_{S_R(y)}(Ry')+\max(Ry)-(Ry)_i\right)-\left(\max(Ry)-(Ry)_i\right).\]
Note that $I(y,y',0)_i<0$, thus for all $\epsilon\in[0,r_i(y')]$, $I(y,y',\epsilon)_i\leq 0$ where
\[r_i(y')=\frac{\max(Ry)-(Ry)_i}{\max\left\{\max(Ry)-(Ry)_i,(Ry')_i-\max_{S_R(y)}(Ry')+\max(Ry)-(Ry)_i\right\}}\]
is a continuous function of $y'$.  Note that $I(y,y',0)_i<0$ for all $(x',y')\in W$ and $i\in\bar{S_R}(y)$. Since $I(y,y',\epsilon))_i$ is continuous in $\epsilon$, there must be some $\epsilon_0>0$, such that $I(y,y',\epsilon_0))_i<0$. Furthermore,
\[\max_{\bar{S_R}(y)}(I(y,y',\epsilon))\leq 0\iff \epsilon\in[0,\min_{i\in\bar{S_R}(y)} r_i(y')].\]
Therefore $\min_{i\in\bar{S_R}(y)} r_i(y')>0$ for all $(x',y')\in W$.

Define $r(x')$ similarly and the parallel results hold for $r(x')$. It is not hard to see that $\epsilon$ satisfies $\max_{\bar{S_R}(y)}(I(y,y',\epsilon))\leq 0$ and $\max_{\bar{S_C}(x)}(I(x,x',\epsilon))\leq 0$ if and only if \[\epsilon\in\left[0,\min\left\{\min_{i\in\bar{S_R}(y)} r_i(y'),\min_{i\in\bar{S_C}(x)} r_i(x')\right\}\right]=[0,r(x',y')].\]

By previous arguments, 
\[r(x',y')=\min\left\{\min_{i\in\bar{S_R}(y)} r_i(y'),\min_{i\in\bar{S_C}(x)} r_i(x')\right\}>0\]
for all $(x',y')\in W$ and $r$ is continuous in $(x',y')$. By compactness of $W$, 
\[\min_{(x',y')\in W}r(x',y')>0.\] 
Take $\gamma$ as $\min_{(x',y')\in W}r(x',y')$. By \Cref{lemma:exact-form-Df}, for all $\epsilon\in[0,\gamma]$ and any $(x',y')\in W$, \[f(x+\epsilon(x'-x),y+\epsilon(y'-y))=f(x,y)+\max\left\{\epsilon Df_R-\epsilon^2Hf_R,\epsilon Df_C-\epsilon^2Hf_C\right\}.\]
\end{proof}

Having proved \Cref{lemma:exact-form-circle}, we can calculate directional derivatives towards stationary point $(x^*,y^*)$ of any point $(x,y)$ near $(x^*,y^*)$. To prove every $\delta$-stationary point $(x,y)$ in a neighborhood of $(x^*,y^*)$ is close to $(x^*,y^*)$, we can calculate $Df(x,y,x^*,y^*)$ and show that it is negative and its absolute value cannot be too small. So we have the following technical lemma.
\begin{lemma}\label{lemma:quadratic-max-increase}
Let $g_i(x)=a_ix-b_ix^2$, $i=1,2$, $x\geq 0$. Suppose $a_1\geq a_2$ and $a_1>0$. Then on the interval $I=\begin{cases}
(0,\frac{a_1}{4b_1}),&\text{if }b_1>0,\\
(0,+\infty),&\text{if }b_1\leq 0
\end{cases}$, function $g(x)=\max\{g_1(x),g_2(x)\}$ is strictly increasing and the left derivatives of $g(x)$ are at least $a_1/2$.
\end{lemma}

\begin{proof}
As $a_1>0$, and $g(x)$ is continuous, we only need to prove that 
\[g'_{-}(x)=\lim_{x'\to x^{-}}\frac{g(x')-g(x)}{x'-x}\geq \frac{a_1}2\] 
holds for every $x\in I$. Then $g$ is strictly increasing on $I$.

One can check that $g_1'(x)=a_1-2b_1x\geq a_1/2>0$
holds for any $x\in I$ in both cases of $b_1>0$ and $b_1\leq 0$. If $b_2\geq b_1$, then $g_2(x)=a_2x-b_2x^2\leq a_1x-b_1x^2=g_1(x)$, so $g(x)=g_1(x)$ holds on $[0,+\infty)$, and thus for $x\in I,~g'_{-}(x)=g'(x)=g_1'(x)\geq a_1/2$.

Now we suppose $b_1>b_2$. If $a_1=a_2$, then $g(x)=g_2(x)$ and $g'_{-}(x)=g_2'(x)=a_2-2b_2x > g_1'(x)\geq a_1/2>0$ for any $x\in I$.

If $a_1>a_2$, then $g_1$ and $g_2$ intersect at $x_1=0$ and $x_2=(a_1-a_2)(b_1-b_2)$. For $x\in(0,x_2]$, $g(x)=g_1(x),~g'_{-}(x)=g_1'(x)$. For $x\in (x_2,+\infty)$,$g(x)=g_2(x)>g_1(x),~g'_{-}(x)=g_2'(x)=a_2-2b_2x=(g_2(x)-b_2x^2)/x>(g_1(x)-b_1x^2)/x=g_1'(x)$. Therefore for any $x\in I,g'_{-}(x)\geq g_1'(x)\geq a_1/2$.

In summary, $g'_{-}(x)\geq a_1/2>0$ holds for every $x\in I$ in all these cases.
\end{proof}

Now we are ready to prove the proposition.

By \Cref{lemma:exact-form-circle}, there exists $\gamma_0>0$ such that for every $(x',y')\in W$ and any $\epsilon\in[0,\gamma_0)$,
\[f(x^*+\epsilon(x'-x^*),y^*+\epsilon(y'-y^*))=f(x^*,y^*)+\max\left\{\epsilon Df_R-\epsilon^2Hf_R,\epsilon Df_C-\epsilon^2Hf_C\right\}.\]


Let
\[a(x',y')=Df(x^*,y^*,x',y')=\max\left\{Df_R(x^*,y^*,x',y'),Df_C(x^*,y^*,x',y')\right\},\]
\[b(x',y')=\max\{(x'-x^*)^T R (y'-y^*),(x'-x^*)^T C (y'-y^*)\}.\] 
Let $x_\theta$ denote $x^*+\theta(x'-x^*)$, $y_\theta$ denote $y^*+\theta(y'-y^*)$, and $g(\theta)=f(x_\theta,y_\theta)-f(x^*,y^*)=\max\{\theta Df_R-\theta^2Hf_R,\theta Df_C-\theta^2Hf_C\}$ for $\theta\in[0,\gamma_0)$.

As $(x^*,y^*)$ is a strict stationary point, we have $a(x',y')> 0$ holds for every $(x',y')\in W$.

By \Cref{lemma:quadratic-max-increase}, for any $\theta\in\left(0,\Theta(x',y')\right)$, we have $g'_{-}(\theta)\geq a(x',y')/2$, where 
\[\Theta(x',y')=\begin{cases}
\min\left\{\gamma_0,\frac{a(x',y')}{4b(x',y')}\right\},& \begin{aligned}
    \text{if }&(Df_R\geq Df_C\text{ and } Hf_R>0)\text{ or}\\
    &(Df_C\geq Df_R\text{ and } Hf_C>0),
\end{aligned}\\
\gamma_0,& \text{otherwise.}
\end{cases}\]
Note that
\begin{align*}
    \min_{x'',y''}Df(x_\theta,y_\theta,x'',y'')&\leq Df(x_\theta,y_\theta,x^*,y^*)\\
    &=\lim_{\epsilon\to 0+}\frac1{\epsilon}(f(x_\theta+\epsilon(x^*-x_\theta),y_\theta+\epsilon(y^*-y_\theta))-f(x_\theta,y_\theta))\\
    &=\lim_{\epsilon\to 0+}\frac1{\epsilon}(g(\theta-\epsilon\theta)-g(\theta))\\
    &=-\theta\cdot g'_{-}(\theta).
\end{align*}
 Assume $(x_\theta,y_\theta)$ is a $\delta$-stationary point, then
 \[-\delta\leq\min_{x'',y''}Df(x_\theta,y_\theta,x'',y'')\leq Df(x_\theta,y_\theta,x^*,y^*)=-\theta\cdot g'_{-}(\theta),\] thus $\theta\leq \delta/g'_{-}(\theta)\leq 2\delta/a(x',y')$.

Note that $a(x',y')$ and $b(x',y')$ are continuous in $(x',y')$. In addition, for any $(x',y')\in W$, $a(x',y')>0$. As $W$ is compact, we have $$\bar a:=\min\{a(x',y'):(x',y')\in W\}>0.$$ Let $\gamma_1=\min\{\gamma_0,\bar a/16\}$, then $\gamma_1\leq \Theta(x',y')$ holds for every $(x',y')\in W$ because $b(x',y')\leq 4$.

Let 
\[U_1=\{(x^*+\theta(x'-x^*),y^*+\theta(y'-y^*)):\theta\in [0,\gamma_1),(x',y')\in W\},\] we can take sufficiently small $r_1>0$ so that $B((x^*,y^*),r_1)\subseteq U_1$. 

For any $\varepsilon>0$, there exists $\theta_0\in(0,1)$ such that for $(x',y')\in W,\theta\in[0,\theta_0]$,
\[\norm{(x^*+\theta(x'-x^*),y^*+\theta(y'-y^*))-(x^*,y^*)}<\varepsilon.\] 
Let $\delta_0=\bar{a}\theta_0/2$.

For any $\delta<\delta_0$ and any $\delta$-stationary point $(x'',y'')$ in $B((x^*,y^*),r_1)$, suppose $x''=x^*+\theta(x'-x^*),y''=y^*+\theta(y'-y^*)$, where $\theta\in [0,\bar a/16),(x',y')\in W$. Then $\theta\leq 2\delta/a(x',y')<2\delta_0/\bar{a}=\theta_0$. Hence we have  $(x'',y'')\in B((x^*,y^*),\varepsilon)$.

For every $(x',y')\in W$, it is clear that $g(\theta)>0$ holds for any $\theta\in(0,\Theta(x',y'))$, so $f(x^*,y^*)<f(x',y')$ holds for every $(x',y')\neq (x^*,y^*)$ in $B((x^*,y^*),r_1)$.

\subsection{Proof of \Cref{prop:strict-tight-SP-exists}}\label{app:proof-strict-tight-SP-exists}
We will prove the proposition by presenting an example and verifying it satisfies the condition. So the proof is naturally divided into two parts: 
\begin{enumerate}
    \item Find a stable tight instance $(R_0,C_0)$ with a tight stationary point $(x^*,y^*)$.
    \item Verify that $(x^*,y^*)$ is a strict stationary point.
\end{enumerate}

Luckily, \cite{DBLP:journals/corr/abs-2107-01471} proposes a tight instance generator and empirically show that there exist stable tight instances, especially for $3\times 3$ game instances. Thus we can find them by experiments! We present the following tight instance. Suppose $m=n=3$, Let
\begin{equation}\label{ex:stable-TS-tight-instance}
    R=\begin{pmatrix}0.6607 & 0 & 0\\
1&1&1\\
1&0.1703&0.1703
\end{pmatrix},\quad C=\begin{pmatrix}0.6607 & 1 & 1\\
0&1&0.0988\\
0&1&0.0988
\end{pmatrix}
\end{equation}
Then \eqref{ex:stable-TS-tight-instance} is a tight instance with strict stationary point $x^*=y^*=\T{(1,0,0)}$ and dual solution $\rho^*=0.582522$ $w^*=\T{(0,0.7744,0.2256)}, z^*=\T{(0,0.5367,0.4633)}$.

After finding an empirically stable tight instance, we have to verify: \begin{enumerate*}[label=(\roman*)]
    \item \eqref{ex:stable-TS-tight-instance} is indeed a tight instance;
    \item $(x^*,y^*)$ is indeed a strict stationary point.
\end{enumerate*}

For the first verification, \cite{DBLP:journals/corr/abs-2107-01471} presents an equivalent condition. We summarize this result in \Cref{thm:appTS:tight-iff}. Recall that we use $S_R(y):=\suppmax(Ry)$ and $S_C(x):=\suppmax(\T{C}x)$. Let
\begin{align*}
    A(\rho,y,z)&:=-\rho Ry+(1-\rho)C(z - y),\\
    B(\rho,x,w)&:=\rho \T{R}(w - x)-(1-\rho)\T{C}x.
\end{align*}
And define quantities
\begin{align*}
\lambda^* &:= \T{(w^*-x^*)}Rz^*,\\
  \mu^* &:=\T{w^*}C(z^*-y^*).
\end{align*}
\begin{theorem}[\cite{DBLP:journals/corr/abs-2107-01471}, Section 4]\label{thm:appTS:tight-iff}
Let $(R,C)$ be a game. Let $b\approx0.3393$ be the tight bound of the TS algorithm. Let $\partial\Lambda$ be the boundary of $[0,1]\times [0,1]$. The following statements characterize a tight instance:
\begin{itemize}
    \item $(x^*,y^*)$ is a stationary point and the corresponding dual solution is $(\rho^*,w^*,z^*)$,
    \item $f(\alpha x^*+(1-\alpha)w^*,\beta y^*+(1-\beta)z^*)\geq b$ for all $(\alpha,\beta)\in\partial\Lambda$.
\end{itemize}
The above statements hold if and only if all below statements hold:
\begin{itemize}
    \item $\displaystyle b=\frac{\lambda^*\mu^*}{\lambda^*+\mu^*}=\frac{1-\min\{\lambda^*,\mu^*\}}{1+\max\{\lambda^*,\mu^*\}-\min\{\lambda^*,\mu^*\}}$.
    \item $\displaystyle\rho^*=\frac{\mu^*}{\lambda^*+\mu^*}$.
    \item $f_R(x^*,y^*)=f_C(x^*,y^*)=b$.
    \item $f_C(w^*,y^*)=f_R(x^*,z^*)=1$.
    \item $\supp(x^*)\subseteq\suppmin(A(\rho^*,y^*,z^*))$ and $\supp(y^*)\subseteq\suppmin(B(\rho^*,x^*,w^*))$.
    \item $\begin{cases}
S_C(x^*)\cap S_C(w^*)\neq\varnothing,& \text{if } f_C(w^*,z^*)> f_R(w^*,z^*),\\
S_R(y^*)\cap S_R(z^*)\neq\varnothing,& \text{if } f_C(w^*,z^*)<f_R(w^*,z^*).
\end{cases}$
\end{itemize}
\end{theorem}
All statements can be verified by direct calculations. To avoid verbose calculations, we omit them. Now we make the second verification, i.e., to prove that $(x^*,y^*)$ is a strict stationary point.

A strict stationary point $(x^*,y^*)$ should have
\[Df(x^*,y^*,x',y')>0,\quad (x',y')\in\Delta_m\times\Delta_n\setminus\{(x^*,y^*)\}.\]
Note further that $Df(x^*,y^*,x^*,y^*)=0$. So equivalently, $Df(x^*,y^*,x',y')=0$ if and only if $(x',y')=(x^*,y^*)$. By \cite{DBLP:journals/corr/abs-2107-01471} and \cite{DBLP:journals/im/TsaknakisS08} (one can also see \Cref{app:TSAlgo}), calculating the steepest descent direction is formalized as dual linear programmings. So the condition becomes $(x^*,y^*)$ is the \emph{unique} optimal solution of the  primal linear programming.

We follow \cite{UniquenessofLPSolution} to check such uniqueness. In general, we consider the following primal LP of finding an $\bar{x}\in\mathbb{R}^n$:
\begin{gather}
    \text{Minimize}\quad \T{p}x \notag\\
    \begin{split}
    \text{s.t.}\quad Ax&=b,\\
    Cx&\geq d,
    \end{split}\label{LP:primal}
\end{gather}
where $p,b,d$ are given vectors in $\mathbb{R}^n,\mathbb{R}^m,\mathbb{R}^k$ respectively, and $A$ and $C$ are given $m\times n$ and $k\times n$ matrices respectively. The corresponding dual LP of finding $\bar{u}\in\mathbb{R}^m$ and $\bar{v}\in\mathbb{R}^n$ is
\begin{gather}
    \text{Maximize}\quad \T{b}u+\T{d}v\notag\\
    \begin{split}
        \text{s.t.}\quad \T{A}u+\T{C}v&=p,\\
        v&\geq 0.
    \end{split}\label{LP:dual}
\end{gather}
Let $C_i$ denote the $i$th row of $C$. Define
\[J=\{i:C_i\bar{x}=d_i\},\]
\[K=\{i:\bar{v}_i>0\},\]
\[L=\{i:C_i\bar{x}=d_i, \bar{v}_i=0\}.\]
Let $C_J$, $C_K$ and $C_L$ be matrices whose rows are $C_i$, $i\in J$, $i\in K$ and $i\in L$ respectively. We have the following result:
\begin{theorem}[\cite{UniquenessofLPSolution}, Theorem 2, (v)]\label{thm:appTS:LP-unique}
Let $\bar{x}$ be the optimal solution of LP \Cref{LP:primal}. Then $\bar{x}$ is the unique optimal solution if and only if both statements below hold:
\begin{itemize}
    \item The rows of $[\T{A}\ \T{C_K}\ \T{C_L}]$ are linear independent.
    \item There is no $x$ satisfying 
    \[Ax=0,\quad C_Kx=0,\quad C_Lx\geq 0.\]
\end{itemize}
\end{theorem}

Using this theorem, we can directly verify by pure calculations that $(x^*,y^*)$ is the unique optimal solution of the primal LP, i.e., $(x^*,y^*)$ is a strict stationary point. Again, we omit the calculation details.
\section{Discussion}\label{sec:discussion}
In this paper, we have studied parametric and convergence stability in nonconvex and nonsmooth optimization. Our conclusion is qualitative and applies to a wide range of functions and optimization algorithms. We summarize them as the following table.
\begin{table}[ht]
    \centering
    \captionsetup{singlelinecheck=off}
    \caption{Summary of our results on two parametric and convergence stability. \checkmark means that such stability always holds. $\times$ means that such stability never holds. ? means that there exist counterexamples so that such stability does not hold. }

\begin{threeparttable}[ht]
    \renewcommand{\arraystretch}{1.25}
    \begin{tabular}{|c|c|c|c|}
    \hline
    Class of Function $f$ &Condition on $x^*$& Parametric Stability&Convergence Stability\\
    \hline
        D &A&\checkmark(\Cref{thm:unique-local-min-is-parametric-stable})&?\footnotemark[26]\\
    \hline
        D&C&?\footnotemark&?\footnotemark[26]\\
    \hline
        $\mathscr{S}$&A and B&\checkmark(\Cref{thm:unique-local-min-is-parametric-stable})&\checkmark(\Cref{thm:max-type-is-convergence-stable})\\
    \hline
        $\mathscr{A}$&A&\checkmark(\Cref{thm:unique-local-min-is-parametric-stable})&\checkmark(\Cref{thm:analytic-max-is-convergence-stable})\\
    \hline
        E&A&\checkmark(\Cref{thm:unique-local-min-is-parametric-stable})&?\footnotemark\\
    \hline
        $\mathscr{S}$&C&?\footnotemark[25]&$\times$(\Cref{thm:max-type-is-convergence-stable})\\
    \hline
    $\mathscr{A}$&C&?\footnotemark[25]&$\times$(\Cref{thm:analytic-max-is-convergence-stable})\\
    \hline
    \end{tabular}
    \begin{tablenotes}
    \item A: $x^*$ is the unique local minimum point in some neighborhood.
    \item B: $x^*$ is the unique stationary point in some neighborhood.
    \item C: $x^*$ is not a local minimum point.
    \item D: $f$ is continuous on $\Omega\times\mathscr{R}$.
    \item E: Finite compositions of smooth mappings and the maximum function.
    \end{tablenotes}
    \label{tab:my_label}
\end{threeparttable}
\end{table}
\footnotetext[25]{For counterexamples, see \cref{construct:saddle-point-not-parametric-stable} and \Cref{app:construct-parametric-stable}.}
\footnotetext[26]{For counterexamples, see \Cref{ex:nonex-for-converge-stable}.}

This work studies some fundamental problems in the algorithmic issues in optimization. Our results reveal some understanding on 
Nash equilibrium computation, nonconvex and nonsmooth optimization, as well as the new optimization methodology of deep neural networks. Questions raise themselves for future explorations. The next immediate open issues include the following future topics.
\begin{enumerate}[fullwidth, listparindent=\parindent]
    \item Further \emph{geometric} and \emph{qualitative} understandings on deep neural networks.
    
   The main methods used to study deep neural networks are probability theory and analysis theory. Deep techniques are developed to do convergence and performance analyses on various neural network models and optimizers \cite{DBLP:conf/icml/Allen-ZhuLS19,DBLP:conf/icml/DuLL0Z19,DBLP:conf/iclr/DuZPS19,DBLP:conf/nips/JacotHG18}. Many results are proved by heavy calculations and delicate inequalities. So these results are very quantitative and able to resolve problems for a specific kind of neural networks and a specific kind of optimizers. Cooper's works \cite{DBLP:journals/corr/abs-2005-04210, DBLP:journals/corr/abs-1804-10200} and our work provide another extremity: a fully qualitative understanding on very general neural networks via geometry theory. There are far fewer calculations. Even though very inaccurate, the results can apply to general neural networks. For instance, from our work we know that it is analyticity rather than merely smoothness making a deep neural network special in convergence stability. Combine understandings from both extremities and we will get a better understanding of the ``black magic'' behind deep neural networks.
   
   \item More \emph{quantitative} analysis on parametric stability and convergence stability.
   
   Our work applies to very general object functions and optimization algorithms; the price for such generality is a lack of quantitative results. Our results are a double-edged sword: It gives a better understanding of stability for general settings while providing less quantitative guidance. For example, if we wish to guarantee a $(\lambda,r)$-convergence stable stationary point, how should we choose the step size and how large can $r$ be? \Cref{prop:quan-result-of-convergence-stable} presents a calculation method, but for most functions, such calculation is almost impossible. If we restrict to a specific class of functions, it is easier to do practical calculations.
   
   \item \emph{Probabilistic considerations} on parametric stability and convergence stability.
   
   Most previous theories, such as learning theory \cite{DBLP:conf/icml/Allen-ZhuLS19,DBLP:conf/icml/DuLL0Z19,DBLP:conf/iclr/DuZPS19,DBLP:conf/nips/JacotHG18} and smoothed analysis \cite{SpielmanCommunACM}, involve the probability into their settings. For instance, a minor perturbation means a mean-zero Gaussian with a small variance in these theories, rather than a small positive number in our work. Involving probability into settings would make further comparisons between these theories possible and thus help us more completely understand stability. Furthermore, probability is conducive to improving local results to global results.
   
   \item Relationships between \emph{trembling-hand-perfect} Nash equilibria and convergence stability.
   
   Trembling-hand-perfection is a typical refinement of Nash equilibrium due to Selten \cite{SeltenPerfectness}. If there is a small uncertainty in the prediction about the opponent's play, weakly dominated strategies are suboptimal and are never chosen by players. Under such uncertainty, only a subset of Nash equilibria are possible. If the TS algorithm finds a Nash equilibrium, a slight perturbation of strategies (possibly due to numerical inaccuracy) is equivalent to a small uncertainty faced by both players. If the Nash equilibrium is convergence-stable, the descent procedure will return to the same Nash equilibrium. Thus the descent direction is still the optimal strategy for both players. So this Nash equilibrium is trembling-hand-perfect. On the other hand, if a Nash equilibrium is trembling-hand-perfect, any tiny uncertainty does not influence both players' choices. Hence the Nash equilibrium is insensitive to small changes in players' strategies and thus convergence-stable. Therefore, it is reasonable to conjecture that every Nash equilibrium found by the TS algorithm is almost surely trembling-hand-perfect.
   
 
\end{enumerate}


We hope our work makes a stepstone toward a full understanding of convergence stability and parametric stability of optimization algorithms in general, especially of deep neural networks. 
\section{Applications}\label{sec:app}
In this section, we present three applications of the general theory. The first is a stability analysis of the state-of-the-art approximate algorithm for finding Nash equilibria in two-player games, the Tsaknakis-Spirakis algorithm. The object function involved in the TS algorithm is much simpler than a general deep neural network — thus, it can be quantitatively analyzed — but is complicated enough to show how to do a nontrivial analysis using our theory. The second one provides a qualitative conclusion on the non-escapable solutions in optimization theory. Finally, the third one presents the connection between general theory and optimization of deep neural networks. 


\section{Constructing Counterexamples}
\subsection{Constructing a Counterexample for Parametric Stability from Smooth Curve}\label{app:construct-parametric-stable}
Suppose local minimum points form a smooth simple curve $\gamma$ diffeomorphic to $[0,1]$. We construct a function $f(x;\lambda)$ such that
\begin{itemize}
    \item when $\lambda=0$, its stationary points form curve $\gamma$,
    \item and every stationary point $x^*\in\gamma$ of $f(\cdot;0)$ is not parametric-stable.
\end{itemize}
We first introduce some concepts from smooth manifolds. See \cite{lee2013sard} for details.

For each $x\in\mathbb{R}^n$, the \emph{tangent space} $ T_x\mathbb{R}^n$ is canonically identified with $\mathbb{R}^n$ itself. It can be viewed as the space of tangent vectors at $x$. Similarly, suppose $M$ is an embedded $m$-dimensional smooth submanifold, its tangent space at $x$ is denoted by $T_xM$. Note that as a Euclidean space, tangent spaces inherit a Euclidean dot product. For each $x\in M$, we define the \emph{normal space to $M$ at $x$} to be the $(n-m)$-dimensional subspace $ N_xM\subseteq  T_x\mathbb{R}^n$ consisting of all vectors that are orthogonal to $ T_xM$ with respect to the Euclidean dot product. The \emph{normal bundle of $M$}, denoted by $NM$, is defined as:
\[NM=\{(x,v)\in\mathbb{R}^n\times\mathbb{R}^n:v\in N_xM\}.\]

Thinking of $NM$ as a submanifold of $\mathbb{R}^n\times\mathbb{R}^n$, we define $E:NM\to\mathbb{R}^n$ by
\[E(x,v)=x+v.\]
This maps each normal space $N_xM$ affinely onto the affine subspace through $x$ and orthogonal to $ T_xM$ . Clearly, $E$ is smooth, because it is the restriction to $NM$ of the addition map $\mathbb{R}^n\times\mathbb{R}^n\to\mathbb{R}^n$. A \emph{tubular neighborhood of $M$} is a neighborhood $U$ of $M$ in $\mathbb{R}^n$ that is the diffeomorphic image under $E$ of an open subset $V\subseteq NM$ of the form
\[V=\{(x,v)\in NM:\norm{v}<\delta(x)\},\]
for some positive continuous function $\delta:M\to\mathbb{R}$.

Our construction relies on the following essential result:
\begin{theorem}[Tubular Neighborhood Theorem, \cite{lee2013sard}, Theorem 6.24]\label{thm:tabular-for-manifold}
Every embedded smooth submanifold of $\mathbb{R}^n$ has a tubular neighborhood.
\end{theorem}

Note that smooth curve $\gamma$ is a one-dimensional smooth submanifold \emph{with boundary} (not smooth submanifold). Thus we need to modify $\gamma$. To do so, we can extend the diffeomorphism at $0$ and $1$ to an open neighborhood of $0$ and $1$ due to smooth invariance of the boundary (one can refer to Theorem 1.46 in \cite{Lee2012}). Then $\gamma$ is extended to $\Tilde{\gamma}:(-\epsilon_1,1+\epsilon_2)\to\mathbb{R}^n$ being a one-dimensional smooth submanifold.

Take a tubular neighborhood $U$ of $\Tilde{\gamma}$, which is diffeomorphic to $V\subseteq N\Tilde{\gamma}$ via a diffeomorphism $\varphi$. Specifically, $\varphi(x)=(u(x),v(x))$, where $u(x)\in \Tilde{\gamma}$ and $v(x)\in N_{u(x)}\Tilde{\gamma}$ for all $x\in U$.  Construct a smooth function $f(x;\lambda)$ on $U$ as follows:
\[f(x;\lambda):= \norm{v(x)}^2+\lambda\Tilde{\gamma}^{-1}(u(x)).\]

Consider the areas
\begin{align*}
    G_1&:=\varphi^{-1}\left(\left\{(x,v)\in NM:x\in\gamma,\norm{v}\leq\frac13\delta(x)\right\}\right),\\
    G_2&:=\varphi^{-1}\left(\left\{(x,v)\in NM:x\in\gamma,\norm{v}\leq\frac23\delta(x)\right\}\right).
\end{align*}
Let $G$ be a polytope (possibly non-convex) such that $G_1\subseteq G\subseteq G_2$. $G$ contains the whole $\gamma$ and its boundary contains the end points $\gamma(\{0,1\})$. Moreover, the restriction of $f(x;\lambda)$ on $G$ is in $\mathscr{F}$. 


Since $U$ is diffeomorphic to $V$, $x\in U$ is a stationary point of $f(\cdot;\lambda)$ if and only if $\varphi(x)\in V$ is a stationary point of $f(\varphi^{-1}(\cdot);\lambda)$. Let $h(u,v)$ denote $f(\varphi^{-1}(u,v);\lambda)$, then $h(u,v)=\norm{v}^2+\lambda\Tilde{\gamma}^{-1}(u)$. It is clear that for any $(u,v)\in V$ such that $v\neq 0_n$, $(u,v)$ is not a stationary point of $f(\cdot;\lambda)$. For all $t\in[0,1]$, $\partial h(\gamma(t),0_n)/\partial t=\lambda$.


Now consider the restriction of $f(x;\lambda)$ in $G$. The set of stationary points of $f(\cdot;0)$ in $G$ is exactly $\gamma$. However, for any $\lambda<0$, the only stationary point of $f(\cdot,\lambda)$ in $G$ is $\gamma(1)$. Meanwhile, for any $\lambda>0$, the only stationary point of $f(\cdot,\lambda)$ in $G$ is $\gamma(0)$. Therefore, for any $x^*\in \gamma$, $x^*$ is not parametric-stable for $\lambda=0$.
\begin{comment}
Define a function $g:\mathbb{R}^{n-1}\to\mathbb{R}$ as 
\[g(x)=\begin{cases}
\norm{x}^2\log\left(1+\norm{x}^{-2}\right),&\norm{x}\neq 0,\\
0,&\norm{x}=0.
\end{cases}\]
$g(x)$ is differentiable everywhere and the unique stationary point (also the unique local minimum point) is $x=0$. Furthermore, $0\leq g(x)<1$. Take a tubular neighborhood $U$ of $\Tilde{\gamma}$, which is diffeomorphic to $V\subseteq N\Tilde{\gamma}$ via a diffeomorphism $\varphi$. Define $h:V\to\mathbb{R}$ as $h(x,v)=g(v)$. Restrict the the definition of $h(\varphi(x))$ to a closed area $G:=\varphi^{-1}\left(\{(x,v)\in NM:x\in\gamma,\norm{v}\leq\delta(x)/2\}\right)$. $G$ contains the whole $\gamma$ and its boundary contains end points $\gamma(\{0,1\})$. Construct $f(x;\lambda)$ on $G$ as follows:
\[f(x;\lambda)=h(\varphi(x))-\lambda\norm{x-\gamma(0)}.\]
Replace $G$ with a smaller compact area (for instance, replace $\delta$ with piecewise-linear function in the definition of $G$) so that $f(x;\lambda)\in\mathscr{F}$. Note that the only stationary points (only local minimum points as well) in $G$ of $f(\cdot;0)$ exactly form the curve $\gamma$. If $\lambda>0$, the minimum points form the set $\Psi=\argmax_{x'\in\gamma}\norm{\gamma(0)-x'}$, whose distance to $\gamma(0)$ is a constant $r>0$ due to closeness of $\gamma$. Now, a minute change on $\lambda$ leads to a giant slide from $\gamma(0)$ to $\Psi$. So $\gamma(0)$ is not parametric-stable for $f(\cdot;0)$.
\end{comment}

\subsection{Verifying (Counter-)\Cref{ex:nonex-for-converge-stable} for Convergence Stability}\label{app:nonexample-convergence-stable}
We do the verification step by step.
\begin{enumerate}
    \item \emph{If $x>0$,
    \begin{align*}
    g'(x)&=x^{-2}e^{-1/x}\sin\left(1/x\right)-x^{-2}e^{-1/x},\\
    h'(x)&=x^{-2}e^{-1/x}\sin\left(1/x\right)-x^{-2}e^{-1/x}\cos\left(1/x\right).
\end{align*}}

    It follows by direct calculations.
    
    \item \emph{For both $g(x)$ and $h(x)$, the differential of every order at $x=0$ is $0$.}
    
    Note that for $x\neq 0$, the $n$th derivative of $g(x)$, denoted by $g^{(n)}(x)$, can be expressed in $\mathrm{poly}(1/x,\sin(1/x),\cos(1/x))e^{-1/|x|}$. Clearly $|g^{(n)}(x)/x|\to 0$ as $x\to 0$. Thus $g^{(n)}(0)=0$ for all $n$. The same fact about $h(x)$ can be proved by a similar argument.
    
    \item \emph{$g(x)$ and $h(x)$ are smooth on $\mathbb{R}$.}
    
    It is clear that $g(x)$ and $h(x)$ are smooth when $x\neq 0$. The case when $x=0$ is verified previously.
    
    \item \emph{$f(x)$ is finite compositions of smooth mappings and the maximum function.}
    
    Note that $\min\{a,b\}=-\max\{-a,-b\}$, $|x|=\max\{x,-x\}$. In addition, $g(x)$ and $h(x)$ are smooth. The result follows.
    
    \item \emph{$f'_{-}(x)>0$ for $\in(0,1)$, so $f(x)$ is strictly increasing on $(0,1)$. These hold symmetrically on $(-1,0)$. So $x^*=0$ is the only stationary point in $B(0,1)$, and is the unique local minimum point.}
    
    Suppose $x>0$. $\sin(1/x)<0$ if and only if $x\in\left(\frac{1}{2k\pi}, \frac{1}{(2k-1)\pi}\right)$, $k=1,2,\ldots$ Thus
    \[f_-'(x)=\begin{cases}
    h'(x)-g'(x),& x\in\left(\frac{1}{2k\pi}, \frac{1}{(2k-1)\pi}\right),k=1,2,\ldots,\\
    -g'(x),& x\in\left(\frac{1}{2(k+1)\pi}, \frac{1}{2k\pi}\right), k=1,2,\ldots,\\
    \max\{0,h'(x)\}-g'(x),& x=\frac{1}{k\pi}, k=1,2,\ldots \end{cases}\]
    
    We discuss case by case to show that $f_-'(x)>0$.
    
    If $x\in\left(\frac{1}{2k\pi}, \frac{1}{(2k-1)\pi}\right)$, $f_-'(x)=h'(x)-g'(x)=x^{-2}e^{-1/x}(1-\cos(1/x))>0$.
    
    If $x\in\left(\frac{1}{2(k+1)\pi}, \frac{1}{2k\pi}\right)$, $f_-'(x)=-g'(x)=x^{-2}e^{-x}(1-\sin(1/x))>0$.
    
    If $x=\frac{1}{(2k-1)\pi}$, $f_-'(x)=h'(x)-g'(x)=2x^{-2}e^{-1/x}>0$.
    
    If $x=\frac{1}{2k\pi}$, $f_-'(x)=-g'(x)=x^{-2}e^{-1/x}>0$.
    
    Note that $f(x)$ is a symmetric function, so $f_+'(x)<0$ for $x\in(-1,0)$.
    
    \item \emph{$x^*=0$ is not convergence-stable for $f(x)$.}
    
    Observe that for $k=1,2\ldots$ and $t\in(0,\pi)$, 
\begin{align*}
    f'\left(\frac{1}{2k\pi-t}\right)&=(2k\pi-t)^2e^{2k\pi-t}(1-\cos(2k\pi-t))\\
    &=(2k\pi-t)^2e^{2k\pi-t}(1-\cos(t))\\
    &\leq \frac{1}{2}(2k\pi)^2e^{2k\pi}t^2=C(k)t^2.
\end{align*}
Take any $r_0>0,\lambda>0$. Let $k=\left\lceil\frac1{2\pi r_0}\right\rceil+1$ and $\epsilon=\frac1{2k\pi+1}$. For any $\delta>0$, let $t=\min\{\sqrt{C^{-1}(k)\delta},1\}$. Suppose an algorithm $\mathcal{A}_{\lambda,\delta}$ starts at $x_0=\frac1{2k\pi-t}\in B(0,r_0)$. It may stop at $x_0$ permanently because $x_0$ is a $\delta$-stationary point. However, $x_0=\frac1{2k\pi-t}\notin B(0,\epsilon)$. This violates the definition of $(r_0,\lambda)$-convergence stability, so $x^*=0$ is not convergence-stable.
\end{enumerate}

\subsection{A General Conclusion on Non-Escapable Solutions to Continuous Optimization}\label{subsec:appGeneralOptimization}
In this application, we consider optimization algorithms for finding a local minimum point of the objective function in $\mathscr{C}$ or $\mathscr{A}$. Such algorithms should be classified into nonsmooth optimization and classical convex optimization. Examples of these algorithms are gradient descents, Newton methods, bundle methods, subgradient methods, gradient sampling methods \cite{bagirov2014introduction, boyd2004convex}. Unlike optimizers in deep learning, their convergence theories have been well-developed.

We say a solution is non-escapable, if when the optimization procedure gets into a neighborhood of the solution, it cannot escape and will eventually converge to the solution. A non-escapable solution is less sensitive to the precision parameter and numerical errors. Thus such solutions are commonly found solutions in practice. For most algorithms, a solution is equivalent to an approximate stationary point. In our theory, a non-escapable stationary point is by definition a convergence-stable stationary point.

To apply our theory, we need to verify that these algorithms are ultimately decreasing, $\delta$-result-stationary and $\lambda$-path-bounded. To provide a provable convergence guarantee, most algorithms are designed to be \emph{descent}, i.e., the point sequence $\{x_k\}_{k\geq 0}$ iteratively generated by the algorithm satisfies $f(x_{k+1})\leq f(x_k)$. Descent algorithms are undoubtedly ultimately decreasing. Also, most algorithms are proved to be $\delta$-result-stationary, i.e., they will find $\delta$-stationary points as all the cluster points of $\{x_k\}$.\footnote{$f$ may have a minimum of $-\infty$. So rigorously speaking, $\{x_k\}$ is $\delta$-result-stationary or tends to $-\infty$. Since we consider stationary points on a compact area, such infinity will not appear.} We have to show they are $\lambda$-path-bounded. The basic idea is to use sufficiently small step sizes. By Cantor's uniform continuity theorem, $f$ on the compact domain $\Omega$ is uniformly continuous, i.e., for all $\epsilon>0$, there exists a $\delta(\epsilon)$ such that for all $x,y\in \Omega$
\[\norm{x-y}<\delta\Longrightarrow |f(x)-f(y)|<\epsilon.\]
We can ensure the algorithm to be $\lambda$-path-bounded by setting the step length $\norm{x_{k+1}-x_k}\leq\delta(\lambda)$.

Our main stability results for these algorithms on functions $\mathscr{C}$ or $\mathscr{A}$ are presented in \Cref{thm:max-type-is-convergence-stable} and \Cref{thm:analytic-max-is-convergence-stable}. They have clear geometric meanings: An isolated local minimum point is stabler than a saddle point or a non-isolated local minimum point. It implies that with small enough step sizes and precision parameters, non-escapability is mainly determined by geometric structures of the stationary point found by the algorithm.
\section{Intrinsic Geometries behind Stability}\label{sec:prop-of-stability}
In this section, we show the deep connection between stability and geometry theory. We utilize such connection to prove some necessary and sufficient conditions for the two notions of stability. These conditions help to understand how it looks like if a stationary point is parametric-stable or convergence-stable. We will see that these two notions of stability are the inherent nature of \emph{objective functions} and depend little on optimization algorithms. Most results are intuitive. To make reading more accessible, we postpone technical proofs to the appendix and only present essential proofs capturing the main idea.

\subsection{Parametric Stability and Deformations}
As examples of parametric stability suggest, a parametric-stable stationary point should be insensitive to minor perturbations of parameters. From the geometric perspective, changes in parameters cause the graph of the function to deform. A parametric-stable stationary point cannot slide to faraway points due to deformations. Thus it should be trapped in the interior of a bowl-shaped area. We first introduce that concept of \emph{equicontinuity}.

\begin{definition}
Let $F:X\times Y\to\mathbb{R}$ be a function, where $X$ and $Y$ are subsets of Euclidean spaces. We say $F(x,y)$ is \emph{equicontinuous} in $y$ if
for every $\epsilon>0$, there exists $\delta>0$ such that for any $x\in X$, $y_1,y_2\in Y$,
\[\norm{y_1-y_2}<\delta\Longrightarrow |F(x,y_1)-F(x,y_2)|<\epsilon.\]
\end{definition}

Then we are able to formalize the previous ``bowl-shaped area'' idea into \Cref{thm:unique-local-min-is-parametric-stable}.

\begin{theorem}\label{thm:unique-local-min-is-parametric-stable}
Let $f(x;R)\in\mathscr{F}$ with domain $\Omega$ and parameter space $\mathscr{R}$. $\Omega$ is a closed area. $f(x;R)$ is equicontinuous in $R$. Suppose stationary point $x^*$ of $f(\cdot;R)$ satisfies that for all $\epsilon>0$, there exists an open neighborhood $U_\epsilon$ of $x^*$, such that $U_\epsilon\subseteq B(x^*,\epsilon)$, and that for any $x'$ in $\Omega\cap\partial U_\epsilon$, $f(x';R)>f(x^*;R)$. Here $\partial U_\epsilon$ denotes the boundary of $U_\epsilon$. Then $x^*$ is parametric-stable for parameter $R$.
\end{theorem}
We have the following direct and intuitive corollary:
\begin{corollary}
Let $f(x;R)\in\mathscr{F}$ with domain $\Omega$ and parameter space $\mathscr{R}$. $\Omega$ is a closed area. $f(x;R)$ is equicontinuous in $R$. Suppose stationary point $x^*$ of $f(\cdot;R)$ satisfies is the unique local minimum point of $f(\cdot;R)$ in $B(x^*,r)\cap \Omega$ for some radius $r>0$, then $x^*$ is parametric-stable for parameter $R$.
\end{corollary}

Before the proof of \Cref{thm:unique-local-min-is-parametric-stable}, we comment a little about this theorem. First, domain $\Omega$ of function $f$ only need to be a closed area. But it can be unbounded and have many holes. Second, equicontinuity is a reasonable condition: Lots of functions in $\mathscr{F}$ satisfy this condition due to the following lemma\footnote{This is a standard result from mathematical analysis. For completeness,  we prove it in \Cref{app:proof-equicontinuity}.}:
\begin{lemma}\label{lemma:equicontinuity}
Let $f(x;R)\in\mathscr{F}$ with domain $\Omega$ and parameter space $\mathscr{R}$. $\Omega\times\mathscr{R}$ is a compact set. $f$ is continuous on $\Omega\times\mathscr{R}$. Then $f$ is equicontinuous in parameter $R$.
\end{lemma}
Last, it is a sufficient condition and in fact rules out many other kinds of stationary points, such as saddle points and local minimum points forming a curve. Although it is not a necessary condition, it is likely to be the \emph{best} condition we can get on a general function. More precisely, there exists a function $f(\cdot;R)$ with a saddle point failing to be parametric-stable. Likewise, such existence of non-parametric-stable stationary points holds for  functions with local minimum points forming a simple smooth curve\footnote{A simple smooth curve is a diffeomorphic function from $[0,1]$ to $\mathbb{R}^n$. Its image is a smooth submanifold in $\mathbb{R}^n$.}.

To see a saddle point would not be parametric-stable, we make the following simple construction:
\begin{construction}\label{construct:saddle-point-not-parametric-stable}
Write $x^*$ as $(x_1^*,\ldots,x_n^*)$. Let 
\[f(x;\lambda):=\sum_{k=1}^{n-1}(x_k-x_k^*)^2+(x_n-x^*_n)^3+\lambda (x_n-x^*_n),\]
$\lambda\in[-1,1]$. It resembles a saddle surface. $x^*=0$ is the \emph{unique} stationary point of $f(\cdot;0)$ and is indeed a saddle point. Similar to \Cref{ex:saddle-point-is-not-parametric-stable}, stationary point $x^*=0$ is not parametric-stable for $f(\cdot;0)$.
\end{construction}

If local minimum points form a simple smooth curve $\gamma$, we can construct function $f(x;\lambda)$ such that:
\begin{enumerate*}[label=(\arabic*)]
    \item All stationary points of $f(\cdot;0)$ are local minimum points and the form curve $\gamma$;
    \item Every stationary point $x^*$ of $f(\cdot;0)$ is not parametric-stable.
\end{enumerate*}
Such a construction need techniques from smooth manifolds. For interested readers, see \Cref{app:construct-parametric-stable}.

Now we prove \Cref{thm:unique-local-min-is-parametric-stable}. The main idea of the proof is to pick a shell $S_{r_0}$ surrounding $x^*$ but at a distance away from $x^*$ and ensure all values in set $f(S_{r_0};R)$ are strictly greater than $f(x^*;R)$. Then by continuity, this property still holds if the graph of $f$ on shell $S_{r_0}$ is deformed under slight perturbations on $R$. 
\begin{proof}[Proof of \Cref{thm:unique-local-min-is-parametric-stable}]
A geometric illustration of the proof is presented in \Cref{fig:illustration-proof-parametric-stable}.

\begin{figure}[ht]
    \centering
    \input{Figures/illustration-proof-parametric-stable}
    \caption{Illustration of the proof of \protect{\Cref{thm:unique-local-min-is-parametric-stable}}.}
    \label{fig:illustration-proof-parametric-stable}
\end{figure}

Below we interpret the geometric proof and fill the missing details. For any $\epsilon>0$, by continuity of $f(\cdot;R)$, there exists $r_0>0$ such that
\[|f(x;R)-f(x^*;R)|<\epsilon/2\]
whenever $\norm{x-x^*}<r_0$. Replace $r_0$ with $\min\{r_0,0.9\epsilon\}$. By assumption of the theorem, choose $U_{r_0}\subseteq B(x^*,r_0)$, an open neighborhood of $x^*$ so that $S_{r_0}:=\partial U_{r_0}\cap \Omega$ satisfies that for all $x'$ in $S_{r_0}$, $f(x';R)>f(x^*;R)$. Note that $S_{r_0}$ is compact, then by Weierstrass' Minimum Theorem, 
\[M_{r_0}:=\min_{y\in S_{r_0}} f(y;R)>f(x^*;R).\] 
Let $\epsilon_0:=\min\{\epsilon/2,(M_{r_0}-f(x^*;R))/2\}$, which is $>0$.
By equicontinuity, there exists a $\delta_0>0$ such that for any $x\in \Omega$, $R'\in\mathscr{R}$ and $\norm{R-R'}<\delta_0$,
\[ |f(x;R)-f(x;R')|<\epsilon_0.\]

Consider $\bar{U}_{r_0}\cap \Omega=(U_{r_0}\cup S_{r_0})\cap \Omega$. $\bar{U}_{r_0}\cap \Omega$ is a compact set and thus $f(\cdot;R')$ has a minimum point $x_0$ on $\bar{U}_{r_0}\cap \Omega$. For all $x'\in S_{r_0}$, we have
\[f(x_0;R')\leq f(x^*;R')< f(x^*;R)+\epsilon_0\leq M_{r_0}-\epsilon_0\leq f(x';R)-\epsilon_0< f(x';R').\]
Every $x'\in S_{r_0}$ satisfies $f(x_0;R')<f(x';R')$, so $x_0\notin S_{r_0}$. Therefore $x_0$ is in open set $U_{r_0}$, and $x_0$ is a stationary point.

Finally, note that $|f(x_0;R')-f(x_0;R)|<\epsilon/2$, and $|f(x_0;R)-f(x^*;R)|<\epsilon/2$, so we have $|f(x_0;R')-f(x^*;R)|<\epsilon$. Since $U_{r_0}\subseteq B(x^*,r_0)$, $\norm{x^*-x_0}<r_0\leq0.9\epsilon<\epsilon$. So $x^*$ is parametric-stable.
\end{proof}

Note that if $f$ is differentiable of second order at stationary point $x^*$, then one can verify $x^*$ is parametric-stable by showing that the Hessian at $x^*$ is positive definite. However, for nonsmooth functions, it would be hard to verify the conditions in \Cref{thm:unique-local-min-is-parametric-stable}. For many problems, we can calculate the approximate or exact value of $\partial f$. It would be useful if we could utilize the differential information to show that a stationary point is parametric-stable. We have the following proposition:

\begin{proposition}\label{prop:strict-sp-is-parametric-stable}
Let $f(x;R)\in\mathscr{F}$ with domain $\Omega$ and parameter space $R$. $f(x;R)$ satisfies Lipschitz condition in $x$ and is equicontinuous in $R$. $\Omega$ is a compact area. If $x^*$ is a strict stationary point of $f(\cdot;R)$, i.e., $\partial f(x^*;R)/\partial s>0$ for every valid unit vector $s$, then $x^*$ is parametric-stable.
\end{proposition}

It is again a very intuitive result: a strict stationary point should be a unique local minimum point. However, it is not trivially true. To understand the subtlety here, recall the definition of $\partial f(x^*)/\partial s$. $\partial f(x^*)/\partial s>0$ means that, $f(x^*+r s)-f(x^*)=r\partial f(x^*)/\partial s+\alpha(r)>\alpha(r)$, where $\alpha(r)/r$ tends to $0$ as $r\to 0$. So $f(x^*+r s)-f(x^*)>0$ only holds for small positive $r$. What we need is a \emph{uniform} $r_0$ such that $f(x^*+r s)-f(x^*)>0$ holds for any $r<r_0$ and any valid $s$. If $\Omega$ is a two-dimensional manifold, there are infinitely many directions $s$ to take into account. Thus such $r_0$ would not exist at all. Fortunately, if $f$ satisfies Lipschitz condition, such $r_0$ can be found. This result is stated as the following technical lemma. We leave the proof to \Cref{app:proof-Df>0circle}.

\begin{lemma}\label{lemma:Df>0circle}
Let $x^*$ be a strict stationary point. There exists $r>0$ and $c>0$ such that for any unit direction vector $s$ and $\theta\in(0,r)$, if $x^*+\theta s\in \Omega$ , then $ f(x^*+\theta s)> f(x^*)+c\theta>f(x^*)$.
\end{lemma}

Now we are ready to prove \Cref{prop:strict-sp-is-parametric-stable}.
\begin{proof}[Proof of \Cref{prop:strict-sp-is-parametric-stable}]
By \Cref{lemma:Df>0circle}, it follows immediately that there is $r>0$ such that for all $x'\in B(x^*,r)\cap \Omega$, $f(x';R)>f(x^*;R)$, i.e., $x^*$ is an isolated local minimum point, as desired.
\end{proof}

We can get a more quantitative result if function $f$ satisfies more conditions.
\begin{proposition}\label{prop:delta-parametric-stable}
Let $f(x;R)\in\mathscr{F}$ with domain $\Omega$ and parameter space $\mathscr{R}$. $\Omega$ is a closed area. Suppose $f(x;R)$ is $M$-Lipschitz in $x$ and uniformly $N$-Lipschitz in $R$, i.e.,
\begin{align*}
    |f(x;R)-f(y;R)|&\leq M\norm{x-y},\quad x,y\in \Omega,\\
    |f(x;R_1)-f(x;R_2)|&\leq N\norm{R_1-R_2},\quad x\in \Omega.
\end{align*}
If $x^*$ is the unique minimum point in $B(x^*,r)$ of $f(\cdot;R)$, then $x^*$ is $\min\{r, \epsilon, (M_{r_0}-f(x^*;R))\}/(2N)$-parametric-stable for $f(\cdot;R)$ where
\begin{gather*}
    r_0:=0.9\min\{r,r/(2M),\epsilon,\epsilon/(2M)\}\quad\text{and}\quad
    M_{r_0}:=\min_{\norm{x'-x^*}=r_0}f(x';R).
\end{gather*}
\end{proposition}
\begin{proof}
We refine the proof of \Cref{thm:unique-local-min-is-parametric-stable}. 

Suppose $\epsilon<r$. Let $r_0:=0.9\min\{\epsilon, \epsilon/(2M)\}$, by $M$-Lipschitz condition in $x$, $U_{r_0}=B(x^*,r_0)$ satisfies the condition in the proof of \Cref{thm:unique-local-min-is-parametric-stable}. Define $S_{r_0}$, $M_{r_0}$ and $\epsilon_0$ as in the proof of \Cref{thm:unique-local-min-is-parametric-stable}. By uniform Lipschitz condition in $R$, there exists a $\delta_0=\epsilon_0/N$ such that for any $x\in \Omega$ and $R'\in\mathscr{R}$ and $\norm{R-R'}<\delta_0$,
\[ |f(x;R)-f(x;R')|\leq N\norm{R-R'}<\epsilon_0.\]
Repeat the argument in the proof of \Cref{thm:unique-local-min-is-parametric-stable}. Then $x^*$ is parametric-stable with $\delta_0=\min\{\epsilon/2, (M_{r_0}-f(x^*;R))/2\}/N$. 

For arbitrary $\epsilon>0$, restrict it to $\min\{r,\epsilon\}$ and the desired result follows.
\end{proof}

\subsection{Convergence Stability and Area Partitioning}
Now we turn to convergence stability. So what kinds of stationary points could be convergence-stable? A saddle point cannot be stable: The optimization procedure will never come back if we perturb it to a descent direction of the saddle point. So the stationary point must be a local minimum point. If a local minimum point is not isolated, i.e., encompassed by infinitely many stationary points, the optimization procedure may get stuck to another distant stationary point. The two scenarios describe a necessary condition for convergence stability. Surprisingly, this is also a sufficient condition for many kinds of functions. Let $\Omega\subseteq\mathbb{R}^n$ be a compact area which is the domain for some $f\in\mathscr{F}$. We have the following result.

\begin{theorem}\label{thm:max-type-is-convergence-stable}
Let function class $\mathscr{S}:=\{\max_{i=1}^s\{f_i\}:f_i\in C^1(\Omega)\}$, i.e., class of max-type functions of $C^1$ functions. Let $f\in\mathscr{S}$ and $x^*$ is a stationary point of $f$. $x^*$ is convergence-stable if and only if $x^*$ is the unique stationary point as well as the unique local minimum point in $B(x^*,r)\cap\Omega$ for some radius $r$.
\end{theorem}

We say a function $f(x)$ defined on $\Omega\subseteq \mathbb{R}^n$ is \emph{finite compositions of analytic mappings and the maximum function}, if 
\[f(x)=f_{m}(f_{m-1}(\cdots f_{2}(f_{1}(x)))),\]
where $f_i:\Omega_{i-1}\to\Omega_{i}$, $\Omega_{i}\subseteq \mathbb{R}^{d_{i}}$, $d_{0}=n,d_{m}=1,\Omega_{0}=\Omega$. $f_i(z_1,z_2,\cdots,z_{d_{i-1}})$ is either an analytic mapping on $\Omega_{i-1}$ , or in the form of 
\[f_i(z_1,z_2,\cdots,z_{d_{i-1}})=\left(\max_{j\in S_{i,1}} z_j,\max_{j\in S_{i,2}} z_j,\cdots,\max_{j\in S_{i,d_{i}}} z_j\right),\] 
where $S_{i,k}\subseteq [d_{i-1}]$ for $k=1,\cdots,d_{i}$.

Then, a more universal version of \Cref{thm:max-type-is-convergence-stable} is

\begin{theorem}\label{thm:analytic-max-is-convergence-stable}
Let function class 
\begin{align*}
    \mathscr{A}:=\{f&\text{ is finite compositions of analytic mappings}\\
    &\text{and the maximum function on cube }\Omega=[-B,B]^n\}\footnotemark.
\end{align*}

Let $f\in\mathscr{A}$ and $x^*$ is a stationary point of $f$. $x^*$ is convergence-stable if and only if $x^*$ is
the unique local minimum point in $B(x^*,r)\Omega$ for some radius $r$.
\end{theorem}

\footnotetext{In fact, we only need the area $\Omega$ to be a compact \emph{semi-analytic set}. Such statement are too technical to follow, so we assume $\Omega$ is a cube for simplicity. But we will prove the semi-analytic case in \Cref{app:proof-analytic-max-convergence-stable}.}

Note that the condition in \Cref{thm:analytic-max-is-convergence-stable} only needs a unique local minimum point, while the condition in  \Cref{thm:max-type-is-convergence-stable} in addition requires a unique stationary point. Such results are, once again, quite intuitive. However, it is by no means trivial. To understand the subtlety, we present the following example:

\begin{example}
Let 
\begin{align*}
g(x)&:=\int_0^x 2u\sin(1/u) du,\\
h(x)&:=\begin{cases}x^2\sin(1/x),&x\neq 0,\\
0, x=0.\end{cases}\\
f(x)&:=\min\{0,h(x)\}-g(x)+|x|.
\end{align*}

We show that $x^*=0$ is the unique stationary point as well as the unique local minimum point in $B\left(0,\frac1\pi\right)$. In the mean time, however, it is not convergence-stable. 

Note that $g(x)$ and $h(x)$ are differentiable on $\mathbb{R}$:
\begin{align*}
    g'(x)&=\begin{cases}2x\sin(1/x),&x\neq 0,\\
    0,&x=0.
    \end{cases}\\
    h'(x)&=\begin{cases}2x\sin(1/x)-\cos(1/x),&x\neq 0,\\
    0,&x=0.
    \end{cases}
\end{align*} 
$f(x)$ can be expressed as finite compositions of differentiable mappings and the maximum function, but has infinitely many non-differentiable points in any neighborhood of $x=0$.

First, we show that if $x\in (0,1)$, then
\[f'_{-}(x)=\lim_{x'\to x-0}\frac{f(x')-f(x)}{x'-x}>0.\] 
Observe that for $x>0$, $h(x)=0$ if and only if $x=\frac1{k\pi}$ for some integer $k>0$.

For $x\in \left(\frac1{(2k+1)\pi},\frac1{(2k-1)\pi}\right]$, $k=1,2,\cdots$, we discuss the left and right derivative of $f$ at $x$.

If $x\in\left(\frac1{2k\pi},\frac1{(2k-1)\pi}\right)$, $f(x)=h(x)-g(x)+x$, so $f^{\prime}_{-}(x)=f^{\prime}_{+}(x)=f'(x)=h'(x)-g'(x)+1=1-\cos(1/x)>0$.

If $x\in\left(\frac1{(2k+1)\pi},\frac1{2k\pi}\right)$, $f(x)=0-g(x)+x$, so $f^{\prime}_{-}(x)=f^{\prime}_{+}(x)=f'(x)=-g'(x)+1=1-2x\sin(1/x)>1-2x>0$.

If $x=\frac1{(2k-1)\pi}$, $f^{\prime}_{-}(x)=1-\cos(1/x)=2,f^{\prime}_{+}(x)=1-2x\sin(1/x)=1$.

For $x=\frac1{2k\pi}$, $f^{\prime}_{-}(x)=1-2x\sin(1/x)=1,f^{\prime}_{+}(x)=1-\cos(1/x)=0$.

Therefore, for any $x\in(0,\frac1\pi)$, $f^{\prime}_{-}(x)>0$, so $x$ is not a stationary point. Similarly we can prove that for any $x\in\left(-\frac1\pi,0\right)$, $f^{\prime}_{+}(x)<0$, so $x$ cannot be a stationary point. It can be checked that $f^{\prime}_{+}(0)=1$ and $f^{\prime}_{-}(0)=-1$, so $x^*=0$ is a stationary point, and that $x^*=0$ is also a strict local minimum point.

However, we can show that $x^*=0$ is not convergence-stable. Observe that for $k=1,2,\cdots$, and $t\in(0,\pi)$, \[f'\left(\frac1{2k\pi-t}\right)=1-\cos\left(2k\pi-t\right)=1-\cos\left(t\right)\leq\frac{1}{2}t^2.\] 
Take any $r_0>0,\lambda>0$. Let $k=\left\lceil\frac1{2\pi r_0}\right\rceil+1$, and $\epsilon=\frac1{2k\pi+1}$. For any $\delta>0$, let $t=\min\{\sqrt{2\delta},1\}$. Suppose an algorithm $\mathcal{A}_{\lambda,\delta}$ starts at $x_0=\frac1{2k\pi-t}\in B(0,r_0)$. It may stop at $x_0$ permanently because $x_0$ is a $\delta$-stationary point. However, $x_0=\frac1{2k\pi-t}\notin B(0,\epsilon)$. This violates the definition of $(r_0,\lambda)$-convergence stability, so $x^*=0$ is not convergence-stable.
\end{example}

We further construct the following similar example which is finite compositions of \emph{smooth} mappings and the maximum function. Such example indicates that the analytic restriction is in some way necessary.
\begin{example}\label{ex:nonex-for-converge-stable}
Let 
\begin{align*}
g(x)&:=\begin{cases}\displaystyle\int_0^x \left(u^{-2}e^{-1/\left|u\right|}\sin\left(1/u\right)-u^{-2}e^{-1/\left|u\right|}\right) du,&x\neq 0,\\
0, &x=0.\end{cases}\\
h(x)&:=\begin{cases}e^{-1/\left|x\right|}\sin(1/x),&x\neq 0,\\
0, &x=0.\end{cases}\\
f(x)&:=\min\{0,h(|x|)\}-g(|x|).
\end{align*}
One can verify the following facts:
\begin{enumerate}
    \item For $x>0$,
    \begin{align*}
    g'(x)&=x^{-2}e^{-1/x}\sin\left(1/x\right)-x^{-2}e^{-1/x},\\
    h'(x)&=x^{-2}e^{-1/x}\sin\left(1/x\right)-x^{-2}e^{-1/x}\cos\left(1/x\right).
    \end{align*}
    \item For both $g(x)$ and $h(x)$, the differential of every order at $x=0$ is $0$.
    \item $g(x)$ and $h(x)$ are smooth on $\mathbb{R}$.
    \item $f(x)$ can be expressed as finite compositions of \emph{smooth} mappings and the maximum function.
    \item $f'_{-}(x)>0$ for $x\in(0,1)$, so $f(x)$ is strictly increasing on $(0,1)$. These hold symmetrically on $(-1,0)$. So $x^*=0$ is the only stationary point in $B(0,1)$, and is the unique local minimum point.
    \item However, $x^*=0$ is not convergence-stable for $f(x)$.
\end{enumerate}
The proof is similar to the previous example, so we put the detailed proof in \Cref{app:nonexample-convergence-stable}.
\end{example}

These examples show that if we do not make proper assumptions on functions to get rid of pathological functions, our primitive intuitions can be unreliable. Structures of functions in class $\mathscr{S}$ or $\mathscr{A}$ are more regular than the arbitrarily picked one in the general class $\mathscr{F}$. Utilizing these structures, we can prove theorems for them.

Now we state the proof. The geometric concept concerned in the proof is \emph{area partitioning}: Partition the domain so that in every part, approximate stationary points own good continuity and convergence properties. The proof is based on a technical lemma, which describes that the locations of approximate stationary points are near an exact stationary point. We prove \Cref{lemma:max-type-delta-SP-distance} for \Cref{thm:max-type-is-convergence-stable} below. The proof involves fewer technical details but captures the essential idea of the proof. The corresponding lemma for the proof of \Cref{thm:analytic-max-is-convergence-stable} is left to \Cref{app:proof-analytic-max-convergence-stable}, which deeply involves techniques in real analytic manifolds and semi-analytic sets. The key property distinguishing analyticity from smoothness is that there exists a delicate area partitioning which guarantees the isolated local minimum point to be an isolated stationary point.

\begin{lemma}\label{lemma:max-type-delta-SP-distance}
For any max-type function $f(x)=\max_{i=1}^{s}f_i(x)$, where $f_i(x)\in C^1(\Omega)\cap\mathscr{F}$ and domain $\Omega$ is a compact set.
If a stationary point $x^*$ satisfies that there are some $r>0$ such that there is no other stationary point in $B(x^*,r)\cap\Omega$, then there exists $r_1>0$ such that for every $\varepsilon>0$, there exists $\delta_0>0$, such that for $\delta\in(0,\delta_0)$, every $\delta$-stationary point in $B(x^*,r_1)\cap\Omega$ is in $B(x^*,\varepsilon)\cap\Omega$.
\end{lemma}
\begin{proof}
To avoid verbose notations, every ball $B(x,r)$ in the proof are replaced by $B(x,r)\cap\Omega$. Such replacement does not affect the correctness of the proof. Let \[Q_i:=\{x\in\Omega:f(x)=f_i(x)\}.\]
Then $\Omega$ is partitioned into closed sets $Q_1,\cdots,Q_s$. See \Cref{fig:illustration-proof-convergence-stable} for a demonstration. 
\begin{figure}[ht]
    \centering
    \input{Figures/illustration-proof-convergence-stable}
    \caption{Illustration of partitioning in \protect{\Cref{lemma:max-type-delta-SP-distance}}. Each colored area represents a closed set $Q_i$.}
    \label{fig:illustration-proof-convergence-stable}
\end{figure}

For every $x\in \Omega$, define
\[
    S_{x}:=\{i\in[s]:f_i(x)=f(x)\}\quad\text{and}\quad S^c_{x}:=\{i\in[s]:f_i(x)<f(x)\}.
\]
The existence of Dini directional derivatives is equivalent to the existence of directional derivatives. Thus below we use Dini derivatives instead. For $x,x'\in\Omega$, we have
\[Df(x,x')=\max_{i\in S_{x}}Df_i(x,x').\]

We prove the lemma by contradiction. Take $r_1=\frac{1}{2}r$. Suppose that there exists $\varepsilon_0>0$ such that for any $\delta_0>0$, there exists $\delta\in(0,\delta_0)$, $x\in B(x^*,r_1)\setminus B(x^*,\varepsilon_0)$, such that $x$ is a $\delta$-stationary point. Equivalently, there is an infinite point sequence $\{x_n\}_{n\geq 1}\subseteq B(x^*,r_1)\setminus B(x^*,\varepsilon_0)$, and $\lim_{n\to\infty}Df(x_n)\geq0$.\footnote{Recall that $Df(x)=\inf_{x'\neq x}Df(x,x')$. $Df(x)\geq 0$ is equivalent to $\partial f(x)\geq 0$.}

Take $U_{\left[\frac12\varepsilon_0,\frac23r\right]}:=\left\{x\in\Omega:\frac12\varepsilon_0\leq\lVert x-x^* \rVert\leq\frac23r\right\}$, then for all $n>0,x_n\in U_{\left[\frac12\varepsilon_0,\frac23r\right]}$. By its compactness, there is a convergent subsequence of $\{x_n\}$ that converges to a point $x_0\in U_{\left[\frac12\varepsilon_0,\frac23r\right]}$. So without loss of generality we assume $\{x_n\}$ is convergent, and $\lim_{n\to\infty}x_n=x_0$.


Since $f_i\in C^1(\Omega)$, $Df_i(x,x')$ is continuous in $x$. Thus for $x'\in\Omega,$
\[\lim_{n\to\infty}Df_{i}(x_n,x')=Df_{i}(x_0,x').\]

Let $r_2:=\min_{i\in S^c_{x_0}}\min_{x\in Q_i}\norm{x-x_0}$ (if $\cup_{i\in S^c_{x_0}}Q_i=\varnothing$, we set $r_2=1$). $r_2$ is the minimum radius that $f(x)$ is determined by $f_i$, $i\in\mathcal{S}_{x_0}$ in $B(x_0,r_2)$. Since $x_0\notin Q_i$ for $i\in S_{x_0}^c$, $\norm{x-x_0}>0$ whenever $x\in Q_i$. By compactness of $Q_i$, $r_2>0$. By convergence of $\{x_n\}$, For some large $N>0$, any $n>N$ satisfies $\norm{x_n-x_0}<r_2$. Therefore for $n>N,\mathcal S_{x_n}\subseteq S_{x_0}$. It follows that for all $x'\in\Omega$,
\begin{align*}
&Df(x_0,x')\\
=&\max_{i\in S_{x_0}}Df_i(x_0,x')\\
=&\max_{i\in S_{x_0}}\lim_{n\to\infty}Df_i(x_n,x')\\
=&\lim_{n\to\infty}\max_{i\in S_{x_0}}Df_i(x_n,x')\\
\geq&\limsup_{n\to\infty}\max_{i\in S_{x_n}}Df_i(x_n,x')\\
=&\limsup_{n\to\infty}Df(x_n,x')\\
\geq&\lim_{n\to\infty}Df(x_n)\\
\geq&0.
\end{align*}
The third equality is due to the continuity of the maximum function. Therefore $x_0$ is a stationary point distinct from $x^*$. Thus $x_0$ is another stationary point in $B(x^*,r)$, a contradiction.
\end{proof}

Now we prove \Cref{thm:max-type-is-convergence-stable} using \Cref{lemma:max-type-delta-SP-distance}. Note that this proof also applies to \Cref{thm:analytic-max-is-convergence-stable} if we instead use function class $\mathscr{A}$'s version of \Cref{lemma:max-type-delta-SP-distance}  (in a stronger form).
\begin{proof}[Proof of \Cref{thm:max-type-is-convergence-stable}]
To avoid verbose notations, every ball $B(x,r)$ in the proof are replaced by $B(x,r)\cap\Omega$. Such replacement does not affect the correctness of the proof.

\emph{Necessity}: Assume that $x^*$ is convergence-stable. Then there is some $r_0>0$ such that for any $x_0\in B(x^*,r_0)$ and $\varepsilon>0$, there exists sufficiently small $\delta,\lambda$ such that the sequence $\{x_k\}_{k\geq 0}$ generated by $\mathcal{A}_{\delta,\lambda}$ with initial point $x_0$ satisfies that every cluster point $\bar{x}$ of $\{x_k\}$ is in $B(x^*,\varepsilon)$. As $\mathcal A_{\delta,\lambda}$ is 
ultimately decreasing, we have $f(\bar{x})\leq f(x_0)$. Since $\varepsilon$ can be arbitrarily small, by the continuity of $f$, we have $f(x^*)\leq f(x_0)$. Furthermore, $f(x^*)=\min_{x\in B(x^*,r_0)}f(x)$ since $x_0$ is arbitrarily chosen. So $x_0$ is a minimum point in $B(x^*,r_0)$.

To prove the necessity, we claim that for any $x'\neq x^*$ in $B(x^*,r_0)$, $f(x')>f(x^*)$. Suppose otherwise that $x'$ is a local minimum point. Then $\partial f(x')/\partial s\geq0$ holds for any valid direction $s$. Thus $x'$ is also a stationary point. But $x'$ cannot be a stationary point, because otherwise the descent procedure starting from $x'$ will stay at $x'$, which violates the assumption.

\emph{Sufficiency}: Assume that there is some $r_0>0$ such that $f(x^*)=\min_{x\in B(x^*,r_0)}f(x)$, and for $x'\neq x^*$ in $B(x^*,r_0)$, $f(x')>f(x^*)$, and $\partial f(x')<0$.

By \Cref{lemma:max-type-delta-SP-distance}, without loss of generality, we may assume that $r_0$ is already sufficiently small, so that for any $\varepsilon>0$, there exists $\delta_0>0$, such that every $\delta_0$-stationary point in $B(x^*,r_0)$ is in $B(x^*,\varepsilon)$. For algorithm $\mathcal{A}_{\delta,\lambda}$ with parameter $\delta<\delta_0$, every cluster point of the generated point sequence $\{x_k\}_{k\geq 0}$ is almost surely a $\delta_0$-stationary point. Thus we only need to prove that $x_k\in B(x^*,r_0)$ holds for all $k\geq 0$.

Take a subset of $B(x^*,r_0)$:
\[U_{[0.8r_0,0.9r_0]}:=\{x\in B(x^*,r_0):\lVert x-x^*\rVert\in[0.8r_0,0.9r_0]\}.\]

As $U_{[0.8r_0,0.9r_0]}$ is a compact set, $\bar{f}:=\min_{x\in U_{[0.8r_0,0.9r_0]}}f(x)>f(x^*)$. Let $\lambda_0:=\frac{\bar{f}-f(x^*)}2$, and let $\mathcal{A}_{\delta,\lambda}$ has parameter $\lambda<\lambda_0$.

By continuity of $f$, there exists $r_1\in(0,0.8r_0)$, such that for $x\in B(x^*,r_1)$, $f(x)\leq f(x^*)+\lambda_0$.
When the starting point $x_0$ is in $B(x^*,r_1)$, we claim that $x_k\in B(x^*,0.8r_0)$ for every $k>0$. That is because otherwise there exists $k_0>0$ and $\theta\in[0,1)$ such that $(1-\theta)x_{k_0}+\theta x_{k_0+1}\in U_{[0.8r_0,0.9r_0]}$. Therefore $f((1-\theta)x_{k_0}+\theta x_{k_0+1})\geq f(x_0)+\lambda_0>f(x_0)+\lambda$, contradicting the fact that $\mathcal{A}_{\delta,\lambda}$ is $\lambda$-path-bounded.

In summary, we prove that for any $\varepsilon>0$, there exists $\delta_0>0$ such that for every $\delta\in(0,\delta_0)$, $\lambda\in(0,\lambda_0)$ and $x_0\in B(x^*,r_1)$, for $\mathcal{A}_{\delta,\lambda}$ starting from $x_0$, every cluster point of the generated point sequence $\{x_k\}_{k\geq 0}$ is almost surely a $\delta_0$-stationary point in $B(x^*,r_0)$, so we have $x_k\in B(x^*,\varepsilon)$ for infinitely many $k$. Hence $x^*$ is $(r_1,\lambda_0)$-convergence-stable by definition.
\end{proof}

In the sufficiency proof of \Cref{thm:max-type-is-convergence-stable}, we actually present a quantitative approach to proving a stationary point to be $(r_1,\lambda_0)$-convergence stable. We summarize it as the following proposition.
\begin{proposition}\label{prop:quan-result-of-convergence-stable}
Let $f\in\mathscr{S}$ or $f\in\mathscr{A}$. $x^*$ is the unique stationary point as well as the unique local minimum point in $B(x^*,r)$. Let
\[\lambda_0:=\frac{1}{2}\left(\min_{\norm{x-x^*}\in[0.8r,0.9r]} f(x)-f(x^*)\right),\]
and 
\[r_1:=\max\{r'\in(0,0.8r]:\text{For every } x\in B(x^*,r'), f(x)\leq f(x^*)+\lambda_0\}.\]
Then $x^*$ is $(r_1, \lambda_0)$-convergence-stable.
\end{proposition}
A direct corollary follows by \Cref{thm:unique-local-min-is-parametric-stable}, \Cref{thm:max-type-is-convergence-stable} and \Cref{thm:analytic-max-is-convergence-stable}:
\begin{corollary}
Let $f(x;R)$ be a function in $\mathscr{S}$ or $\mathscr{A}$ with domain $\Omega$ and parameter space $\mathscr{R}$. Suppose further that $f$ is equicontinuous in $R$. $x^*$ is a stationary point for $f(\cdot;R)$. If $x^*$ is convergence-stable, then $x^*$ is also parametric-stable for parameter $R$.
\end{corollary}

That is an interesting corollary. Parametric stability has nothing to do with optimization algorithms but convergence stability involves optimization algorithms heavily. The one implies the other is somewhat unexpected merely from the definitions. From a geometric perspective, however, this implication is natural: Both of them characterize a bowl-shaped graph of the function.

Using \Cref{lemma:x0-in-A*} and \Cref{lemma:constant-on-hatA} in the proof of \Cref{thm:analytic-max-is-convergence-stable}, for function $f\in\mathscr{A}$, we are able to characterize the geometric properties of the locus of all stationary points of $f$. For interested readers, see \Cref{app:proof-analytic-max-convergence-stable} and \Cref{app:proof-geometric-char-of-sp}. Formally, we have the following proposition.
\begin{proposition}[geometric characterizations on stationary points]\label{prop:geometric-char-of-sp}
Stationary points of $f\in\mathscr{A}$ form a finite union of connected subsets $A_i$ satisfying all the following conditions:
\begin{itemize}
    \item Each $A_i$ is a sub-analytic set as well as an analytic manifold (possibly a single point).
    \item Each point in $A_i$ is a stationary point.
    \item $f$ is a constant on $\bar{A}_i$, the closure of $A_i$.
    \item If $x^*$ is a non-isolated local minimum point, then there exists a neighborhood $U$ of $x^*$, such that for every $A_i$ and every $x\in A_i\cap U$, $x$ is also a non-isolated local minimum point.
\end{itemize}
\end{proposition}

At last, we mention that all results in this section hold for Dini stationary points. So this theory can be applied to more situations than gradient-based settings.
\subsection{Stability Analysis of the Tsaknakis-Spirakis Algorithm for Approximate Nash Equilibria}\label{subsec:appTSAlgo}
Nash equilibrium is the key concept in game theory. Nash proved the existence of such equilibria in normal form games \cite{nash1951non}. For computer science, a major concern is to \emph{compute} a Nash equilibrium. Papadimitriou \cite{DBLP:journals/jcss/Papadimitriou94} introduced the complexity class PPAD. Finding an approximate solution in general was proved to be PPAD-complete for 3NASH by Daskalakis, Goldberg and Papadimitriou \cite{DBLP:journals/im/TsaknakisS08} and for 2NASH by Chen, Deng and Teng \cite{DBLP:journals/jacm/ChenDT09}. It is conjectured that PPAD-hard problem cannot be solved efficiently. These results lead to attempts to find an $\epsilon$-approximate Nash equilibrium in polynomial time for some small constant $\epsilon>0$. The state-of-the-art algorithm for two-player normal form games\footnote{We refer it as TS algorithm for short.} is given by Tsaknakis and Spirakis \cite{DBLP:journals/im/TsaknakisS08} with an approximation ratio of at most $b\approx 0.3393$. Improving the approximation ratio has been unsuccessful ever since their work.

Until recently, little understanding of the TS algorithm has been given. Chen \emph{et al}. \cite{DBLP:journals/corr/abs-2107-01471} gives a complete worst-case analysis and shows that the TS algorithm has a lower bound of $0.3393$ as well. They profile all tight instances as feasible solutions of a linear program. Their experiments indicate that most tight instances are unstable:  A small perturbation near a $0.3393$ solution may help the TS algorithm find another faraway solution with a much better approximation ratio. How to understand such instability remains blank. In this application, we propose a thorough stability analysis for the TS algorithm based on the theory presented above. The analysis naturally provides some computable sufficient conditions to judge whether a tight instance is stable and shows that $0.3393$ is still a lower bound even when stability is taken into consideration, which coincides with the empirical results by \cite{DBLP:journals/corr/abs-2107-01471}. Our results may help further understand the TS algorithm and improve the approximation ratio finally.

\subsubsection{Tsaknakis-Spirakis Algorithm}
We first introduce some concepts. We use $R_{m\times n}$ and $C_{m\times n}$ to denote the payoff matrices of the row player and column player, where the row player and the column player have $m$ and $n$ strategies, respectively. Without loss of generality, assume all entries of $R$ and $C$ are normalized to $[0, 1]$. Probability vectors are used to describe players' behavior. The $i$th entry in a probability vector indicates the probability that a player plays pure strategy  $i$. In other words, row player's strategy and column player's strategy lie in $\Delta_m$ and $\Delta_n$ respectively, where
\begin{align*}
    \Delta_k &:=\{x\in\mathbb{R}^k:x_i\geq 0, \sum_i x_i = 1\}.
\end{align*}
A strategy pair $(x, y) \in \Delta_m \times \Delta_n$ is called an \emph{$\epsilon$-approximate Nash equilibrium}, if for any $x'\in\Delta_m$, $y'\in\Delta_n$, the following inequalities hold:
\begin{align*}
    \T{x'}Ry\leq \T{x}Ry+\epsilon,\\
    \T{x}Cy'\leq \T{x}Cy+\epsilon.
\end{align*}
A Nash equilibrium is an $\epsilon$-approximate Nash equilibrium with $\epsilon = 0$.

To simplify our further discussion, for any probability vector $u$, we use the following notations:
\begin{align*}
    \max(u)&:=\max\{u_i\},\\
    \max_S(u)&:=\max\{u_i:i\in S\},\\
    \supp(u)&:=\{i:u_i>0\},\\
    \suppmax(u) &:= \{i:\forall j,\ u_i\geq u_j\},\\
    \suppmin(u) &:= \{i:\forall j,\ u_i\leq u_j\},
\end{align*}

The TS algorithm tries to find the approximate Nash equilibrium problem by directly minimizing the approximation ratio itself. Specifically, we define the following functions:
\begin{align*}
    f_R(x,y) &:= \max(Ry)-\T{x}Ry, \\
    f_C(x,y) &:= \max(\T{C}x)-\T{x}Cy, \\
    f(x, y) &:= \max\left\{f_R(x, y), f_C(x, y)\right\}.
\end{align*}
The goal is to minimize $f(x,y)$ over $\Delta_m\times\Delta_n$. Clearly, $f$ is the approximation ratio itself.

The idea of the TS algorithm is to find a Dini stationary point\footnote{In this application, we will always use Dini stationary point, so we simply call it a stationary point thereafter.} of the objective function $f$ by a descent algorithm starting at some certain initial point and make a further adjustment by making a convex combination over the stationary point and corresponding best response. We only care about stationary point searching since the adjustment step is only related to the stationary point itself, not the searching process. For readers who want to refer to the full algorithm, see Appendix~\ref{app:TSAlgo}. 

\subsubsection{A Robust Tight Approximation Ratio}
We say a stationary point $(x^*,y^*)$ is \emph{tight} if the TS algorithm finds it and the output approximation ratio is $b\approx 0.3393$. Our main goal is to prove the following two theorems:

\begin{theorem}\label{thm:appTS:strict-SP-exists}
There exists a tight instance $(R,C)$ such that the corresponding tight stationary point $(x^*,y^*)$ is a strict one being both parametric-stable and convergence-stable.
\end{theorem}

\begin{theorem}\label{thm:appTS:tight-in-practice}
Let $\mathcal{F}$ be a distribution over $\Delta_m\times\Delta_n$. Given $(x_0,y_0)\sim\mathcal{F}$, let $\mathsf{out}(x_0,y_0)$ be a random variable equal to the approximation ratio of the output of the TS algorithm given initial point $(x_0,y_0)$. Let $b\approx 0.3393$ is the tight bound of the TS algorithm. Then there exists a tight instance $(R,C)$ such that for any distribution $\mathcal{F}$ with positive density function,
\[\Pr_{(x_0,y_0)\sim\mathcal{F}}\left[\mathsf{out}(x_0,y_0)=b\right]>0.\]
\end{theorem}

Theorem~\ref{thm:appTS:tight-in-practice} is a natural corollary of Theorem~\ref{thm:appTS:strict-SP-exists} and we will prove it at the end of this subsection. It shows that there exists a game that one can \emph{indeed} find a $0.3393$ solution by doing experiments. The lower bound by \cite{DBLP:journals/corr/abs-2107-01471} is thus a robust result from this angle.

To prove Theorem~\ref{thm:appTS:strict-SP-exists}, we need to show
\begin{enumerate}
    \item The descent algorithm $\mathcal{A}$ is decreasing (thus ultimately decreasing), $\delta$-result-stationary for any $\delta>0$ and $\lambda$-path-bounded with $\lambda=0$.
    \item $f\in\mathscr{S}$, i.e., $f$ is a max-type function of $C^1$ functions with a compact polyhedron domain.
    \item $f(x,y;R,C)$ is equicontinuous in $(R,C)$.
    \item For any strict stationary point $(x^*,y^*)$, there exists a radius $r>0$ such that $(x^*,y^*)$ is the unique local minimum point as well as the unique stationary point in $B((x^*,y^*),r)$.
    \item There exists a tight instance $(R,C)$ with a strict tight stationary point $(x^*,y^*)$.
\end{enumerate}

By \cite{DBLP:journals/im/TsaknakisS08}, $\mathcal{A}$ is decreasing (thus ultimately decreasing) and $\delta$-result-stationary. To prove $\mathcal{A}$ is $\lambda$-path-bounded with $\lambda=0$, we have the following lemma:
\begin{lemma}\label{lemma:descent-monotone}
Any iteratively generated point sequence $\{x_k\}_{k\geq 0}$ by the descent algorithm $\mathcal{A}$ starting from any initial point $x_0$ satisfies that for any $\theta\in[0,1]$, the 
\[f(x_k,y_k)\geq f((1-\theta)x_k+\theta x_{k+1},(1-\theta)y_k+\theta y_{k+1}).\]
\end{lemma}
Since the proof of \Cref{lemma:descent-monotone} involves many details of the descent algorithm $\mathcal{A}$, we prove the lemma in \Cref{app:proof-descent-monotone}. For readers interested in the proof, see \Cref{app:TSAlgo} for a complete descent algorithm first.

By writing every item in $\max(Ry)$ and $\max(C^Tx)$ explicitly, we see that $f(x,y)$ is a max-type function of bilinear functions, which are undoubtedly $C^1$ functions. The domain of $f$ is $\Delta_m\times\Delta_n$, which is clearly a compact polyhedron domain. Note that $\Delta_m\times\Delta_n\times[0,1]^m\times[0,1]^n$ is a compact set, thus by \Cref{lemma:equicontinuity}, $f(x,y;R,C)$ is equicontinuous in $(R,C)$.

The only two tough conditions are the existence of radius $r$ and the existence of tight stationary point. We state them as the following two propositions.
\begin{proposition}\label{prop:strict-delta-SP-distance}
If $(x^*,y^*)$ is a strict stationary point, then there exists $r_1>0$ such that for any $\varepsilon>0$, there are some $\delta_0>0$, such that for $\delta\in(0,\delta_0)$, every $\delta$-stationary point in $B((x^*,y^*),r_1)$ is in $B((x^*,y^*),\varepsilon)$. Moreover, $f(x^*,y^*)<f(x',y')$ holds for every $(x',y')\neq (x^*,y^*)$ in $B((x^*,y^*),r_1)$.
\end{proposition}

\begin{proposition}\label{prop:strict-tight-SP-exists}
There exists a tight instance $(R_0,C_0)$ with a strict tight stationary point $(x^*,y^*)$. In other words, $f(x^*,y^*;R_0,C_0)=b\approx 0.3393$ and the output of the TS algorithm reaches a approximation ratio of $b\approx 0.3393$.
\end{proposition}
The proofs of both propositions are every technical and heavily rely on work \cite{DBLP:journals/corr/abs-2107-01471} and \cite{DBLP:journals/im/TsaknakisS08}. We postpone them to \Cref{app:proof-strict-delta-SP-distance} and \Cref{app:proof-strict-tight-SP-exists}, respectively. It is worth noting that our proofs also show a detailed approach to calculating $r_1$ and finding a tight stationary point. Thus even though we can prove \Cref{thm:appTS:strict-SP-exists} using the combination of \Cref{lemma:Df>0circle}, \Cref{thm:analytic-max-is-convergence-stable} and \Cref{prop:delta-parametric-stable} rather than \Cref{prop:strict-delta-SP-distance}, we here demonstrate the possibility to use a primitive but more quantitative approach.

Now we clear the ``todo list'' and thus \Cref{thm:appTS:strict-SP-exists} is proved. Below we prove \Cref{thm:appTS:tight-in-practice} for completeness.
\begin{proof}[Proof of \Cref{thm:appTS:tight-in-practice}]
Take game $(R,C)$ and strict stationary point $(x^*,y^*)$ in \Cref{thm:appTS:strict-SP-exists}. By \Cref{prop:strict-delta-SP-distance} and \Cref{thm:max-type-is-convergence-stable}, there exists a radius $r>0$ such that $(x^*,y^*)$ is $(r,0)$-convergence stable. It suffices to show $Q(\mathcal{F})=\Pr_{(x_0,y_0)\sim\mathcal{F}}\left[(x_0,y_0)\in B((x^*,y^*),r))\right]>0$ for all $\mathcal{F}$. First, consider the probability measure $P$ induced by the Lebesgue measure on $\Delta_m\times\Delta_n$ (view as a product of simplexes with a positive measure on $\mathbb{R}^{m-1}\times\mathbb{R}^{n-1}$). Clearly, under Euclidean norm, $Q(P)>0$. Now consider $\mathcal{F}$ in the theorem. Since $\mathcal{F}$ has positive density function, $\mathrm{d}\mathcal{F}/\mathrm{d} P>0$, $P$-a.s. Then by Radon-Nikodym Theorem, 
\begin{align*}
    Q(\mathcal{F})=\int_{(x_0,y_0)\in B((x^*,y^*),r))} \mathrm{d}\mathcal{F}=\int_{(x_0,y_0)\in B((x^*,y^*),r))}\frac{\mathrm{d}\mathcal{F}}{\mathrm{d} P}\mathrm{d} P.
\end{align*}
By the property of Lebesgue measure, $Q(\mathcal{F})=0$ if and only if $\mathrm{d} \mathcal{F}/\mathrm{d} P=0$, $P$-a.s. on $B((x^*,y^*),r)$. However, this contradicts the assumption. Thus $Q(\mathcal{F})>0$.
\end{proof}
\section{Parametric Stability and Convergence Stability}\label{sec:notion-of-stability}
In this section, we introduce two different kinds of stability involved in optimization problems: \emph{parametric stability} and \emph{convergence stability}. They arise due to numerical precision considerations but from different perspectives.

\subsection{Parametric Stability}
To introduce parametric stability, we focus on parameterized function $f(x;r)\in\mathscr{F}$. In computer, numerical calculations for problems about real numbers are mostly inaccurate, implying that the input parameter $r$ could be slightly perturbed to some other parameter $r'$ with $\norm{r-r'}<\epsilon$. Such a practical consideration elicits a natural question: Does such tiny inaccuracy lead to drastic changes to the stationary point $x^*$? In some situations, which point is found is non-negligible. For instance, in machine learning, if we view finding an optimal model as an optimization problem, different solutions with the same loss could have very different test accuracies. Such observation motivates the following formal definition of \emph{parametric stability}:

\begin{definition}
Let $f(x;r)$ be a function in $\mathscr{F}$, where $x\in \Omega$ and $r\in\mathscr{R}$. A stationary point $x^*\in \Omega$ is $\delta(\epsilon)$-\emph{parametric-stable} for parameter $r_0$ if for all $\epsilon>0$, there exists $\delta(\epsilon)>0$ such that for every $r'\in\mathscr{R}$ satisfying $\norm{r'-r_0}<\delta$ there exists a stationary point $x'$ for parameter $r'$ such that
\[\norm{x'-x^*}<\epsilon\text{ and } |f(x^*;r_0)-f(x';r')|<\epsilon.\]
If we simply say $x^*\in \Omega$ is \emph{parametric-stable}, it means that such $\delta$ exists but has no explicit form.
\end{definition}
In other words, a parametric-stable stationary point is less sensitive to small perturbations in the sense of its position and $f$ value. Let us see several examples of one-dimensional functions about parametric stability.

\begin{example}\label{ex:parametric-stable}
Let $f(x;\lambda)=|x|^3-\lambda x, x\in[-1,1]$, $\lambda\in[-1,1]$. The unique stationary point of $f$ is
\[x^*(\lambda)=\begin{cases}\sqrt{\lambda/3},&0<\lambda\leq 1,\\
0,&\lambda=0,\\
-\sqrt{-\lambda/3},&-1\leq \lambda<0.
\end{cases}\]
The corresponding $f$ value is $f(x^*(\lambda))=-\frac{2}{3\sqrt{3}}|\lambda|^{3/2}$. By definition, $x^*(0)=0$ is $\sqrt{\epsilon/3}$-parametric-stable for $\lambda=0$. See \Cref{fig:parametric-stable-ex} for a geometric demonstration.

\begin{figure}[ht]
    \centering
    \input{Figures/parametric-stable-ex}
    \caption{Demonstration of \protect{\Cref{ex:parametric-stable}}. The location and the $f$ value of the stationary point change mildly as $\lambda$ changes.}
    \label{fig:parametric-stable-ex}
\end{figure}

Indeed, we will see later that $x^*$ is parametric-stable for all $\lambda$. 
\end{example}

A saddle point is usually not parametric-stable:
\begin{example}\label{ex:saddle-point-is-not-parametric-stable}
Let $f(x;\lambda)=x^3-\lambda x$, $x\in[-1,1]$, $\lambda\in[-1,1]$. When $\lambda=0$, $x^*=0$ is a stationary point (saddle point). Any slight perturbation on $\lambda$ towards the negative direction makes $f$ a strictly increasing function, and then $x^*$ slides to $-1$. Thus $x^*=0$ is not parametric-stable for $\lambda=0$. See \Cref{fig:saddle-point-not-parametric-stable} for a geometric demonstration.

\begin{figure}[ht]
    \centering
    \input{Figures/saddle-point-not-parametric-stable}
    \caption{Demonstration of \protect{\Cref{ex:saddle-point-is-not-parametric-stable}}. Stationary point $x^*$ slides to $-1$ no matter how slightly $\lambda$ changes.}
    \label{fig:saddle-point-not-parametric-stable}
\end{figure}
\end{example}

Even a local minimum point may not be parametric-stable:
\begin{example}\label{ex:local-min-point-is-not-parametric-stable}
Let $f(x;\lambda)=|x+1|+|x-1|-\lambda x$, $x\in[-2,2]$, $\lambda\in[-1,1]$.  When $\lambda=0$, $x^*=0$ is a stationary point (local minimum point). Any slight perturbation on $\lambda$ towards the negative direction makes $x^*$ slides to $-1$; any slight perturbation on $\lambda$ towards the positive direction makes $x^*$ slides to $1$. Thus $x^*=0$ is not parametric-stable for $\lambda=0$. See \Cref{fig:local-min-point-not-parametric-stable} for a geometric demonstration.

\begin{figure}[ht]
    \centering
    \input{Figures/local-min-point-not-parametric-stable}
    \caption{Demonstration of \protect{\Cref{ex:local-min-point-is-not-parametric-stable}}. Stationary point $x^*$ slides to $-1$ no matter how slightly $\lambda$ changes.}
    \label{fig:local-min-point-not-parametric-stable}
\end{figure}
\end{example}

\subsection{Convergence Stability}
While the parametric stability is intuitive, the other concept, convergence stability, is much subtler. We consider a specific optimization algorithm $\mathcal{A}$ for function $f$. $\mathcal{A}$ will choose some initial point $x_0$ and do many steps of an optimization procedure. Each step yields a point $x_k$. When certain conditions are satisfied, $\mathcal{A}$ will terminate; otherwise, $\mathcal{A}$ never stops. Although a perfect algorithm should find an exact stationary point, it is impossible due to limited precision and time. Thus, most algorithms aim to find an approximate stationary point. Precisely, given a precision $\delta$ as the input parameter, a reasonable algorithm should
\begin{itemize}
    \item either terminate at a $\delta$-stationary point,
    \item or generate an infinite series $\{x_k\}$ with $\delta$-stationary points as its cluster\footnote{A cluster point $x'$ of set $\{x_k\}$ is defined as for any $r>0$, for infinitely many $k$, $x_k\in B(x',r)\setminus\{x'\}$.} points.
\end{itemize}

Our goal is to define a stability concept about approximate stationary points. Naturally, stability means a $\delta$-stationary point should be close to an exact stationary point when $\delta$ is close to $0$; otherwise, a slight difference on $\delta$ may lead to a significant change of the found solution. Such changes are ubiquitous, especially in nonconvex optimization problems. Moreover, the changes may result in a worse $f$ value or other worse properties (recall the example of test accuracies in the part of parametric stability). 

We consider the behavior of the optimization algorithm $\mathcal{A}$ when the precision $\delta$ is positive but sufficiently small. Intuitively, a stable stationary point $x^*$ should be accompanied with $\delta$-stationary points that converge to $x^*$ as $\delta\to 0$. Since there could be many exact stationary points scattered in the domain, approximate stationary points could also be scattered. To focus on the nearby approximate stationary points, we only inspect a neighborhood of the exact stationary point $x^*$. When the optimization algorithm $\mathcal{A}$ starts from a point $x_0$ near the exact stationary point $x^*$, the following two statements both describe the stability of a stationary point $x^*$:
\begin{itemize}
    \item For any sufficiently small precision parameter $\delta$, $\mathcal{A}$ will find a $\delta$-stationary point \emph{near} $x^*$  regardless of the choices of the initial point in the neighborhood of $x^*$.
    \item any $\delta$-stationary point in the neighborhood gets close to $x^*$ when $\delta\to 0$.
\end{itemize}

This observation motivates the formal definition of convergence stability. Specifically, for algorithm $\mathcal{A}$, we denote the starting point by $x_0$, and the point after the $k$-th iteration by $x_k$. For convenience, we assume that if $\mathcal{A}$ finds a $\delta$-stationary point $x_{k_0}$, then for all $k>k_0$, it sets $x_k=x_{k_0}$, so that $\{x_k\}$ is always an infinite sequence. A \emph{limit} point $x'$ of $\{x_k\}$ is then defined as: for any $r>0$, for infinitely many $k$, $x_k\in B(x',r)$.\(\)

We need to restrict types of optimization algorithms to obtain nontrivial results. In many practical scenarios\footnote{See \Cref{subsec:appGeneralOptimization} for a further discussion.}, we observe the following three properties:

\begin{enumerate}
    \item The objective function $f(x)$ decreases after each iteration, i.e., for all $k=0,1,2,\cdots$, 
    \[f(x_{k+1})\leq f(x_{k}).\]
    \item Every limit point of $\{x_k\}$ is a $\delta$-stationary point, where $\delta$ can be controlled by $\mathcal A$'s input parameters.
    \item The step size $\norm{x_k-x_{k+1}}$ is limited, $i.e.$, there is $\lambda>0$ such that for all $k=0,1,2,\cdots$, $\norm{x_k-x_{k+1}}\leq\lambda$. If we assume $f(x)$ is $L$-Lipschitz, then this implies that for any $\theta\in[0,1]$,
    \[f(\theta x_{k+1}+(1-\theta)x_{k})\leq f(x_{k})+L\lambda.\]
\end{enumerate}

Our analysis needs $\{x_k\}$ to have the following more relaxed properties.

\begin{definition}\label{def:commonalgorithmproperties}
We say $\{x_k\}_{k\geq 0}$ is 
\begin{itemize}
    \item \emph{ultimately decreasing}, if every limit point $\bar{x}$ of $\{x_k\}$ satisfies $f(\bar{x})\leq f(x_0)$.
    \item $\delta$-\emph{result-stationary}, if every limit point of $\{x_k\}$ is $\delta$-stationary point.
    \item $\lambda$-\emph{path-bounded}, if for any $ k=0,1,\cdots$ and any $\theta\in[0,1]$, $f(\theta x_{k+1}+(1-\theta)x_{k})\leq f(x_0)+\lambda$.
\end{itemize}
\end{definition}


We say an algorithm $\mathcal A$ has these properties, if the point sequence $\{x_k\}$ generated by $\mathcal A$ is guaranteed to always (or almost surely\footnote{Throughout this paper below, the term "almost surely" (with respect to a probability measure) only accounts for the randomness inside the algorithm. Choices of initial points are not taken into consideration.}) satisfy these properties, regardless of the starting point $x_0$. The ultimate decrease property is very reasonable as long as an optimization algorithm does try to find a local minimum point. The $\delta$-result-stationary property is usually theoretically proven for many effective optimization algorithms. The $\lambda$-path-bounded property is usually satisfied in practice\footnote{For instance, we consider the descent optimization process on a continuous function $f$. By setting the step size sufficiently small, we will have path-bounded property. See \Cref{subsec:appGeneralOptimization} for detail discussions.}, and can be easily guaranteed for the descent optimization process on an $L$-Lipschitz function $f$ by limiting the step size $\lVert x_k-x_{k+1}\rVert< \lambda/L$.

Now we present the definition of convergence stability. We use notation $A_{\delta, \lambda}$ to denote an optimization algorithm satisfying ultimate decreasing, $\delta$-result-stationary and $\lambda$-path-bounded. 

\begin{definition}\label{def:converge-stab}
A stationary point $x^*\in \Omega$ is called $(r_0,\lambda)$-\emph{convergence-stable} if for every $\varepsilon>0$, there exists $\delta>0$,
such that for any $\delta_1\in(0,\delta),\lambda_1\in(0,\lambda)$ and any initial point $x_0\in B(x^*,r_0)$, the point sequence $\{x_k\}_{k\geq 0}$ generated by an \emph{arbitrary} algorithm $\mathcal A_{\delta_1,\lambda_1}$ almost surely satisfies
\[\limsup_{k\to\infty}\norm{x_k-x^*}< \varepsilon.\]
If we simply say $x^*\in \Omega$ is \emph{convergence-stable}, it means such $r_0$ and $\lambda$ exist but it has no explicit form.
\end{definition}

The definition of convergence stability is less intuitive, so we give a few simple examples to demonstrate it.

\begin{example}
Let $f(x)=x^2, x\in[-1,1]$. We want to show that stationary point $x^*=0$ is $(r_0,\lambda)$-convergence-stable for some $r_0$ and $\lambda$. $Df(x,x')=2x(x'-x)$. $Df(x)=\min_{x'} Df(x,x')=-2|x|(1+|x|)$. Thus the only $\delta$-stationary points are
\[\Tilde{x}_\pm(\delta)=\pm\frac{1}{2}\left(\sqrt{2\delta+1}-1\right).\]
By definition, any algorithm $\mathcal{A}_{\delta_1,\lambda_1}$ with $\delta_1\in(0,\delta)$, $\lambda_1>0$ and $x_0\in[-1,1]$ can find $\Tilde{x}_+(\delta)$ and $\Tilde{x}_-(\delta)$ as all possible limit points. As $\delta\to 0$, $\norm{\Tilde{x}_\pm(\delta)}\to 0$. Thus for any $r_0\in(0,1]$ and $\lambda>0$, $x^*=0$ is convergence-stable.
\end{example}

We also show a nonexample:
\begin{example}\label{ex:nonex-for-convergence-stable}
Let $f(x)=|x+1|+|x-1|, x\in[-2,2]$. We prove that $x^*=0$ is not convergence-stable. It suffices to note that any $\delta$-stationary point with $\delta<2$ lies in $(-1,1)$. Thus for any $r>0$ and $\lambda>0$, set $\epsilon=\min\{1,r/2\}$. An initial point $x_0=0.9r$ will make some algorithms get stuck and terminate. Thus $x^*$ fails the definition of convergence stability.
\end{example}

\Cref{ex:nonex-for-convergence-stable} implies that convergence stability is a property mainly determined by the stationary point $x^*$ itself. In other words, it makes the weakest assumptions on optimization algorithms and guarantees for \emph{any} ultimately decreasing, $\delta$-result-stationary and $\lambda$-path-bounded algorithm, the optimization solution cannot escape from a neighborhood of $x^*$ and will always find a solution near $x^*$ once it goes into that neighborhood.
\section{Preliminaries}\label{sec:preliminaries}
\subsection{Basic Notations and Definitions}\label{subsec:notations-and-defs}
\begin{itemize}
    \item For two sets $A$ and $B$, the intersection is $A\cap B$; the union is $A\cup B$; the difference is $A\setminus B=\{x\in A:x\notin B\}$.
    \item The closure of a set $A$ is denoted by $\bar{A}$.
    \item For a positive integer $n$, define $[n]:=\{1,2,\ldots, n\}$.
    \item Denote by $1_k$ a $k$-dimension column vector with all entries equal to $1$. Denote by $0_k$ a $k$-dimension column vector with all entries equal to $0_k$.
    \item Denote the real number set by $\mathbb{R}$ and $n$-dimensional Euclidean space by $\mathbb{R}^n$. 
    \item The Euclidean norm in $\mathbb{R}^n$ is $\norm{\cdot}$.
    \item The open ball centered at $x$ with radius $r$,  $\{x'\in\mathbb{R}^n:\norm{x-x'}<r\}$, is denoted by $B(x,r)$. A closed ball can be expressed by $\bar{B}(x,r)$.
    \item An open area on $\mathbb{R}^n$ is a connected open set. A closed area on $\mathbb{R}^n$ is the closure of an open area, which is a connected closed set. Not every open set is an open area and not every closed set is a closed area. We only consider open or closed area, so we simply call a set an area if it is an open or closed area. The definition of area is only used for function domains.
    \item Let $\Omega\subseteq \mathbb{R}^n$. $s$ has a valid direction at $x\in \Omega$ if there exists $x'\in \Omega\setminus\{x\}$ such that $s=\theta_0(x'-x)$ for some $\theta_0>0$ and $x+\theta(x'-x)\in \Omega$ for all $\theta\in(0,1)$. In other words, although $s$ can be arbitrarily long, its direction is the same as one parallel segment
    starting at $x$ and totally contained in $\Omega$.
\end{itemize}
Let $F$ be a mapping. If $F$ maps to $\mathbb{R}$, we call $F$ a function. Below we suppose $F:\Omega\to\mathbb{R}^n$ is a mapping. The following notations and definitions apply directly to functions as well.
\begin{itemize}
    \item The image of $S\subseteq \Omega$ is $F(S):=\{F(x):x\in S\}$.
    \item The pre-image of $S\subseteq\mathbb{R}$ is $F^{-1}(S):=\{x\in \Omega:F(x)\in S\}$.
    \item We say $F$ is continuous on set $A\subseteq \Omega$ if for every open set $U\subseteq \mathbb{R}^n$, $F^{-1}(U)\cap A$ is an open set in $A$ (under the subspace topology of $A$).
    \item If $F$ is univariate, denote its derivative at $x$ by $F'(x)$ or $\mathrm{d} F(x)/\mathrm{d} x$.
    \item Suppose $F$ is multi-variate. Its directional derivative in direction $s$ at $x$ is $\partial F(x)/\partial s$. Its total differential at $x$ is $\mathrm{d} F(x)$. Its gradient at $x$ is $\mathrm{grad} F(x)$.
    \item We say $F$ is $L$-Lipschitz if for all $x$ and $y$ in the domain, $\norm{F(x)-F(y)}\leq L\norm{x-y}$. 
    \item $F\in C^n(\Omega)$ if it is continuously differentiable of the $n$th order on $\Omega$. $F\in C^\infty(\Omega)$ if $F\in C^n(\Omega)$ for all positive integer $n$. If $F\in C^\infty(\Omega)$, we call $F$ a smooth mapping.
    \item We say $F$ is analytic if  for every $x\in \Omega$, $F$ is equal to its Taylor expansion at $x$ in some neighborhood of $x$. If $F$ is analytic, then it is smooth. However, the converse is not true in general.
\end{itemize}
For less frequently used definitions and notations, we will define them when using them.

\subsection{Function class $\mathscr{F}$}
We consider the optimization problems of minimizing function $f: \Omega\to\mathbb{R}$, where $\Omega$ is an area on $\mathbb{R}^n$. It is worth noting that $f$ can be nonconvex and nonsmooth. We require that at every point in $\Omega$, \emph{valid unit} direction vector $s$ form a closed set. This requirement restricts the possible domain of $f$. For example, $\Omega$ can be an open ball or a closed polyhedron but not a closed ball\footnote{The surface of a closed ball has tangent vectors which cannot be described by a segment contained in the ball.}. We make such an assumption on function domains so that the minimum value of directional derivatives can be attained. Furthermore, $f$ should be continuous and have directional derivatives in every valid direction at every point in $\Omega$. We also demand that the directional derivative $\partial f(x)/\partial s$, as a function of unit direction vector $s$, be continuous in $s$.\footnote{Note that the domain of $s$ could be disconnected. For instance, directions in $\mathbb{R}$ are $\{-1,1\}$, which is a discrete set. Our continuity concept involves a topological meaning. See \Cref{subsec:notations-and-defs} for details. For the case of $\{-1,1\}$, its topology is discrete. Thus any function defined on it is continuous.} Function $\partial f(x)/\partial s$ has a minimum in $s$, denoted by $\partial f(x)$. Nevertheless, $\partial f(x)/\partial s$ needs not to be continuous in $x$. Denote the set of functions satisfying these properties by $\mathscr{F}$. Let us examine some examples of $\mathscr{F}$. In the following examples, unless otherwise specified, the domains of functions are an open convex area $\Omega\subseteq\mathbb{R}^n$, and hence all directions in $\mathbb{R}^n$ are valid. Also, we assume direction vector $s$ is a unit vector in the examples.
\begin{example}\label{ex:differentiable}
Differentiable function $f$. All smooth functions are included. It is by definition that $f$ is directional differentiable on every possible direction at every point in $\Omega$. $\partial f(x)/\partial s=\mathrm{grad} f(x)\cdot s$ is linear in $s$, so is $\partial f(x)/\partial s$ is continuous in $s$.
\end{example}

\begin{example}\label{ex:max-type}
Max-type functions of differentiable functions. Specifically, $f=\max_{i=1}^s\{f_i\}$, where $f_i$ is differentiable on $\Omega$. Below we show that $f\in\mathscr{F}$.

Let $x\in\Omega$. We first prove that $\partial f(x)/\partial s$ exists. Let $K$ be the index set so that $f(x)=f_i(x)$ for $i\in K$. By the definition of differentiablity, $f_i(x+\alpha s)=f_i(x)+\alpha\partial f_i(x)/\partial s+\beta_i(\alpha)$, where $\beta_i(\alpha)/\alpha\to 0$ as $\alpha\to 0+$. For sufficiently small $\alpha$, by continuity of $f_i$,
\begin{align*}
    f(x+\alpha s)&=\max_{i\in K}\left\{f_i(x)+\alpha\partial f_i(x)/\partial s+\beta_i(\alpha)\right\}\\
    &=f(x)+\alpha\max_{i\in K}\{\partial f_i(x)/\partial s\}+\gamma(\alpha),
\end{align*}
where $\gamma(\alpha)=\max_{i\in K}\{\beta_i(\alpha)\}$. Note that $\gamma(\alpha)/\alpha\to 0$ as $\alpha\to 0+$. Thus we have $\partial f(x)/\partial s$ exists and is equal to $\max_{i\in K} \partial f_i(x)/\partial s$.

To show that $\partial f(x)/\partial s$ is continuous in $s$, note that by \Cref{ex:differentiable}, $\partial f_i(x)/\partial s$ is continuous in $s$. The result follows by noting that the maximum function preserves continuity.
\end{example}

As \Cref{ex:max-type} shows, there are nonsmooth functions in $\mathscr{F}$. We extend the function types in $\mathscr{F}$ by the following example:
\begin{example}\label{ex:composite}
Let $f_i$'s be functions in \Cref{ex:max-type}. Let $g$ be a differentiable function on $\mathbb{R}$. We show that $h(x)=g\left(\max_{i=1}^s\{f_i(x)\}\right)=g(f(x))\in\mathscr{F}$. Let $H(\alpha)=h(x+\alpha s)$. We first prove that $H$ is differentiable. By \Cref{ex:max-type}, $f(x)=\max_i \{f_i(x)\}$ is directional differentiable. Thus $f(x+\alpha s)$ is differentiable about $\alpha$. The proof follows from the composition law of the differentials. $\partial h(x)/\partial s$ can be calculated as:
\begin{align*}
    \partial h(x)/\partial s&=H'(0)\\
    &=g'(f(x))\cdot\left.\frac{\mathrm{d} f(x+\alpha s)}{\mathrm{d}\alpha}\right|_{\alpha=0}\\
    &=g'(f(x))\partial f(x)/\partial s.
\end{align*}
By \Cref{ex:max-type}, $\partial f(x)/\partial s$ is continuous in $s$. Thus $\partial h(x)/\partial s$ is continuous in $s$.
\end{example}

Note that proofs in \Cref{ex:max-type} and \Cref{ex:composite} apply on cases where $f_i$ and $g$ simply belong to $\mathscr{F}$ (no need to be differentiable), so we have the proposition below:
\begin{proposition}\label{prop:F-is-closed}
$\mathscr{F}$ is closed under operations of arithmetic, composition and maximum function. In other words, 
\begin{itemize}
    \item If $f,g\in\mathscr{F}$ and they have the same domain, then $f\pm g\in\mathscr{F}$ and $fg\in\mathscr{F}$. If $g\neq 0$, $f/g\in\mathscr{F}$.
    \item If $f_i(x)\in\mathscr{F}$, $i=1,2,\ldots, s$ and they have the same domain, then $\max_{i=1}^s\{f_i(x)\}\in\mathscr{F}$.
    \item If $f:\mathbb{R}\to\mathbb{R}$ is in $\mathscr{F}$ and $g\in\mathscr{F}$, then $f(g(x))\in\mathscr{F}$.
\end{itemize}
\end{proposition}
\begin{proof}
The second and third statements can be proven by the same arguments in \Cref{ex:max-type} and \Cref{ex:composite}, respectively. We only prove the first statement.

By \Cref{lemma:two-derivative} and basic properties of differentials,
\begin{align*}
\frac{\partial(f\pm g)}{\partial s}&=\frac{\partial f}{\partial s}\pm \frac{\partial g}{\partial s},\\ \frac{\partial (fg)}{\partial s}&=g(x)\frac{\partial f}{\partial s}+f(x)\frac{\partial g}{\partial s},\text{ and}\\ \frac{\partial (f/g)}{\partial s}&=\frac{1}{g(x)}\frac{\partial f}{\partial s}-\frac{f(x)}{(g(x))^2}\frac{\partial g}{\partial s}.
\end{align*}
The result follows directly by that $\partial f/\partial s$ and $\partial g/\partial s$ are both continuous in $s$.
\end{proof}
These examples and \Cref{prop:F-is-closed} suggest that $\mathscr{F}$ covers a majority of object functions discussed in optimization, machine learning and other application fields. We will consider subsets of $\mathscr{F}$ (i.e., put more restrictions on function types)  sharing enough good properties and derive nontrivial conclusions for them. Particularly, we are interested in:
\begin{itemize}
    \item Functions with parameters: $f$ can be written as $f(x;r)$ for some parameter $r$.
    \item Functions satisfying Lipschitz condition: there exists some constant $M>0$ such that $|f(x)-f(y)|\leq M\norm{x-y}$.
    \item Max-type functions of $C^1$ functions: $f=\max_{i=1}^s\{f_i\}$, where $f_i$'s have continuous gradients.
    \item Finite compositions of analytic mappings and the maximum function.
\end{itemize}

It sometimes simplifies notations to consider \emph{Dini directional differentials} \cite{Demyanov2010}:
\begin{definition}
Suppose $x,x'\in \Omega$. The \emph{Dini directional differential} in direction $x'-x$ at point $x$, denoted by $Df(x,x')$, is defined as the following limit (if the limit exists)
\[\lim_{\alpha\to 0+}\frac{f(x+\alpha (x'-x))-f(x)}{\alpha}.\]
\end{definition}
Note that the usual derivatives should be defined as the following limit:
\[\frac{\partial f(x)}{\partial s}:=\lim_{\alpha\to 0+}\frac{f(x+\alpha s)-f(x)}{\norm{\alpha s}}.\]
The relation between the two limits is direct:
\begin{lemma}\label{lemma:two-derivative}
Let $s=x'-x$. Then $\frac{\partial f(x)}{\partial s}$ exists if and only if $Df(x,x')$ exists. Moreover,
\[Df(x,x')=\frac{\partial f(x)}{\partial s}\norm{s}.\]
\end{lemma}
Similar to $\partial f$, we define $Df(x):=\inf_{x'\in \Omega\setminus\{x\}}Df(x,x')$ for Dini directional derivatives.

\subsection{Terminologies in Optimizations}
Now we involve some concepts in optimization theory. From the optimization perspective, an optimization algorithm is designed to search for a minimum point (or locally minimum point) of an objective function. In most settings, the optimization procedure is equivalent to finding a \emph{stationary point (or critical point)}, which can be defined as follows for functions in $\mathscr{F}$:
\begin{definition}
Let $f\in\mathscr{F}$ with domain $\Omega$. $x\in \Omega$ is a \emph{stationary point (or critical point)} if $\partial f(x)/\partial s\geq 0$ for every valid direction $s$, or equivalently, $\partial f(x)\geq 0$.
\end{definition}
Note that a stationary point could be a saddle point or a local minimum point. To ensure a local minimum point, we have the following definition:
\begin{definition}
Let $f\in\mathscr{F}$ with domain $\Omega$. $x\in \Omega$ is a \emph{strict stationary point} if $\partial f(x)/\partial s>0$ for every valid direction $s$, or equivalently, $\partial f(x)>0$.
\end{definition}

Due to precision restrictions, it is almost impossible for most optimization algorithms to find an exact stationary point. Thus \emph{approximate stationary points} are proposed:
\begin{definition}
Let $f\in\mathscr{F}$ with domain $\Omega$. Let $\delta>0$. $x\in \Omega$ is a $\delta$-\emph{stationary point} if $\partial f(x)/\partial s\geq -\delta$ for every valid direction $s$ or equivalently, $\partial f(x)\geq -\delta$.
\end{definition}
We can similarly define different concepts about \emph{Dini stationary points} by replacing $\partial f(x)/\partial s$ with $Df(x,x')$ in the definitions above.
\subsection{Connections to Optimization of Deep Neural Networks}\label{sect:connection-NN}
Deep neural networks play a fundamental role in deep learning. Ever since the triumph of AlexNet \cite{DBLP:conf/nips/KrizhevskySH12} on the classification problem of ImageNet \cite{DBLP:conf/cvpr/DengDSLL009}, empirical studies of deep learning thrive and have great impacts on traditional optimization, computer vision, natural language processing and almost every study field \cite{Alom2018TheHB}. There are also many theoretical results on various topics trying to reveal the ``black magic'' behind deep neural networks. A well-developed example is NTK theory \cite{DBLP:conf/nips/JacotHG18}, which is applied to prove several global convergence results \cite{DBLP:conf/icml/Allen-ZhuLS19,DBLP:conf/icml/DuLL0Z19,DBLP:conf/iclr/DuZPS19}.

In this part, we show the connections between our theory and optimization of neural networks. It turns out that our theory is tightly linked to the last steps in the training process when the learning rate (or step size in optimization jargon) is sufficiently small and the traning loss is nearly stationary. Also, we characterize the geometric properties of the set of stationary points. Such results may give a better understanding of difficulties in fine-tuning and initialization.

A neural network can be described as a prediction function or a score function
\[f:\mathcal{W}\times\mathcal{X}\to\mathbb{R}^d,\]
where $\mathcal{W}\subseteq\mathbb{R}^p$ is the parameter space (or weight space) and $\mathcal{X}\subseteq\mathbb{R}^m$ is the input space. For regression problems, $d=1$ and $f$ is a prediction function. For classification problems, $d$ is the total number of classes and $f$ is a score function.

For almost all neural network structures, including multi-layer perceptrons (MLP), convolutional neural networks (CNN),  ResNet, recurrent neural networks (RNN), and transformers, $f$ is commonly composed of several modules, such as matrix multiplication, local convolution, normalization (such as BatchNorm, LayerNorm and WeightNorm), , pooling operation (such as max-pooling and average-pooling), residual connection, self-attention, and activation functions (such as ReLU, sigmoid and GeLU) \cite{DLSurveyACMComputSurv,attentionisallyouneed}.

In the training stage, we are given $N$ training inputs $\{x_i\}_{i=1}^N$ and the corresponding labels (or values) $\{y_i\}_{i=1}^N$. If we are doing classification tasks, assume $y_i\in[n]$, our goal is to minimize the loss function 
\[L(W)=\mathcal{L}(W,f, \{x_i\},\{y_i\})+R(W),\]
where $\mathcal{L}(W,f, \{x_i\},\{y_i\})$ is the loss of the inputs, and $R(W)$ is the regularization term such as $\lambda\lVert W\rVert_2^2$ or $\lambda\lVert W\rVert_1$, or simply $0$.

The most popular loss function is the average loss of all inputs, i.e., 
\[\mathcal{L}(W,f, \{x_i\},\{y_i\})=\frac{1}{N}\sum_{i=1}^N L_i.\]

For simplicity, denote $f(x_i)$ by $s$, which is a $d$-dimensional column vector. The most common forms of $L_i$ are
\begin{itemize}
    \item quadratic loss (square loss) for $n=1$: $L_i=(s-y_i)^2/2$,
    \item SVM loss (Hinge loss) for $n>1$: $L_i=\sum_{j\neq y_i}\max\left\{0,s_j-s_{y_i}+1\right\}$, and 
    \item softmax loss (cross-entropy loss) for $n>1$: $L_i=-\log\left(e^{s_{y_i}}/\sum_j e^{s_j}\right)$. 
\end{itemize}


Even though $f$ and loss function vary, the modules and functions involved in a neural network are composed of the most basic operators, such as linear operators, maximum operators and exponential operators. Most of these operators can be expressed as finite compositions of analytic mappings and the maximum function. However, operators like square root, multiplicative inverse and logarithm may lead to numerical instability or infinity in gradient calculations. The common approaches to avoiding such numerical problems are to restrict the domain of inputs of the operator\footnote{It is usually accomplished by previous operators producing these inputs. For example, $\log x$ in the cross-entropy loss never encounters $x=0$, since $x=e^{s_{y_i}}/\sum_j e^{s_j}>0$.} or replace the operators with numerically stable substitutes\footnote{For instance, $\sqrt{x}$ is often replaced with $\sqrt{x+\epsilon}$, where $\epsilon>0$ is a small constant.}. Taking these factors into consideration, we have the following observation:

\begin{observation}
We can assume that the parameter space $\mathcal{W}$ is restricted to a compact area, and that $L(W)$ belongs to function class $\mathscr{A}$.
\end{observation}

That is exactly what \Cref{thm:analytic-max-is-convergence-stable} needs.

In the training stage, we use an optimizer to iteratively generate a sequence of parameters $\{W_i\}_{i\geq 1}$. In each iteration, the optimizer finds a direction to move and steps a specific size (called \emph{learning rate}) towards that direction. When the procedure nearly meets its end, i.e., loss $L(W)$ only has minor changes and the learning rate is sufficiently small, our theory arises. That leads to our second observation:

\begin{observation}
If we can guarantee that an optimizer $\mathcal{A}$ converges to an approximate stationary point $\Tilde{W}$ and the loss is $\lambda$-path-bounded with small $\lambda>0$ in a small neighbourhood of $\Tilde{W}$, then $\mathcal{A}$ satisfies the requirement of optimization algorithms in the definition of convergence stability.
\end{observation}

In this application, we consider the most common optimizer in deep learning, \emph{gradient descent (GD)}. So does it satisfy these conditions? In 2019, Du \emph{et al}. \cite{DBLP:conf/icml/DuLL0Z19} showed that with high probability, GD could find a global minimum of a smooth deep residual convolutional neural network with a Gaussian initialization when the network is over-parameterized (or wide enough). In the same year, Zhu \emph{et al}. \cite{DBLP:conf/icml/Allen-ZhuLS19} proved that with high probability, GD could find a global minimum of a deep ResNet or CNN with a Gaussian initialization and ReLU activation functions when the network is over-parameterized. Thus in a reasonable sense, the GD optimizer does converge and find an approximate stationary point. Also, the convergence result implies ultimately decreasing. How about $\lambda$-path-bounded? Since GD guarantees a descent loss, it is trivially true for GD as long as the learning rate is sufficiently small.

Now, combining these two observations, we can interpret \Cref{thm:analytic-max-is-convergence-stable} as the following thoerm:
\begin{theorem}\label{thm:appNN:iff-convergence-stable}
Let $f:\mathcal{W}\times\mathcal{X}\to\mathbb{R}^d$ be a neural network and $L(W)$ be the loss function. Suppose further that $L\in\mathscr{A}$, i.e., $L$ can be expressed in finite compositions of analytic mappings and the maximum function on compact area $\mathcal{W}$. Suppose that with appropriate parameters, optimizer $\mathcal{A}_{\delta,\lambda}$ is ultimately decreasing, $\delta$-result-stationary and $\lambda$-path-bounded for any $\delta,\lambda>0$. Then stationary point $W^*$ is convergence-stable about $L$ and $\mathcal{A}_{\delta,\lambda}$ if and only if $W^*$ is the unique local minimum point in some neighborhood of $W^*$.
\end{theorem}

Cooper \cite{DBLP:journals/corr/abs-1804-10200} proved the following theorem showing the landscape of the loss of overparameterized smooth neural networks.
\begin{theorem}[\cite{DBLP:journals/corr/abs-1804-10200}, Theorem 2.1]
Let $f:\mathcal{W}\times\mathcal{X}\to\mathbb{R}$ be a smooth neural network and $L(W)$ be the quadratic loss function with no regularization term. Suppose that the number of parameters of the neural network is greater than the number of data points that it is training on, i.e., $\dim\mathcal{W}>N$. The set $M=L^{-1}(0)$ is generically (that is, possibly after an arbitrarily small change to the data set) a smooth $(\dim\mathcal{W}-N)$ dimensional submanifold (possibly empty) of $\mathbb{R}^p$.
\end{theorem}
By combining these results, we get the following negative result, indicating that even a well-trained neural network can be not stable.
\begin{corollary}\label{cor:appNN:neither-nor-stable}
Let $f:\mathcal{W}\times\mathcal{X}\to\mathbb{R}^d$ be a smooth neural network and $L(W)$ be the quadratic loss function with no regularization term. Suppose that the number of parameters of the neural network is greater than the number of data points that it is training on, i.e., $\dim\mathcal{W}>N$. Suppose further that with appropriate parameters, optimizer $\mathcal{A}_{\delta,\lambda}$ is ultimately decreasing, $\delta$-result-stationary and $\lambda$-path-bounded for any $\delta,\lambda>0$. Then all zero-training loss points of $L$ are not convergence-stable about $L$ and $\mathcal{A}_{\delta,\lambda}$.
\end{corollary}

Any procedure of training a neural network involves initialization. \Cref{cor:appNN:neither-nor-stable} implies that even the initialization is close to a global minimum point, there is no guarantee to find the same global minimum point by the gradient descent method.

In \Cref{sec:prop-of-stability}, we have already characterized the locus of all stationary points in a deep neural network: They form finitely many connected subsets $A_i$ satisfying all the following conditions:
\begin{itemize}
    \item Each $A_i$ is a sub-analytic set as well as an analytic manifold (possibly a single point).
    \item Each point in $A_i$ is a stationary point.
    \item $f$ is a constant on $\bar{A}_i$, the closure of $A_i$.
    \item If $x^*$ is a non-isolated local minimum point, then there exists a neighborhood $U$ of $x^*$, such that for every $A_i$ and every $x\in A_i\cap U$, $x$ is a non-isolated local minimum point.
\end{itemize}

The last statement allows us to generalize \Cref{cor:appNN:neither-nor-stable} to any non-isolated local minimum point in most deep neural networks. An interesting corollary following from these properties is that there are only \emph{finite} possible loss values of stationary points. Thus the final training loss can only take \emph{finite} possible values. 

