Using Negative Curvature in Solving Nonlinear Programs

TL;DR

This paper extends negative curvature methods from unconstrained to equality constrained nonlinear optimization, proving convergence to second-order stationary points and discussing potential generalizations.

Contribution

It introduces a negative curvature approach for constrained nonlinear programs and proves its convergence to second-order stationary points.

Findings

Proves convergence to second-order stationary points.

Extends negative curvature methods to equality constrained problems.

Discusses potential for broader nonlinear program applications.

Abstract

Minimization methods that search along a curvilinear path composed of a non-ascent nega- tive curvature direction in addition to the direction of steepest descent, dating back to the late 1970s, have been an effective approach to finding a stationary point of a function at which its Hessian is positive semidefinite. For constrained nonlinear programs arising from recent appli- cations, the primary goal is to find a stationary point that satisfies the second-order necessary optimality conditions. Motivated by this, we generalize the approach of using negative curvature directions from unconstrained optimization to nonlinear ones. We focus on equality constrained problems and prove that our proposed negative curvature method is guaranteed to converge to a stationary point satisfying second-order necessary conditions. A possible way to extend our proposed negative curvature method to general nonlinear programs is also briefly discussed.

Tables

Table 1. Table 1: Examples of Ω Ω \Omega and Π Ω ​ ( ⋅ ) subscript Π Ω ⋅ \Pi_{\Omega}(\cdot)

Constraint Sets	$\Omega$	${\Pi }_{\Omega }(\cdot )$
spherical	$\left\{𝑿\in {ℝ}^{n\times m}\right|{\Vert 𝑿\Vert }_F=1\}$	$𝑿/{\Vert 𝑿\Vert }_F$
multiple spherical	${\bigcup }_i\left\{{𝒙}_i\in {ℝ}^{n_i}\right|{\Vert {𝒙}_i\Vert }_2=1\}$	${𝒙}_i/{\Vert {𝒙}_i\Vert }_2,\forall i$
orthogonality	$\left\{𝑿\in {ℝ}^{n\times m}\right|{𝑿}^T𝑿=I\}$	$𝑼{𝑽}^{\top },\text{with}𝑿=𝑼𝚺{𝑽}^{\top }\text{ as SVD}$

Equations

$$\text{minimize}f(\mathbf{x})\text{subject to}{c}_i(\mathbf{x})=0,i\in \left\{1,2,\dots ,m\right\},$$

$$\mathcal{L}(\mathbf{x},\mathbf{\lambda }):=f(\mathbf{x})-\sum _{i=1}^m{\lambda }_i{\mathbf{c}}_i(\mathbf{x}),$$

$$\mathcal{G}({\mathbf{x}}^{\star }):={\nabla }_{\mathbf{x}}\mathcal{L}\left({\mathbf{x}}^{\star },{\mathbf{\lambda }}^{\star }\left({\mathbf{x}}^{\star }\right)\right)=0$$

$$\mathcal{H}({\mathbf{x}}^{\star }):={\mathcal{P}}_{{\mathbf{x}}^{\star }}^{\top }{\nabla }_{\mathbf{x}\mathbf{x}}^2\mathcal{L}\left({\mathbf{x}}^{\star },{\mathbf{\lambda }}^{\star }\left({\mathbf{x}}^{\star }\right)\right){\mathcal{P}}_{{\mathbf{x}}^{\star }}⪰0$$

$$\mathcal{T}=\sum _{i=1}^n\underset{ptimes}{\underbrace{{\mathbf{v}}_i\otimes {\mathbf{v}}_i\otimes \cdots \otimes {\mathbf{v}}_i}}\in {\mathbb{R}}^{\stackrel{ptimes}{\overbrace{n\times n\times \cdots \times n}}},$$

$$\underset{\mathbf{x}\in {\mathbb{R}}^n}{\min }T(\mathbf{x},\mathbf{x},\mathbf{x},\mathbf{x})=\sum _{i=1}^n\sum _{j=1}^n\sum _{k=1}^n\sum _{l=1}^n{\mathcal{T}}_{ijkl}{x}_i{x}_j{x}_k{x}_l\text{mbox}s.t.{\Vert \mathbf{x}\Vert }_2=1.$$

$$\left\{\mathbf{x}\in {\mathbb{R}}^n|G(\mathbf{x})=0,H(\mathbf{x})⪰0\right\}=\left\{\pm {\mathbf{v}}_1,\dots ,\pm {\mathbf{v}}_n\right\}.$$

$$\underset{[{\mathbf{x}}_1,{\mathbf{x}}_2,\dots ,{\mathbf{x}}_n]\in {\mathbb{R}}^{n\times n}}{\min }\sum _{i=1}^n\sum _{j\ne i}\mathcal{T}({\mathbf{x}}_i,{\mathbf{x}}_i,{\mathbf{x}}_j,{\mathbf{x}}_j)\text{mbox}s.t.{\Vert {\mathbf{x}}_i\Vert }_2=1,i=1,\dots ,n,$$

$$\underset{\mathbf{X}\in {\mathcal{S}}_{+}^{n\times n}}{\min }\text{mbox}Tr(\mathbf{C}\mathbf{X})\text{mbox}s.t.\text{mbox}Tr({\mathbf{A}}_i\mathbf{X})=b_{i}i=1,\dots ,m,$$

$$\underset{\mathbf{L}\in {\mathbb{R}}^{n\times p}}{\min }\text{mbox}Tr({\mathbf{L}}^{\top }\mathbf{C}\mathbf{L})\text{mbox}s.t.\text{mbox}Tr({\mathbf{L}}^{\top }{\mathbf{A}}_i\mathbf{L})=b_{i}i=1,\dots ,m.$$

$$\underset{\mathbf{X}\in {\mathcal{S}}_{+}^{n\times n}}{\min }\text{mbox}Tr(\mathbf{C}\mathbf{X})\text{mbox}s.t.\mathrm{diag}\left(\mathbf{X}\right)=1,$$

$$\underset{\mathbf{L}\in {\mathbb{R}}^{n\times p}}{\min }\text{mbox}Tr({\mathbf{L}}^{\top }\mathbf{C}\mathbf{L})\text{mbox}s.t.\Vert {\mathbf{e}}_i^{\top }\mathbf{L}\Vert =1i=1,\dots ,n.$$

$${\mathbf{\lambda }}^{\star }(\mathbf{x})\in \arg \underset{\mathbf{\lambda }}{\min }\Vert {\nabla }_{\mathbf{x}}\mathcal{L}(\mathbf{x},\mathbf{\lambda })\Vert .$$

$$\mathcal{G}(\mathbf{x}):={\nabla }_{\mathbf{x}}\mathcal{L}\left(\mathbf{x},{\mathbf{\lambda }}^{\star }\left(\mathbf{x}\right)\right)$$

$$\mathcal{H}(\mathbf{x}):={\mathcal{P}}_{\mathbf{x}}^{\top }{\nabla }_{\mathbf{x}\mathbf{x}}^2\mathcal{L}\left(\mathbf{x},{\mathbf{\lambda }}^{\star }\left(\mathbf{x}\right)\right){\mathcal{P}}_{\mathbf{x}}.$$

$$\left|f({\Pi }_{\Omega }({\mathbf{x}}_0+\mathbf{\delta }))-f({\mathbf{x}}_0)-\mathcal{G}({\mathbf{x}}_0{)}^{\top }\mathbf{\delta }\right|\le C_0{\Vert \mathbf{\delta }\Vert }^2,\text{and}$$

$$\left|f({\Pi }_{\Omega }({\mathbf{x}}_0+\mathbf{\delta }))-f({\mathbf{x}}_0)-\mathcal{G}({\mathbf{x}}_0{)}^{\top }\mathbf{\delta }-\frac{1}{2}{\mathbf{\delta }}^{\top }\mathcal{H}({\mathbf{x}}_0)\mathbf{\delta }\right|\le C_5{\Vert \mathbf{\delta }\Vert }^3$$

