A gradient descent akin method for inequality constrained optimization

TL;DR

This paper introduces a simplified gradient descent-like method for inequality constrained optimization, analyzing its convergence and demonstrating its effectiveness through numerical examples including shape optimization.

Contribution

It presents a new computationally efficient first-order method derived from a previous search direction approach, with a rigorous mathematical analysis of its convergence properties.

Findings

The method converges to local solutions under certain conditions.

Numerical tests show the method's effectiveness on standard and shape optimization problems.

The approach simplifies previous algorithms while maintaining convergence guarantees.

Abstract

We propose a first-order method for solving inequality constrained optimization problems. The method is derived from our previous work [12], a modified search direction method (MSDM) that applies the singular-value decomposition of normalized gradients. In this work, we simplify its computational framework to a "gradient descent akin" method, i.e., the search direction is computed using a linear combination of the negative and normalized objective and constraint gradient. The main focus of this work is to provide a mathematical aspect to the method. We analyze the global behavior and convergence of the method using a dynamical systems approach. We then prove that the resulting trajectories find local solutions by asymptotically converging to the central path(s) for the logarithmic barrier interior-point method under the so-called relative convex condition. Numerical examples are reported, which include both common test examples and applications in shape optimization.

Figures

Tables

Table 1. Table 1 : Results with the parameter ζ = 0.98 𝜁 0.98 \zeta=0.98 and tuned step sizes

Prob.	step size	iters.	f(${𝐱}^{\star }$)	f(${𝐱}_{\zeta }$)	abs. error
G01	0.002	2362	-15	-14.7215	1.86e-2
G04	0.2	136	-3.0665e+4	-3.0657e+4	2.85e-4
G06	0.002	4826	-6.9618e+3	-6.8371e+3	1.79e-2
G07	0.0027	3009	24.3062	24.7876	1.98e-2
G08	0.01	66	-9.5825e-2	-9.5063e-2	0.80e-2
G09	0.05	120	6.8063e+2	6.9238e+2	1.73e-2
G10	0.35	5319	7.0492e+3	7.1898e+3	1.99e-2
G18	0.01	257	-0.8660	-0.8546	1.32e-2
G19	0.05	294	32.6556	2.7120e+2	7.30
G24	0.02	268	-5.5080	-5.4147	1.69e-2

Equations

$$\begin{array}{rl}& \text{minimize}f(x),\\ & \text{subject to}{g}_i(x)\le 0,i=1,...,m,\end{array}$$

$$\underset{x\to \infty }{\lim }f(x)=+\infty ;$$

$$\frac{dx}{dt}=-\frac{\nabla f}{|\nabla f|}-\zeta \frac{\nabla g}{|\nabla g|},\zeta \in [0,1),$$

$$\Phi (x)=-\sum _{i=1}^m\log (-g_i(x)),i=1,...,m,$$

