DEFT-FUNNEL: an open-source global optimization solver for constrained grey-box and
black-box problems

The local-search algorithm starts by building an initial interpolation set either from a simplex or by drawing samples from an uniform distribution. The construction of the interpolation set is done in the function deft_funnel_build_initial_sample_set, which is called only once during a local search within deft_funnel. The choice between random sampling and simplex is done in deft_funnel when calling the function deft_funnel_set_parameters, which defines the majority of parameters of the local search. If random sampling is chosen, it checks if the resulting interpolation set is well poised and, if not, it is updated using the Algorithm 6.3 described in Chapter 6 in [12], which is implemented in deft_funnel_repair_Y.

For the sake of simplicity, we assume henceforth that the objective function is also a black box. Given a poised set of sample points ${\cal Y}_{0}=\{y^{0},y^{1},\ldots,y^{p}\}$ with an initial point $x_{0}\in{\cal Y}_{0}$, the next step of our algorithm is to replace the objective function $f(x)$ and the black-box constraint functions $c(x)=(c_{1}(x),c_{2}(x),\ldots,c_{q}(x))$ by surrogate models $m^{f}(x)$ and $m^{c}(x)=(m^{c_{1}}(x),m^{c_{2}}(x),\ldots,m^{c_{q}}(x))$ built from the solution of the interpolation system

$$M(\phi,{\cal Y})\alpha_{\phi}=\Upsilon({\cal Y}),$$

(25)

where

$$M(\phi,{\cal Y})=\begin{pmatrix}\phi_{0}(y^{0})&\phi_{1}(y^{0})&\cdots&\phi_{b}(y^{0})\\ \phi_{0}(y^{1})&\phi_{1}(y^{1})&\cdots&\phi_{b}(y^{1})\\ \vdots&\vdots&\ddots&\vdots\\ \phi_{0}(y^{p})&\phi_{1}(y^{p})&\cdots&\phi_{b}(y^{p})\end{pmatrix},\;\;\Upsilon({\cal Y})=\begin{pmatrix}\Upsilon(y^{0})\\ \Upsilon(y^{1})\\ \vdots\\ \Upsilon(y^{p})\end{pmatrix},\;\;p\leq b,$$

where $\phi$ is the basis of monomials and $\Upsilon(x)$ is replaced by the objective function $f(x)$ or some black-box constraint function $c_{j}(x)$.

We consider underdetermined quadratic interpolation models that are fully linear and that are enhanced with curvature information along the optimization process. Since the linear system (25) is potentially underdetermined, the resulting interpolating polynomials will be no longer unique and so we provide to the user four approaches to construct the models $m^{f}(x)$ and $m^{c_{j}}(x)$ that can be chosen by passing a number from 1 to 4 to the input ‘whichmodel’: 1 - subbasis selection approach; 2 - minimum $\ell_{2}$-norm models (by defatul); 3 - minimum Frobenius norm models; and 4 - regression mnodels (recommended for noisy functions). Further details about each one can be found in [12].

If $p_{0}=|{\cal Y}_{0}|=n+1$, a linear model rather than an underdetermined quadratic model is built for each function. The reason is that, despite both having error bounds that are linear in $\Delta$ for the first derivatives, the error bound for the latter includes also the norm of the model Hessian, as stated in Lemma 2.2 in [46], which makes it worse than the former.

Whenever $n+1<p_{k}\leq(n+1)(n+2)/2=p^{\max}$, the algorithm builds underdetermined quadratic models based on the choice of the user between the approaches described above. If regression models are considered instead, we set $p^{\max}=(n+1)(n+2)$, which means that the sample set is allowed to have twice the number of sample points required for fully quadratic interpolation models. Notice that having a number of sample points much larger than the required for quadratic interpolation can also worsen the quality of the interpolation models as the sample set could contain points that are too far from the iterate, which is not ideal for models built for local approximation.

It is also possible to choose the initial degree of the models between fully linear, quadratic with a diagonal Hessian or fully quadratic. This is done within deft_funnel by setting the input argument cur_degree in the call to deft_funnel_set_parameters to one of the following options: model_size.plin, model_size.pdiag or model_size.pquad.