$${\mathbf{x}}_{k+1}\leftarrow {\Pi }_{\Omega }\left({\mathbf{x}}_k-t\mathcal{G}({\mathbf{x}}_k)\right)$$

$${\mathbf{x}}_{k+1}\leftarrow {\Pi }_{\Omega }\left({\mathbf{x}}_k-t_1\mathcal{G}({\mathbf{x}}_k)-t_2\mathbf{d}\right)$$

$$\hat{\mathbf{x}}\leftarrow {\Pi }_{\Omega }({\mathbf{x}}_k-t{\mathcal{G}}_k).$$

$$\hat{\mathbf{x}}\leftarrow {\Pi }_{\Omega }({\mathbf{x}}_k-t{\mathcal{G}}_k+t^{\alpha }{\mathbf{d}}_k).$$

$$f(\mathbf{x}(t))-f({\mathbf{x}}_k)\le -\Vert {\mathcal{G}}_k{\Vert }^{2}t+C_0\Vert {\mathcal{G}}_k{\Vert }^2{t}^2;$$

$$f(\mathbf{x}(t))-f({\mathbf{x}}_k)\le -\Vert {\mathcal{G}}_k{\Vert }^{2}t+\frac{1}{2}\Vert {\mathcal{H}}_k\Vert \Vert {\mathcal{G}}_k{\Vert }^2{t}^2-\frac{1}{2}|{\lambda }_k{|}^3{t}^{2\alpha }+8C_5\Vert {\mathcal{G}}_k{\Vert }^3{t}^3+8C_5|{\lambda }_k{|}^3{t}^{3\alpha }.$$

$$f(\mathbf{x}(t))-f({\mathbf{x}}_k)\le ⟨{\mathcal{G}}_k,-t{\mathcal{G}}_k⟩+C_0\Vert {\mathcal{G}}_k{\Vert }^2{t}^2=-\Vert {\mathcal{G}}_k{\Vert }^{2}t+C_0\Vert {\mathcal{G}}_k{\Vert }^2{t}^2.$$

$$f(\mathbf{x}(t))-f({\mathbf{x}}_k)$$

$$\le ⟨{\mathcal{G}}_k,-t{\mathcal{G}}_k+t^{\alpha }{\mathbf{d}}_k⟩+\frac{1}{2}(-t{\mathcal{G}}_k+t^{\alpha }{\mathbf{d}}_k{)}^{\top }{\mathcal{H}}_k(-t{\mathcal{G}}_k+t^{\alpha }{\mathbf{d}}_k)+C_5\Vert -t{\mathcal{G}}_k+t^{\alpha }{\mathbf{d}}_k{\Vert }^3$$

$$\le -\Vert {\mathcal{G}}_k{\Vert }^{2}t+⟨{\mathcal{G}}_k,{\mathbf{d}}_k⟩t^{\alpha }+\frac{1}{2}{\mathcal{G}}_k^{\top }{\mathcal{H}}_k{\mathcal{G}}_k{t}^2-{\mathcal{G}}_k^{\top }{\mathcal{H}}_k{\mathbf{d}}_k{t}^{\alpha +1}+\frac{1}{2}{\mathbf{d}}_k^{\top }{\mathcal{H}}_k{\mathbf{d}}_k{t}^{2\alpha }+C_5\Vert -t{\mathcal{G}}_k+t^{\alpha }{\mathbf{d}}_k{\Vert }^3.$$

$$⟨{\mathcal{G}}_k,{\mathbf{d}}_k⟩t^{\alpha }\le |{\lambda }_k|\text{mbox}sign(-{\mathbf{v}}_k^{\top }{\mathcal{G}}_k){\mathbf{v}}_k^{\top }{\mathcal{G}}_k{t}^{\alpha }\le 0$$

$$\frac{1}{2}{\mathcal{G}}_k^{\top }{\mathcal{H}}_k{\mathcal{G}}_k{t}^2\le \frac{1}{2}\Vert {\mathcal{H}}_k\Vert {\Vert {\mathcal{G}}_k\Vert }^2{t}^2$$

$$-{\mathcal{G}}_k^{\top }{\mathcal{H}}_k{\mathbf{d}}_k{t}^{\alpha +1}=-{\mathcal{G}}_k^{\top }|{\lambda }_k|\text{mbox}sign(-{\mathbf{v}}_k^{\top }{\mathcal{G}}_k){\mathcal{H}}_k{\mathbf{v}}_k{t}^{\alpha +1}={\lambda }_k|{\lambda }_k|\text{mbox}sign(-{\mathbf{v}}_k^{\top }{\mathcal{G}}_k)(-{\mathcal{G}}_k^{\top }{\mathbf{v}}_k)t^{\alpha +1}\le 0$$

$$\frac{1}{2}{\mathbf{d}}_k^{\top }{\mathcal{H}}_k{\mathbf{d}}_k{t}^{2\alpha }=\frac{1}{2}|{\lambda }_k{|}^2{\mathbf{v}}_k^{\top }{\mathcal{H}}_k{\mathbf{v}}_k{t}^{2\alpha }=-\frac{1}{2}|{\lambda }_k{|}^3{t}^{2\alpha }$$

$$\Vert -t{\mathcal{G}}_k+t^{\alpha }{\mathbf{d}}_k{\Vert }^3\le {\left(t\Vert {\mathcal{G}}_k\Vert +t^{\alpha }\Vert {\mathbf{d}}_k\Vert \right)}^3\le 8\max (t\Vert {\mathcal{G}}_k\Vert ,t^{\alpha }\Vert {\mathbf{d}}_k\Vert {)}^3\le 8t^3{\Vert {\mathcal{G}}_k\Vert }^3+8t^{3\alpha }|{\lambda }_k{|}^3.$$

$$\Vert \mathcal{H}(\mathbf{x})\Vert$$

$$\le \Vert {\nabla }_{\mathbf{x}\mathbf{x}}^2\mathcal{L}(\mathbf{x},{\mathbf{\lambda }}^{\star }(\mathbf{x}))\Vert$$

$$\le \Vert {\nabla }^{2}f(\mathbf{x})-\sum _{i\in [m]}{\lambda }_i^{\star }(\mathbf{x}){\nabla }^2{c}_i(\mathbf{x})\Vert$$

Full text

1 Introduction

This paper is concerned with solving general smooth nonlinear optimization problems. Considering the difficulty of obtaining a global optimal solution, our goal is to develop an efficient method to locate a stationary point that satisfies the second-order necessary conditions. This goal is motivated by recent applications and developments, where stationary points satisfying second-order necessary conditions are of primary interest. Specifically, in many nonlinear programs arising from dictionary learning SQW152 , tensor decomposition GHJY15 , phase retrieval SQW16 , and semidefinite programming BM03 , stationary points satisfying the second-order necessary conditions are the targets of the application and can often be proven to be global optimizers.

Even though the scope of our concern is more general than unconstrained nonlinear programs, studying such problems helps us to better understand how to achieve our goal. To find a stationary point of an unconstrained function satisfying second-order necessary conditions, apparently, using gradient information alone in a deterministic method is not sufficient as such methods can be trapped by saddle points. Therefore, second-order information should be utilized in (hopefully) some computationally efficient manner. Methods, e.g. Newton’s method and trust-region methods, that involve computing the inverse, Cholesky factorization or eigen-decomposition of the Hessian or some modifications of it, might become impractical when the number of variables grow to tens of thousands. As a remedy, methods that utilize negative curvature directions were proposed in MC77 ; MS79 ; GOLD80 and later further actively developed in grandinetti1984nonlinear ; martinez1990algorithm ; ferris1996nonmonotone ; lucidi1998curvilinear ; gould2000exploiting ; conforti2001curvilinear ; apostolopoulou2010curvilinear based on inexpensive computations involving the Hessian. Specifically, these methods imitate the steepest descent method but also incorporate a negative curvature direction besides the negative of the gradient when computing the step to take at each iteration. Consequently, saddle points are usually avoided and these methods can be guaranteed to converge to a stationary point satisfying the second-order necessary conditions without requiring much additional computational cost.