$$\frac{dx}{dt}=\left\{\begin{split}&-\frac{\nabla f}{|\nabla f|}-\zeta\frac{\nabla\Phi}{|\nabla\Phi|},~{}~{}~{}~{}~{}&\text{if}~{}~{}\nabla\Phi\neq 0;\\ &-\nabla f,&\text{if}~{}~{}\nabla\Phi=0.\end{split}\right.$$

$$\frac{dx}{dt}=-\nabla f.$$

$$f(x+d)\approx f(x)+\nabla f(x)d.$$

$$\begin{array}{rl}& \text{minimize}\nabla f(x)d,\\ & \text{subject to}|d|=1.\end{array}$$

$$d=-\frac{\nabla f}{|\nabla f|},$$

$$\text{minimize}f(x)+\eta \Phi (x),$$

$$0=\nabla f(x^{\ast })+\eta \nabla \Phi (x^{\ast }),x^{\ast }\in {\Omega }_{-},$$

$$\frac{\nabla f}{|\nabla f|}+\frac{\nabla g}{|\nabla g|}=0.$$

$$\frac{\nabla f}{|\nabla f|}+\frac{\nabla \Phi }{|\nabla \Phi |}=0.$$

$$\frac{dx}{dt}=-\frac{\nabla f}{|\nabla f|}(1-\zeta )+\left(-\frac{\nabla f}{|\nabla f|}-\frac{\nabla \Phi }{|\nabla \Phi |}\right)\zeta ,\zeta \in [0,1).$$

$$\mathbf{m}=\left(\begin{array}{c}\frac{\nabla f}{|\nabla f|}\\ \frac{\nabla g}{|\nabla g|}\end{array}\right).$$

$$\begin{array}{r}df=\nabla fdx,\\ dg=\nabla gdx.\end{array}$$

$$\left(\begin{array}{c}\frac{df}{|\nabla f|}\\ \frac{dg}{|\nabla g|}\end{array}\right)=\mathbf{m}dx.$$

$$\mathbf{m}=\mathbf{U\Sigma }{\mathbf{V}}^T=\sum _{i=1}^{min(2,n)}{\sigma }_i{\mathbf{u}}_i{\mathbf{v}}_i^T.$$

$$-\frac{\nabla f}{|\nabla f|}=\cos {\alpha }_1{\mathbf{v}}_1^T+\cos {\alpha }_2{\mathbf{v}}_2^T,$$

$${\mathbf{s}}_c=\cos {\alpha }_1{\mathbf{v}}_1^T+c\cdot \cos {\alpha }_2{\mathbf{v}}_2^T,$$

$$\begin{array}{rl}{\mathbf{v}}_1^T& =\frac{1}{\sqrt{2-2\cos \theta }}\left(-\frac{\nabla f}{|\nabla f|}+\frac{\nabla g}{|\nabla g|}\right),\\ {\mathbf{v}}_2^T& =\frac{1}{\sqrt{2+2\cos \theta }}\left(-\frac{\nabla f}{|\nabla f|}-\frac{\nabla g}{|\nabla g|}\right),\end{array}$$

$$\begin{array}{r}\cos {\alpha }_1=-<\frac{\nabla f}{|\nabla f|},{\mathbf{v}}_1^T>=\frac{\sqrt{1-\cos \theta }}{\sqrt{2}},\\ \cos {\alpha }_2=-<\frac{\nabla f}{|\nabla f|},{\mathbf{v}}_2^T>=\frac{\sqrt{1+\cos \theta }}{\sqrt{2}}.\end{array}$$

$${\mathbf{s}}_c=-\frac{\nabla f}{|\nabla f|}-\frac{(c-1)}{2}\left(\frac{\nabla f}{|\nabla f|}+\frac{\nabla g}{|\nabla g|}\right).$$

$${\mathbf{s}}_{\zeta }=-\frac{\nabla f}{|\nabla f|}-\zeta \frac{\nabla g}{|\nabla g|},$$

$$\begin{array}{rl}& \text{minimize}f(x_1,x_2)=\frac{1}{2}(x_1^2+x_2^2),\\ & \text{subject to}g(x_1,x_2)=-x_2+10\le 0.\end{array}$$

$${\mathbf{s}}_{\zeta }=\frac{-1}{\sqrt{x_1^2+x_2^2}}\left\{x_1,x_2-\zeta \sqrt{x_1^2+x_2^2}\right\}.$$

$$x_2+\sqrt{x_1^2+x_2^2}=2{\bar{x}}_2{\left|\frac{x_1}{x_1^0}\right|}^{1-\zeta }.$$

$$(x_{2,\zeta }{)}^{\zeta }=\frac{2}{1+\zeta }\frac{{\bar{x}}_2}{|x_1^0{|}^{1-\zeta }}{\left(\frac{\sqrt{1-{\zeta }^2}}{\zeta }\right)}^{1-\zeta },$$

$$|x_{1,\zeta }|=\frac{{\bar{x}}_2}{\zeta }\sqrt{1-{\zeta }^2}.$$

$$x_1^2=4{\bar{x}}_2^2-4x_2{\bar{x}}_2,x_1\in (0,x_1^0)(\text{mbox}or(x_1^0,0),)$$

Full text

1 Introduction

In this paper, we propose a first-order method for solving inequality constrained optimization problems. The problem of interest is

[TABLE]

where $f,g_{1},...,g_{m}:\mathbb{R}^{n}\rightarrow\mathbb{R}$ are twice differentiable.

The present method is derived from our previous work Chen , a modified search direction method (MSDM) for inequality constrained optimization problems. MSDM computes a descent direction of the objective function using the singular-value decomposition that exploits the normalized gradient information of a constrained optimization problem. Numerical experiments show that the MSDM finds local solutions by traversing along the central path(s) for the logarithmic barrier method. However, there is a lack of the mathematical theory of the method. The intrinsic optimization parameter in MSDM has remained heuristic. In this work, we simplify the computational framework of MSDM to a “gradient descent akin” method, i.e., we compute the search direction using a linear combination of the negative and normalized objective and constraint gradient. We analyze the global behavior and convergence using a dynamical system perspective. We prove that the resulting trajectories find local solutions by asymptotically converging to the central path of the logarithmic barrier interior-point method. We should note that in this paper, we do not address the design and analysis of a practical implementation of the method, which is currently an active area of research.

1.1 Main results

For problem (1), we make the following assumptions:

(A1) Coercive condition for the objective function $f(x)$,

[TABLE]

(A2) $\nabla f(x)\neq 0$ in the feasible set $\Omega=\{x:g_{i}(x)\leq 0,~{}i=1,...,m\}$;

(A3) $\nabla g_{i}(x)\neq 0,~{}i=1,...,m$ in the feasible set $\Omega$;

(A4) $f$ and $g_{i}$ are twice continuously differentiable functions.

Under assumptions (A1)$-$(A4), we propose a dynamical system for solving the problem (1) that includes a single inequality constraint $g(x)$,

[TABLE]

where “$|\cdot|$” denotes the Euclidean norm. $\nabla f$ and $\nabla g$ are the gradient row vectors of the objective function and constraint function, respectively.

We show analysis of the global behavior of the method. Among others, we propose an error measure $\epsilon$ for the optimization solution, and show that the time in which the present system finds a local solution is $\mathcal{O}(1/\epsilon)$, which is the ergodic convergence rate for first-order methods. We prove global and local convergence of the method. The present system results in a trajectory, which, under the so-called relative convex condition, converges to the central path of the interior-point methods as $\zeta\rightarrow 1^{-}$.

Provided the well-known logarithmic barrier function

[TABLE]

where $\log(\cdot)$ denotes the natural logarithm, we propose a generalization of the system (2) for the multiple constrained problem with $\zeta\in[0,1)$ as

[TABLE]

Compared to our previous work Chen , we simplified the framework to compute the search direction so that the implementation effort and computational cost are reduced significantly. We report computational experiences with various test examples that include both common benchmark problems and applications in shape optimization. Finally, we provide a large-scale real world shape optimization, which, to the best of our knowledge, is unlike any other case previously presented in the literature.

1.2 Organization of paper

In section 2, we first review related works with the focus on theoretical aspects. In section 3, we derive the method for single inequality constrained problems. Preliminary studies are then presented in section 4. Our main results are presented in section 5 - 7. We report numerical examples in section 8, and finally give a conclusion in section 9.

2 Related works

We review related works to share our view of the present method from different well-established perspectives. Here, the emphasis is on the theoretical connections/differences as the main focus of this paper is to provide a theoretical foundation for the present method. We also briefly review the shape optimization as it has motivated the development of the method in this work.

2.1 Gradient descent method

The gradient descent method, originally proposed by Cauchy, is a first-order method for unconstrained optimization that uses the negative gradient of the objective function for the design update. In a dynamical system, it writes

[TABLE]

The direction of the gradient descent is the steepest descent direction in the Euclidean norm. A first-order Taylor expansion at current iterate $x$ of the objective function reads

[TABLE]

The steepest descent direction $d$ is found by the optimization problem

[TABLE]

From the Cauchy-Schwarz inequality, we obtain

[TABLE]

which is the negative and normalized objective gradient. Comparing (5) and (2), which linearly combines the negative and normalized objective and constraint gradient, we consider the present system: a gradient descent akin method for inequality constrained optimization problems. We especially note that the present method results in trajectories that are homotopic with the gradient descent trajectory (see Remark 1).

2.2 Dynamical systems approaches

Dynamical systems approaches have been used to study optimization methods in many works of literature. Extensive studies on the connections between interior-point flows with linear programming methods can be found in (helmke2012optimization, , Chapter 4) and the references therein. For quadratic programming problems, dorr2012smooth proposes a dynamical system, which results in trajectories that converge to the saddle point of the associated Lagrangian function. su2014differential studies the celebrated Nesterov’s accelerated gradient method using a dynamical system as the analysis tool. In lessard2016analysis , a framework based on dynamical systems is proposed to analyze and design first-order unconstrained optimization methods. In this work, we study the behavior of the present method using the dynamical systems approach. The ODE interpretation of the method allows us to give a rigorous analysis of its global and local convergence.

In some literature, optimization methods that use dynamical systems are called trajectory methods. These methods construct optimization paths in a way so that one or all solutions to the optimization problem are a priori known to lie on these paths diener1995trajectory . Typically, these optimization paths are solution trajectories to ODE of first or second-order. Trajectory methods are mainly studied for unconstrained optimizations for finding local solutions behrman1998efficient botsaris1978differential , and global solutions griewank1981generalized snyman1987multi . Studies for constrained optimization are, however, very limited, see ali2018trajectory and the references therein. In this work, we propose a new dynamical system for inequality constrained optimization problems. The resulting trajectories find local solutions by the limiting behavior as $\zeta\rightarrow 1^{-}$.

2.3 Interior-point methods

Interior-point methods (IPMs), which are based on the Newton method, are among the most competitive methods for constrained optimization problems. The signature of IPM is the existence of continuously parameterized families of approximate solutions that converge to the exact solution asymptotically Forsgren . IPMs find a wide variety of applications of convex and nonconvex optimizations in broad fields. There are a vast amount of excellent works that have been devoted to IPM. A comprehensive review of this class of methods is certainly beyond the scope of the present paper, however we refer the interested reader to Forsgren Potra , more recently, in gondzio2012interior , and many other excellent optimization books.

The present method results in trajectories that asymptotically converge to the central path for a particular IPM, the logarithmic barrier method fiacco1990nonlinear . We introduce its connections/differences to the present method in the next. Consider a convex optimization problem of the form (1), we start with its approximated unconstrained problem using the logarithmic barrier function $\Phi$,

[TABLE]

where the barrier parameter $\eta$ is a positive parameter. As $\eta\rightarrow 0$, the solution of the approximated problem converges to the original one. The central path is characterized by the set of points that satisfy the necessary and sufficient conditions boyd :

[TABLE]

where $\Omega_{-}=\{x:g_{i}(x)<0,~{}i=1,...,m\}$. The conditions (7) are interpreted as a modified KKT system in the literature boyd byrd2000trust Forsgren . The barrier method finds an approximated solution for the original problem by 1) iteratively decreasing the barrier parameter $\eta$, and 2) in each iteration, solving the subproblem (the modified KKT system) defined by (7) using the Newton method. Therefore, the barrier parameter $\eta$ can be considered a central path parameter.

In the present method, we do not parameterize the central path. Instead, we normalized the objective and constraint gradients. For optimization problems with a single inequality constraint, we propose the normalized central path condition as

[TABLE]

The condition (8) is used in section 5 and 6 for the analysis. A generalization of the normalized central path condition from single inequality to multiple inequalities is introduced by the use of the logarithmic barrier function (see section 7) as

[TABLE]

Compared with the barrier method, the path parameter $\eta$ has vanished. The condition (9) characterizes the central path, which differs from (7), which instead characterizes a point on the central path.

To illustrate the behavior of the present method in a very rough way, we rewrite the first ODE of system (4) as

[TABLE]