The interpolation system (25) is solved using a QR factorization of the matrix $M(\phi,{\cal Y})$ within the function deft_funnel_computeP, which is called by deft_funnel_build_models.

In order to evaluate the error of the interpolation models and their derivatives with respect to the original functions $f$ and $c$, we make use of the measure of well poisedness of ${\cal Y}$ given below.

As it is shown in [12], the error bounds for at most fully quadratic models depend linearly on the constant $\Lambda$; the smaller it is, the better the interpolation models approximate the original functions. We also note that the error bounds for undetermined quadratic interpolation models are linear in the diameter of the smallest ball containing ${\cal Y}$ for the first derivatives and quadratic for the function values.

Finally, since the constraint functions $c(x)$ are replaced by surrogate models $m^{c}(x)$ in the algorithm, we define the models $m^{z}(x)\stackrel{{\scriptstyle\rm def}}{{=}}(m^{c}(x),h(x))$ and $m^{z}(x,s)\stackrel{{\scriptstyle\rm def}}{{=}}m^{z}(x)-s$, which are those used for computing new directions.

In this subsection, we explain how the subspace minimization is employed in our algorithm. We define the subspace ${\cal S}_{k}$ at iteration $k$ as

$${\cal S}_{k}\stackrel{{\scriptstyle\rm def}}{{=}}\{x\in\mathbb{R}^{n}\;|\;\left[x\right]_{i}=\left[l^{x}\right]_{i}\;\textrm{ for }\;i\in{\cal L}_{k}\;\textrm{ and }\;\left[x\right]_{i}=\left[u^{x}\right]_{i}\;\textrm{ for }\;i\in{\cal U}_{k}\},$$

where ${\cal L}_{k}\stackrel{{\scriptstyle\rm def}}{{=}}\{i\;|\;\left[x_{k}\right]_{i}-\left[l^{x}\right]_{i}\leq\epsilon_{b}\}$ and ${\cal U}_{k}\stackrel{{\scriptstyle\rm def}}{{=}}\{i\;|\;\left[u^{x}\right]_{i}-\left[x_{k}\right]_{i}\leq\epsilon_{b}\}$ define the index sets of (nearly) active variables at their bounds, for some small constant $\epsilon_{b}>0$. If ${\cal L}_{k}\neq\emptyset$ or ${\cal U}_{k}\neq\emptyset$, a new well-poised interpolation set ${\cal Z}_{k}$ is built and a recursive call is made in order to solve the problem in the new subspace ${\cal S}_{k}$.

If the algorithm converges in a subspace ${\cal S}_{k}$ with an optimal solution $(\tilde{x},\tilde{s})$, it checks if the latter is also optimal for the full-space problem, in which case the algorithm stops. If not, the algorithm continues by attempting to compute a new direction in the full space. This procedure is illustrated in Figure 2.

The dimensionality reduction of the problem mitigates the chances of degeneration of the interpolation set when the sample points become too close to each other and thus affinely dependent. Figure 3 gives an example of this scenario as the optimal solution is approached.

In order to check the criticality in the full-space problem, a full-space interpolation set of degree $n+1$ is built in an $\epsilon$-neighborhood around the point $x^{*}_{{\cal S}}$, which is obtained by assembling the subspace solution $\tilde{x}$ and the $|{\cal L}_{k}\,\cup\,{\cal U}_{k}|$ fixed components $\left[x_{k}\right]_{i}$. The models $m^{f}_{k}$ and $m^{c}_{k}$ are then updated and the criticality step is entered.

The complete subspace minimization step is described in Algorithm 3.2.2 and it is implemented in the function deft_funnel_subspace_min, which is called inside deft_funnel_main.

 Algorithm 3.2: SubspaceMinimization(${\cal S}_{k-1}$, ${\cal Y}_{k}$, $x_{k}$, $s_{k}$, $\Delta_{k}^{f}$, $\Delta_{k}^{z}$, $v_{k}^{\max}$)  

1:

Check for (nearly) active bounds at $x_{k}$ and define ${\cal S}_{k}$. If there is no (nearly) active bound or if ${\cal S}_{k}$ has already been explored, go to Step 5. If all bounds are active, go to Step 4.

2:

Build a new interpolation set ${\cal Z}_{k}$ in ${\cal S}_{k}$.

3:

Call recursively LocalSearch(${\cal S}_{k}$, ${\cal Z}_{k}$, $x_{k}$, $s_{k}$, $\Delta_{k}^{f}$, $\Delta_{k}^{z}$, $v_{k}^{\max}$) and let $(x^{*}_{{\cal S}},s^{*}_{{\cal S}})$ be the solution of the subspace problem after adding the fixed components.

4:

If $dim({\cal S}_{k-1})<n$, return $(x^{*}_{{\cal S}},s^{*}_{{\cal S}})$. Otherwise, reset $(x_{k},s_{k})=(x^{*}_{{\cal S}},s^{*}_{{\cal S}})$, construct new set ${\cal Y}_{k}$ around $x_{k}$, build $m_{k}^{f}$ and $m_{k}^{c}$ and recompute $\pi^{f}_{k-1}$ (optimality measure).

5:

If ${\cal S}_{k}$ has already been explored, set $(x_{k+1},s_{k+1})=(x_{k},s_{k})$, reduce the trust regions radii $\Delta_{k+1}^{f}=\gamma\Delta_{k}^{f}$ and $\Delta_{k+1}^{z}=\gamma\Delta_{k}^{z}$, set $\Delta_{k+1}=\min[\Delta_{k+1}^{f},\Delta_{k+1}^{z}]$ and build a new poised set ${\cal Y}_{k+1}$ in ${\cal B}(x_{k+1};\Delta_{k+1})$.

 

The normal step aims at reducing the constraint violation at given point $(x,s)$ defined by

$\displaystyle v(x,s)\stackrel{{\scriptstyle\rm def}}{{=}}{\scriptstyle\frac{1}{2}}\|m^{z}(x,s)\|^{2}.$

(26)

To ensure that the step $n_{k}$ is normal to the linearized constraint $m^{z}(x_{k},s_{k})+J(x_{k},s_{k})n=0$, where $J(x,s)\stackrel{{\scriptstyle\rm def}}{{=}}(J(x)\;-I)$ is the Jacobian of $m^{z}(x,s)$ with respect to $(x,s)$ and $J(x)$ is the Jacobian of $m^{z}(x)$ with respect to $x$, we require that

$\displaystyle\|n_{k}\|_{\infty}\leq\kappa_{n}\|m^{z}(x_{k},s_{k})\|,$

(27)

for some $\kappa_{n}\geq 1$.

The computation of $n_{k}$ is done by solving the constrained linear least-squares subproblem