Inspired by the unconstrained case, we generalize the negative curvature approach to equality constrained nonlinear programs. In order to accomplish this, one natural question is how to generalize the notions of the gradient and the Hessian when constraints are imposed. For unconstrained problems, we can regard the gradient and the Hessian as quantities that characterize local optimality. Specifically, according to the first and second-order necessary optimality conditions, a local minimizer must have a zero gradient and a positive semidefinite Hessian. For constrained nonlinear programs, it seems plausible that a generalized gradient and Hessian could be developed based on the first-order and second-order optimality conditions. Essentially, this is the approach we adopt here. Our generalized gradient and Hessian are closely related to the Riemannnian gradient and Hessian in the context of optimization on Riemannian manifolds petersen2006riemannian ; AMS09 ; absil2013extrinsic ; absil2009all . This explicit connection is discussed in the paper. Our approach of deriving generalized gradient and Hessian from the classical optimality point of view is quite well-motivated and self-contained, and moreover keeps the technical difficulty to a minimum level.

Organization.

The rest of the paper is organized as follows. In Section 2, we review some basic results on equality constrained optimization. In Section 3, we present our negative curvature line search method for equality-constrained nonlinear programs in particular. The notions of gradient and Hessian are generalized based on local optimality conditions. Using these generalizations, we propose and study a negative curvature method for equality constrained problem. Last, we briefly discuss how to extend the negative curvature method proposed for the equality constrained problem to general nonlinear programming problems.

2 Preliminaries

In this section, we review several fundamental results regarding the equality-constrained problem:

[TABLE]

where the objective function $f:\mathbb{R}^{n}\to\mathbb{R}$ and constraint functions $c_{i}:\mathbb{R}^{n}\to\mathbb{R}$ are smooth so that it is possible to characterize local optimality conditions based on their derivatives. The results and notation in this section are quite classical nocedal2006numerical . Regarding the feasible set $\Omega:=\left\{\bm{x}\in\mathbb{R}^{n}\;|\;c_{i}(\bm{x})=0\right\}$, we first impose a regularity condition which is known as LICQ in optimization literature.

Definition 1** (LICQ).**

We say that the linear independence constraint qualification (LICQ) holds at $\bm{x}\in\Omega$, if the set of constraint gradients $\left\{\nabla c_{1}(\bm{x}),\nabla c_{2}{(\bm{x})},\ldots,\nabla c_{m}{(\bm{x})}\right\}$ are linearly independent, i.e. the matrix $\nabla\bm{c}({\bm{x}}):=[\nabla c_{1}(\bm{x}),\nabla c_{2}{(\bm{x})},\ldots,\nabla c_{m}{(\bm{x})}]\in\mathbb{R}^{n\times m}$ has full column rank.

From a geometric point of view, the LICQ condition guarantees ${\Omega}$ to be a differentiable manifold of dimension $n-m$. Now we are ready to state the definitions of tangent and normal subspaces, which are particularly useful geometry concepts to characterize the variational properties of $f(\cdot)$.

Definition 2** (Tangent and Normal Subspaces).**

The normal subspace $\mathcal{N}_{\Omega}(\bm{x})$ at a point $\bm{x}\in\Omega$ is defined as the subspace spanned by the set of constraint gradients $\left\{\nabla c_{1}(\bm{x}),\ldots,\nabla c_{m}{(\bm{x})}\right\}$, i.e. $\mathcal{N}_{\Omega}(\bm{x}):=\mathfrak{R}(\nabla\bm{c}({\bm{x}}))\subseteq\mathbb{R}^{n}$ the range of the matrix $\nabla\bm{c}({\bm{x}})$. We denote by $\mathcal{P}_{\mathcal{N}_{\Omega}(\bm{x})}$ projection onto the normal subspace $\mathcal{N}_{\Omega}(\bm{x})$.

The tangent subspace $\mathcal{T}_{\Omega}(\bm{x})$ at a point $\bm{x}\in\Omega$ complements the normal subspace in $\mathbb{R}^{n}$, i.e. $\mathcal{T}_{\Omega}(\bm{x})=\mathcal{N}^{\perp}_{\Omega}(\bm{x})=\mathfrak{R}(\nabla\bm{c}({\bm{x}}))^{\perp}=\mathfrak{N}(\nabla\bm{c}({\bm{x}})^{\top})$, is the nullspace of the matrix $\nabla\bm{c}({\bm{x}})^{\top}$. We denote by $\mathcal{P}_{\bm{x}}$ the projection onto the tangent subspace $\mathcal{T}_{\Omega}(\bm{x})$.

Remark 1**.**

We may drop the subscripts $\Omega$ and $\bm{x}$, respectively, from $\mathcal{N}_{\Omega}(\bm{x})$ and $\mathcal{T}_{\Omega}(\bm{x})$ if they are clear from the context.

We are now fully prepared to state the first-order and second-order necessary optimality conditions for problem (2.1) based on the Lagrangian function,

[TABLE]

where $\lambda_{i}$ is the Lagrange multiplier associated with the $i$-th constraint $c_{i}(\bm{x})=0$.

Theorem 1** (Necessary Optimality Conditions nocedal2006numerical ).**

Suppose that $\bm{x}^{\star}$ is a local minimizer of problem (2.1) and satisfies the LICQ condition. Then $\bm{x}^{\star}$ satisfies the following first-order and second-order conditions

[TABLE]

where $\bm{\lambda}^{\star}(\bm{x}^{\star})\in\arg\min_{\bm{\lambda}}\;\left\|\nabla_{\bm{x}}\mathcal{L}(\bm{x}^{\star},\bm{\lambda})\right\|.$

Feasible points satisfying conditions (2.3) and (2.4) are typically referred as second-order critical points. The gist of our paper is to design iterative numerical methods to locate them. We next briefly describe examples of optimization problems that arise in signal processing, machine learning and statistics where such solutions are essentially sought.

Symmetric Orthogonal Tensor Decomposition (SOTD).

Tensor is a multidimensional array and the symmetric orthogonal tensor decomposition (SOTD) naturally generalizes the spectral decomposition of a symmetric matrix. Here, we focus on the $p$-way $n$-dimensional symmetric orthogonal (SOD mu2015successive ; wang2017tensor ; mu2017 ) tensor

[TABLE]

where $\bm{V}=[\bm{v}_{1},\bm{v}_{2},\ldots,\bm{v}_{n}]\in\mathbb{R}^{n\times n}$ is an orthogonal matrix and $\otimes$ denotes the usual outer product, so the $(i_{1},i_{2},\ldots,i_{p})$-th entry of $\bm{v}\otimes\bm{v}\otimes\cdots\otimes\bm{v}$ is the scalar $v_{i_{1}}v_{i_{2}}\cdots v_{i_{p}}$. The problem addressed by SOTD, is to find (up to sign and permutation) the components $\bm{v}_{i}$’s given $\bm{\mathcal{T}}$, has many applications, including higher-order statistical estimation mccullagh1987tensor , independent component analysis comon1994independent ; comon2010handbook , and parameter estimation for latent variable models anandkumar2014tensor ).

Ge et al. GHJY15 consider the SOTD specifically for the $p=4$ case and analyze the geometry of the following minimization problem

[TABLE]

They prove that for problem (2.6), the points satisfying the second-order necessary condition coincide with the component $\bm{v}_{i}$’s. Specifically, they show that

[TABLE]

Therefore, solving problem (2.6) yields one component $\bm{v}_{i}$, and after that, one can apply standard deflation procedures to obtain the others one by one.