While $-\frac{\nabla f}{|\nabla f|}$ is the steepest descent direction of the objective function, the term $\left(-\frac{\nabla f}{|\nabla f|}-\frac{\nabla\Phi}{|\nabla\Phi|}\right)$ contributes to the centering behavior (referred to Theorem 3). There are three major differences compared to the barrier method:

While $\eta$ is a path parameter that controls the asymptotic convergence progress for the barrier method, $\zeta$ is a homotopy parameter that determines the shape of the optimization trajectory in the present method;

• 2)

There is no subproblem (modified KKT system) defined and solved in the present method. Instead, we solve directly for the trajectory of the proposed dynamical system. The optimization solutions are known a priori to lie on the resulting trajectory;

• 3)

The local and global convergence behavior of the present method are markedly different from the barrier method.

The asymptotic convergence for the barrier method along the central path has been extensively studied in fiacco1990nonlinear . For general inequality-constrained problems (1) under mild assumptions (A1)-(A4), to force global convergence, IPM usually implements a merit function associated line-search method or trust region framework byrd2000trust vanderbei1999interior wachter2006implementation . To provide a global and local convergence theory for the present method, we prove the following:

The resulting trajectory converges to a critical point of the objective function as $\zeta\in[0,1)$ (see Theorem 35);

• -

The trajectory finds a KKT solution upon reaching the boundary of the feasible set $\Omega$ as $\zeta\rightarrow 1^{-}$, provided that there is no critical point of the objective function in $\Omega$ (see Theorem 1, 5, and Remark 4);

• -

As $\zeta\rightarrow 1^{-}$, the second-order optimality conditions are automatically satisfied at KKT solutions (see section 6);

• -

The trajectory is able to switch between central paths to remain a descent direction of the objective function without using an additional framework (referred to Theorem 8 and Lemma 5.2).

From a more abstract point of view, the difference between the present method and the barrier method may be analogous to the difference between the gradient descent method and the Newton-based method for unconstrained minimizations. While the former seek critical points of the objective function or KKT solutions that satisfy the respective second-order optimality conditions using the negative gradient information, the latter seek the respective solutions using the Newton method.

2.4 Feasible direction methods

The method of feasible directions (MFD) dates back to the 1960’s by the work of Zoutendijk zoutendijk1960methods and has enjoyed fruitful developments for decades. MFDs have been especially popular in the engineering community because of the importance of ending up with a design that satisfies the hard specifications expressed by a set of inequalities chen2000methods . The general idea behind the MFD is to move from one feasible design to an improved feasible design iteratively so that a local solution can be found Arora . In the present method, remaining feasibility is not a mandatory mechanism. The idea behind the method design is to find a search direction that approaches the central path while maintaining a descent direction of the objective function. Due to this major difference, we do not categorize our method as a method of feasible directions.

In the present method, we compute a search direction that uses normalized gradients. In stander1993new , the authors also present a feasible direction method that applies normalized gradients. Their work is then continued and further developed in de1994feasible and stander1995robustness . In these works, the active-set strategy is used. The common idea is to formulate a linear system under given input criteria on a chosen working set of (active) constraints, and a feasible descent search direction is obtained by solving the linear system. In the present method, we use a barrier function based formulation to treat multiple inequality constraints. The search direction is computed as a linear combination of the negative and normalized objective and barrier gradient.

2.5 Shape optimization

As a subset of structural optimization, shape optimization is characterized by a very large or even infinite number of design variables that describe the varying boundary in the optimization process. Introductions to shape optimization are given in haslinger2003introduction sokolowski1992introduction . Shape optimization is distinct from another well-known problem in structural engineering: topology optimization bendsoe2013topology (sometimes referred to as the homogenization method allaire2012shape ). The main difference is that the topology optimization method removes smoothness and topological constraints in shape optimization, which results in different optimization formulations. Many topology optimization problems can be formulated in an (equivalent) convex optimization problem, while shape optimizations are typically nonlinear, and often nonconvex hoppe2007adaptive . This difference partially contributes to the fact that there are successful implementations of IPM for large-scale topology optimization problems jarre1998optimal kocvara2016primal maar2000interior , but only a few works have presented a shape optimization that uses an IPM as the optimizer antil2007path herskovits2000shape . In the latter works, the size of the shape optimization problem is only moderate so that the power of IPM is not fully exploited. One of the most successful methods for nonlinear topology optimization is the method of moving asymptotes (MMA) that was introduced by Svanberg in 1987 svanberg1987method . In each iteration, MMA generates and solves an approximated convex problem related to the original one. For shape optimization, however, there is as yet no literature that discusses a large-scale problem using MMA.

A notable difficulty for shape optimization is the computation of shape Hessians, which are complex objects even for moderate problems. Analysis of aerodynamic optimization in arian1999analysis shows that shape Hessians are ill-conditioned for three-dimensional problems. Recently, several works compute approximated shape Hessians and use a Newton-based method for the design optimization schillings2011efficient schmidt2013three . For some disciplines, such as computational fluid dynamics or transient coupled problems, even the computation of shape gradient can be a challenge. See for example albring2016efficient korelc2009automation reuther1999constrained , which are actively undergoing investigation.

In general, large-scale shape optimization is mainly performed using gradient descent type method so far schulz2015 . In engineering practice, a large number of constraints may be considered. The lack of literature in this regard has motivated our development of MSDM in Chen . In the present work, we further simplify the computational framework of MSDM to a gradient descent akin method. We provide a mathematical basis for the method’s optimization behavior. As a result, the implementation effort and computational cost are reduced significantly. It opens the possibility of shape optimization to a wider range of applications.

3 Deriving the search direction for single inequality constrained optimizations

In this section, we show the consistent derivation of the system (2) from our previous work Chen considering a single inequality constraint. First, we review the basic ideas of the modified search direction method and then show the derivation.

In MSDM, at each iteration, we construct a sensitivity matrix

[TABLE]

The change in the objective function $df$ and the constraint function $dg$ resulting from an arbitrary design change $dx\in\mathbb{R}^{n}$ reads

[TABLE]

A perspective from an input-output system established by the sensitivity matrix $\mathbf{m}$ gives

[TABLE]

Applying singular-value decomposition to the sensitivity matrix $\mathbf{m}$,

[TABLE]

Thus, an orthonormal bases set $\mathbf{v}_{i},~{}i=1,2$ is obtained. Each $\mathbf{v}_{i}$ can be used as a base search direction for the design update, and $\mathbf{v}_{1}$ and $\mathbf{v}_{2}$ are defined as follows:

$\mathbf{v}_{1}$: by taking $\delta\mathbf{v}_{1}$ as the design change, we obtain a change in objective as well as in constraint function $[\frac{df}{|\nabla f|},\frac{dg}{|\nabla g|}]^{T}=\sigma_{1}\delta\mathbf{u}_{1}$, which is a decrease in the objective function and an increase in the constraint function.

• -

$\mathbf{v}_{2}$: by taking $\delta\mathbf{v}_{2}$ as the design change, we obtain a change in the objective as well as in the constraint function $[\frac{df}{|\nabla f|},\frac{dg}{|\nabla g|}]^{T}=\sigma_{2}\delta\mathbf{u}_{2}$, which is a decrease in the objective function and a decrease in the constraint function.

It is worth mentioning that the design vector $\delta\mathbf{v}_{1}$ provides a similar result as the filter approach presented in fletcher2002nonlinear , which tries to minimize the so-called bi-objective optimization problem with two goals of minimizing the objective function $f$ and the constraint violation $|g|$ (with the difference that $g$ takes different signs in both cases).

With $\mathbf{v}_{1}$ and $\mathbf{v}_{2}$, we can then rewrite the normalized steepest descent direction $-\frac{\nabla f}{|\nabla f|}$ as

[TABLE]

where $\alpha_{1}$ is the angle between $\mathbf{v}^{T}_{1}$ and $-\frac{\nabla f}{|\nabla f|}$, and $\alpha_{2}$ is the angle between $\mathbf{v}^{T}_{2}$ and $-\frac{\nabla f}{|\nabla f|}$. The modified search direction proposed in Chen reads

