On the global convergence of randomized coordinate gradient descent for non-convex optimization

Coordinate gradient descent is a popular approach in optimization, see e.g., review articles [Wright-15, shi2016primer]. Advantages of coordinate gradient method include that compared with the full gradient descent, it allows larger stepsize [Nesterov-12] and enjoys faster convergence [Saha-13], and it is also friendly for parallelization [Liu-15, Richtarik-16].

The convergence of coordinate gradient descent has been analyzed in several settings on the property of the objective function and on the strategy of coordinate selection. The understanding of convergence for convex problems is quite complete: For methods with cyclic choice of coordinates, the convergence has been established in [Beck-13, Saha-13, Sun-15] and the worst-case complexity is investigated when the objective function is convex and quadratic in [Sun-19]. For methods with random choice of coordinates, it is shown in [Nesterov-12] that $\mathbb{E}f(x_{t})$ converges to $f^{*}=\min_{x\in\mathbb{R}^{d}}f(x)$ sublinearly in the convex case and linearly in the strongly convex case. Convergence of objective function in high probability has also been established in [Nesterov-12]. We also refer to [Richtarik-14, Liu-14, Liu-15, Wright-15] for further convergence results for random coordinate selection for convex problems. More recently, convergence of methods with random permutation of coordinates (i.e., a random permutation of the $d$ coordinates is used for every $d$ step of the algorithm) have been analyzed, mostly for the case of quadratic objective functions [Lee-19, Oswald-17, Gurbuzbalaban-20, Wright-20]. It has been an ongoing research direction to compare various coordinate selection strategies in various settings. In addition, in the non-convex and non-smooth setting, the convergence of coordinate/alternating descent methods can be analyzed for tame/semi-algebraic functions with Kurdyka-Łojasiewicz property (see e.g. [attouch2013convergence, attouch2010proximal, bolte2014proximal, boct2020inertial]).

For non-convex objective functions, the global convergence analysis is less developed, as the situation becomes more complicated. Escaping strict saddle points has been a focused research topic in non-convex optimization, motivated by applications in machine learning. It has been established that various first-order algorithms with gradient noise or added randomness to iterates would escape strict saddle points, see e.g., [Ge-15, Levy-16, Jin-17, Jin-18, Jin-19, Guo-20] for works in this direction.

Among previous works for escaping saddle points, perhaps the closest in spirit to our current result are [Lee-16, ONeill-19, Lee-2019, leadeigen-19], where algorithms without gradient or iterate randomness are studied. It is proved in [Lee-16] that for almost every initial guess, the trajectory of the gradient descent algorithm (without any randomness) with constant stepsize would not converge to a strict saddle point. The result has been extended in [Lee-2019] to a broader class of deterministic first-order algorithms, including coordinate gradient descent with cyclic choice of coordinate. The global convergence result for cyclic coordinate gradient descent is also proved in [leadeigen-19] under slightly more relaxed conditions. Similar convergence result is also obtained for heavy-ball method in [ONeill-19]. Let us emphasize that in the case of coordinate algorithms, it is not merely a technical question whether the algorithm can escape the strict saddle points without randomly perturbing gradients or iterates. In fact, one simply cannot employ such random perturbations, e.g., adding a random Gaussian vector to the iterate, since doing so would destroy the coordinate nature of the algorithm.

The analysis in works [Lee-16, Lee-2019, ONeill-19, leadeigen-19] is based on viewing the algorithm as a deterministic dynamical system, and applying center-stable manifold theorem for deterministic dynamical system [Shub-87], which characterizes the local behavior near a stationary point of nonlinear dynamical systems. Such a framework obviously does not work for randomized algorithms. To some extent, our analysis can be understood as a natural generalization to the framework of random dynamical systems, which allows us to analyze the long time behavior of randomized algorithms, in particular coordinate gradient descent with random coordinate selection.

Let us mention that various stable, unstable, and center manifolds theorems have been established in the literature of random dynamical systems, see e.g., [Arnold_RDS, Ruelle-1979, Ruelle-82, Boxler-89-center, Liu-06]. These sample-dependent random manifolds also characterize the local behavior of random dynamical systems. However, as far as we can tell, one cannot simply apply such “off-the-shelf results” for the analysis of Algorithm 1. Instead, for study of the algorithm, we have to carry out a quantitative finite block analysis for the random dynamical system near the stationary points. Our proof technique is inspired by stability analysis of Lyapunov exponent of random dynamical systems, as in [Ledrappier-91, Froyland-15].

