An Alternative Globalization Strategy for Unconstrained Optimization

TL;DR

This paper introduces a novel globalization strategy for unconstrained optimization that employs multiple points per iteration to enhance convergence speed and robustness, supported by theoretical guarantees and numerical experiments.

Contribution

It presents a new multi-point based globalization approach for unconstrained optimization with proven convergence and practical parallel implementation.

Findings

Supports rapid convergence from remote starting points

Demonstrates promising numerical results

Provides parallel implementation details

Abstract

We propose a new globalization strategy that can be used in unconstrained optimization algorithms to support rapid convergence from remote starting points. Our approach is based on using multiple points at each iteration to build a representative model of the objective function. Using the new information gathered from those multiple points, a local step is gradually improved by updating its direction as well as its length. We give a global convergence result and also provide parallel implementation details accompanied with a numerical study. Our numerical study shows that the proposed algorithm is a promising alternative as a globalization strategy.

Figures

Tables

Table 1. Table 1: Run times (in seconds) for varying values of n 𝑛 n and p 𝑝 p

Problem	$n$	$p=1$	$p=2$	$p=4$	$p=8$	$p=16$
	5M	8.80	4.41	2.24	1.29	0.89
COSINE	25M	43.98	22.07	11.19	6.33	4.35
	125M	219.77	111.74	57.98	31.62	21.47
	5M	25.19	12.57	6.49	4.08	2.88
NONCVXUN	25M	118.59	58.63	30.44	18.39	13.29
	125M	496.45	249.43	127.38	75.43	55.07
	5M	45.93	22.69	11.89	7.91	6.11
ROSENBR	25M	219.52	111.16	57.02	35.74	28.31
	125M	1113.83	554.07	283.83	178.41	140.60

Equations

$$\underset{x\in {\Re }^n}{\min }f(x),$$

$$\begin{array}[]{lll}f_{k}:=f(x_{k}),&g_{k}:=\nabla f(x_{k}),&x_{k}^{t}:=x_{k}+s^{t}_{k},\\[5.69054pt] f^{t}_{k}:=f(x_{k}^{t}),&g^{t}_{k}:=\nabla f(x_{k}^{t}),&y^{t}_{k}:=g^{t}_{k}-g_{k}.\end{array}$$

$$f_k^t-f_k\le \rho g_k^{⊺}{s}_k^t\text{mbox}forsome\rho \in (0,1).$$

$$m_k^t(s):={\alpha }_k^0(s)l_k^0(s)+{\alpha }_k^t(s)l_k^t(s-s_k^t),$$

$$l_k^0(s):=f_k+g_k^{⊺}s\text{mbox}and{l}_k^t(s-s_k^t):=f_k^t+(g_k^t{)}^{⊺}(s-s_k^t).$$

$${\alpha }_k^0(s)=\frac{(s-s_k^t{)}^{⊺}(-s_k^t)}{(-s_k^t{)}^{⊺}(-s_k^t)}\text{mbox}and{\alpha }_k^t(s)=\frac{s^{⊺}{s}_k^t}{(s_k^t{)}^{⊺}{s}_k^t},$$

$$\Vert s{\Vert }^2+\Vert s-s_k^t{\Vert }^2\le \Vert s_k^t{\Vert }^2.$$

$$\begin{array}[]{ll}\mbox{minimize\hskip 11.38109pt }&m_{k}^{t}(s)=\alpha_{k}^{0}(s)l_{k}^{0}(s)+\alpha_{k}^{t}(s)l_{k}^{t}(s-s^{t}_{k})\\ \mbox{subject to\hskip 11.38109pt}&\|s\|^{2}+\|s-s_{k}^{t}\|^{2}\leq\|s_{k}^{t}\|^{2},\end{array}$$

$$m_k^t(s)=f_k+({\bar{g}}_k^t{)}^{⊺}s+s^{⊺}{B}_k^{t}s,$$

$${\bar{g}}_k^t:=g_k+\frac{1}{\Vert s_k^t{\Vert }^2}(f_k^t-(g_k^t{)}^{⊺}{s}_k^t-f_k)s_k^t$$

$$B_k^t:=\frac{1}{\Vert s_k^t{\Vert }^2}{s}_k^t(g_k^t-g_k{)}^{⊺}.$$

$$\nabla m_k^t(s)=g_k+\frac{1}{\Vert s_k^t{\Vert }^2}(f_k^t-(g_k^t{)}^{⊺}{s}_k^t-f_k)s_k^t+\frac{1}{\Vert s_k^t{\Vert }^2}[s_k^t(y_k^t{)}^{⊺}+y_k^t(s_k^t{)}^{⊺}]s.$$

$${\alpha }^{\ast }:=\underset{0\le \alpha \le 1}{\mathrm{argmin}}{m}_k^t(\alpha s_k^t)=\min \left\{\frac{(y_k^t{)}^{⊺}{s}_k^t-f_k^t+f_k}{2(y_k^t{)}^{⊺}{s}_k^t},1\right\}.$$