[TABLE]

where $c\geq 1$ is introduced to enlarge the contribution of the design mode $\mathbf{v}_{2}$.

In the following, we show the derivation of (2) from (16). According to SVD and the definition of $\mathbf{v}_{1}$ and $\mathbf{v}_{2}$, we have

[TABLE]

where $\theta$ is the angle between the objective function gradient $\nabla f$ and the constraint function gradient $\nabla g$. With $\theta$ we also have

[TABLE]

Inserting (17), (18) into (16) we have

[TABLE]

As we are mainly interested in the direction of the vector field $\mathbf{s}_{c}$, we can rewrite it as

[TABLE]

with $\zeta=\frac{c-1}{c+1}$. With $c\in[1,+\infty)$ we have $\zeta\in[0,1)$ and thus we get the present dynamical system (2).

4 Preliminary studies and intuition

In this section, we demonstrate the present method on two simple examples and conjecture the behavior of the resulting optimization trajectory with intuition.

4.1 An analytical 2D optimization example

We first show a 2D optimization problem and solve it analytically. The optimization problem reads,

[TABLE]

The present search direction field $\mathbf{s}_{\zeta}$ reads

[TABLE]

We define an initial design as $(x_{1}^{0},x_{2}^{0})$. Let $\bar{x}_{2}=\frac{1}{2}\left(x_{2}^{0}+\sqrt{(x_{1}^{0})^{2}+(x_{2}^{0})^{2}}\right)$, then, the trajectory $\Gamma^{\zeta}$ of the present search direction field $\mathbf{s}_{\zeta}$ is

[TABLE]

Let $(x_{1,\zeta},x_{2,\zeta})$ be a point on the trajectory $\Gamma^{\zeta}$ with a maximal $x_{2}$ component, then we have

[TABLE]

and

[TABLE]

Let $\zeta\rightarrow 1^{-}$, then, $x_{1,\zeta}\rightarrow 0$, $x_{2,\zeta}\rightarrow\bar{x}_{2}$, the trajectory $\Gamma^{\zeta}$ will converge to the curve $\Gamma$ that is a union of the parabola

[TABLE]

and the interval $(0,\bar{x}_{2})$ on $x_{2}$-axis as is shown in figure 1. As $\zeta\rightarrow 1^{-}$, a part of the trajectory $\Gamma^{\zeta}$ converges to the central path. For $\zeta\in[0,1)$, the resulting trajectories are homotopic relative to their endpoints, which are the initial design and the critical point of the objective function.

By $\bar{x}_{2}=\frac{1}{2}\left(x_{2}^{0}+\sqrt{(x_{1}^{0})^{2}+(x_{2}^{0})^{2}}\right)$, we have $\bar{x}_{2}\geq x_{2}^{0}$. This means for any feasible initial design, as $\zeta\rightarrow 1^{-}$, the resulting optimization trajectory always reaches first a close neighborhood of the central path (at point $(0,\bar{x}_{2})$) , where it is at a larger distance to the boundary of the feasible set compared to the initial design. It then follows the central path and reaches the optimal solution.

4.2 A nonconvex optimization example

We show a second example that includes a nonconvex constraint,

[TABLE]

In figure 2, we plot the optimization trajectory with $\zeta=0.9999$ together with a few depicted contours of both the objective function and the constraint function. We plot three points A, B, and C on the central path, where the plotted contours of the objective and constraint function are tangent to each other. We choose an initial design $\mathbf{x}^{0}$ that is close to the point A. It can be observed: instead of heading to the left side of the central path, the optimization trajectory finds its way to the right side. It then follows the central path, but leaves at point C and reaches the other central path. It eventually finds the optimal solution by following this central path. Along the optimization trajectory, the objective function value decreases steadily. The question is now, why does the optimization trajectory choose one side (point B side) of the same central path over another (point A side)? The answer may lie in the difference in the curvatures of the contours of the objective and constraint function between point A and B. We observe the following fact:

Let $\kappa_{f}$ and $\kappa_{g}$ be the curvature of the contours of objective and constraint function at the central path, respectively, both with the normal vector $\frac{\nabla f}{|\nabla f|}$, then

at point A: $\kappa_{f}-\kappa_{g}>0$;

• -

at point B: $\kappa_{f}-\kappa_{g}<0$.

Based on this observation, we conjecture the behavior of the optimization trajectory: As $\zeta\rightarrow 1^{-}$, the optimization trajectory is able to approach and follow a central path, on which the central point satisfies the condition

[TABLE]

We call condition (28) the relative convex condition. This behavior can also be used to explain why the optimization trajectory leaves the central path at point C, where $\kappa_{f}=\kappa_{g}$, and heads towards another central path. In section 6, we give a local convergence analysis starting with this conjecture.

5 Global behavior and convergence under coercive condition

We show analysis for the global behavior for the present dynamical system (2) under assumptions (A1)$-$(A4). To this purpose, we partly use the geometric analysis in jost2011riemannian .

Consider a $\mu-$neighborhood of a central path with $\mu\in[0,1]$,

[TABLE]

Obviously, $\Theta_{\mu}$ shrinks to the central path as $\mu\rightarrow 1^{-}.$

Recall the present search direction field $\textbf{s}_{\zeta}$,

[TABLE]

Let $x(t;\zeta,x_{0})$ be the solution of the system:

[TABLE]

with $x_{0}$ as the initial design in the feasible set and $T_{\zeta,x_{0}}$ as the maximal existence interval of $x(t;\zeta,x_{0})$ in the whole $\mathbb{R}^{n}$. We show some properties of the trajectory $x(t;\zeta,x_{0})$.

Lemma 5.1**.**

Let $f$ and $g$ satisfy assumptions (A1)(A3) and (A4), then

(i) The trajectory $x(t;\zeta,x_{0})$ always stays within a bounded domain, i.e.,

[TABLE]

(ii) The constraint function $g$ decreases along the trajectory $x(t;\zeta,x_{0})$ out of the cone neighborhood $\Theta_{\zeta}$, and increases in $\Theta_{\zeta}$.

Proof.

Based on the deformation of the objective function $f$ and constraint function $g$ along the trajectory $x(t;\zeta,x_{0})$:

[TABLE]

We get a proof directly. ∎

Theorem 1**.**

Suppose that assumptions (A1) - (A4) hold. Then

(i) the trajectory $x(t;\zeta,x_{0})$ must go out of the feasible set with $\zeta\in[0,1)$;

(ii) the minimum time in which the trajectory reaches the boundary of the feasible set is at most $\frac{C}{1-\zeta}$ with $C$ independent of $\zeta$.

Proof.

We prove (i) by contradiction. Suppose that $x(t;\zeta,x_{0})$ stays in the feasible set $\Omega$ for any $t<T_{\zeta,x_{0}}$, then the vector field $\textbf{s}_{\zeta}$ keeps $C^{1}$ continuous in a neighborhood of trajectory $x(t;\zeta,x_{0})|_{[0,T_{\zeta,x_{0}})}$. Lemma 1 ensures that

[TABLE]

with some positive numbers $A$ and $B$. On the other hand, Picard’s existence theorem implies that $T_{\zeta,x_{0}}=+\infty.$ The integral of the first formula of (30) shows

[TABLE]

This is to say

[TABLE]

Notice that

[TABLE]

So

[TABLE]

By assumption (A4) and Lemma 5.1, there is a constant $l$ so that

[TABLE]

Hence, we have a Lipschitz continuity for $|\nabla f(x(t;\zeta,x_{0}))|$:

[TABLE]

Now, we claim

[TABLE]

Otherwise, there is a sequence of $t_{j}\rightarrow+\infty$ and a positive constant $b$ so that

[TABLE]

Choosing $\delta=\frac{b}{2l}$, then

[TABLE]

for any $|t_{j}-t|\leq\delta.$ Therefore,

[TABLE]

This is a contradiction to (32). Hence