The rest of this paper will be organized as follows. In Section 2, we review the preliminaries of random dynamical system, for convenience of readers. Our main result is stated in Section 3. The proofs can be found in Section 4.

Let us first specify the random dynamical system corresponding to the Algorithm 1.

• •

  Probability space. For each $t\in\mathbb{N}$, denote $(\Omega_{t},\Sigma_{t},\mathbb{P}_{t})$ the usual probability space for the distribution $\mathcal{U}_{\{1,2,\dots,n\}}\times\mathcal{U}_{[\alpha_{\min},\alpha_{\max}]}$, where $\mathcal{U}_{\{1,2,\dots,n\}}$ and $\mathcal{U}_{[\alpha_{\min},\alpha_{\max}]}$ are the uniform distributions on the set $\{1,2,\dots,n\}$ and interval $[\alpha_{\min},\alpha_{\max}]$ respectively. Let $(\Omega,\mathcal{F},\mathbb{P})$ be the product probability space of all $(\Omega_{t},\Sigma_{t},\mathbb{P}_{t}),\ t\in\mathbb{N}$. Denote $\pi_{t}$ as the projection from $(\Omega,\mathcal{F},\mathbb{P})$ onto $(\Omega_{t},\Sigma_{t},\mathbb{P}_{t}),\ t\in\mathbb{N}$. Thus, a sample $\omega\in\Omega$ can be represented as a sequence $((i_{0},\alpha_{0}),(i_{1},\alpha_{1}),\dots)$, where $(i_{t},\alpha_{t})=\pi_{t}(\omega),\ t\in\mathbb{N}$. Let $\{\mathcal{F}_{t}\}_{t\in\mathbb{N}}$ be the filtration defined by

  $$\mathcal{F}_{t}=\sigma\left\{(B_{0}\times\cdots\times B_{t})\times\left(\prod_{j>t}\Omega_{j}\right):B_{i}\in\Sigma_{i},\ i=0,1,\dots,t\right\}.$$

• •

  Metric dynamical system. The metric dynamical system on $\Omega$ is constructed by the (left) shifting operator $\tau:\Omega\rightarrow\Omega$ defined as

  $$\tau(\omega)=\tau(\pi_{0}(\omega),\pi_{1}(\omega),\cdots):=(\pi_{1}(\omega),\pi_{2}(\omega),\cdots),$$

  which is clearly measurable and $\mathbb{P}$-preserving. The metric dynamical system is then given by $\theta(t)=\tau^{t}$ for $t\in\mathbb{N}$.

• •

  Random dynamical system. For any $\omega\in\Omega$ and $t\in\mathbb{N}$, we define $\phi(\omega)$ to be a (nonlinear) map on $\mathbb{R}^{d}$ as

  $$\begin{split}\phi(\omega):\mathbb{R}^{d}&\rightarrow\quad\quad\quad\mathbb{R}^{d}\\ x\ &\mapsto x-\alpha e_{i}e_{i}^{\top}\nabla f(x),\end{split}$$

  where $(i,\alpha)=\pi_{0}(\omega)$ is the first pair/element in the sequence $\omega$, and we define the map $\varphi(t,\omega)$ via

  $$\varphi(t,\omega)=\phi(\tau^{t-1}\omega)\circ\cdots\circ\phi(\tau\omega)\circ\phi(\omega),\qquad\text{for }t\geq 1,$$

  while $\varphi(0,\omega)$ is the identity operator. It is clear that $\varphi(t,\omega)$ satisfies the cocycle property (2.1) and hence defines a random dynamical system on $X=\mathbb{R}^{d}$ over $\{\tau^{t}\}_{t\in\mathbb{N}}$. The iterate of Algorithm 1 follows the random dynamical system as

  $$x_{t}=\phi(\tau^{t-1}\omega)x_{t-1}=\cdots=\phi(\tau^{t-1}\omega)\circ\cdots\circ\phi(\tau\omega)\circ\phi(\omega)x_{0}=\varphi(t,\omega)x_{0}.$$

  It can be seen that $\{x_{t}\}_{t\in\mathbb{N}}$ is $\{\mathcal{F}_{t}\}$-predictable, i.e., $x_{t}$ is $\mathcal{F}_{t-1}$-measurable for any $t\in\mathbb{N}_{+}$, since $x_{t}$ is determined by samples $(i_{0},\alpha_{0}),(i_{1},\alpha_{1}),\dots,(i_{t-1},\alpha_{t-1})$.