$$\begin{array}[]{ll}\mbox{minimize\hskip 11.38109pt }&[(f_{k}^{t}-(g_{k}^{t})^{\intercal}s_{k}^{t}-f_{k})s_{k}^{t}+\|s_{k}^{t}\|^{2}g_{k}]^{\intercal}s+s^{\intercal}[s_{k}^{t}(g_{k}^{t}-g_{k})^{\intercal}]s\\ \mbox{subject to\hskip 11.38109pt}&s^{\intercal}(s-s_{k}^{t})\leq 0.\end{array}$$

$$(f_k^t-(g_k^t{)}^{⊺}{s}_k^t-f_k-\beta )s_k^t+\Vert s_k^t{\Vert }^2{g}_k+s_k^t(y_k^t{)}^{⊺}{s}_k^{t+1}+y_k^t(s_k^t{)}^{⊺}{s}_k^{t+1}+2\beta s_k^{t+1}$$

$$\beta (s_k^{t+1}-s_k^t{)}^{⊺}{s}_k^{t+1}$$

$$s_k^{t+1}=c_g(\beta )g_k+c_s(\beta )s_k^t+c_y(\beta )y_k^t,$$

$$c_g(\beta )=-\frac{\Vert s_k^t{\Vert }^2}{2\beta },c_s(\beta )=-\frac{\Vert s_k^t{\Vert }^2}{2\beta }\left(\frac{\delta }{\Vert s_k^t{\Vert }^2}-\frac{(y_k^t{)}^{⊺}{s}_k^t+2\beta }{\theta }((y_k^t{)}^{⊺}{g}_k+\frac{\delta }{\Vert s_k^t{\Vert }^2}(y_k^t{)}^{⊺}{s}_k^t)+\frac{\Vert y_k^t{\Vert }^2}{\theta }((s_k^t{)}^{⊺}{g}_k+\delta )\right),$$

$$c_y(\beta )=-\frac{\Vert s_k^t{\Vert }^2}{2\beta }\left(-\frac{(y_k^t{)}^{⊺}{s}_k^t+2\beta }{\theta }((s_k^t{)}^{⊺}{g}_k+\delta )+\frac{\Vert s_k^t{\Vert }^2}{\theta }((y_k^t{)}^{⊺}{g}_k+\frac{\delta }{\Vert s_k^t{\Vert }^2}(y_k^t{)}^{⊺}{s}_k^t)\right),$$

$$\delta =f_k^t-f_k-(g_k^t{)}^{⊺}{s}_k^t-\beta \text{mbox}and\theta ={\left((y_k^t{)}^{⊺}{s}_k^t+2\beta \right)}^2-\Vert s_k^t{\Vert }^2\Vert y_k^t{\Vert }^{.}$$

$$\frac{1}{4\beta {\theta }^2}P(\beta )=0,$$

$$\underset{s\in {\mathbb{R}}^n}{\min }\{m_k^t(s)+{\sigma }_1(s-s_k^t{)}^{⊺}s\}.$$

$${\sigma }_1=\frac{1}{\Vert s_k^t{\Vert }^2}(f_k^t-(g_k^t{)}^{⊺}{s}_k^t-f_k),$$

$$\underset{s\in {\mathbb{R}}^n}{\min }\{m_k^t(s)+{\sigma }_1(s-s_k^t{)}^{⊺}s+{\sigma }_2{s}^{⊺}s\}.$$

$$\underset{s\in {\mathbb{R}}^n}{\min }\{g_k^{T}s+s^{⊺}\frac{1}{\Vert s_k^t{\Vert }^2}[\sigma I+s_k^t(y_k^t{)}^{⊺}]s\},$$

$$\sigma =({\sigma }_1+{\sigma }_2)\Vert s_k^t{\Vert }^2=(f_k^t-(g_k^t{)}^{⊺}{s}_k^t-f_k)+{\sigma }_2\Vert s_k^t{\Vert }^2.$$

$$\Vert s_k^{t+1}\Vert \le \eta \Vert s_k^t\Vert$$

$$s_k^{t+1}=-\Vert s_k^t{\Vert }^2({\bar{B}}_k^t{)}^{-1}{g}_k,$$

$${\bar{B}}_k^t=2\sigma I+s_k^t(y_k^t{)}^{⊺}+y_k^t(s_k^t{)}^{⊺}.$$

$$\Vert s_k^{t+1}\Vert$$

$$=\left(\Vert s_k^t\Vert \Vert g_k\Vert {\lambda }_{\min }({\bar{B}}_k^t{)}^{-1}\right)\Vert s_k^t\Vert ,$$

Full text

1 Introduction.

In unconstrained optimization, frequently used algorithms, like quasi-Newton or trust-region methods, need to involve mechanisms that ensure convergence to local solutions from remote starting points. Roughly speaking, these globalization strategies guarantee that the improvement obtained by the algorithm is comparable to the improvement obtained with a gradient step [NW06].

In this paper we present a new globalization strategy that can be used in unconstrained optimization methods for solving problems of the form

[TABLE]