[TABLE]

Notice that this convergence contradicts (31). This completes the proof for (i).

To prove (ii), let $T_{\zeta,x_{0}}^{\sharp}$ be the time in which the trajectory first reaches the boundary of feasible set, then (31) holds for $0<t<T_{\zeta,x_{0}}^{\sharp}$. So

[TABLE]

Hence

[TABLE]

which implies

[TABLE]

so that $T_{\zeta,x_{0}}^{\sharp}\leq\frac{C}{1-\zeta}$ holds. This ends the proof. ∎

Theorem 2**.**

Suppose that assumptions (A1),(A3), (A4) hold and $\nabla g\neq 0$ in the whole $\mathbb{R}^{n}$. The trajectory $x(t;\zeta,x_{0})$ converges to a connected subset of critical points of the objective function by $\zeta\in[0,1)$. Especially, if the critical points of objective function is isolated, then

[TABLE]

Proof.

First notice that the system (29) is not well defined when $\nabla f=0$. The trajectory $x(t;\zeta,x_{0})$ will terminate at these points upon finding them. Therefore, the maximal interval $T_{\zeta,x_{0}}$ may be finite. To overcome this difficulty, we use an equivalent Y$-$system:

[TABLE]

The system (36) has the same orbit as (29) but different parameterization. It has the solution $y(\tau;\zeta,x_{0})$ with infinite existence interval in the whole $\mathbb{R}^{n}$ by Picard’s existence theorem.

For $\zeta\in[0,1)$, consider an integral along $y(\tau;\zeta,x_{0})$,

[TABLE]

Lemma 5.1 ensures

[TABLE]

with some positive number $M$ independent of $T$. Hence

[TABLE]

With $1+\zeta\cos\theta\geq 1-\zeta$, a similar method to the proof for Theorem 1(i) shows that

[TABLE]

This proves that $y(\tau,\zeta,x_{0})$ approaches the connected subset of the critical points. An isolated condition makes sure that

[TABLE]

for some critical point $x_{c}$ of the objective function. Notice that along the trajectory,

[TABLE]

hence

[TABLE]

∎

Remark 1*.*

The global behavior of the system (29) looks like a gradient flow. Obviously, as $\zeta=0$, the system reduces to the normalized and negative gradient flow of the objective function. Especially, the present trajectory is homotopic with the gradient descent trajectory for $\zeta\in(0,1)$ by Theorem 35 and the continuous dependence of the solutions to differential equations on parameters.

Remark 2*.*

The analytical trajectory may find a critical point of constraint $g$ without condition (A3). In practical implementations, we suggest using the gradient descent $-\nabla f$ to escape these critical points.

Lemma 5.2**.**

Under assumptions (A1)$-$(A4), the trajectory $x(t;\zeta=1,x_{0})$ converges, as $t\rightarrow+\infty$, into any $\mu-$neighborhood $\Theta_{\mu}$ of the central path $L$ with $\mu\in[0,1)$, given that the initial design $x_{0}\in\Omega$.

Proof.

Along the trajectory $x(t;1,x_{0})$, we derive the constraint function $g$ with respect to $t$ and have

[TABLE]

By Picard’s existence theorem, the vector field $\textbf{s}_{1}$ keeps $C^{1}$ continuity in a neighborhood of the trajectory $x(t;1,x_{0})$ and $T_{1,x_{0}}=+\infty$. Integrating (37) gives

[TABLE]

Based on Lemma 5.1 and the assumption (A4), we know that the integral function $|\nabla g|(1+\cos\theta)$ is continuously differentiable with respect to $t$. Same procedure as in the proof for Theorem 1 shows that

[TABLE]

resulting with (A3) that $(1+\cos\theta)\rightarrow 0$ as $t\rightarrow+\infty$. In other word, trajectory $x(t;1,x_{0})$ comes into any $\mu-$neighborhood $\Theta_{\mu}$ of the central path $L$, $\mu\in[0,1)$. ∎

Theorem 3**.**

Under assumptions (A1)$-$(A4), the trajectory $x(t;1,x_{0})$ converges to a point $\widehat{x}$ on central path $L$ generally, i.e.,

[TABLE]

provided $x_{0}\in\Omega$.

Proof.

According to (39) we have

[TABLE]

Choose a sequence $t_{j}\rightarrow+\infty$ as $j\rightarrow+\infty$ with

[TABLE]

for some point $\widehat{x}$. Hence

[TABLE]

and

[TABLE]

This says that the point $\widehat{x}\in L$, and the values $f(\widehat{x}),g(\widehat{x})$ depend only on $x_{0}$. In general, the point with these conditions is unique111If the central path and the contour of the constraint function $g$ intersect transversally, then the point $\widehat{x}$ is unique. By Sard’s theorem, transversal intersection occurs with the probability of 1.. Therefore

[TABLE]

∎

Theorem 4**.**

Under assumptions (A1)$-$(A4), the trajectory $x(t;\zeta,x_{0})$ is convergent to the trajectory $x(t;1,x_{0})$ uniformly for $t\in(0,\alpha\log\frac{1}{1-\zeta})$ for some $\alpha>0$.

Proof.

By system (29), we have the integral equation:

[TABLE]

So

[TABLE]

In the bounded domain $\{x\in\mathbb{R}^{n}:f(x)\leq f(x_{0}),g(x)\leq 0\}$, there is a constant $M$ such that

[TABLE]

Set

[TABLE]

Then

[TABLE]

Hence

[TABLE]

Or

[TABLE]

Integrating this inequality gives

[TABLE]

It follows that for $T\leq-\alpha\log(1-\zeta),$ we have

[TABLE]

Therefore, for $t\in(0,-\alpha\log(1-\zeta))$

[TABLE]

A choice of $\alpha>0$ with $1-2M\alpha>0$ completes the proof of the theorem. ∎

Remark 3*.*

Theorem 4 indicates that the time in which the trajectory needs to arrive at the boundary of the feasible set must be larger than $\alpha\log\frac{1}{1-\zeta}$.

Theorem 5**.**

Let $x_{\zeta}^{\sharp}$ be the first point where the trajectory reaches the boundary of the feasible set, then

[TABLE]

This is to say that the point $x_{\zeta}^{\sharp}$ belongs to the closure of the $\zeta$-neighborhood of the central path. Especially, the limit of $x_{\zeta}^{\sharp}$ as $\zeta\rightarrow 1^{-}$ is at the intersection of the central path and the boundary of the feasible set.

Proof.

Let $x_{\zeta}^{\sharp}=x(t^{\sharp};\zeta,x_{0})$, by the choice of the point $x_{\zeta}^{\sharp}$,

[TABLE]

Hence

[TABLE]

By (30),

[TABLE]

Therefore,

[TABLE]

This ends the proof. ∎

Remark 4*.*

According to Theorem 5, obviously, we have $\cos\theta\rightarrow-1^{+}$ as $\zeta\rightarrow 1^{-}$. The normalized centrality condition (8)

[TABLE]

can be satisfied. Given that $x_{\zeta}^{\sharp}$ is a point on the boundary of the feasible set, the first-order necessary conditions (KKT conditions) are therefore satisfied. This is straightforward as the Lagrange multiplier $\lambda^{\star}$ associated with the Lagrangian function $\mathcal{L}(x,\lambda)$ for the considered problem is

[TABLE]

and

[TABLE]

Under assumption (A2) and (A3), we have $\lambda^{\star}>0$. With $g(x_{\zeta}^{\sharp})=0$ strict complementarity holds. To find out whether the point $x_{\zeta}^{\sharp}$ is a local solution, one needs to check the second-order sufficient conditions. We discuss this in the next section.

Remark 5*.*

Based on Theorem 5, we propose an error measure $\epsilon>0$ for the optimization solutions

[TABLE]

According (52), we then have

[TABLE]