$\displaystyle\left\{\begin{array}[]{rc}\displaystyle\min_{n=(n^{x},n^{s})}&{\scriptstyle\frac{1}{2}}\|m^{z}(x_{k},s_{k})+J(x_{k},s_{k})n\|^{2}\\ \textrm{s.t.:}&l^{s}\leq s_{k}+n^{s}\leq u^{s},\\ &l^{x}\leq x_{k}+n^{x}\leq u^{x},\\ &x_{k}+n^{x}\in{\cal S}_{k},\\ &n\in{\cal N}_{k},\end{array}\right.$

(33)

where

$\displaystyle{\cal N}_{k}\stackrel{{\scriptstyle\rm def}}{{=}}\{n\in\mathbb{R}^{n+m}\mid\|n\|_{\infty}\leq\min\left[\,\Delta^{z}_{k},\,\kappa_{n}\,\|m^{z}(x_{k},s_{k})\|\,\right]\,\},$

(34)

for some trust-region radius $\Delta_{k}^{z}>0$. Finally, a funnel bound $v_{k}^{\max}$ is imposed on the constraint violation $v_{k}\stackrel{{\scriptstyle\rm def}}{{=}}v(x_{k},s_{k})$ for the acceptance of new iterates to ensure the convergence towards feasibility.

We notice that, although a linear approximation of the constraints is used for calculating the normal step, the second-order information of the quadratic interpolation model $m^{z}(x,s)$ is still used in the SQO model employed in the tangent step problem as it is shown next.

The subproblem (33) is solved within the function deft_funnel_normal_step, which makes use of an original active-set algorithm where the unconstrained problem is solved at each iteration in the subspace defined by the currently active bounds, themselves being determined by a projected Cauchy step. Each subspace solution is then computed using a SVD decomposition of the reduced matrix. This algorithm is implemented in the function deft_funnel_blls and is intended for small-scale bound-constrained linear least-squares problems.

The tangent step is a direction that improves optimality and it is computed by using a SQO model for the problem (14) after the normal step calculation. The quadratic model for the function objective function is defined as

$\displaystyle\psi_{k}((x_{k},s_{k})+d)\stackrel{{\scriptstyle\rm def}}{{=}}m^{f}(x_{k},s_{k})+\langle g_{k},d\rangle+{\scriptstyle\frac{1}{2}}\langle d,B_{k}d\rangle,$

(35)

where $m^{f}(x_{k},s_{k})\stackrel{{\scriptstyle\rm def}}{{=}}m^{f}(x_{k})$, $g_{k}\stackrel{{\scriptstyle\rm def}}{{=}}\nabla_{(x,s)}m^{f}(x_{k},s_{k})$, and $B_{k}$ is the approximate Hessian of the Lagrangian function

	$\displaystyle{\cal L}(x,s,\mu,\xi^{s},\tau^{s},\xi^{x},\tau^{x})$	$\displaystyle=m^{f}(x,s)+\langle\mu,m^{z}(x,s))\rangle+\langle\tau^{s},s-u^{s}\rangle+\langle\xi^{s},l^{s}-s\rangle$
		$\displaystyle\quad\;+\langle\tau^{x},x-u^{x}\rangle+\langle\xi^{x},l^{x}-x\rangle$

with respect to $(x,s)$, given by

$\displaystyle B_{k}=\begin{pmatrix}H_{k}+\sum_{i=1}^{m}[\hat{\mu}_{k}]_{i}Z_{ik}&0\\ 0&0\end{pmatrix},$

(36)

where $\xi^{s}$ and $\tau^{s}$ are the Lagrange multipliers associated to the lower and upper bounds, respectively, on the slack variables $s$, and $\xi^{x}$ and $\tau^{x}$, the Lagrange multipliers associated to the lower and upper bounds on the $x$ variables. In (36), $H_{k}=\nabla^{2}_{xx}m^{f}(x_{k},s_{k})$, $Z_{ik}=\nabla^{2}_{xx}m^{z}_{ik}(x_{k},s_{k})$ and the vector $\hat{\mu}_{k}$ may be viewed as a local approximation of the Lagrange multipliers with respect to the equality constraints $m^{z}(x,s)=0$.

By applying (24) into (35), we obtain

$$\begin{array}[]{rcl}\psi_{k}((x_{k},s_{k})+n_{k}+t)&=&\psi_{k}((x_{k},s_{k})+n_{k})+\langle g_{k}^{N},t\rangle+{\scriptstyle\frac{1}{2}}\langle t,B_{k}t\rangle,\end{array}$$

(37)

where

$\displaystyle g_{k}^{N}\stackrel{{\scriptstyle\rm def}}{{=}}g_{k}+B_{k}\,n_{k}.$

(38)

Since (37) is a local approximation for the function $m^{f}((x_{k},s_{k})+n_{k}+t)$, a trust region with radius $\Delta_{k}^{f}$ is used for the complete step $d=n_{k}+t$:

$\displaystyle{\cal T}_{k}\stackrel{{\scriptstyle\rm def}}{{=}}\{d\in\mathbb{R}^{n+m}\mid\|d\|_{\infty}\leq\Delta^{f}_{k}\}.$

(39)

Moreover, given that the normal step was also calculated using local models, it makes sense to remain in the intersection of both trust regions, which implies that