Ge et al. GHJY15 further propose a larger optimization problem

[TABLE]

to find all the components at one shot. Although problem (2.8) is substantially more complicated than problem (2.6), similar geometrical phenomenon and the advantageous property of (2.6) are preserved. Specifically, Ge et al. GHJY15 proves that any solution of problem (2.8) that satisfies the second-order necessary conditions in Theorem 1 is a signed and permuted version of $[\bm{v}_{1},\bm{v}_{2},\ldots,\bm{v}_{n}]\in\mathbb{R}^{n\times n}$.

Semidefinite Programming (SDP).

SDP is one of the most exciting developments in mathematical optimization and has been successfully applied to model and solve problems in traditional convex constrained optimization, control theory, and combinatorial optimization. In general, an SDP problem can be defined as

[TABLE]

where $\mathcal{S}_{+}^{n\times n}$ is the cone of symmetric positive semidefinite matrices.

Although interior-point methods can solve SDP problems in polynomial time, scalability is a problem and in practice interior-point methods run out of memory and time once $n$ is greater than $1e3$. To address this issue, Burer and Monteiro burer2003nonlinear ; burer2005local , take advantage of the low-rank structure of SDP optimal solutions characterized by Pataki pataki1998rank and Barvinok barvinok1995problems . Explicitly they replaced $\bm{X}=\bm{L}\bm{L}^{\top}$ where $\bm{L}$ has size $n\times p$, which leads to the following non-convex but low-dimensional surrogate:

[TABLE]

Burer and Monteiro burer2003nonlinear ; burer2005local also applied local optimization methods to solve the above general nonlinear program (2.10) with surprisingly good performance, even though the theoretical justifications between problems (2.9) and (2.10) were unclear.

Recently, Boumal et al. boumal2016non proved that whenever $p$ satisfies $p(p+1)>2m$ (and some other technical but mild conditions hold), problem (2.10) is equivalent to problem (2.9) in the sense that for any second-order critical point of (2.10) $\bm{L}^{\star}$, the square matrix $\bm{L}^{\star}(\bm{L}^{\star})^{\top}$ is optimal to problem (2.9).

One straightforward but powerful application of the above discussion is to solve the SDP problem derived as relaxation for the max-cut problem. In their celebrated work goemans1995improved , Goemans and Williamson tackled the NP-hard max-cut problem by the following SDP

[TABLE]

whose optimal solution (after rounding) will yield an approximate solution to max-cut within a ratio of .878. The corresponding Burer-Monteiro reformulation can be written as

[TABLE]

Based on the above discussion, when we choose $p=\lceil\sqrt{2n}\rceil$, even though problem (2.12) has substantially smaller dimension than problem (2.11), any second-order critical point $\bm{L}^{\star}$ of (2.12) is able to recover a global optimal solution of (2.11) with $\bm{X}^{\star}=\bm{L}^{\star}(\bm{L}^{\star})^{\top}$.

There is a large body of recent literature, especially in signal processing and learning theory, where finding second-order critical points is the fundamental target for the application at hand. Interested readers can find examples in complete dictionary recovery SQW152 , generalized phase retrieval SQW16 , matrix completion ge2016matrix , phase synchronization and community detection boumal2016nonconvex ; bandeira2016low .

3 Equality-Constrained Problem

In this section, we extend the definitions (2.3) and (2.4) of $\mathcal{G}(\bm{x}^{\star})$ and $\mathcal{H}(\bm{x}^{\star})$ at a local minimizer $\bm{x}^{\star}$ to a general feasible point $\bm{x}\in\Omega$, and verify that $\mathcal{G}(\bm{x})$ and $\mathcal{H}(\bm{x})$ can be regarded as natural generalizations of the conventional gradient and Hessian. Based on this, we extend the classical negative curvature algorithm from the unconstrained problem to the equality constrained one (2.1), and prove this method converges to a second-order critical point of (2.1).

3.1 Generalized gradient and Hessian

Based on Theorem 1, we note that $\mathcal{G}(\bm{x}^{\star})$ and $\mathcal{H}(\bm{x}^{\star})$ provide an elegant characterization of local optimality for the equality constrained problem in analogy with the roles of the gradient and Hessian for the unconstrained one. Inspired by this, we extrapolate the definitions of (2.3) and (2.4) from constrained local minimizers $\bm{x}^{\star}$ to any feasible point $\bm{x}\in\Omega,$ with the intention of using these generalized quantities as one uses the gradient and Hessian in the algorithmic design for unconstrained optimizations.

Definition 3** (Generalized Gradient and Hessian).**

For any feasible $\bm{x}\in\Omega$, let

[TABLE]

Then we define the generalized gradient at $\bm{x}$ as

[TABLE]

and the generalized Hessian at $\bm{x}$ as

[TABLE]

In the following, we verify that $\mathcal{G}(\cdot)$ and $\mathcal{H}(\cdot)$ behave like the gradient and Hessian not only for a local optimal solution but also for any feasible point $\bm{x}\in\Omega$ in an approximate sense, under the mild conditions described below:

Assumption 2**.**

(a)

$\nabla f$, $\nabla^{2}f$, $\nabla c_{i}$ and $\nabla^{2}c_{i}$ are Lipschitz continuous over $\Omega$ with Lipschitz constants $L_{f,1}$, $L_{f,2}$, $L_{c_{i},1}$ and $L_{c_{i},2}$.

(b)

$\sup_{{\bm{x}}\in\Omega}\left\|\nabla f({\bm{x}})\right\|\leq\gamma_{f,1}$, $\sup_{{\bm{x}}\in\Omega}\left\|\nabla^{2}f({\bm{x}})\right\|\leq\gamma_{f,2}$, $\sup_{{\bm{x}}\in\Omega}\left\|\nabla c_{i}({\bm{x}})\right\|\leq\gamma_{c_{i},1}$ and $\\ \sup_{{\bm{x}}\in\Omega}\left\|\nabla^{2}c_{i}({\bm{x}})\right\|\leq\gamma_{c_{i},2}$.

(c)

[$\sigma_{0}$-LICQ]     For every ${\bm{x}}\in\Omega$, $\sigma_{\text{min}}(\nabla\bm{c}({\bm{x}}))\geq\sigma_{0}$ for some $\sigma_{0}>0$, where $\sigma_{\text{min}}(\nabla\bm{c}({\bm{x}}))$ is the smallest singular value of the matrix $\nabla\bm{c}(\bm{x})$

Remark 2**.**

The last assumption is slightly stronger that LICQ as LICQ only requires $\sigma_{\text{min}}(\nabla\bm{c}({\bm{x}}))>0.$

We next prove a key lemma regarding $\mathcal{G}(\cdot)$ and $\mathcal{H}(\cdot)$. Loosely speaking, we show that for every $\bm{x}\in\Omega$, given a small perturbation $\bm{\delta}\in\mathcal{T}_{\Omega}(\bm{x})\subset\mathbb{R}^{n}$, $f({\bm{x}})+\mathcal{G}({\bm{x}})^{\top}{\bm{\delta}}$ and $f({\bm{x}})+\mathcal{G}({\bm{x}})^{\top}{\bm{\delta}}+\frac{1}{2}{\bm{\delta}}^{\top}\mathcal{H}({\bm{x}}){\bm{\delta}}$, respectively approximates $f(\bm{x}+\bm{\delta})$ up to the first order and the second order. This result is crucial for our generalization of negative curvature methods. But one issue we need to fix in advance is that $\bm{x}+\bm{\delta}$ may possibly lie outside of $\Omega$ and thus the objective $f$ might not even be defined at $\bm{x}+\bm{\delta}$. In order to resolve this infeasibility issue, we introduce the following projection operator:

Definition 4**.**