The KKT conditions are satisfied for $x_{\zeta}^{\sharp}$ when $\cos\theta=-1$ as mentioned in Remark 4. A small value $\epsilon>0$ seems to be a natural error measure for the present method. Additionally, $\epsilon$ is defined using the intrinsic parameter $\zeta$ that determines the shape of the optimization trajectory. With the error measure $\epsilon$ we can interpret Theorem 1(ii) and Theorem 4 as follows. Theorem 1(ii) implies the time for the trajectory to find a first-order $\epsilon$-optimal solution is $\mathcal{O}(1/\epsilon)$, which is the ergodic rate of convergence for first-order methods; Theorem 4 implies the trajectory is convergent to the trajectory $x(t;1,x_{0})$ uniformly for $t\in(0,\alpha\log(1/\epsilon))$.

6 Local convergence analysis

Consider optimization problem (1) with a single inequality constraint. For simplicity, let the origin be a point on the central path, and $x_{n}$ lies in the direction of $\frac{\nabla f}{|\nabla f|}$, i.e.,

[TABLE]

By the implicit function theorem, we have a function $x_{n}=\phi(x_{1},...,x_{n-1})$, which satisfies

[TABLE]

Furthermore:

• ①

At the origin,

[TABLE]

• ②

In a neighborhood of the origin,

[TABLE]

Thus we have

[TABLE]

The contours of the objective function is the graph of the implicit function $x_{n}=\phi(x_{1},...,x_{n-1})$. For $n=2$, the curvature formula gives

[TABLE]

For $n>2$, we choose a direction $\omega=(\omega_{1},\cdots,\omega_{n-1})$ with $|\omega|=1$. Then, in the two dimensional subspace defined by $\mathbb{R}^{2}=\{\omega t,x_{n}\}=\{\omega_{1}t,\cdots,\omega_{n-1}t,x_{n}\}$, the function $x_{n}=\phi(\omega t)$ has curvature in direction $\omega$,

[TABLE]

Similar analysis on $\psi$, which is deduced from the implicit function $g=\text{constant}$, gives:

[TABLE]

Now the conjectured relative convex condition (28) in the two dimensional subspace $\mathbb{R}^{2}=\{\omega t,x_{n}\}$ would be $\kappa_{f}^{\omega}<\kappa_{g}^{\omega}$, i.e.,

[TABLE]

for any direction $\omega$. Let $H_{\phi}(\mathbf{0})$ and $H_{\psi}(\mathbf{0})$ be the Hessian matrix of $\phi$ and $\psi$ about variable $\tilde{x}$, respectively. Then, in the sense of positive definite matrix, we have

[TABLE]

At $(x_{1},...,x_{n-1})=(0,...,0)$, we have

[TABLE]

Notice that

[TABLE]

Thus, we have

[TABLE]

Similarly, we have for $H_{\psi}(\mathbf{0})$

[TABLE]

Due to the positive definiteness (62), we have

[TABLE]

Recall at the origin we have

[TABLE]

Therefore we have

[TABLE]

Set

[TABLE]

Definition 1**.**

A point $x$ on the central path is called nondegenerate if $\tilde{C}(x)$ is invertible; it is relative convex if $\tilde{C}(x)$ is a positive definite matrix.

Remark 6*.*

Let a point $x^{\star}$ satisfy the KKT conditions. The relative convex condition $\tilde{C}(x^{\star})\succ 0$ is equivalent to the second-order sufficient conditions for constrained optimization. This is straightforward as the matrix

[TABLE]

is equivalent to the projected Hessian being positive definite (Nocedal, , p. 348), with a Lagrange multiplier $\lambda^{\star}$ satisfying the KKT conditions and strictly complementarity holding,

[TABLE]

Here, we deduce the relative convex condition from the perspective of the difference in the curvatures of the function contours in the feasible set. It is defined on the central path and may be seen as a perturbed version of the second-order sufficient conditions.

At the origin, the Jacobian matrix for $\frac{\nabla f}{|\nabla f|}$ reads

[TABLE]

For $1\leq i\leq n-1$, $\partial_{x_{i}}f=0$ at the origin, thus

[TABLE]

Recall also at the origin we have

[TABLE]

At the origin, $\frac{\partial_{x_{n}}f}{|\nabla f|}$ has a maximum value $1$, resulting in its derivatives being zero:

[TABLE]

Thus, the Jacobian matrix $J_{f}^{n}$ for $\frac{\nabla f}{|\nabla f|}$ at the origin reads

[TABLE]

where $J_{f}^{n-1}$ is an $(n-1)\times(n-1)$ matrix:

[TABLE]

Similarly, we have the Jacobian matrix $J_{g}^{n}$ for $\frac{\nabla g}{|\nabla g|}$ at the origin that reads

[TABLE]

where $J_{g}^{n-1}$ is also an $(n-1)\times(n-1)$ matrix:

[TABLE]

So we get matrix $\tilde{C}=J_{f}^{n-1}+J_{g}^{n-1}$ once more again.

By choosing a new coordinate system, say $\tilde{x}=(x_{1},...,x_{n-1})$ again, whose basis vectors are the eigenvectors $\mathbf{v}$ of $\tilde{C}$, the matrix $\tilde{C}$ transforms to a diagonal matrix $\tilde{C}_{\lambda}$ in the new coordinate system:

[TABLE]

with $\lambda_{1},...,\lambda_{n-1}$, which are the eigenvalues of the symmetric matrix $\tilde{C}$. Hence

[TABLE]

where $\mu_{i}=\frac{\partial}{\partial x_{n}}\left(\frac{\partial x_{i}f}{|\nabla f|}\right)+\frac{\partial}{\partial x_{n}}\left(\frac{\partial x_{i}g}{|\nabla g|}\right)$.

Assume that at the origin, the tangent of the central path $\mathcal{L}:\frac{\nabla f}{|\nabla f|}+\frac{\nabla g}{|\nabla g|}=0$ writes as

[TABLE]

The directional derivative of $\frac{\nabla f}{|\nabla f|}+\frac{\nabla g}{|\nabla g|}$ along the tangent of cental path $L$ is zero, we have

[TABLE]

We denote at origin:

[TABLE]

The Taylor formula of the negative search direction $\frac{\nabla f}{|\nabla f|}+\zeta\frac{\nabla g}{|\nabla g|}$ reads

[TABLE]

Summarizing the above analysis, we have

Lemma 6.1**.**

If a point on the central path $L$ is nondegenerate, then the system (29) may locally be considered as a perturbation of the linear equations:

[TABLE]

In the following, we use the matrix representations

[TABLE]

where $\delta_{ij}$ is the Kronecker-delta.

Lemma 6.2**.**

Given an initial point $(x_{1}^{0},x_{2}^{0},\cdots,x_{n}^{0}),~{}x_{n}^{0}>0$, the linearized equation system (84) has the solution

[TABLE]

with

[TABLE]

Proof.

Given an initial point $(x_{1}^{0},x_{2}^{0},\cdots,x_{n}^{0}),~{}x_{n}^{0}>0$, we have

[TABLE]

Therefore, for $1\leq i\leq n-1$ we have

[TABLE]

We write the ODE (87) in matrix form,

[TABLE]

(88) has the solution

[TABLE]

where $C_{0}$ is a vector of constants depending on the initial point $(x_{1}^{0},\cdots,x_{n}^{0})$. ∎

Theorem 6**.**

If $\tilde{C}_{\lambda}$ is positive definite, then there is a positive $\zeta_{0}<1$ such that the solution of the system (84) has the asymptotic line

[TABLE]

as $t\rightarrow+\infty$ for any $\zeta_{0}<\zeta<1$.

Proof.

The distance of a point on the trajectory (89) to the asymptotic line (90) is at most $|C_{0}e^{-\Lambda_{\zeta}t}|$. Choose a $\zeta_{0}$ such that $\Lambda_{\zeta}$ is a positive definite matrix, then the distance goes to zero as $t\rightarrow+\infty$ for any $\zeta_{0}<\zeta<1$. This is a proof of the theorem. ∎