where $f$ is a first order differentiable function. Conventional methods for solving this problem mostly use local approximations that are very powerful once the algorithm arrives at the close proximity of a stationary point. The main idea behind the proposed strategy is based on using additional information collected from multiple points to construct an adaptive approximation of the function $f$ in (1). In particular, we update our approximate model constructed around the current iterate and a sequence of trial points. Our objective is to come up with a better local model than the one obtained by the current iterate only. We observe that the collection of this additional information has a profound effect on the performance of the method as well. Before moving to the next iterate, the step is improved by updating its direction as well as its length simultaneously. This step computation involves only the inner products of vectors. Therefore, each iteration of the algorithm is amenable to a parallel implementation. Though acquiring the additional information from multiple points may add an extra burden on the algorithm, this burden can be alleviated by using the readily available parallel processors. Furthermore, our numerical experiments demonstrate that the additional computations at each step may reduce the total number of iterations, since we learn more about the function structure.

Parallel execution of linear algebra operations is common in the parallel optimization literature; generally in designing parallel implementations of existing methods. Among the earliest work is [BRS88], where parallelization of linear algebra steps and function evaluations in the BFGS method is discussed. Around the same time, [NS89] propose to use truncated-Newton methods combined with computation of the search direction via block iterative methods. In particular, they focus on the block conjugate-direction method. By using block iterative methods, they parallelize some steps of the algorithm, in particular, the linear algebra operations. For extensive reviews of the topic, we refer to [Zen89, MTK03, DW03].

The use of multiple trial points or directions is also a recurring idea in the field of parallel nonlinear optimization. However, unlike our case, these mostly depend on concurrent information collection / generation. [Pat84] proposes an algorithm, which is based on evaluating multiple points in parallel to determine the next iterate. In the parallel line search implementation of [PFZ98], multiple step sizes are tried in parallel. [vL85] discusses three parallel quasi-Newton algorithms that are also based on considering multiple directions to extend the rank-one updates. To solve the original problem by solving a series of independent subproblems, [Han86] suggests using conjugate subspaces for quasi-Newton updates. Thus, the resulting method can be applied in a parallel computing environment.

The rest of this paper is organized as follows. In Section 2, we discuss the derivation and the theoretical properties of the proposed strategy. In Section 3, the success of the main idea of the method is numerically tested, and its parallel implementation is discussed in Section 4. Conlusion of the study is presented in Section 5.

2 Proposed Globalization Strategy.

At the core of the proposed strategy lies an extended model function, which is constructed by using two linear models at two reference points. The first of these points is the current iterate and the second one is a trial point. When one of these trials leads to a successful step, then the algorithm moves to a new iterate. To make our discussion more concrete, we start by describing an iteration of the algorithm.

Let $x_{k}$ and $s_{k}$ denote the current iterate and the $k^{th}$ step of the algorithm, respectively. Then, the next iterate becomes $x_{k+1}=x_{k}+s_{k}$. These are the outer iterations. The vector $s_{k}$ is obtained by executing a series of inner iterations and computing trial steps $s_{k}^{t}$, $t=0,1,\cdots$. We can further simplify our exposition by defining

[TABLE]

The trial step $s_{k}^{t}$ is accepted as new $s_{k}$, if it provides sufficient decrease in the sense

[TABLE]

This is nothing but the well-known Armijo condition [NW06].

Up to this point, the algorithm behaves like a typical unconstrained minimization procedure. However, to come up with the proposed globalization strategy, we design a new subprocedure for computing the trial steps, $s_{k}^{t}$. The inner step computation is done by solving a subproblem that consists of an extended model function updated at every inner iteration and a constraint restricting the length of the steps. The extended model function is an approximation to the objective function, $f$. It is constructed by using the information gathered around the region bordered by $x_{k}$ and $x^{t}_{k}$. Actually, the extended model function is a combination of the linear models of $f$ at $x_{k}$ and $x^{t}_{k}$. That is,

[TABLE]

where

[TABLE]

This construction is illustrated in Figure 1. Note that both weights, $\alpha_{k}^{0}(s)$ and $\alpha_{k}^{t}(s)$ are functions of $s$. We want to make sure that if $s$ is closer to $s_{k}^{t}$, then the weight $\alpha_{k}^{t}(s)$ of the linear model around $x_{k}$ increases. Similarly, $\alpha_{k}^{t}(0)$ should increase as $(s-s_{k}^{t})$ gets closer to $(-s_{k}^{t})$. This is achieved by measuring the length of the projections of $(s-s_{k}^{t})$ and $s$ on $s_{k}^{t}$ by

[TABLE]

respectively.

Next, we introduce a constraint to handle the possibly hard quadratic function $m_{k}^{t}$ as well as control the step length. Recall that our model function is constructed around the previous trial step $s_{k}^{t}$ by considering the area lying between $x_{k}$ and $x_{k}+s_{k}^{t}$. Thus, we impose the constraint

[TABLE]

Clearly, this constraint controls not only the length of $s$ but also its deviation from the previous trial point $s_{k}^{t}$. In addition, this choice of the constraint ensures that both weight functions $\alpha_{k}^{0}(s)$ and $\alpha_{k}^{t}(s)$ are in $[0,1]$ and they add up to

1. Hence, our model function $m_{k}^{t}$ is the convex combination of the linear functions $l_{k}^{0}$ and $l_{k}^{t}$. This observation is formally presented in Lemma 2.1 and its proof is given in Appendix A.