In our analysis, we will use linearization of the dynamical system $\varphi(t,\omega)$ at a critical point $x^{*}$ of $f$. Without loss of generality, we assume $x^{*}=0$; otherwise we consider the system with state being $x-x^{*}$. The resulting linear system, which depends on $H=\nabla^{2}f(x^{*})=(H_{ij})_{1\leq i,j\leq d}$, is given by (here and in the sequel, we use the superscript $H$ to indicate dependence on the matrix)

where

Note that $A^{H}(\cdot)$ is bounded in $\Omega$. We know that $\left(\log\left\|A^{H}(\cdot)\right\|\right)_{+}$ is integrable. When $\alpha<1/|H_{ii}|$, the matrix $A^{H}(\omega)=I-\alpha e_{i}e_{i}^{\top}H$ is invertable, and the inverse is given explicitly by applying the Sherman-Morrison formula:

In particular, we have

Thus, if we take the maximal stepsize $\alpha_{\max}$ such that $\alpha_{\max}<1/\max_{1\leq i\leq d}|H_{ii}|$, $\left\|A^{H}(\cdot)^{-1}\right\|$ is bounded in $\Omega$, and as a result $\left(\log\left\|A^{H}(\cdot)^{-1}\right\|\right)_{+}$ is also integrable. Therefore, the assumptions of Theorem 2.3 hold. The shifting operator $\tau$ is ergodic on $(\Omega,\mathcal{F},\mathbb{P})$ by Kolmogorov’s $0$–$1$ law. Then Theorem 2.3 applies for $\theta=\tau$ with $p^{H}(\cdot)$, $\lambda_{i}^{H}(\cdot)$, and $d_{i}^{H}(\cdot)$ all being a.e. constant. For any $\omega\in\widetilde{\Omega}$ that is the set in Theorem 2.3 satisfying $\mathbb{P}(\widetilde{\Omega})=1$, we denote

Then the following invariant property holds

Note that $W_{-}^{H}(\omega)$ works as a center-stable subspace. That is, for any $x\in W_{-}^{H}(\omega)$ and any $\epsilon>0$, it holds that $\left\|\Phi(t,\omega)x\right\|\leq e^{t\epsilon}$, for sufficiently large $t$, and for $x\notin W_{-}^{H}(\omega)$, $\left\|\Phi(t,\omega)x\right\|$ grows exponentially as $t\rightarrow\infty$ with rate greater than $\min_{\lambda_{i}>0}\lambda_{i}-\epsilon$ for any $\epsilon>0$.

In this section, we specify the assumptions of the objective function $f$ in this paper. The first is a standard smoothness assumption of $f$.

An optimization algorithm is expected to converge, under some reasonable assumptions, to a critical point of $f$ where the gradient vanishes. Our aim is to further characterize the possible limits of the algorithm iterates. For this purpose, we distinguish $\mathrm{Crit}_{s}(f)$, the set of all strict saddle points (including local maxima with non-degenerate Hessian) of $f$:

where we use the subscript $s$ to emphasize that it is strict. Due to the presence of negative eigenvalue of Hessian, if we were considering the gradient dynamics near the critical point, the saddle point would be an unstable equilibrium. Our first result is that this instability also occurs in the linear random dynamical system $\Phi^{H}(t,\omega)$ where $H=\nabla^{2}f(x^{*})$. In other words, the dimension of $W^{H}_{+}(\omega)$ defined in (3.5) is greater than $0$. While this would mainly serve as a preliminary step for our analysis of the nonlinear dynamics, the conclusion by itself might be of interest, and is stated as follows. The proof will be deferred to Section 4.1.

Our goal is to generalize such results to the nonlinear dynamics near strict saddle points of $f$, for which, we would require two additional assumptions, as the following.

Assumption 2 is also a standard assumption, which in particular guarantees that each strict saddle point is isolated, due to the non-degenerate Hessian. For each strict saddle point, Proposition 3.1 guarantees that the corresponding unstable subspace $W^{H}_{+}(\omega)$ is non-trivial (has dimension at least $1$). We would in fact require a stronger technical assumption on its structure.