For any $\bm{y}\in\mathbb{R}^{n}$, we denote by $\Pi_{\Omega}(\bm{y})\in\Omega$ a point that is closest to $\bm{y}$, i.e. $\Pi_{\Omega}(\bm{y})\in\arg\min_{\bm{x}\in\Omega}\;\left\|\bm{x}-\bm{y}\right\|$, where $\left\|\cdot\right\|$ denotes the Euclidean norm.

We are now ready to state the key lemma.

Lemma 1**.**

Under Assumptions 2(b), 2(c) and 2(d), for any ${\bm{x}}_{0}\in\Omega$ and ${\bm{\delta}}\in\mathcal{T}_{\Omega}(\bm{x}_{0})$, we have

[TABLE]

where $C_{0}=\gamma_{f,1}\{1+\frac{\Gamma_{1}}{\sigma_{0}^{2}}\}\frac{4}{R^{2}}+4(L_{f,1}+\frac{\gamma_{f,1}\sqrt{\Gamma_{1}\Lambda_{1}}}{\sigma_{0}^{2}})$, $C_{5}=8C_{1}+C_{2}+C_{3}+C_{4},C_{1}=\frac{L_{f,2}}{2}+\frac{\gamma_{f,1}}{2\sigma_{0}^{2}}\sqrt{\Gamma_{1}\Lambda_{2}}$, $C_{2}=\{\gamma_{f,2}+\frac{\gamma_{f,1}}{2\sigma_{0}^{2}}\sqrt{\Gamma_{1}\Gamma_{2}}\}\frac{4}{R}$, $C_{3}=\{\gamma_{f,2}+\frac{\gamma_{f,1}}{2\sigma_{0}^{2}}\sqrt{\Gamma_{1}\Gamma_{2}}\}\frac{2}{R}$, $C_{4}=\{\gamma_{f,1}+\frac{\gamma_{f,1}}{\sigma_{0}^{2}}\Gamma_{1}\}(\frac{2\sqrt{\Gamma_{1}\Lambda_{1}}}{\sigma_{0}^{2}}+\frac{2\sqrt{\Gamma_{1}^{3}\Lambda_{1}}}{\sigma_{0}^{4}})\frac{8}{R}$, $R=\sigma_{0}/\sqrt{\Lambda_{1}}$, $\Gamma_{1}=\sum_{i=1}^{m}\gamma_{c_{i},1}^{2}$, $\Gamma_{2}=\sum_{i=1}^{m}\gamma_{c_{i},2}^{2}$, $\Lambda_{1}=\sum_{i=1}^{m}L_{c_{i},1}^{2}$, and $\Lambda_{2}=\sum_{i=1}^{m}L_{c_{i},2}^{2}$.

Proof.

See the appendix. ∎

Remark 3**.**

If $\Omega$ is a differentiable manifold, our $\mathcal{G}(\cdot)$ and $\mathcal{H}(\cdot)$ defined in (3.2) and (3.3) can be proven to be the Riemannian gradient and Riemannian Hessian of $f$ over $\Omega$ (see subsection 4.2 for details). The expressions (3.2) and (3.3) concretize these abstract geometrical concepts by providing an explicit algebraic way to compute them.

3.2 Negative curvature method

In this subsection, we will present a general framework for using negative curvature of the generalized Hessian $\mathcal{H}(\bm{x})$ to solve (2.1) and show that the sequence generated by this framework converges to a second-order critical point.

It is clear that using the simple projected gradient method

[TABLE]

might be trapped by the saddle points of (2.1). To overcome this difficulty, it is natural to consider using second-order information, that is, information about the Hessian of $f$ and the $c_{i}$’s. In particular, we can use a negative curvature direction of $\mathcal{H}(\bm{x})$ defined in subsection 3.1. We say $\bm{d}$ is a negative curvature direction of $\mathcal{H}(\bm{x})$ if it has the property that $\bm{d}^{\top}\mathcal{H}(\bm{x})\bm{d}<0$. Based on the second-order approximation result revealed in Lemma 1, it is intuitive that moving in a negative curvature direction will enable an algorithm to escape from a saddle point. Therefore, it is tempting to move along a direction that combines the negative gradient direction with a descent direction of negative curvature:

[TABLE]

especially for the regions close to saddle points.

Moreover, negative-curvature directions can be obtained at a relatively small cost. For example, when a Hessian matrix $H(\bm{x})$ is indefinite, the eigenvector corresponding to its algebraically smallest eigenvalue is a negative curvature direction. It can be obtained by executing the power iteration method on $H(\bm{x})$ to obtain a eigenvalue, eigenvector pair $(\lambda_{\text{dom}},\bm{d}_{\text{dom}})$ corresponding to the eigenvalue of largest magnitude, the so-called dominant eigenvalue $\lambda_{\text{dom}}$. If $\lambda_{\text{dom}}<0$, $\bm{d}_{\text{dom}}$ is the direction of most negative curvature. Otherwise, we can perform the method again on $\mathcal{H}-\lambda_{\text{dom}}\mathcal{I}$. More efficient and robust methods using variants of Lanzos algorithm to compute the algebraically smallest eigenpair can be found in larsen2004propack ; stathopoulos2010primme .

In Algorithm 1, we provide a general framework for using negative curvature directions to solve problem (2.1). This algorithm integrates the first-order and second-order methods. Specifically, when the iterate $\bm{x}_{k}$ is far from any saddle point, we only use the first-order information to make progress:

[TABLE]

When $\bm{x}_{k}$ is near saddle points, we combine the negative gradient and the negative curvature direction:

[TABLE]

If $\alpha$ is chosen to be $2$ and $\Omega=\mathbb{R}^{n}$ , (3.8) is reduced to the one leveraged in GOLD80 ; if $\alpha$ is set to $1/2$ and $\Omega=\mathbb{R}^{n}$, (3.8) is equivalent to the one used in MS79 .

In the rest of this section, we will confirm our intuition that Algorithm 1 is capable of escaping saddle points by proving the following theorem:

Theorem 3**.**

For the sequences $\{\bm{\mathcal{G}}_{k}\}$ and $\{\lambda_{k}\}$ generated by Algorithm 1, one has $\bm{\mathcal{G}}_{k}\rightarrow\bm{0}$ and $\lambda_{k}\rightarrow 0$ as $k\rightarrow\infty$.

Remark 4**.**

Based on Theorem 3, any cluster point of $\left\{\bm{x}_{k}\right\}$ generated by Algorithm 1 is a second-order critical point of problem (2.1).

To prove Theorem 3, we need several lemmas.

Lemma 2**.**

Consider the $k$-th iteration.

(a)

For the parametrized curve ${\bm{x}}(t)=\Pi_{\Omega}({\bm{x}}_{k}-t\bm{\mathcal{G}}_{k})$, we have

[TABLE]

(b)

For the parametrized curve ${\bm{x}}(t)=\Pi_{\Omega}({\bm{x}}_{k}-t\bm{\mathcal{G}}_{k}+t^{\alpha}\bm{d}_{k})$, we have

[TABLE]

Proof.

For part (a), we can directly apply (3.4) in Lemma 1:

[TABLE]

Now let us focus on part (b).

For part (b), based on (3.5) of Lemma 1, we first have

[TABLE]

Next, we will bound each term in (3.11):

[TABLE]

With (3.12)-(3.16) plugged in (3.11), we reach (3.10).

Lemma 3**.**

There exists a constant $\gamma_{h}\geq 0$ such that $\sup_{\bm{x}\in\Omega}\left\|\mathcal{H}(\bm{x})\right\|\leq\gamma_{h}.$

Proof.

For any $\bm{x}\in\Omega$, one has

[TABLE]

Based on Assumption 2, we have

[TABLE]

From the definition (3.1). we can derive that

[TABLE]

Together with Assumption 2, it can be obtained that

[TABLE]