To observe the behavior of the linearized equation (84) when $\zeta\rightarrow 1^{-}$, we eliminate the variable $t$ with $t=\frac{x_{n}^{0}-x_{n}}{1-\zeta}$, and rewrite the solution (89) as

[TABLE]

Notice that $x_{n}^{\prime}(t)<0$, and we have $x_{n}(t)=x_{n}^{0}-(1-\zeta)t<x_{n}^{0}.$

Theorem 7**.**

Let $\zeta\rightarrow 1^{-}$ and the point on the central path be relative convex, then the solution of system (84) converges to the tangent of the central path $\mathcal{L}$ defined in (80).

Proof.

Following the observation described above, $x_{n}<x_{n}^{0}.$ Let $\zeta\rightarrow 1^{-}$, then $(x_{n}^{0}-x_{n})/(1-\zeta)\rightarrow+\infty$. With the relative convex condition, we have $\lambda_{i}>0$. If $\zeta\rightarrow 1^{-}$, the eigenvalues of $\Lambda_{\zeta}$ $>\frac{1}{2}\min(\lambda_{i})>0$. Resulting in $\Lambda_{1}\succ 0$ by denoting $\Lambda_{\zeta=1}=\Lambda_{1}$. Therefore,

[TABLE]

The solution (91) converges to

[TABLE]

which is the tangent of the central path $\mathcal{L}$ as is given in (80). ∎

Similar method gives:

Theorem 8**.**

Suppose a point $x$ on the central path is nondegenerate and $\tilde{C}_{\lambda}(x)$ has at least one negative eigenvalue. Let $\zeta\rightarrow 1^{-}$, then the solution of system (84) will leave a neighborhood of the central path.

Remark 7*.*

The present trajectory converges to a central path with $\zeta\rightarrow 1^{-}$ under the relative convex condition. In Remark 4, we state that the first-order necessary conditions can be obtained as the trajectory reaches the boundary of the feasible set as $\zeta\rightarrow 1^{-}$. In Remark 6, we state that the second-order sufficient conditions is equivalent to the relative convex condition at the boundary of the feasible set. Combining both statements we conclude: by choosing $\zeta\rightarrow 1^{-}$, the present trajectory approaches a local solution by traversing along a central path.

7 The method for multiple constraints: a formulation based on the logarithmic barrier function

Recall the logarithmic barrier function for the problem (1),

[TABLE]

We first notice that the barrier function grows without bound if $g_{i}(x)\rightarrow 0^{-}$ and it is not differentiable at the boundary of the feasible set. To overcome this difficulty, we consider a subset of the feasible set,

[TABLE]

where M is a sufficiently large positive number so that $\Omega_{M}$ approximates the original feasible set $\Omega$. The barrier function is twice continuously differentiable in $\Omega_{M}$ and thus fulfills the assumption (A4).

Set

[TABLE]

The original problem (1) can be approximately reformulated as

[TABLE]

so that the present system (2) may be applied.

Notice that the global behavior shown in section 5 is based on the assumption that $\nabla g\neq 0$ in the feasible set $\Omega$. In general, this may not always be the case for the function $G(x),\forall x\in\Omega_{M}$. To escape the points where $\nabla\Phi(x)=0$, we suggest using the gradient descent $-\nabla f(x)$. Thus we get the system (4) for solving multiple inequality constrained optimization with $\zeta\in[0,1)$ as

[TABLE]

Note that in computational practice, $\nabla\Phi=0$ rarely occurs. The second ODE of the above system serves as a safeguard for the present method.

Corollary 1**.**

Let $x_{\zeta}^{\sharp}$ be the first point where the resulting trajectory of the system (4) reaches the boundary of $\Omega_{M}$, then

[TABLE]

This is to say that the point $x_{\zeta}^{\sharp}$ belongs to the closure of $\zeta$-neighborhood of the central path. Especially, the limit of $x_{\zeta}^{\sharp}$ as $\zeta\rightarrow 1^{-}$ is at the intersection of the central path and the boundary of the subset $\Omega_{M}$.

Proof.

A similar method as used in Theorem 5 gives the proof. ∎

Let $\mathbf{x}^{\star}\in\Omega_{M}$ be a point on the central path of the barrier function method, set

[TABLE]

Definition 2**.**

A point $\mathbf{x}^{\star}\in\Omega_{M}$ on the central path is relative convex if $\tilde{C}_{\Phi}(\mathbf{x}^{\star})$ is a positive definite matrix.

Corollary 2**.**

Let $\zeta\rightarrow 1^{-}$, and the point $\mathbf{x}^{\star}$ on the central path be relative convex. Then, the solution of the linearized system of (4) converges to the tangent line of the central path $\mathcal{L}$ at $\mathbf{x}^{\star}$.

Proof.

In the subset $\Omega_{M}$, $\nabla\Phi$ is Lipschitz continuous. Under the assumption (A2), we have $\nabla\Phi\neq 0$ along a central path. A uniformly continuous argument shows that $\nabla\Phi\neq 0$ in a close neighborhood of the central path. A similar method shown in section 6 gives a proof. ∎

To conclude, using the barrier function formulation for problem (1), the resulting trajectory achieves an approximated local solution that locates on the boundary of the subset $\Omega_{M}$. Notice, too, that the system (4) does not depend on the choice of $M$, and $\Omega_{M}$ exhaust $\Omega$, i.e., $\Omega=\cup_{M>0}\Omega_{M}$. This means that the resulting trajectory keeps approaching the boundary of the original feasible set $\Omega$ by crossing the boundary of any subset $\Omega_{M}$ dependent on $M$. Therefore, it eases the practical implementation of the method since no extra parameter $M$ needs to be defined for the stopping criterion, but rather, a check for each constraint violation would be sufficient.

8 Numerical experiments

When implementing the present method numerically, the canonical first-order optimization procedure is defined by the differential equation

[TABLE]

where

[TABLE]

In practical applications, this procedure might result in poor performance. We first notice that, according to Theorem 5, the accuracy of the solutions increases as the parameter $\zeta\rightarrow 1^{-}$. However, as $\zeta\rightarrow 1^{-}$, $|\mathbf{s}_{\zeta}|\downarrow 0$ at the central path, thus resulting in slow convergence rate. Another difficulty results from the potential poor scaling of the logarithmic barrier function, which results in ill-conditioned Hessian near the boundary of the feasible set (Nocedal, , p. 500-502). We will show test examples that suffer from this problem in the following.

To overcome these difficulties, a step size rule that considers and adapts the parameter $\zeta$ may be needed. A good practical self-adaptive $\zeta$ is expected to result in an optimization trajectory that behaves similar to the Mehrotra’s practical implementation for IPM Mehrotra , which is well-known for its very good performance and its difficulty in deriving a convergence theory. Design and analyze such a self-adaptive $\zeta$ require studies and analysis that exceeds the scope of this manuscript.

Another way to improve the performance might be the implementation of a momentum method or the Nesterov Accelerated Gradient nesterov1983method for the present system. Although our first implementations have shown promising performance on some test problems, we leave a systematical study together with the design of a step size rule (that considers self-adaptive $\zeta$) for future work.

In this work, rather, we propose a more fundamental framework so that the focus is on the numerical behavior of the present trajectory, neglecting any advanced modification. In this framework, we use the fixed parameter $\zeta$. In doing so, we can only get solutions that are $(1-\zeta)$-suboptimal solutions (referred in Theorem 5 and Remark 5). We let $\mathcal{X}_{\zeta}$ denote the set of all $(1-\zeta)$-suboptimal solutions,

[TABLE]

given that the solutions exist only on the boundary of the feasible set. A different way of defining the suboptimal solutions is given in skaf2010techniques and some useful applications using suboptimal solutions are shown.

To get better performance compared to the procedure (96), we suggest the following modified dynamical system:

[TABLE]