Lemma 2.1

If constraint (5) holds at inner iteration $t+1$ of iteration $k$, then both $\alpha_{k}^{0}(s)$ and $\alpha_{k}^{t}(s)$ are nonnegative, and they satisfy $\alpha_{k}^{0}(s)+\alpha^{t}(s)=1$.

We also make the following observations about the region defined by constraint (5). First, this region is never empty, since any step $s=\gamma s_{k}^{t}$ for $\gamma\in(0,1]$ is in the region. This implies in a sense that a backtracking operation on the previous trial step, $s_{k}^{t}$ would give a new feasible trial step, $s_{k}^{t+1}$. Second, there are always feasible steps in the steepest descent direction whenever $g_{k}^{\intercal}s_{k}^{t}\leq 0$ holds. That is, the steps $s=-\xi g_{k}$ are feasible for $0\leq\xi\leq-\frac{g_{k}^{\intercal}s_{k}^{t}}{\|g_{k}\|^{2}}$. Thus, as long as the first trial step, $s_{k}^{0}$ is gradient related, each feasible region constructed around a subsequent trial step involves a feasible point that is also gradient related.

Algorithm 1 gives the outline of the proposed method. The stopping condition is the usual check whether the current iterate, $x_{k}$ is a stationary point. We start the inner iterations by computing the first trial step, $s_{k}^{0}$. We shall elaborate on how to select this initial trial step in Section 4. As long as the trial steps $s_{k}^{t}$, $t=0,1,2,\ldots$ are not acceptable, we carry on with constructing and minimizing the model functions. Once an acceptable step is obtained by the inner iteration, we evaluate the next iterate $x_{k+1}$ and continue with solving the overall problem.

In Algorithm 1, the crucial part is solving subproblem (6), which has a quadratic objective function and a quadratic constraint. After some derivation, the objective function of this subproblem simplifies to

[TABLE]

where

[TABLE]

and

[TABLE]

The Hessian of this quadratic objective function is the rank-one matrix $B^{t}_{k}$. Its only nonzero eigenvalue is $(s_{k}^{t})^{\intercal}y_{k}^{t}$, and the eigenvectors corresponding to this eigenvalue are $s_{k}^{t}$ and $-s_{k}^{t}$. Furthermore, its eigenvectors corresponding to the zero eigenvalues are perpendicular to $y_{k}^{t}$. Let us note at this point that the gradient of the model $m_{k}^{t}(s)$ is given by

[TABLE]

2.1 Convex Objective Function.

Before considering the general case, we discuss in this section how subproblem (6) can be solved when the objective function is convex. We shall soon observe that it is not straightforward to obtain trial steps that yield sufficient descent with this initial design. Thus, we will propose modifications over this first attempt in Section 2.2 to guarantee convergence. However, we still present here our observations with the convex case as it constitutes the basis of our approach for general functions given in the next section.

We first check whether the subproblem always provides a nonzero step at every inner iteration. It turns out that if the most recent trial step, $s_{k}^{t}$ is gradient related, then the next trial step, $s_{k}^{t+1}$ cannot be the zero vector. Furthermore, if the objective function value of the recent trial step is strictly greater than the current objective function value, then the optimal solution of the subproblem is in the interior of the region defined by the constraint (5). Lemma 2.2 summarizes this discussion. The proof of this lemma is given in Appendix A.

Lemma 2.2

Let $f$ be a convex function and consider the inner iteration $t+1$ of iteration $k$ with $\|s_{k}^{t}\|\neq 0$. If $g_{k}^{\intercal}s_{k}^{t}\leq 0$, then $\|s_{k}^{t+1}\|\neq 0$ and

[TABLE]

If we additionally have $f_{k}^{t}>f_{k}$, then $\alpha^{*}<1$.

Note that we require condition $g_{k}^{\intercal}s_{k}^{t}\leq 0$ hold for all trial steps to guarantee nonzero trial steps. So, a convergence argument for the proposed algorithm is not clear unless it ensures this condition. We shall visit this issue later in this section.

We next show how subproblem (6) can be solved, if the conditions given in Lemma 2.2 hold. This subproblem can be rewritten as

[TABLE]

Note that this is a convex program as $(s_{k}^{t})^{\intercal}(g_{k}^{t}-g_{k})\geq 0$ by the convexity assumption on $f$ in this section. Let $s_{k}^{t+1}$ denote the optimal solution of the above problem and $\beta\geq 0$ be the Lagrange multiplier corresponding to the constraint. The Karush-Kuhn-Tucker optimality conditions imply

[TABLE]

Simplifying (8) leads to

[TABLE]

where

[TABLE]

[TABLE]

with

[TABLE]

Consider first the case when $\beta=0$. Then, using Lemma 2.2, we have $s_{k}^{t+1}=\alpha^{*}s_{k}^{t}$. When $\beta>0$, the optimal solution is at the boundary of the constraint. If we now use (10) with (9), then after some derivation we obtain

[TABLE]