The lemma can be established by substituting (3.18) and (3.20) into (3.17). Hence $\gamma_{h}=\sum_{i\in[m]}\gamma_{c_{i},2}+\frac{1}{\sigma_{0}^{2}}\sqrt{\sum_{i\in[m]}\gamma^{2}_{c_{i},1}}\cdot\gamma_{f,1}\sum_{i\in[m]}\gamma_{c_{i},2}$ ∎

∎

Lemma 4**.**

$\{t_{k}\}$* is uniformly bounded from below, i.e., there exists $\beta>0$ such that $t_{k}\geq\beta$ for every $k\in\mathbb{N}$.*

Proof.

Consider the first case $\|\bm{\mathcal{G}}_{k}\|\geq\epsilon$. For ${\bm{x}}(t)=\Pi_{\Omega}({\bm{x}}_{k}-t\bm{\mathcal{G}}_{k})$, as shown in Lemma 2, one has

[TABLE]

Therefore, whenever $t\leq(1-\sigma)/C_{0}$, we have

[TABLE]

Thus, $t_{k}\geq t_{0}\rho^{\lceil\log_{\rho}\{(1-\sigma)/(C_{0}t_{0})\}\rceil}$ for this case.

Now let us consider the other case $\|\bm{\mathcal{G}}_{k}\|<\epsilon$. For ${\bm{x}}(t)=\Pi_{\Omega}({\bm{x}}_{k}-t\bm{\mathcal{G}}_{k}+t^{\alpha}\bm{d}_{k})$, it follows from part (b) of Lemma 2 that

[TABLE]

When $t\leq\underline{t}:=\min\left\{\left(\frac{1-\sigma}{16C_{5}}\right)^{1/\alpha},\;\frac{2-2\sigma}{\gamma_{h}+16C_{5}\epsilon},\;1\right\}$, it can be verified that

[TABLE]

Combining (3.24) and (3.25), we have

[TABLE]

Therefore, $t_{k}\geq t_{0}\rho^{\lceil\log_{\rho}\underline{t}/t_{0}\rceil}$ for this case.

Taking both cases into consideration, we have proved this lemma.

∎

Now we are ready to prove Theorem 3.

Proof of Theorem 3.

Based on Lemma 4, we have

[TABLE]

Since $\lim_{k\rightarrow\infty}f({\bm{x}}_{k})>-\infty$, we must have

[TABLE]

∎

Remark 5**.**

In Algorithm 1, we can relax the requirement of finding the algebraically smallest eigenpair of $\bm{\mathcal{H}}_{k}$. Instead, it is sufficient to find $(\bar{\lambda},\bar{\bm{v}})$ satisfying

[TABLE]

where $\delta>0$ is a prescribed constant.

3.3 Examples of $\Omega$ and $\Pi_{\Omega}(\cdot)$

For Algorithm 1 to be practical, the projection of a point onto the feasible set $\Omega$ must be affordable. In this subsection, we enumerate a number feasible sets $\Omega$, defined by equality constraints, that are frequently encountered in optimization problems and have computationally tractable projections $\Pi_{\Omega}(\cdot)$:

There is a wide range of applications in which one is interested in solving optimization problems with constraints as listed in Table 1. The paper wen2013feasible provides an extensive list of such problems and references to particular applications. These include eigenvalue and subspace tracking problems arising in signal processing; low-rank matrix optimization problems such as those that arise in SDP relaxation of combinatorial problems (e.g., the max-cut problem described in section 2); $p$-harmonic flows and other problems involving normal preserving constraints such as those that arise in $1$-bit compressive sensing, color image denoising, micromagnetics, liquid crystal theory, and directional diffusion; homogeneous polynomial optimization with spherical constraints arising in tensor eigenvalue problems signal processing, MRI, data training, approximation theory, portfolio selection and computation of the stability number of a graph; sparse principal component analysis, electronic structures computation, etc.

3.4 Discussion: extension to general constrained problems

We briefly discuss possibilities of generalizing the proposed algorithm to optimization problems with inequality constraints, i.e., problems of the form:

[TABLE]

By adding a squared slack variable to each inequality constraint, since $c_{i}({\bm{x}})\leq 0\Leftrightarrow c_{i}({\bm{x}})+z_{i}^{2}=0$, one can transform a problem with inequality constraints to one that has only equality constraints. Moreover, it can be verified that the conditions imposed in Assumption 2 on the constraints $c_{i}({\bm{x}})$, carry over to the transformed constraints. Therefore, algorithm 1 can be applied again. But one caveat of this reformulation is that additional second-order critical points might be introduced. We leave it as future work to investigate better approaches to handling general inequality constraints.

4 Appendix

4.1 Proof of Lemma 1

We need the proposition below in the proof.

Proposition 1**.**

((GHJY15, , Lemma 33)) Assume that Assumption 2 holds and define $R=\sqrt{1/(\sum_{i=1}^{m}L_{c_{i},1}^{2}/\sigma_{0}^{2})}$. For any ${\bm{x}}_{0}\in\Omega$, any $\bm{v}\in\mathbb{R}^{m}$, let ${\bm{x}}_{1}={\bm{x}}_{0}+\bm{v}$ and ${\bm{x}}_{2}={\bm{x}}_{0}+\mathcal{P}_{{\bm{x}}_{0}}(\bm{v})$. Then we have

[TABLE]

Now we are ready to prove Lemma 1.

Proof of Lemma 1.

To prove (3.4), for any $\bm{x}\in\Omega$, let $\bm{y}=\Pi_{\Omega}(\bm{x}+\bm{\delta})$ and consider the Lagrangian $\mathcal{L}(\bm{y},\bm{\lambda}^{*}(\bm{x}))=f(\bm{y})-\sum_{i=1}^{m}\lambda^{*}_{i}(\bm{x})c_{i}(\bm{y})$. Since $\bm{x},\bm{y}\in\Omega$, $\bm{c}(\bm{x})=\bm{c}(\bm{y})=\bm{0}$ and hence $\mathcal{L}(\bm{y},\bm{\lambda}^{*}(\bm{x}))=f(\bm{y})$ and $\mathcal{L}(\bm{x},\bm{\lambda}^{*}(\bm{x}))=f(\bm{x})$. Therefore by Taylor expansion, we have for some $s\in(0,1)$

[TABLE]

Hence using Assumption 2(b), Cauchy-Schwartz and the fact that $s<1$, we have

[TABLE]

Both $\|\bm{\lambda}^{*}(\bm{x})\|$ and $\|\mathcal{G}(\bm{x})\|$ can be bounded by constants. From (3.19), use of the assumptions that $\inf_{\bm{x}\in\Omega}\sigma_{\text{min}}(\nabla\bm{c}({\bm{x}}))\geq\sigma_{0}$, $\sup_{\bm{x}\in\Omega}\|\nabla f(\bm{x})\|\leq\gamma_{f,1}$, $\sup_{\bm{x}\in\Omega}\|\nabla c_{i}(\bm{x})\|\leq\gamma_{c_{i},1}$ and the definition of Frobenius norm,

[TABLE]

While

[TABLE]

To bound $\|\bm{y}-\bm{x}\|$, by definition of $\bm{y}$ and $\Pi_{\Omega}(\cdot)$ we have

[TABLE]

As above $\|\bm{y}-(\bm{x}+\bm{\delta})\|\leq\|\bm{\delta}\|$. However, we also need to bound $\|\bm{y}-(\bm{x}+\bm{\delta})\|$ in terms of $\|\bm{\delta}\|^{2}$ to facilitate our analysis. Specifically, since $\bm{\delta}\in\mathcal{T}_{\Omega}(\bm{x})$ and $\bm{x}\in\Omega$, by Proposition 1, we have

[TABLE]

where $R=\sigma_{0}/\sqrt{\Lambda_{1}}$. Now plugging (4.5), (4.4), (4.3), and (4.2) into (4.1) we obtain the desired result (3.4), that is,