Referring to the work murray2019revisiting , we can say that the systems (98) and (96) are topologically equivalent and solutions of (98) are merely arc-length reparameterizations of (96) solutions. As the system (96) may move slowly at the close neighborhood of a central path, the system (98) moves in the same orbit with constant speed. In murray2019revisiting , the authors show that the normalized gradient descent method escapes saddle points “quickly”. This might be beneficial when solving nonconvex optimization problems, where the gradient of the function $\nabla f(x)$ vanishes at saddle points. In fact, this phenomenon may be similar to the present system (96), where $|\mathbf{s}_{\zeta}|\downarrow 0$ at the central path as $\zeta\rightarrow 1^{-}$. In the following, we show numerical experiments applying the system (98) with appropriate constant step sizes. Note, by using the logarithmic barrier function formulation for multiple constraints in the system (98), we may still suffer from potential poor scaling behavior near feasible set boundaries.

8.1 Experiments with common benchmarks

We first show numerical experiments for the inequality constrained problems in the EA competition at the 2006 IEEE Congress on Evolutionary Computation liang2006problem . These benchmarks are widely used among the community of evolutionary algorithms. We choose them as test examples for three reasons. First, they are well defined constrained optimization problems and have different characteristics rao2016jaya . Second, they are nontrivial to solve with first-order methods. Third, the relatively simple formulation of the optimization problems allows us to gain deeper insight into the numerical behavior of the present method. We choose the inequality constrained optimization problems for the experimentation. The problem $G12$ is excluded because it has a feasible set consisting of $9^{3}$ disjointed spheres. For the problem $G24$, which has a feasible set consisting of two disconnected sub-regions, we choose initial designs in the sub-region that contains the reported optimal solution.

We conduct numerical experiments with the parameter $\zeta=0.98$ and tuned fixed-step sizes. All bound constraints are treated as inequalities. Random initializations, that are away from the reported optimal solution, are selected in the feasible sets. As shown in table 1, apart from the problem G19, we find solutions for the test problems that have an absolute error less than 2e-2 compared to the reported optima. More accurate results can be obtained when we choose shorter step sizes and larger parameter $\zeta$.

For the problem G19, the reported optimal solution $\mathbf{x}^{\star}=$ (1.6699e-17, 3.9637e-16, 3.9459, 1.060e-16, 3.2831, 9.9999, 1.1283e-17, 1.2026e-17, 2.5071e-15, 2.2462e-15, 0.3708, 0.2785, 0.5238, 0.3886, 0.2982) is not achievable with a fixed-step size that has a Euclidean norm of $5e-2$.

For the problems G01, G04, G06, G07, G08, G24, we find sub-optimal solutions that are close to the reported optimal designs. For the problems G09, G10, G18, G19, sub-optimal solutions close to local minima are found. It appears that the “no free lunch theorem” wolpert1997no may apply to the present method implementation: while we solve some of the test problems in no more than a few hundred of iterations, a much larger number of iterations are needed for the remaining problems.

Still, to get more insight into the numerical behavior of the method implementation, we plot the centrality measure $\cos\theta$ over the optimization process for each test problem in figure 3. The dashed-lines indicate the $\zeta$-neighborhood. The results may be summarized in three categories:

Category 1. G04 and G09: optimization traverses within the $\zeta$-neighborhood;

Category 2. G01, G07, G08, G19, G19, and G24: optimization traverses to the $\zeta$-neighborhood but zigzags when it gets close to the solutions;

Category 3. G06 and G10: optimization zigzags around the $\zeta$-neighborhood during the optimization process.

The reasons for the zigzagging may be two-folds. First, it may be due to the poor scaling of the logarithmic barrier function near the boundary of the feasible set. The second reason may be the “overshooting”: a large fixed-step size is unable to achieve a small $\zeta$-neighborhood when close to an optimal solution. In problem G06 and G10, the central paths locate closely to the boundaries of the respective feasible sets, thus resulting in zigzagging throughout the whole optimization process. In figure 4, we show the narrow feasible set of the test problem G06. The central path traverses close to the two boundaries of the feasible set. A similar phenomenon can be observed in problem G10. In problem G08, the zigzagging disappears when sufficiently short step sizes are chosen. It thus supports our argument of the “overshooting”.

8.2 Test example in shape optimization

Although it is shown that the potential poor scaling of the logarithmic barrier function may result in increasing computational effort, we want to point out that, for the problems of categories 1 and 2, the logarithmic barrier function formulation can be very efficient in treating a large number of constraints. We support our argument with a shape optimization problem.

Here, we show an academic convex problem. The objective is to maximize the volume of a small sphere. This small sphere is located inside a bigger sphere that acts as the geometric constraint. The shape of both spheres are represented with finite element meshes zienkiewicz1977finite . The optimization problem writes

[TABLE]

where $V(x)$ is the volume function, $g_{i}(x)$ is a point-wise defined geometric constraint for the $i$-th design node, $m$ is the number of nodes of the design mesh, and $x\in\mathbb{R}^{3m}$ is the field of nodal coordinates of the design sphere mesh. The number of nodes of the small sphere (design sphere) is 19897. Thus, the total number of design variables is 59691 and the total number of constraints is 19897. We use the logarithmic barrier function for multiple constraints and choose the parameter $\zeta=0.95$. In figure 5, the shape variation process with depicted iterations is shown. Initially, the design sphere is located close to the boundary of the constraint sphere. During the shape variation process, it moves towards the center of the constraint sphere, while adapting its shape at each iteration. One could recognize easily that the central path of this optimization problem is being approached and followed until the solution is found when constraints become active.

8.3 Real-world application to shape optimization

We consider a real-world application to shape optimization. The present method is implemented in ShapeModule, which is a flexible solver-agnostic optimization platform and provides optimization algorithms as well as shape control methods, such as Vertex Morphing Hojjat . The optimization problem is to minimize the mass of a frame structure under load-displacement constraint (i.e., the displacement of every surface node is bounded). The optimization problem writes

[TABLE]

where $M(x)$ is the function for the mass, $g_{i}(x)$ is a point-wise formulated displacement constraint for the $i$-th node, $m$ is the number of nodes of the design surface mesh, and $x\in\mathbb{R}^{3m}$ is the field of nodal coordinates of the design surface mesh. The number of design variables is 144423, and the number of constraints is 48141. Note that for multiple constraints, we can use the logarithmic barrier function formulation as in the previous test examples. Each single displacement constraint gradient can be efficiently computed using the adjoint sensitivity analysis. In this application, we use the load-displacement sensitivity provided by the software OptiStruct to conform with a standard industrial design chain. We choose the parameter $\zeta=0.95$.

In figure 6 we show the initial frame design and the shape optimized design after 194 iterations. The mass of the structure is reduced by $41\%$ as shown in figure 7. In figure 8, we show a plot of maximum constraint value $g=max\{g_{i}\}$ of each iteration. In figure 9, we show that the optimization is able to approach and follow a central path within the $\zeta$-neighborhood.

Remark 8*.*

By following a central path, an intermediate design improves not only the objective function but also the constraint function. Take the design of iteration 80 as an example: the mass is reduced by $20.5\%$, and the displacement is reduced by $12.1\%$. These designs may enrich the design options if the original problem is reformulated as a bi-objective optimization problem, in which both mass and the maximum displacement are set as objectives. The resulting intermediate designs alongside a central path are approximated Pareto solutions.

9 Conclusion

This paper proposes a gradient descent akin method for solving inequality constrained nonlinear programming problems. We show the global behavior and convergence of the method and prove local convergence under the introduced relative convex condition. Robustness is shown in various computational test problems. It is also shown that the method exhibits the potential for solving nonconvex nonlinear problems. The method may be especially suited for large-scale problems due to the very cheap computational cost in each iteration. Practical implementations based on the present method are of interest in future work.

Acknowledgments.

The work by the author LC was done during his work as one of the coordinators at the Bavarian Graduate School of Computational Engineering (BGCE). The working experiences and financial support are gratefully acknowledged. The authors are grateful to the ShapeModule team at the BMW Group for providing a real-world model and their framework. The authors also thank Jian Cui from the Helmholtz Pioneer Campus for his generosity in proofreading this manuscript.