We expect that Assumption 3 holds generically. We also show in Appendix A, Assumption 3 can be verified when $H$ has no zero off-diagonal elements (and $W^{H}_{+}(\omega)$ is not trivial). However, there exist cases that Assumption 3 might not hold. One example is $H=\nabla^{2}f(x^{*})=\text{diag}(H_{1},H_{2})$ where $H_{1}\in\mathbb{R}^{d_{1}\times d_{1}}$ only has positive eigenvalues and $H_{2}\in\mathbb{R}^{d_{2}\times d_{2}}$ only has negative eigenvalues, which implies that $W_{-}^{H}(\omega)=\text{span}\{e_{i}:1\leq i\leq d_{1}\}$ and $W_{+}^{H}(\omega)=\text{span}\{e_{i}:d_{1}+1\leq i\leq d\}$.

We also remark that Assumption 1 is essential in our framework, since the linearized system is defined using the Hessian matrix. Analysis of randomized coordinate method for non-smooth optimization problems requires new techniques and deserves future research.

Given an initial guess $x_{0}$, the behavior of the algorithm, in particular, the limit of $x_{t}$, depends on the particular sample $\omega\in\Omega$. For any $x^{*}\in\mathrm{Crit}_{s}(f)$, we denote the set of all $\omega$ such that the algorithm starting at $x_{0}$ would converge to $x^{*}$:

We further define the set $\Omega_{s}(\mathrm{Crit}_{s}(f),x_{0})$ as the union of all $\Omega_{s}(x^{*},x_{0})$ over $x^{*}\in\mathrm{Crit}_{s}(f)$,

Thus, if $\omega\notin\Omega_{s}(\mathrm{Crit}_{s}(f),x_{0})$, the limit $\lim_{t\to\infty}x_{t}$, if exists, will not be one of the strict saddle points. Our main result in this paper proves that the set is of measure zero, i.e., for any initial guess $x_{0}$ that is not a strict saddle point, with probability $1$, Algorithm 1 will not converge to a strict saddle point.

The intuition behind the proof of Theorem 1 is to compare the nonlinear dynamics around a strict saddle point $x^{*}\in\mathrm{Crit}_{s}(f)$ with its linearization, as the linear dynamics has non-trivial unstable subspace, thanks to Proposition 3.1. Ideally, one would hope that the nonlinear dynamics would closely follow the linear dynamics, and thus leave the neighborhood of $x^{*}$ eventually; the obstacle for such argument is however that the approximation of the linearization is only valid for a finite time interval. Therefore, to establish the instability behavior of the nonlinear dynamics, we would need a much more refined and quantitative argument using the instability of the linear system. In particular, we would need to show that over a finite interval, with high probability, the linear system, and hence the nonlinear system, would drive $x_{t}$ away from the strict saddle point with quantitative bounds; see Theorem 4.4 in Section 4.2. Theorem 1 then follows from an argument with a similar spirit as law of large numbers; see Section 4.3.

As an application of our main result Theorem 1, we can obtain the global convergence to stationary points with no negative Hessian eigenvalues for Algorithm 1. More specifically, denote by

the set of all critical points of $f$, we have the following corollary, which will also be proved in Section 4.3.

We will first study the linear dynamical system, for which we assume the objective function is given by

where $H$ is a symmetric matrix in $\mathbb{R}^{d\times d}$ with at least one negative eigenvalue. In this case, the coordinate descent algorithm is given by

which corresponds to the linear dynamical system $\Phi^{H}(t,\omega)$ with single step map $A^{H}(\omega)$, defined in (3.1) and (3.2) respectively.

Our main goal in this subsection is to prove Proposition 3.1 for this linear dynamical system which states that at least one Lyapunov exponent of $\Phi^{H}(t,\omega)$ is positive. It suffices to show that there exists some $x_{0}$, such that $\left\|x_{t}\right\|$ grows exponentially to infinity, which will follow from an energy argument, similar to the proof of [Lee-2019, Proposition 5]. Although we consider a randomized coordinate gradient descent algorithm instead of a cyclic one, one step, i.e., Lemma 4.3, in the proof of Proposition 4.1 follows closely the proof in [Lee-2019, Appendix A]. We start from $x_{0}$ with $f^{H}(x_{0})<0$ and consider a finite time interval with length $m\geq d$. Proposition 4.1 establishes a quantitative decay estimate for $f^{H}(x_{t+m})$ compared with $f^{H}(x_{t})$, which leads to our desired result Proposition 3.1.

We finish the proof by establishing Lemma 4.3 below.