[TABLE]

where $C_{0}=\gamma_{f,1}\{1+\frac{\Gamma_{1}}{\sigma_{0}^{2}}\}\frac{4}{R^{2}}+4(L_{f,1}+\frac{\gamma_{f,1}\sqrt{\Gamma_{1}\Lambda_{1}}}{\sigma_{0}^{2}})$.

To prove (3.5), for every ${\bm{x},\bm{y}}\in{\Omega}$, by the definition of the Lagrangian, we have $f({\bm{x}})=\mathcal{L}({\bm{x}},\bm{\lambda}^{\star}({\bm{x}}))$ and $f(\bm{y})=\mathcal{L}(\bm{y},\bm{\lambda}^{\star}({\bm{x}}))$. Let $\bm{\eta}={\bm{y}}-{\bm{x}}$. Then by Taylor’s theorem, there exists a $t\in(0,1)$ such that

[TABLE]

Then the last term in $(\ref{feasible_expansion_estimate})$ can be bounded as follows

[TABLE]

where C_{1}=\big{\{}\frac{L_{f,2}}{2}+\frac{\gamma_{f,1}}{2\sigma_{0}^{2}}\sqrt{\Gamma_{1}\Lambda_{2}}\big{\}} Combining the above inequality with (4.6), we obtain

[TABLE]

We consider two case: (i) $\bm{x}_{0}+\bm{\delta}\in\Omega$, and (ii) $\bm{x}_{0}+\bm{\delta}\not\in\Omega$.

Case (i): Substituting $\bm{y}=\bm{x}_{0}+\bm{\delta}=\Pi_{\Omega}(\bm{x}_{0}+\bm{\delta})$ and $\bm{x}=\bm{x}_{0}$ into (4.7), and noting that $\bm{\eta}=\bm{y}-\bm{x}=\bm{\delta}$, we have from (4.7) and the facts that $\mathcal{H}(\bm{x}_{0})=\mathcal{P}_{\bm{x}_{0}}^{\top}\nabla_{\bm{x}\bm{x}}^{2}\mathcal{L}(\bm{x}_{0},\bm{\lambda}^{\star}(\bm{x}_{0}))\mathcal{P}_{\bm{x}_{0}}$ and $\bm{\delta}\in\mathcal{T}_{\Omega}(\bm{x}_{0})$, hence that $\mathcal{P}_{\bm{x}_{0}}\bm{\delta}=\bm{\delta}$.

[TABLE]

Case (ii): ${\bm{x}}_{0}+{\bm{\delta}}\not\in\Omega$. Letting ${\bm{y}_{0}}=\Pi_{\Omega}({\bm{x}}_{0}+{\bm{\delta}})$, we have from (4.4) and then from (4.7) that

[TABLE]

Let $\bm{\xi}=(\bm{x}_{0}+\bm{\delta})-\bm{y}_{0}=(\bm{x}_{0}+\bm{\delta})-\Pi_{\Omega}({\bm{x}}_{0}+{\bm{\delta}})$. Then clearly $\bm{y}_{0}-\bm{x}_{0}=\bm{\delta}-\bm{\xi}$, and (4.8) can be rewritten as

[TABLE]

We further note that

[TABLE]

Next, we will show that the last three terms in (4.11) satisfy

[TABLE]

The following are helpful in establishing this. First, from the definition of $\bm{\xi}=(\bm{x}_{0}+\bm{\delta})-\Pi_{\Omega}({\bm{x}}_{0}+{\bm{\delta}})$ and Proposition 1, one has $\|{\bm{\xi}}\|\leq 4\|{\bm{\delta}}\|^{2}/R$ as well as the bound $\|\bm{\xi}\|\leq\|\bm{\delta}\|$. Second, since $\mathcal{G}({\bm{x}}_{0})\in\mathcal{T}_{\Omega}(\bm{x}_{0})$, we have $\mathcal{G}({\bm{x}}_{0})^{\top}{\bm{\xi}}=\mathcal{G}({\bm{x}}_{0})^{\top}\mathcal{P}_{{\bm{x}}_{0}}{\bm{\xi}}$. Third, since $\bm{y}_{0}=\Pi_{\Omega}(\bm{x}_{0}+\bm{\delta})$, $\mathcal{P}_{\bm{y}_{0}}\{(\bm{x}_{0}+\bm{\delta})-\Pi_{\Omega}(\bm{x}_{0}+\bm{\delta})\}=\bm{0}$. Using these observations and the fact that $\|{\bm{y}_{0}}-{\bm{x}}_{0}\|\leq 2\|{\bm{\delta}}\|$, we are ready to bound these three terms one by one,

[TABLE]

where $C_{2}=\{\gamma_{f,2}+\frac{\gamma_{f,1}}{2\sigma_{0}^{2}}\sqrt{\Gamma_{1}\Gamma_{2}}\}\frac{4}{R}$. Applying exactly the same method to bound $\frac{1}{2}{\bm{\xi}}^{\top}\{\nabla^{2}_{\bm{x}\bm{x}}\mathcal{L}({\bm{x}_{0}},\bm{\lambda}^{\star}({\bm{x}_{0}}))\}{\bm{\xi}}$ we have

[TABLE]

where $C_{3}=\{\gamma_{f,2}+\frac{\gamma_{f,1}}{2\sigma_{0}^{2}}\sqrt{\Gamma_{1}\Gamma_{2}}\}\frac{2}{R}$. To bound the last term $\mathcal{G}(\bm{x}_{0})^{\top}\bm{\xi}$, we first note that $\mathcal{G}(\bm{x}_{0})^{\top}\bm{\xi}=\mathcal{G}(\bm{x}_{0})^{\top}\mathcal{P}_{\bm{x}_{0}}\bm{\xi}$ since $\mathcal{G}(\bm{x}_{0})\in\mathcal{T}_{\Omega}(\bm{x}_{0})$ and $\mathcal{P}_{\bm{y}_{0}}\bm{\xi}=\bm{0}$ since $\bm{\xi}=\bm{y}_{0}-\Pi_{\Omega}\bm{y}_{0}$. Therefore

[TABLE]

We first derive a bound for $\left\|\mathcal{P}_{\bm{x}_{0}}\bm{\xi}-\mathcal{P}_{\bm{y}_{0}}\bm{\xi}\right\|$. Note that $\mathcal{P}_{\bm{x}}=\nabla\bm{c}({\bm{x}})\{\nabla\bm{c}({\bm{x}})^{\top}\nabla\bm{c}({\bm{x}})\}^{-1}\nabla\bm{c}({\bm{x}})^{\top}$; thus

[TABLE]

Since

[TABLE]

and

[TABLE]

Moreover, by

[TABLE]

Therefore $\|\nabla\bm{c}({\bm{x}_{0}})\{\nabla\bm{c}({\bm{x}_{0}})^{\top}\nabla\bm{c}({\bm{x}_{0}})\}^{-1}\|\|\nabla\bm{c}(\bm{x}_{0})-\nabla\bm{c}(\bm{y}_{0})\|$ in (4.16) can be bounded by

[TABLE]

We need to further simplify $\|\nabla\bm{c}({\bm{x}_{0}})\{\nabla\bm{c}({\bm{x}_{0}})^{\top}\nabla\bm{c}({\bm{x}_{0}})\}^{-1}-\nabla\bm{c}({\bm{y}_{0}})\{\nabla\bm{c}({\bm{y}_{0}})^{\top}\nabla\bm{c}({\bm{y}_{0}})\}^{-1}\|$ in order to obtain an upper bound for (4.16).

[TABLE]

An upper bound of (4.21) can be obtained by combining (4.18) and the $\sigma_{0}$-LICQ condition, that is,

[TABLE]