$\displaystyle d_{k}\in{\cal R}_{k}\stackrel{{\scriptstyle\rm def}}{{=}}{\cal N}_{k}\cap{\cal T}_{k}\stackrel{{\scriptstyle\rm def}}{{=}}\{d\in\mathbb{R}^{n+m}\mid\|d\|_{\infty}\leq\Delta_{k}\},$

(40)

where $\Delta_{k}=\min[\Delta_{k}^{z},\Delta_{k}^{f}]$.

In order to make sure that there is still enough space left for the tangent step within ${\cal R}_{k}$, we first check if the following constraint on the normal step is satisfied:

$\displaystyle\|n_{k}\|_{\infty}\leq\kappa_{{\cal R}}\Delta_{k},$

(41)

for some $\kappa_{{\cal R}}\in(0,1)$. If (41) holds, the tangent step is calculated by solving the following subproblem

$\displaystyle\left\{\begin{array}[]{rc}\displaystyle\min_{t=(t^{x},t^{s})}&\langle g_{k}^{N},t\rangle+{\scriptstyle\frac{1}{2}}\langle t,B_{k}t\rangle\\ \textrm{s.t.:}&J(x_{k},s_{k})t=0,\\ &l^{s}\leq s_{k}+n^{s}_{k}+t^{s}\leq u^{s},\\ &l^{x}\leq x_{k}+n^{x}_{k}+t^{x}\leq u^{x},\\ &x_{k}+n^{x}_{k}+t^{x}\in{\cal S}_{k}.\\ &n_{k}+t\in{\cal R}_{k},\\ \end{array}\right.$

(48)

where we require that the new iterate $x_{k}+d_{k}$ belongs to subspace ${\cal S}_{k}$ and that it satisfies the bound constraints  (16) and (17). In the Matlab code, the tangent step is calculated by the function deft_funnel_tangent_step, which in turn calls either the Matlab solver linprog or our implementation of the nonmonotone spectral projected gradient method [8] in deft_funnel_spg to solve the subproblem (48). The choice between both solvers is based on whether $\|B_{k}\|\leq\epsilon$, for a small $\epsilon>0$, in which case we assume that the problem is linear and therefore linprog is used.

We define our $f$-criticality measure as

$\displaystyle\pi_{k}^{f}\stackrel{{\scriptstyle\rm def}}{{=}}-\langle g_{k}^{N},r_{k}\rangle,$

(49)

where $r_{k}$ is the projected Cauchy direction obtained by solving the linear optimization problem

$\displaystyle\left\{\begin{array}[]{rc}\displaystyle\min_{r=(r^{x},r^{s})}&\langle g_{k}^{N},r\rangle\\ \textrm{s.t.:}&J(x_{k},s_{k})r=0,\\ &l^{s}\leq s_{k}+n^{s}_{k}+r^{s}\leq u^{s},\\ &l^{x}\leq x_{k}+n^{x}_{k}+r^{x}\leq u^{x},\\ &x_{k}+n^{x}_{k}+r^{x}\in{\cal S}_{k}.\\ &\|r\|_{\infty}\leq 1.\end{array}\right.$

(56)

By definition, $\pi_{k}^{f}$ measures how much decrease could be obtained locally along the projection of the negative of the approximate gradient $g_{k}^{N}$ onto the nullspace of $J(x_{k},s_{k})$ intersected with the region delimited by the bounds. This measure is computed in deft_funnel_compute_optimality, which uses lingprog to solve the subproblem (56).

A new local estimate of the Lagrange multipliers $(\mu_{k},\xi^{s}_{k},\tau^{s}_{k},\xi^{x}_{k},\tau^{x}_{k})$ are computed by solving the following problem:

$\displaystyle\left\{\begin{array}[]{rc}\displaystyle\min_{(\mu,\hat{\xi}^{s},\hat{\tau}^{s},\hat{\xi}^{x},\hat{\tau}^{x})}&{\scriptstyle\frac{1}{2}}\|{\cal M}_{k}(\mu,\hat{\xi}^{s},\hat{\tau}^{s},\hat{\xi}^{x},\hat{\tau}^{x})\|^{2}\\ \textrm{s.t.:}&\hat{\xi}^{s},\hat{\tau}^{s},\hat{\xi}^{x},\hat{\tau}^{x}\geq 0,\end{array}\right.$

(59)

where

	$\displaystyle{\cal M}_{k}(\mu,\hat{\xi}^{s},\hat{\tau}^{s},\hat{\xi}^{x},\hat{\tau}^{x})$	$\displaystyle\stackrel{{\scriptstyle\rm def}}{{=}}$	$\displaystyle\left(\begin{array}[]{c}g_{k}^{N}\\ 0\end{array}\right)+\left(\begin{array}[]{c}J(x_{k})^{T}\\ -I\end{array}\right)\mu+\left(\begin{array}[]{c}0\\ I_{\tau}^{s}\end{array}\right)\hat{\tau}^{s}+\left(\begin{array}[]{c}0\\ -I_{\xi}^{s}\end{array}\right)\hat{\xi}^{s}$		(73)
			$\displaystyle+\left(\begin{array}[]{c}I_{\tau}^{x}\\ 0\end{array}\right)\hat{\tau}^{x}+\left(\begin{array}[]{c}-I_{\xi}^{x}\\ 0\end{array}\right)\hat{\xi}^{x},$

the matrix $I$ is the $m\times m$ identity matrix, the matrices $I_{\xi}^{s}$ and $I_{\tau}^{s}$ are obtained from $I$ by removing the columns whose indices are not associated to any active (lower and upper, respectively) bound at $s_{k}+n^{s}_{k}$, the matrices $I_{\xi}^{x}$ and $I_{\tau}^{x}$ are obtained from the $n\times n$ identity matrix by removing the columns whose indices are not associated to any active (lower and upper, respectively) bound at $x_{k}+n^{x}_{k}$, and the Lagrange multipliers $(\hat{\xi}^{s},\hat{\tau}^{s},\hat{\xi}^{x},\hat{\tau}^{x})$ are those in $(\xi^{s},\tau^{s},\xi^{x},\tau^{x})$ associated to active bounds at $s_{k}+n^{s}_{k}$ and $x_{k}+n^{x}_{k}$. All the other Lagrange multipliers are set to zero.

The subproblem (59) is also solved using the active-set algorithm implemented in the function deft_funnel_blls.

The algorithm computes normal and tangent steps depending on the measures of feasibility and optimality at each iteration. Differently from [42, 43], where the computation of the normal steps depends on the measure of optimality, here the normal step is computed whenever the following condition holds

$\displaystyle\|z(x_{k},s_{k})\|>\epsilon,$

(74)

for some small $\epsilon>0$, i.e. preference is always given to feasibility. This choice is based on the fact that, in many real-life problems with expensive functions and a small budget, one seeks to find a feasible solution as fast as possible and that a solution having a smaller objective function value than the current strategy is already enough. If (74) fails, we set $n_{k}=0$.

We define a $v$-criticality measure that indicates how much decrease could be obtained locally along the projection of the negative gradient of the Gauss-Newton model of $v$ at $(x_{k},s_{k})$ onto the region delimited by the bounds as

$$\pi_{k}^{v}\stackrel{{\scriptstyle\rm def}}{{=}}-\langle J(x_{k},s_{k})^{T}z(x_{k},s_{k}),b_{k}\rangle,$$

where the projected Cauchy direction $b_{k}$ is given by the solution of

$\displaystyle\left\{\begin{array}[]{rc}\displaystyle\min_{b=(b^{x},b^{s})}&\langle J(x_{k},s_{k})^{T}z(x_{k},s_{k}),b\rangle\\ \textrm{s.t.:}&l^{s}\leq s_{k}+b^{s}\leq u^{s},\\ &l^{x}\leq x_{k}+b^{x}\leq u^{x},\\ &x_{k}+b^{x}\in{\cal S}_{k},\\ &\|b\|_{\infty}\leq 1.\end{array}\right.$

(80)

We say that $(x_{k},s_{k})$ is an infeasible stationary point if $z(x_{k},s_{k})\neq 0$ and $\pi_{k}^{v}=0$, in which case the algorithm terminates.

The procedure for the calculation of the normal step is given in the algorithm below. In the code, it is implemented in the function deft_funnel_normal_step, which calls the algorithm in deft_funnel_blls in order to solve the normal step subproblem (33).

 Algorithm 3.3: NormalStep($x_{k}$, $s_{k}$, $\pi_{k}^{v}$, $v_{k}$, $v_{k}^{\max}$)  

1:

If $z(x_{k},s_{k})\neq 0$ and $\pi_{k}^{v}=0$, STOP (infeasible stationary point).

2:

If (74), compute a normal step $n_{k}$ by solving (33). Otherwise, set $n_{k}=0$.

 

If the solution of (56) is $r_{k}=0$, then by (49) we have $\pi_{k}^{f}=0$, in which case we set $t_{k}=0$. If the current iterate is farther from feasibility than from optimality, i.e., for a given a monotonic bounding function $\omega_{t}$, the condition

$\displaystyle\pi_{k}^{f}>\omega_{t}(\|z(x_{k},s_{k})\|)$

(81)

fails, then we skip the tangent step computation by setting $t_{k}=0$.

After the computation of the tangent step, the usefulness of the latter is evaluated by checking if the following conditions

$\displaystyle\|t_{k}\|>\kappa_{{\cal Z}{\cal S}}\|n_{k}\|$

(82)

and

$\displaystyle\delta_{k}^{f}\stackrel{{\scriptstyle\rm def}}{{=}}\delta^{f,t}_{k}+\delta^{f,n}_{k}\geq\kappa_{\delta}\delta^{f,t}_{k},$

(83)

where

$\displaystyle\delta^{f,t}_{k}\stackrel{{\scriptstyle\rm def}}{{=}}\psi_{k}((x_{k},s_{k})+n_{k})-\psi_{k}((x_{k},s_{k})+n_{k}+t_{k})$

(84)

and

$\displaystyle\delta^{f,n}_{k}\stackrel{{\scriptstyle\rm def}}{{=}}\psi_{k}(x_{k},s_{k})-\psi_{k}((x_{k},s_{k})+n_{k}),$

(85)

are satisfied for some $\kappa_{{\cal Z}{\cal S}}>1$ and for $\kappa_{\delta}\in(0,1)$. The inequality (83) indicates that the predicted improvement in the objective function obtained in the tangent step is substantial compared to the predicted change in $f$ resulting from the normal step. If (82) holds but (83) fails, the tangent step is not useful in the sense just discussed, and we choose to ignore it by resetting $t_{k}=0$.

The tangent step procedure is stated in Algorithm 3.2.5 and it is implemented in the function deft_funnel_tangent_step.

 Algorithm 3.4: TangentStep($x_{k}$, $s_{k}$, $n_{k}$)  

1:

If (41) holds, then

1.1:

select a vector $\hat{\mu}_{k}$ and define $B_{k}$ as in (36);

1.2:

compute $\mu_{k}$ by solving (59);

1.3:

compute the modified Cauchy direction $r_{k}$ by solving (56) and define $\pi_{k}^{f}$ as (49);

1.4:

if (81) holds, compute a tangent step $t_{k}$ by solving (48).

2:

If (41) fails, set $\mu_{k}=\mu_{k-1}$. In this case, or if (81) fails, or if (82) holds but (83) fails, set $t_{k}=0$ and $d_{k}=n_{k}$.

3:

Define $(x_{k}^{+},s_{k}^{+})=(x_{k},s_{k})+d_{k}$.

 

Depending on the contributions of the current iteration in terms of optimality and feasibility, we classify it into three types: $\mu$-iteration, $f$-iteration and $z$-iteration. This is done by checking if some conditions hold for the trial point defined as

$\displaystyle(x_{k}^{+},s_{k}^{+})\stackrel{{\scriptstyle\rm def}}{{=}}(x_{k},s_{k})+d_{k}.$

(86)