where $P(\beta)$ is a sixth order polynomial function of $\beta$ (see the online supplement111http://people.sabanciuniv.edu/sibirbil/glob_strat/OnlineSupplement.pdf for the details). Thus, the subproblem can be solved by computing the roots of this polynomial.

In fact, it is not hard to show that if the inner iteration step $s_{k}^{t}$ satisfies $s_{k}^{t}=-M^{t}g_{k}$ for some matrix $M^{t}$, then at the next inner iteration we have $s_{k}^{t+1}=-M^{t+1}g_{k}$ However, it is not possible to guarantee that the matrix $M^{t+1}$ will be positive definite given $M^{t}$ is. That is because the gradient of the model function is given as the conical combination of $g_{k}$ and $-s_{k}^{t}$. First of all, the model gradient can be zero for $s=0$ at a point $x_{k}$ with $g_{k}\neq 0$. Also, when the deviation of the objective function $f$ from linearity is relatively large, then the model gradient can point a backward direction; i.e. the reverse direction of $s_{k}^{t}$. Since the constraint of the model does not allow such a backward step, this can cause $s_{k}^{t+1}$ to be zero. Moreover, the model gradient is scaled with a matrix, and the scaled step adds a third component along $y_{k}$ to the linear combination whose coefficient is not necessarily positive. So, the resulting $s_{k}^{t+1}$ might not satisfy $g_{k}^{T}s_{k}^{t+1}\leq 0$ even if $s_{k}^{t}$ does.

2.2 General Objective Function.

Using our observations on the convex case, we next concentrate on a construction where the model function $m_{k}^{t}$ has the same gradient value at $s=0$ as the original objective function, $f$. Then, we also add a regularization term so that the subproblem can be made convex even if $f$ is not.

We start by relaxing the constraint (5) and moving it to the objective function. That is

[TABLE]

We choose

[TABLE]

so that $\nabla m_{k}^{t}(0)=g_{k}$. Observe that the above choice of $\sigma_{1}$ can be negative (for $f$ nonconvex), in which case the constraint operates the reverse way; i.e. it tries to keep $s$ away from $s_{k}^{t}$. Since that causes us to loose our control on the size of $s$, we add a regularization term and the subproblem becomes

[TABLE]

After simplifying the objective function, we obtain our new subproblem for the general case

[TABLE]

where

[TABLE]

Note that the model function (11) is always convex provided that the regularization parameter $\sigma_{2}$ is chosen sufficiently large. Beyond the convexity of $m_{k}^{t}$, we want to guarantee for some $\eta\in(0,1)$ that the steplength condition

[TABLE]

holds. Suppose $s_{k}^{t+1}$ is such a step and it is the optimal solution for problem (11) with a convex objective function. Then, the first order optimality condition implies

[TABLE]

where

[TABLE]

Thus, we obtain

[TABLE]

where $\lambda_{\min}$ denotes the smallest eigenvalue of its matrix parameter. It is not difficult to see that $2\sigma$ is an eigenvalue of $\bar{B}_{k}^{t}$ with multiplicity $n-2$, and the remaining two eigenvalues are the extreme eigenvalues of $\bar{B}_{k}^{t}$ given by

[TABLE]

where $\lambda_{\max}$ denotes the largest eigenvalue of its matrix parameter.

To obtain a convex model, we need

[TABLE]

On the other hand, the steplength condition (12) requires

[TABLE]

Since the latter bound is larger, using relation (15) for choosing $\sigma_{2}$ provides the convexity requirement and satisfies the steplength condition. That is

[TABLE]

Our last step is to state the minimizer $s^{t+1}_{k}$ of (11) with selected $\sigma_{1}$ and $\sigma_{2}$ values. Following similar steps as in Section 2.1, after some derivation we obtain

[TABLE]

where

[TABLE]

[TABLE]

with

[TABLE]

and

[TABLE]

It is important to observe that the solution of the subproblem does not require storing any matrices and the main computational burden is only due to inner products, which could be done in parallel (see Section 4).

When it comes to the convergence of the proposed algorithm, we are saved basically by keeping the directions of its steps gradient related. However, we could not directly refer to the existing convergence results, since the directions in our algorithm change during inner iterations and the steplength is not computed by using a one dimensional function. In the following convergence note, we only assume additionally that the objective function $f$ is bounded below and its gradient $\nabla f$ is Lipschitz continuous with parameter $L$.

Lemma 2.3

Suppose that the first trial step $s_{k}^{0}$ satisfies $m_{0}\|g_{k}\|\leq\|s_{k}^{0}\|\leq M_{0}\|g_{k}\|$ and $0\geq s_{k}^{0}g_{k}\geq-\lambda_{0}\|g_{k}\|^{2}$ for some $m_{0},M_{0},\lambda_{0}\in(0,\infty)$. Then, at any iteration $k$, the optimal solution $s_{k}^{t+1}$ in (16) becomes an acceptable step satisfying (2) in finite number of inner iterations at a nonstationary point $x_{k}$.

Proof.

Consider any iteration $k$. Since $s_{k}^{t+1}$ is obtained by (13), the steps computed at each inner iteration satisfy $-g_{k}^{\intercal}s_{k}^{t+1}\geq\lambda_{t+1}\|g_{k}\|^{2}$, $t=0,1,2,\cdots$ for $\lambda_{t+1}=\|s_{k}^{t}\|^{2}(\lambda_{max}(\bar{B}_{k}^{t}))^{-1}$. Note for $t=0,1,\cdots$ that $\lambda_{t}>0$ by the choice of $\lambda_{0}$ and the relation (14). We have

[TABLE]

Using (12), we have $\|s_{k}^{t}\|\leq\eta_{t}M_{0}\|g_{k}\|$. Then,

[TABLE]

Therefore, we obtain

[TABLE]

Suppose that the acceptance criterion (2) is never satisfied as $t\rightarrow\infty$. Then,

[TABLE]

This implies

[TABLE]

Next the desired inequality is obtained by

[TABLE]

Since $\eta\in(0,1)$, the right hand side of the last inequality increases without a bound as $t\rightarrow\infty$. This gives a contradiction as $L$ is bounded. So, an acceptable $\|s_{k}^{t}\|$ should be obtained in finite number of inner iterations. ∎

Note that the conditions on the initial step length, $s_{k}^{0}$ in Lemma 2.3 are simply satisfied, if we choose $s_{k}^{0}=-\epsilon g_{k}$ for any $\epsilon\in(0,1]$. This implies $m_{0}=M_{0}=\lambda_{0}=\epsilon$.

Theorem 2.1

Let $\{x_{k}\}$ be the sequence of iterates generated by Algorithm 2 and $\{s_{k}^{0}\}$ satisfy the requirements in Lemma 2.3. Then, any limit point of $\{x_{k}\}$ is a stationary point of the objective function $f$.

Proof.

By Lemma 2.3, for any $k$, (2) is satisfied after a finite number of inner iterations, say $t(k)$. Also, note that relation (13) implies $\|s_{k}^{t}\|\geq\lambda_{t}\|g_{k}\|$ with the same $\lambda_{t}$ as in (18). We will now derive a lower bound for $\lambda_{t(k)}$:

[TABLE]

So,

[TABLE]

Recall that for the accepted step $s_{k}$ of iteration $k$ we have $-s_{k}^{\intercal}g_{k}\geq\lambda_{t(k)}\|g_{k}\|^{2}$. Clearly, $\lambda_{t(k)}>0$ is bounded away from zero as $t(k)$ is finite. Consider now any subsequence of $\{x_{k}\}$ with indices $k\in\mathcal{K}$ such that

[TABLE]

Then, for any $k,k^{\prime}\in\mathcal{K}$ with $k^{\prime}>k$ we have

[TABLE]

Since $\{x_{k}\}_{k\in\mathcal{K}}$ converges to $\hat{x}$, the continuity of $f$ implies that $f_{k}\rightarrow\hat{f}$ as $k\in\mathcal{K}$. Therefore $f_{k}-f_{k^{\prime}}\rightarrow 0$ as $k,k^{\prime}\rightarrow\infty$. Thus, we obtain $\nabla f(\hat{x})=0$ by (19) since $\lambda_{t(k)}$ is bounded away from zero. ∎

2.3 Relationship with Quasi-Newton Updates.

Observing the formula in (13), we next question the relationship between the step computation of the proposed algorithm and the quasi-Newton updates, in particular the DFP formula [NW06]. The key to understanding this relationship is to treat the current iterate $x_{k}$ as the previous iterate, whereas our trial step $x_{k}^{t}$ becomes the next iterate in standard quasi-Newton updates.

Let $D_{k}$ denote the quasi-Newton approximation to the Hessian at iterate $x_{k}$. Then, the infamous DFP formula yields the new approximation to the Hessian in the next iterate $x_{k+1}$ as

[TABLE]

where $s_{k}=x_{k+1}-x_{k}$ and $y_{k}=\nabla f(x_{k+1})-\nabla f(x_{k})$. We first substitute $s_{k}^{t}$ for $s_{k}$ and $y_{k}^{t}$ for $y_{k}$ in the above equation. Then, by setting

[TABLE]

we obtain

[TABLE]

Let us now denote the Hessian of the objective our model function given in relation (11) by $H_{k}^{t}$. Then, we have

[TABLE]

Next we plug in the value of $\sigma$ from (17) and obtain

[TABLE]

The last term above, in a sense, serves the purpose of obtaining a positive definite matrix. Thus, the objective function becomes convex. Looking at our approach from the quasi-Newton approximation point of view, we again confirm that the additional information collected from the trial iterate $x_{k}^{t}$ is used to come up with a better model around the current iterate $x_{k}$. Our derivation above also shows that the construction of the model function $m^{t+1}(s)$ from one inner iteration to the next is completely independent.

3 Practical Performance.

In this section, we shall conduct a numerical study to demonstrate the novel features of the proposed globalization strategy. We have compiled a set of unconstrained optimization problems from the well-known CUTEst collection [GOL04]. We have also embedded the proposed strategy (PS) into the L-BFGS package [Oka10] as a new globalization alternative to the existing two line search routines; backtracking (BT) and More-Thuente line search (MT).

First, let us give a two-dimensional example to illustrate how the proposed strategy can be beneficial. Figure 2 shows the contours of the Rosenbrock function. The inner iterates of PS are denoted by $+$ markers, the inner iterates of BT are denoted by $\circ$ markers, and the inner iterates of MT are denoted by $\times$ markers. All three globalization strategies start from the same initial point. The numbers next to the trial points denote the inner iterations. The corresponding objective function values are given in the subplot at the southeast corner of the figure.

Figure 2 illustrates the main point of the proposed strategy: PS can change the direction of the trial step as well as its length. It is important to note that PS is not a trust-region method because it builds a new model function by using the current trial step for computing the next trial step. However, in a trust-region method the model function is not updated within the same iteration and the current trial step does not contribute to the computation of the next trial step. For this particular two-dimensional problem, PS performs better than the other two globalization strategies as it attains the lowest objective function value after three inner iterations.

Next, we test the effect of the direction updates introduced by the proposed strategy. We solve the CUTEst problems by employing the backtracking line search with Wolfe conditions using the default reduction factor of 0.5 given in the L-BFGS package. For the proposed strategy, we also set the parameter $\eta$ to 0.5. Therefore, the new strategy requires the same amount of steplength decrease (but it may additionally alter the search direction). In our tests, we start with both the L-BFGS step with memory set to 5 ($m=5$) and the gradient step ($m=0$). The maximum number of iterations is set to 500 for the L-BFGS step, and to 1,000 for the gradient step. The comparison is carried out for those problems, where both methods obtained the same local solution (which is assessed by comparing the objective function values and the first two components of the solution vectors). The results are reported for 30 problems when $m=0$, and for 63 problems when $m=5$.

Figure 3 summarizes our observations. We compare two strategies in terms of the total number of function evaluations required to solve the problems. When $m=0$, solving $60\%$ of the problems required less number of function evaluations than the line search. The same figure becomes $50.94\%$ when $m=5$, and solving $26.42\%$ of the problems took exactly the same number of iterations with both strategies. Overall, these results show that the direction updates of the new strategy can indeed provide improvements.

4 Parallel Implementation.

In this section, we will provide an outline of the subproblem solution operations while we discuss their parallel execution. Then, we conduct some numerical experiments to illustrate the parallel performance of the proposed algorithm.

As relation (16) shows, the basic trial step at inner iteration $t+1$ is in the subspace spanned by $g_{k}$, $g_{k}^{t}$, and the previous trial step $s_{k}^{t}$. We should also emphasize that the computation of $s_{k}^{t+1}$ requires just a few floating point operations, once the necessary inner products are completed. In particular, the computation of (16) requires the following inner products.

[TABLE]

Once the scalar values $v_{i}$, $i=1,\cdots,6$ are available, the multipliers $c_{g}(\sigma)$, $c_{s}(\sigma)$, $c_{y}(\sigma)$ in (16) can be computed in a total of 32 floating point operations as follows:

[TABLE]

Clearly, the total cost of synchronization operations is negligible. Thus, the parallel execution of the above inner products determine the parallel performance as well as the two remaining components of the algorithm: the computation of initial trial step $s_{k}^{0}$ and function/gradient evaluations. We require that $s_{k}^{0}$ is computed without introducing too much sequential work to the overall algorithm, and satisfy the conditions in Section 2. A simple choice could be the gradient step, which we also preferred in our own implementation.

The parallel implementation given in Algorithm 2 is considered for a shared memory environment with $p$ physical cores (threads). Note that the parallel implementation of function/gradient computations is clearly problem dependent. However –at least partial– separability is likely to occur for large-scale problems. This issue, we assume, is taken care of by the user.

We coded the algorithm in C using OpenMP. We have selected three test problems from the CUTEst collection [GOL04], whose dimensions can be varied to obtain small- to large-scale problems, namely COSINE, NONCVXUN, and ROSENBR. The objective functions of the problems, and the initial points used in our tests are as follows:

[TABLE]

where $x_{i}$ denotes the $i$th component of point $x$ and $x^{0}$ is the initial point. We note that the function and gradient evaluations of all three test problems are in the order of $O(n)$. So, function evaluations are not dominant to the computations of the new strategy. Moreover, the evaluations of all three functions can be done in parallel with one synchronization in computing the objective function value.

To test the parallelization performance of the algorithm, we set the dimension of the problems to $5$ million, $25$ million and $125$ million, and solve each of the resulting instances for various values of the total number of threads, $p\in\{1,2,4,8,16\}$. We run the algorithm by setting $\eta=0.5$, and the step acceptability is checked using condition (2) with $\rho=0.1$. We set both the maximum number of inner iterations and the maximum number of outer iterations to $100$. If the algorithm cannot reach a solution with the desired accuracy

[TABLE]

then we report the clock time elapsed until the maximum number of outer iterations is reached. We summarize these results in Table 1 We observe that the speed-up is close to linear for up to 8 cores. Since there is no data reuse, if we further increase the number of cores, the main memory bandwidth becomes the bottleneck. This is a common issue when memory-bounded computations, like inner product evaluations, are parallelized.

Finally, we consider the load balance issue. To observe the usage of capacity and the workload distribution, we plot the CPU usage during the solution of 1M-dimensional instance of the problem COSINE (see Figure 4). The plots are obtained by recording the CPU usage information per second via the mpstat command-line tool (see Linux manual pages). Figure 4(a) reveals the idle resources when the program runs sequentially. In fact, during the sequential run, the average resource usage stays at a level of $11.12\%$. This value is computed by taking the average usage of eight cores. When eight threads are used, this average raises up to $78.65\%$. Figure 4(c) shows the usage of all available resources as well as the distribution of the workload.

5 Conclusion and Future Research.

The current study is a part of our research efforts that aim at harnessing parallel processing power for solving optimization problems. Our main motivation here was to present the design process that we went through to come up with an alternate globalization strategy for unconstrained optimization that can be implemented on a shared memory parallel environment.

The proposed algorithm works with a model function and considers its optimization over an elliptical region. However, it does not employ a conventional trust-region approach. Nor does it match with any one of the well-known quasi-Newton updates. The basic idea is to use multiple trial points for learning the function structure until a good point is obtained. These trials constitute the inner iterations. Fortunately, all these computations of the algorithm are domain-separable with two $O(1)$ synchronization points at each inner iteration.

Our numerical experience verifies that the direction updates of the resulting algorithm can reduce the total number of function evaluations required, and it is scalable to a degree allowed by the inner iterations with a balanced distribution of the workload among parallel processors. Like any parallel optimization method, the parallel performance of the proposed algorithm can be compromised if function evaluations are computationally intensive and unfit for parallelization.

We only solved a set of examples with our algorithm. It is yet to be seen whether the algorithm is apt for solving other large-scale problems, especially those ones arising in different applications. In our implementation, we chose the initial trial point simply by using the negative gradient step. One may try to integrate other, potentially more powerful but still parallelizable, steps to the algorithm. An example could be using a truncated Newton step; perhaps computed by a few iterations of the conjugate gradient on the initial quadratic model. Then switching to the globalization strategy, as described here, could adjust the direction and the length of the initial step. Finally, our observations on the relationship of the new strategy with the quasi-Newton update formulas can be studied further in future research.

A Omitted Proofs.

Lemma 2.1* If constraint (5) holds at inner iteration $t+1$ of iteration $k$, then both $\alpha_{k}^{0}(s)$ and $\alpha_{k}^{t}(s)$ are nonnegative, and they satisfy $\alpha_{k}^{0}(s)+\alpha^{t}(s)=1$. *

Proof.

We have

[TABLE]

Since constraint (5) implies $s^{\intercal}s\leq s^{\intercal}s_{k}^{t}$, we obtain

[TABLE]

The result follows from the definitions of $\alpha^{0}(s)$ and $\alpha^{t}(s)$. ∎

Lemma 2.2* Let $f$ be a convex function and consider the inner iteration $t+1$ of iteration $k$ with $\|s_{k}^{t}\|\neq 0$. If $g_{k}^{\intercal}s_{k}^{t}\leq 0$, then $\|s_{k}^{t+1}\|\neq 0$ and*

[TABLE]

*If we additionally have $f_{k}^{t}>f_{k}$, then $\alpha^{*}<1$. *

Proof.

By using relation (7), we have

[TABLE]

Since $f$ is a convex function, we know that $f_{k}^{t}-(g_{k}^{t})^{\intercal}s_{k}^{t}-f_{k}\leq 0$. If $f_{k}^{t}-(g_{k}^{t})^{\intercal}s_{k}^{t}-f_{k}=0$, then $\|\nabla m_{k}^{t}(0)\|=\|g_{k}\|\neq 0$ because $x_{k}$ is not a stationary point. Now, consider the case $f_{k}^{t}-(g_{k}^{t})^{\intercal}s_{k}^{t}-f_{k}<0$, and suppose for contradiction that $\|\nabla m_{k}^{t}(0)\|=0$. Then,

[TABLE]

Multiplying both sides with $g_{k}$, we obtain

[TABLE]

Since $g_{k}^{\intercal}s_{k}^{t}\leq 0$, we obtain a contradiction and hence $\|\nabla m_{k}^{t}(0)\|\neq 0$. This shows that $\|s_{k}^{t+1}\|\neq 0$. Note that

[TABLE]

This implies that $s_{k}^{t}$ is a descent direction for $m_{k}^{t}(s)$ at $s=0$. Then, it is easy to solve the one-dimensional convex optimization problem to obtain

[TABLE]

Again by using relation (7), we have

[TABLE]

which gives

[TABLE]

If $f_{k}^{t}<f_{k}$ along with $g_{k}^{\intercal}s_{k}^{t}\leq 0$, then

[TABLE]

holds. Therefore, we have $-(s^{t}_{k})^{\intercal}\nabla m_{k}^{t}(s^{t}_{k})\leq 0$, which implies that $(-s_{k}^{t})$ is a descent direction for $m_{k}^{t}(s)$ at $s=s_{k}^{t}$. In this case, the minimizer along $s_{k}^{t}$ should be in the interior. Thus, $\alpha^{*}<1$. ∎