To upper bound $\|\{\nabla\bm{c}({\bm{x}_{0}})^{\top}\nabla\bm{c}({\bm{x}_{0}})\}^{-1}-\{\nabla\bm{c}({\bm{y}_{0}})^{\top}\nabla\bm{c}({\bm{y}_{0}})\}^{-1}\|$ in (4.22), we need to utilize the fact that for any invertible matrices $A,B$, $\left\|A^{-1}-B^{-1}\right\|\leq\left\|A^{-1}\right\|\left\|A-B\right\|\left\|B^{-1}\right\|$ and the $\sigma_{0}$-LICQ condition. More specifically,

[TABLE]

Therefore we can bound $\|\nabla\bm{c}({\bm{x}_{0}})\{\nabla\bm{c}({\bm{x}_{0}})^{\top}\nabla\bm{c}({\bm{x}_{0}})\}^{-1}-\nabla\bm{c}({\bm{y}_{0}})\{\nabla\bm{c}({\bm{y}_{0}})^{\top}\nabla\bm{c}({\bm{y}_{0}})\}^{-1}\|$ in (4.17) by

[TABLE]

With the upper bounds for (4.16) and (4.17), we bound $\left\|\mathcal{P}_{\bm{x}_{0}}\bm{\xi}-\mathcal{P}_{\bm{y}_{0}}\bm{\xi}\right\|$ by

[TABLE]

Let us go back to the task of bounding $\mathcal{G}(\bm{x}_{0})^{\top}\bm{\xi}$. By (4.15) and the inequality above, we have

[TABLE]

where $C_{4}=\{\gamma_{f,1}+\frac{\gamma_{f,1}}{\sigma_{0}^{2}}\Gamma_{1}\}(\frac{2\sqrt{\Gamma_{1}\Lambda_{1}}}{\sigma_{0}^{2}}+\frac{2\sqrt{\Gamma_{1}^{3}\Lambda_{1}}}{\sigma_{0}^{4}})\frac{8}{R}$. Finally, we can bound the last three terms in (4.11) by combining (4.13), (4.14, and (4.26), that is,

[TABLE]

Plugging (4.27) into (4.9) yields the final result, we obtain,

[TABLE]

∎

4.2 Riemannian gradient and Riemannian Hessian

Recall that $\Omega=\{{\bm{x}}\in\mathbb{R}^{n}\;|\;c_{i}({\bm{x}})=0,i=1,\ldots,m\}$, $\nabla\bm{c}({\bm{x}})=[\nabla c_{1}({\bm{x}}),\cdots,\nabla c_{m}({\bm{x}})]$ and $\sigma_{\min}({\nabla\bm{c}}({\bm{x}}))$ the minimum singular value of ${\nabla\bm{c}}({\bm{x}})$. Assume that $\inf\{\sigma_{\min}({\nabla\bm{c}}({\bm{x}}))\ |\ {\bm{x}}\in\Omega\}>\alpha$ for some $\alpha>0$. For second-order differentiable functions $f$ and $c_{i}(\bm{x}),i=1,\ldots,m$, let $\mathcal{L}({\bm{x}},\bm{\lambda})=f({\bm{x}})-\sum_{i=1}^{m}\lambda_{i}c_{i}({\bm{x}})$ and $\bm{\lambda}^{\star}({\bm{x}})=\arg\!\min_{\bm{\lambda}}\|\nabla_{\bm{x}}\mathcal{L}({\bm{x}},\bm{\lambda})\|$. Define $\mathcal{G}({\bm{x}})=\nabla_{\bm{x}}\mathcal{L}({\bm{x}},\bm{\lambda}^{\star})$ and $\mathcal{H}({\bm{x}})=\nabla^{2}_{{\bm{x}\bm{x}}}\mathcal{L}({\bm{x}},\bm{\lambda}^{\star})$. Note that $\Omega$ can be considered as a $n-m$ dimensional Riemannian sub-manifold of $\mathbb{R}^{n}$ under the LICQ assumption. We will show that $\mathcal{G}(\cdot)$ and $\mathcal{H}(\cdot)$ are Riemannian gradient and Riemannian Hessian of $f$ over $\Omega$ in the following lemma. A similar argument can be found in absil2013extrinsic and absil2009all

Lemma 5**.**

$\mathcal{G}(\cdot)$* and $\mathcal{H}(\cdot)$ defined in 3.2 and 3.3 are the Riemannian gradient and Riemannian Hessian of $f(\cdot)$ over $\Omega$.*

Proof.

Let $\text{grad}f({\bm{x}})$ denote the Riemannian gradient of $f$ at ${\bm{x}}\in\Omega$. By definition, for every $\eta\in\mathcal{T}_{\Omega}({{\bm{x}}})$, $Df({\bm{x}})[\eta]=\langle\text{grad}f({\bm{x}}),\eta\rangle_{\bm{x}}$ where $Df({\bm{x}})[\eta]$ denotes the directional derivative of $f$ at ${\bm{x}}$ along the direction $\eta$ and $\langle\cdot,\cdot\rangle_{\bm{x}}$ denotes the Riemannian metric on $\mathcal{T}_{\Omega}({{\bm{x}}})$. Since $\Omega$ is an embedded Riemannian sub-manifold of $\mathbb{R}^{n}$, $\langle\cdot,\cdot\rangle_{\bm{x}}$ coincides with the Euclidean inner product. Therefore $\mathcal{G}({\bm{x}})$ is the projection of $\nabla f({\bm{x}})$ onto $\mathcal{T}_{\Omega}({{\bm{x}}})$. Note that for every $\bm{\eta}\in\mathcal{T}_{\Omega}({{\bm{x}}})$, we have ${\nabla\bm{c}}({\bm{x}})\eta=\bm{0}$. Then, we have

[TABLE]

The last equality follows from the definition of directional derivative. Therefore $\langle\text{grad}f({\bm{x}}),\bm{\eta}\rangle_{\bm{x}}=\langle\mathcal{G}({\bm{x}}),\bm{\eta}\rangle_{\bm{x}}$ for every $\bm{\eta}\in\mathcal{T}_{\Omega}({{\bm{x}}})$, that is, $\mathcal{G}({\bm{x}})=\text{grad}f({\bm{x}})$. Next, we will show that $\mathcal{H}({\bm{x}})$ is the Riemannian Hessian of $f$.

Let $\text{Hess}f({\bm{x}})$ denote the Riemannian Hessian of $f$ at ${\bm{x}}\in\Omega$. By definition, for all ${\bm{\xi}},\bm{\eta}\in\mathcal{T}_{\Omega}({\bm{x}})$, $\text{Hess}f({\bm{x}})[{\bm{\xi}},\bm{\eta}]=\langle\overline{\nabla}_{\xi}\text{grad}f({\bm{x}}),\bm{\eta}\rangle_{\bm{x}}$ where $\overline{\nabla}$ denotes the Riemannian connection on $\Omega$. First note that under the LICQ assumption, ${\nabla\bm{c}}({\bm{x}})^{\top}{\nabla\bm{c}}({\bm{x}})$ is invertible for every ${\bm{x}}\in\Omega$ and straight forward calculation gives $\bm{\lambda}^{\star}({\bm{x}})=({\nabla\bm{c}}({\bm{x}})^{\top}{\nabla\bm{c}}({\bm{x}}))^{-1}{\nabla\bm{c}}({\bm{x}})^{\top}\nabla f({\bm{x}})$. By the inverse function theorem, $\bm{\lambda}^{\star}({\bm{x}})$ is differentiable. By Proposition 5.3.2 in AMS09 and $\mathcal{G}({\bm{x}})=\text{grad}f({\bm{x}})$,

[TABLE]

Therefore $\text{Hess}f({\bm{x}})[{\bm{\xi}},\bm{\eta}]=\mathcal{H}({\bm{x}})({\bm{\xi}},\bm{\eta})$. ∎

Acknowledgment

This research was supported in part by NSF Grant CCF1527809.