Convergence analysis of an infeasible quasi-Newton bundle method for nonsmooth convex programming

Abstract

By utilizing the improvement function, we change the nonsmooth convex constrained optimization into an unconstrained optimization, and construct an infeasible quasi-Newton bundle method with proximal form. It should be noted that the objective function being minimized in unconstrained optimization subproblem may vary along the iterations (it does not change if the null step is made, otherwise it is updated to a new function). It is necessary to make some adjustment in order to obtain the convergence result. We employ the main idea of infeasible bundle method of Sagastizàbal and Solodov, and under the circumstances that each iteration point may be infeasible for primal problem, we prove that each cluster point of the sequence generated by the proposed algorithm is the optimal solution to the original problem. Furthermore, for BFGS quasi-Newton algorithm with strong convex objective function, we obtain the condition which guarantees the boundedness of quasi-Newton matrices and the R-linear convergence of the iteration points.

Keywords

Non-smooth optimization / convex constraint / improvement function / bundle method / quasi-Newton direction

1 Introduction

In the field of operations research, it is often necessary to minimize a function that has no derivatives at some points, and such problems are called non-differentiable optimization (NDO) or nonsmooth optimization (NSO). Many practical problems are NSO problems. Let the nonnegative variable x represent income. In many economic systems, income x corresponds to a tax function T(x), which usually has discontinuous derivatives. Given some threshold and tax rate 0=r_0,r_1,\dots ,r_{m-1}. T is given in the following form: , where T_0=0,T_i=T_{i-1}+a_i(r_{i-1}-r_i). It is easy to verify that if , then T is a continuously increasing convex function. Further, it is easy to verify that T has the following form: T(x)=max\{T_i+r_{i}x\mid i=0,1,\dots ,m-1\}. We usually want to minimize the tax rate, and thus the following optimization problem arises: min\{T(x)\mid -x⩽0\}, where T(x) is a convex function and the set of constraints \{x\in \mathbb{R}\mid -x⩽0\} is a convex set, which is a convex programming, and its more general form is

(1)\begin{array}{r}\text{minimize}f(x)\text{s.t.}c(x)⩽0,\end{array}

where f,c:{\mathbb{R}}^n\to \mathbb{R} are convex function, which is generally not differentiable. Suppose the Slater bound norm of Problem (1) holds, i.e., there exists x\in {\mathbb{R}}^n such that . We also assume that there exists an oracle that for any given x\in {\mathbb{R}}^n can compute the function value f(x),c(x) and one subgradient of each function, i.e., g_f(x)\in \mathrm{\partial }f(x) and g_c(x)\in \mathrm{\partial }c(x).

In general, NSO problems are difficult to solve. Among the many solution methods, the well-known methods are the subgradient method, the cutting plane method, the analytic central tangent plane method, and the bundle method. Among them, the latter two are currently considered to be the most stable and trustworthy methods for solving NSO problems. The NSO problems with constraints are more complicated [4]. For general nonsmooth convex constrained problems, one approach is to equivalently solve the unconstrained optimization problem with an exact penalty function [7]. However, the application of penalty functions inevitably involves the problem of taking values of penalty parameters. If a large parameter value is required to guarantee the accuracy of a given penalty function, which will increase the difficulty in numerical computation. There are also some bundle methods that do not use the idea of penalty functions [12]. However, all descent steps in these methods are required to be feasible, including the initial point. We know that this requirement is likely to be as difficult as solving Problem (1). As a result, the entire computational burden of solving the problem increases substantially. Other work on bundle methods for solving constrained optimization problems can be found in the literature [5, 8, 9, 14]. One of them, [14] gives a dual stable bundle method resulting from combining the proximal bundle method and the horizontal bundle method, which combines the advantages of both, while proposing an effective stopping criterion. When the objective function is nonconvex, [5] constructed a reassignment bundle method. Under the condition of using only the approximate information of the objective and constraint functions, [8] gave a proximal bundle method, which requires the error between the approximate and real values to be fixed but can be unknown, and finally obtains the approximate optimal solution of the problem.

We give the improvement function related to Problem (1): for a given x\in {\mathbb{R}}^n, let

\begin{array}{r}{h}_x(y)=max\{f(y)-f(x),c(y)\},\mathrm{\forall }y\in {\mathbb{R}}^n.\end{array}

It is called the improvement function of Problem (1). \overline{x} is the solution of Problem (1) when and only when \overline{x} is the optimal solution of the unconstrained optimization problem min{h}_{\overline{x}}(\cdot ), see Theorem 2.1 in the next section. As a theoretical tool in the convergence analysis of the bundle method, the use of h_{\overline{x}}(\cdot ) dates back to literature [11]. However, the method requires that the resulting sequence of descent steps is feasible. The piecewise linear model of the improvement function has been used in feasible direction methods for solving smooth problems [15], while these methods still require that the descent steps are feasible and the improvement function itself is not addressed in the algorithm. Infeasible bundle methods are less common, with the literature [6] and the “stage I–stage II” improvement of the feasible methods in the constrained bundle methods before [3]. In these methods, it is necessary to predict the values of the objective function and parameters, which is obviously not feasible for general functions. In this paper, the infeasible bundle method of Sagastizdbal and Solodov [16] and the Moreau-Yosida regularization idea [13] are employed, where the objective function changes with the generation of descent step, but the adjusted model is still an upper and lower approximation of the updated objective function.

The full paper is structured as follows: Section 2 performs the problem transformation and gives the upper and lower approximation models of the objective function and its related properties, where the objective function is constructed by using the bundle method idea. Section 3 gives the specific infeasible quasi-Newtonian bundle method for solving the convex constrained optimization problem. Section 4 presents the global convergence analysis of the given algorithm. Section 5 gives the convergence results of the BFGS bundle method. The last section is the conclusion.

In the whole text, the following standard notation is used: ⟨x,y⟩=\sum _{i=1}^n{x}_i{y}_i denotes the inner product in Euclidean space and the corresponding norm denoted by \Vert \cdot \Vert, x^{+}=max\{x,0\}, for a subset X of {\mathbb{R}}^n, conv X represents its convex hull and \mathrm{\partial }h(x) represents the sub-differential of the convex function h(x) at the point x.

2 Bundle method

Theorem 2.1 [6]　 Assuming that Problem (1) satisfies the Slater constraint specification, the following conclusions are equivalent:

(1) \overline{x} is the solution of Problem (1);

(2) min\{h_{\overline{x}}(y)|y\in {\mathbb{R}}^n\}=h_{\overline{x}}(\overline{x})=0;

(3) 0\in \mathrm{\partial }{h}_{\overline{x}}(\overline{x}).

In view of the conclusion of Theorem 2.1, for a given x\in {\mathbb{R}}^n, consider the following problem:

(2)\begin{array}{r}min\{h_x(y)\mid y\in {\mathbb{R}}^n\}.\end{array}

According to the Moreau-Yosida regularization conclusion, the above problem is equivalent to

(3)\begin{array}{r}min\{H_x^M(x)\mid x\in {\mathbb{R}}^n\},\end{array}

where

(4)\begin{array}{r}{H}_x^M(x)=min\left\{h_x(y)+\frac{1}{2}\Vert y-x{\Vert }_M^2\mid y\in {\mathbb{R}}^n\right\},\end{array}

M is an n\times n symmetric positive definite matrix, \Vert \cdot {\Vert }_M=⟨\cdot ,M\cdot ⟩. (4) is called the Moreau-Yosida regularization of h_x(y) associated with M. It is well known that H_{\cdot }^M(\cdot ) is a continuous differentiable convex function defined on {\mathbb{R}}^n, although h_x(\cdot ) may be non-differentiable. The derivative of H_{\cdot }^M(\cdot ) at point x is G_x^M(x):=▽H_x^M(x)=M(x-p(x))\in \mathrm{\partial }{h}_x(p(x)), where p(x) is the unique optimal solution to Problem (4). Further, G_x^M(x) is a global Lipschitz continuous function which module is \Vert M\Vert. \overline{x} is a minimal point of h_x(\cdot ) when and only when G_x^M(\overline{x})=0 and \overline{x}=p(\overline{x}). Let y=x+d. Then (4) is equivalent to

(5)\begin{array}{r}{H}_x^M(x)=min\left\{h_x(x+d)+\frac{1}{2}\Vert d{\Vert }_M^2\mid d\in {\mathbb{R}}^n\right\}.\end{array}

The following considers the approximation of h_x(x+d) using the bundle method. Assume that the information is already available in the bundle (e_{f_i},e_{c_i},g_f^i\in {\mathrm{\partial }}_{e_{f_i}}f(x),g_c^i\in {\mathrm{\partial }}_{e_{c_i}}c(x)). The linearization errors of f_i and c_i are defined as e_{f_i}=f(x)-f(y^i)-⟨g_f^i,x-y^i⟩ and e_{c_i}=c(x)-c(y^i)-⟨g_c^i,x-y^i⟩, then bundle meets .

Lemma 2.1 [16]　 For each i\in B_i^{\text{oracle}}, define

Then e_i⩾0,g_{h_x}^i\in {\mathrm{\partial }}_{e_i}{h}_x(x).

From Lemma 2.1, we know that h_x(y)⩾h_x(x)+⟨g_{h_x}^i,y-x⟩-e_i,\mathrm{\forall }y\in {\mathbb{R}}^n. Let y=x+d. Then h_x(x+d)⩾c^{+}(x)+max\{⟨g_{h_x}^i,d⟩-e_i\mid i\in B_l^{\text{oracle}}\}=:{\check{h}}_x(x+d). Furthermore, make

(6)\begin{array}{rl}{\check{H}}_x^M(x)& =min\left\{{\check{h}}_x(x+d)+\frac{1}{2}\Vert d{\Vert }_M^2\mid d\in {\mathbb{R}}^n\right\}\\ & =c^{+}(x)+min\left\{max\left\{-e_i+⟨g_{h_x}^i,d⟩\mid i\in B_l^{\text{oracle}}\right\}+\frac{1}{2}{d}^{\text{T}}Md\mid d\in {\mathbb{R}}^n\right\}.\end{array}

Let d(x) be the optimal solution of (6). Let v(x)=max\{-e_i+⟨g_{h_x}^i,d(x)⟩\mid i\in B_i^{\text{oracle}}\}. Then {\check{H}}_x^M(x)=c^{+}(x)+v(x)+\frac{1}{2}d(x{)}^{\text{T}}Md(x). Let a(x)=x+d(x). It is an approximation of the real solution of (5). Then let {\hat{H}}_x^M(x)=h_x(a(x))+\frac{1}{2}d(x{)}^{\text{T}}Md(x),x\in {\mathbb{R}}^n. For the upper and lower approximation functions of H_x^M(x), the following conclusion holds.

Lemma 2.2　 \mathrm{\forall }x\in {\mathbb{R}}^n:

(i) {\hat{H}}_x^M(x)⩾H_x^M(x)⩾{\check{H}}_x^M(x);

(ii) {\hat{H}}_x^M(x)=H_x^M(x)\iff a(x)=p(x).

Proof　(i) Due to the fact that for x\in {\mathbb{R}}^n,h_x(x+d)⩾{\check{h}}_x(x+d),\mathrm{\forall }d\in {\mathbb{R}}^n, then H_x^M(x)⩾{\check{H}}_x^M(x). By the definition of a(x) and the optimal solution d(x) of (6), we have h_x(a(x))+\frac{1}{2}d(x{)}^{\text{T}}Md(x)⩾{\check{h}}_x(a(x))+\frac{1}{2}d(x{)}^{\text{T}}Md(x). Thus, {\hat{H}}_x^M(x)⩾H_x^M(x).

(ii) Apparently holds.□

Let \epsilon (x)={\hat{H}}_x^M(x)-{\check{H}}_x^M(x). Then \epsilon (x)⩾0. Let C be a given positive number and \delta (x) be a given positive number during each inner loop iteration. If

(7)\begin{array}{r}\epsilon (x)⩽\delta (x)min\left\{d(x{)}^{\text{T}}Md(x),C\right\},\end{array}

then a(x) is considered to be an approximation of p(x). If (7) does not hold, let y^{i+1}=x+d(x) and compute e_{f_{i+1}},e_{c_{i+1}},f_{i+1},c_{i+1},g_f^{i+1},g_c^{i+1}, and then e_{i+1},g_{h_x}^{i+1}. According to Lemma 2.1, add the information of the new trial point y^{i+1} to the bundle, construct a new approximate model for h_x(x+d), and then resolve (6) to find the new d(x),\epsilon (x), and check it with (7). If condition (7) never holds, using a similar proof in [13], the following two lemmas can be obtained.

Lemma 2.3　 Assuming that x is not an optimal solution for h_x(\cdot ), and (7) is never satisfied in the sub-algorithm, then \epsilon (x)\to 0.

Lemma 2.4　 Let {\tilde{G}}_x^M(x)=M(x-a(x))=-Md(x). Then

(8)\begin{array}{r}{\Vert G_x^M(x)-{\tilde{G}}_x^M(x)\Vert }_{M^{-1}}=\Vert p(x)-a(x){\Vert }_M⩽\sqrt{2\epsilon (x)},\end{array}

(9)\begin{array}{r}\Vert G_x^M(x)-{\tilde{G}}_x^M(x)\Vert ⩽\sqrt{2\epsilon (x)\Vert M\Vert }.\end{array}

The following lemma shows that each inner loop iteration must stop in a finite number of steps, thus ensuring the smooth execution of the entire algorithm.

Lemma 2.5　 If x is not an optimal solution of h_x(\cdot ), then after a finite number of computational iterations of the sub-problem (6), it is always possible to find a solution d(x) of the sub-problem such that (7) holds.

Proof　Assuming that the conclusion does not hold, we will keep adding information of the new trial point y^{i+1} to the bundle, leading to i\to \mathrm{\infty }. By Lemma 2.3, we have \epsilon (x)\to 0. By Lemma 2.4, \Vert G_x^M(x)-{\tilde{G}}_x^M(x)\Vert \to 0. Since x is not an optimal solution for h_x(\cdot ), there exists a positive {\delta }_0 and a positive integer i_0 such that \Vert {\tilde{G}}_x^M(x)\Vert ⩾{\delta }_0. This holds for all i⩾i_0. By d(x{)}^{\text{T}}Md(x)={\left({\tilde{G}}_x^M(x)\right)}^{\text{T}}{M}^{-1}{\tilde{G}}_x^M(x) and the fact that (7) does not hold, there is \epsilon (x)>\delta (x)min\{\frac{{\delta }_0^2}{\Vert M\Vert },N\},i⩾i_0, which contradicts \epsilon (x)\to 0.□

For the decreasing step iteration point x^{k+1} generated by the algorithm, the following lemma ensures that the new model {\check{h}}_{x^{k+1}}(\cdot ) is still a valid lower approximation to the new objective function h_{x^{k+1}}(\cdot ) after the objective function is updated from {\check{h}}_{x^k}(\cdot ) to {\check{h}}_{x^{k+1}}(\cdot ) in (6), thus ensuring the next iteration to proceed smoothly.

Lemma 2.6 [16]　Assume that the algorithm produces descent step iteration point x^{k+1}, its linearization error associated with x^{k+1} is defined as follows:

\begin{array}{r}\begin{array}{rl}& e_{f_i}^{k+1}=e_{f_i}^k+f\left(x^{k+1}\right)-f\left(x^k\right)+⟨g_f^i,x^k-x^{k+1}⟩,\\ & e_{c_i}^{k+1}=e_{c_i}^k+c\left(x^{k+1}\right)-c\left(x^k\right)+⟨g_c^i,x^k-x^{k+1}⟩.\end{array}\end{array}

Then g_{h_{k+1}}^i\in {\mathrm{\partial }}_{e_i^{k+1}}{h}_{k+1}(x^{k+1}), where e_i^{k+1}⩾0, and g_{h_{k+1}}^i is obtained by converting x to the corresponding x^{k+1} in the result of Lemma 2.1.

3 Specific algorithms

Algorithm 3.1　Infeasible quasi-Newton bundle methods for convex constrained optimization problems.

Step 1 (Initial step)　 \sigma ,\rho ,C are given constants and is a sequence of positive numbers satisfying . x^0 is the initial estimated solution, B_0 is an n\times n symmetric positive definite matrix. Let k=0, choose d^0,{\epsilon }_0 such that {\epsilon }_0⩽{\delta }_{0}min\{(d^0{)}^{\text{T}}M{d}^0,C\}. Start the iterative process of the sub-problem with i=1, u^1=x^0.

Step 2 (Compute the search direction)　If \Vert {\tilde{G}}_{x^k}^M(x^k)\Vert =0, the algorithm stops and x^k is an optimal solution. Otherwise, calculate

(10)\begin{array}{r}{s}^k=-B_k^{-1}{\tilde{G}}_{x^k}^M(x^k).\end{array}

Step 3 (Line search)　Starting from l=0, let j_k be the smallest non-negative integer l such that

(11)\begin{array}{r}{\check{H}}_{x^k+{\rho }^l{s}^k}^M\left(x^k+{\rho }^l{s}^k\right)⩽{\hat{H}}_{x^k}^M\left(x^k\right)+\sigma {\rho }^l{\left(s^k\right)}^{\text{T}}{\tilde{G}}_{x^k}^M\left(x^k\right),\end{array}

where {\check{H}}_{x^k+{\rho }^l{s}^k}^M\left(x^k+{\rho }^l{s}^k\right) satisfies

(12)\begin{array}{rl}& {\hat{H}}_{x^k+{\rho }^l{s}^k}^M\left(x^k+{\rho }^l{s}^k\right)-{\check{H}}_{x^k+{\rho }^l{s}^k}^M\left(x^k+{\rho }^l{s}^k\right)\\ & ⩽{\delta }_{k+1}min\{(d(x^k+{\rho }^l{s}^k){)}^{\text{T}}Md(x^k+{\rho }^l{s}^k),C\}.\end{array}

Setting {\tau }_k:={\rho }^{j_k},x^{k+1}:=x^k+{\tau }_k{s}^k.

Step 4 (Correct the quasi-Newton matrix)　Let \mathrm{\Delta }{x}^k=x^{k+1}-x^k,\mathrm{\Delta }{y}^k={\tilde{G}}_{x^{k+1}}^M(x^{k+1})-{\tilde{G}}_{x^k}^M(x^k). If (\mathrm{\Delta }{x}^k{)}^{\text{T}}\mathrm{\Delta }{y}^k>0. Correct B_k to B_{k+1}, so that B_{k+1} is a symmetric positive definite matrix and satisfies the quasi-Newton equation B_{k+1}\mathrm{\Delta }{x}^k=\mathrm{\Delta }{y}^k. Otherwise, set B_k=M. Let k=k+1, and go back to Step 2.

Note 1　Take not of the stopping criterion in Step 2. If \Vert {\tilde{G}}_{x^k}^M(x^k)\Vert =0, then we can infer that d(x^k)=0 from {\tilde{G}}_{x^k}^M(x^k)=-Md(x^k), and from (7), we know that \epsilon (x^k)=0. Then we get G_{x^k}^M(x^k)=0 and x^k=p(x^k),0\in \mathrm{\partial }{h}_{x^k}(x^k) from (9). From Theorem 2.1, we know that x^k is the optimal solution to the original Problem (1).

Note 2　In Step 3, the process of finding j_k requires repeatedly solving the sub-problem (6) to find a descent step iteration point x^{k+1} that satisfies conditions (11) and (12). During this process, the objective function of the sub-problem (6) does not change (Step 1). Once j_k is found, Step 3 is skipped and the next iteration is performed. At this time, the objective function of (6) is changed (the so-called descent step), and the original approximate model {\check{h}}_{x^k}(\cdot ) for h_{x^k}(\cdot ) becomes the approximate model {\check{h}}_{x^{k+1}}(\cdot ) for h_{x^{k+1}}(\cdot ), i.e., not only the approximate model is changed, but the original objective function that is being approximated is also changed. According to Lemma 2.1 and Lemma 2.6, {\check{h}}_{x^{k+1}}(\cdot ) is still a valid approximation of h_{x^{k+1}}(\cdot ), which ensures the smooth iteration of the algorithm. This is different from the previous bundle methods, in which the objective function in the sub-problem of the bundle method is fixed and does not change with the generation of descent steps.

Note 3　As the iteration proceeds, the number of elements in the bundle B_l^{\text{oracle}}increases, and it is necessary to selectively delete and compress the bundle when the number of elements in the bundle becomes large. Let the optimal solution of the dual problem of Problem (6) be \alpha \in {\mathbb{R}}^{|B_l^{\text{oracle}}|} with {\alpha }_i⩾0,\sum _{i\in B_l^{\text{oracle}}}{\alpha }_i=1. It is well known that the optimal solution of Problem (6) can be expressed according to the optimal solution \alpha of its dual problem. However, in this formulation, only the sub-gradient g_{h_x}^i corresponding to {\alpha }_i>0 takes effect, and this is also in effect for the linearization error e_i. Thus, when |B_l^{\text{oracle}}|=n^{max}, we perform the following selection deletion and aggregation into techniques, where n^{max} represents the upper limit of the number of elements that can be stored in the bundle:

— Select the point pair (e_i,g_{h_x}^i) corresponding to {\alpha }_i=0, which can be removed directly from the bundle.

— If there are still many remaining point pairs corresponding to {\alpha }_i>0, their information is compressed into a point pair in the form of a convex combination (with {\alpha }_i as the convex combination factor). A new affine function is constructed using this point pair, which is inserted into the model. After which, the point pairs corresponding to {\alpha }_i>0 can be arbitrarily selected until the number of elements in bundle B_l^{\text{oracle}} is less than n^{max}. This technique can both preserves the information of the positive point pairs in the bundle and ensures that the number of elements in the bundle is within the specified range.

4 Global convergence analysis

The following theorem is a reasonable explanation for the existence of j_k in each iteration of the algorithm.

Theorem 4.1　 If x^k is not the minimal value point of h_{x^k}(\cdot ), then there exists {\overline{\tau }}_k>0 such that

(13)\begin{array}{r}{\check{H}}_{x^k+\tau s^k}^M\left(x^k+\tau s^k\right)⩽{\hat{H}}_{x^k}^M\left(x^k\right)+\sigma \tau {\left(s^k\right)}^{\text{T}}{\tilde{G}}_{x^k}^M\left(x^k\right)\end{array}

holds for all \tau \in (0,{\overline{\tau }}_k), where {\check{H}}_{x^k+\tau s^k}^M\left(x^k+\tau s^k\right) satisfies

(14)\begin{array}{rl}{\hat{H}}_{x^k+\tau s^k}^M\left(x^k+\tau s^k\right)& -{\breve{H}}_{x^k+\tau s^k}^M\left(x^k+\tau s^k\right)\\ & ⩽{\delta }_{k+1}min\{(d(x^k+\tau s^k){)}^{\text{T}}Md(x^k+\tau s^k),C\}.\end{array}

Proof　Since x^k is not a minimal point of h_{x^k}(\cdot ), there exists a positive number {\tilde{\tau }}_k>0 such that \mathrm{\forall }\tau \in \left(0,\tilde{\tau }\right],x^k+\tau s^k is not a minimal point of H_{\cdot }^M(\cdot ) either. By Lemma 2.5, for each \tau \in \left(0,\tilde{\tau }\right], there is always a d(x^k+\tau s^k) such that (14) holds. If {\hat{H}}_{x^k}^M(x^k)=H_{x^k}^M(x^k), by Lemma 2.2, we have a(x^k)=p(x^k), so {\tilde{G}}_{x^k}^M(x^k)=M(x^k-a(x^k))=M(x^k-p(x^k))=G_{x^k}^M(x^k). Since x^k is not an optimal solution, (10) implies that . Also, since H_{\cdot }^M(\cdot ) is continuously differentiable and , there is a positive number {\overline{\tau }}_k,{\overline{\tau }}_k⩽{\tilde{\tau }}_k such that \mathrm{\forall }\tau \in \left(0,{\overline{\tau }}_k\right]. While we have H_{x^k+\tau s^k}^M\left(x^k+\tau s^k\right)⩽H_{x^k}^M+\sigma \tau {\left(s^k\right)}^{\text{T}}{G}_{x^k}^M\left(x^k\right), which implies that (13) holds. If {\hat{H}}_{x^k}^M(x^k)>H_{x^k}^M(x^k), then there also exists {\overline{\tau }}_k>0 (sufficiently small) such that {\check{H}}_{x^k+\tau s^k}^M\left(x^k+\tau s^k\right)⩽H_{x^k+\tau s^k}^M\left(x^k+\tau s^k\right)\stackrel{\tau \to 0}{⟶}{H}_{x^k}^M\left(x^k\right)⩽H_{x^k}^M\left(x^k\right) + \frac{1}{2}\left({\hat{H}}_{x^k}^M\left(x^k\right)-H_{x^k}^M\left(x^k\right)\right)⩽{\hat{H}}_{x^k}^M\left(x^k\right)+\sigma \tau {\left(s^k\right)}^{\text{T}}{\tilde{G}}_{x^k}^M\left(x^k\right),\mathrm{\forall }\tau \in \left(0,{\overline{\tau }}_k\right], and thus (13) also holds.□

Since it is assumed that , there exists a constant A>0 such that

(15)\begin{array}{r}\sum _{k=0}^{\mathrm{\infty }}{\delta }_k⩽A.\end{array}

Introducing the notation D_k\{x\in {\mathbb{R}}^n\mid H_{x^k}^M(x)⩽H_{x^0}^M(x^0)+2AC\}, if we assume that both f and c are lower bounded functions, then \{D_k{\}}_{k⩾0} is bounded.

Lemma 4.1　 For all k⩾0, we have

(16)\begin{array}{r}{H}_{x^{k+1}}^M\left(x^{k+1}\right)⩽H_{x^k}^M\left(x^k\right)+C\left({\delta }_k+{\delta }_{k+1}\right),x^k\in D_k\text{.}\end{array}

Proof　According to the algorithm criterion, for all k⩾0, there is

\begin{array}{r}\begin{array}{rl}{H}_{x^{k+1}}^M\left(x^{k+1}\right)& ⩽{\check{H}}_{x^{k+1}}^M\left(x^{k+1}\right)+C{\delta }_{k+1}\\ & ⩽{\hat{H}}_{x^k}^M\left(x^k\right)+\sigma {\rho }^{j_k}{\left(s^k\right)}^{\text{T}}{\tilde{G}}_{x^k}^M\left(x^k\right)+C{\delta }_{k+1}\\ & ={\hat{H}}_{x^k}^M\left(x^k\right)-\sigma {\rho }^{j_k}{\tilde{G}}_{x^k}^M{\left(x^k\right)}^{\text{T}}{B}_k^{-1}{\tilde{G}}_{x^k}^M\left(x^k\right)+C{\delta }_{k+1}\\ & ⩽{\hat{H}}_{x^k}^M\left(x^k\right)+C{\delta }_{k+1}\\ & ⩽{\check{H}}_{x^k}^M\left(x^k\right)+C{\delta }_k+C{\delta }_{k+1}.\end{array}\end{array}

Thus, for all k⩾0, then H_{x^{k+1}}^M(x^{k+1})⩽H_{x^0}^M(x^0)+2AC. Since x^0\in D_0, so x^k\in D_k,\mathrm{\forall }k⩾0. The conclusion is confirmed.□

The global convergence results of the algorithm are given below.

Theorem 4.2　 Assuming that f and c are lower bounded and there exist two positive constants a_1 and a_2 such that \Vert B_k\Vert ⩽a_1 and \Vert B_k^{-1}\Vert ⩽a_2 hold for all k⩾0, the algorithm produces any cluster of the sequence \{x^k\} that is a minimal value point for Problem (1).

Proof　By Lemma 4.1 and the assumption that f and c are lower bounded, combined with (15) and (16), there exists H_{\ast }^M such that

\begin{array}{r}\underset{k\to \mathrm{\infty }}{lim}{H}_{x^k}^M(x^k)=H_{\ast }^M.\end{array}

Since {\delta }_k\to 0, according to Lemma 2.1 and the algorithm criterion, we have {\epsilon }_k\to 0 and \underset{k\to \mathrm{\infty }}{lim}{\hat{H}}_{x^k}^M\left(x^k\right)=\underset{k\to \mathrm{\infty }}{lim}{\check{H}}_{x^k}^M\left(x^k\right)=H_{\ast }^M. Thus, according to the assumptions on B_k,

(17)\begin{array}{r}\underset{k\to \mathrm{\infty }}{lim}{\tau }_k\Vert {\tilde{G}}_{x^k}^M(x^k){\Vert }^2=0.\end{array}

Let \overline{x} be any cluster-point of \{x^k\}, and \{x^k{\}}_{k\in K} be a sub-sequence converging to \overline{x}. According to Lemma 2.3,

(18)\begin{array}{r}\underset{k\to \mathrm{\infty },k\in K}{lim}{\tilde{G}}_{x^k}^M(x^k)=G_{x^k}^M(\overline{x}).\end{array}

If lim\underset{k\to \mathrm{\infty },k\in K}{inf}{\tau }_k>0, then according to (17) and (18), we have G_{\overline{x}}^M(\overline{x})=0, which means that \overline{x} is the optimal solution to the original problem. Otherwise, lim\underset{k\to \mathrm{\infty },k\in K}{inf}{\tau }_k=0. Assuming that {\tau }_k\to 0(k\to \mathrm{\infty },k\in K), according to the line search criterion and Lemma 2.1, we have

(19)\begin{array}{r}\frac{H_{x^k+{\rho }^{j_k-1}{s}^k}^M\left(x^k+{\rho }^{j_k-1}{s}^k\right)-H_{x^k}^M\left(x^k\right)}{{\rho }^{j_k-1}}>\sigma {\left(s^k\right)}^{\text{T}}{\tilde{G}}_{x^k}^M\left(x^k\right).\end{array}

Then by (l8), \{{\tilde{G}}_{x^k}^M(x^k){\}}_{k\in K} is bounded and with the assumption that B_k is bounded, we get \{s^k{\}}_{k\in K} is bounded. We may assume again \underset{k\in K,k\to \mathrm{\infty }}{lim}{s}^k=\overline{s}. Since {\rho }^{j_k-1}\to 0,k\in K,k\to \mathrm{\infty }, we get {\overline{s}}^{\text{T}}{G}_{\overline{x}}^M(\overline{x})⩾\sigma {\overline{s}}^{\text{T}}{G}_{\overline{x}}^M(\overline{x}). By and the assumption for B_k, {\overline{s}}^{\text{T}}{G}_{\overline{x}}^M(\overline{x})⩽-\frac{1}{c_2}\Vert \overline{s}{\Vert }^2, which implies that {\overline{s}}^{\text{T}}{G}_{\overline{x}}^M(\overline{x})=0 as well as \overline{s}=0, we finally get G_{\overline{x}}^M(\overline{x})=0.□

Note 4　We will discuss the assumption of boundedness of \{B_k\} and \{B_k^{-1}\} in Theorem 4.2 in the next section.

5 Convergence analysis of infeasible BFGS bundle methods

This section mainly considers the convergence of the infeasible BFGS bundle method.

The infeasible BFGS bundle method is obtained by taking B_0=M in the initial step of Algorithm 3.1 and taking the positive sequence \{{\delta }_k{\}}_{k=1}^{\mathrm{\infty }} to satisfy . The selection of B_k is specified in Step 4, and the infeasible BFGS bundle method is obtained as follows.

Algorithm 5.1　Infeasible BFGS bundle method for convex constrained optimization problem.

Step 4 (Correction of the quasi-Newton matrix)　Let \mathrm{\Delta }{x}^k=x^{k+1}-x^k,\mathrm{\Delta }{y}^k={\tilde{G}}_{x^{k+1}}^M(x^{k+1})-{\tilde{G}}_{x^k}^M(x^k). If there exist two constants a_3>0,a_4\in (0,1) such that the iteration of the algorithm satisfies the following two conditions:

(20)\begin{array}{r}{\Vert \mathrm{\Delta }{x}^k\Vert }_M\left(\sqrt{2{\epsilon }_k}+\sqrt{2{\epsilon }_{k+1}}\right)⩽a_3{\left(\mathrm{\Delta }{x}^k\right)}^{\text{T}}\mathrm{\Delta }{y}^k,\end{array}

(21)\begin{array}{r}2{\Vert \mathrm{\Delta }{y}^k\Vert }_M\left(\sqrt{2{\epsilon }_k}+\sqrt{2{\epsilon }_{k+1}}\right)⩽min\left\{a_4,{\delta }_k^{\frac{1}{3}}+{\delta }_{k+1}^{\frac{1}{3}}\right\}{\Vert \mathrm{\Delta }{y}^k\Vert }^2.\end{array}

Let B_{k+1}=B_k-\frac{B_k\mathrm{\Delta }{x}^k{\left(\mathrm{\Delta }{x}^k\right)}^{\text{T}}{B}_k}{{\left(\mathrm{\Delta }{x}^k\right)}^T{B}_k\mathrm{\Delta }{x}^k}+\frac{\mathrm{\Delta }{y}^k{\left(\mathrm{\Delta }{y}^k\right)}^{\mathbf{T}}}{{\left(\mathrm{\Delta }{x}^k\right)}^{\text{T}}\mathrm{\Delta }{y}^k}, otherwise make B_{k+1}=M. Let k=k+1, go back to Step 2.

From Step 4, if (20) holds, then {\left(\mathrm{\Delta }{x}^k\right)}^{\text{T}}\mathrm{\Delta }{y}^k>0; if B_k is a positive definition matrix, B_{k+1} is also a positive definition matrix and satisfies the quasi-Newton equation

\begin{array}{r}{B}_{k+1}\mathrm{\Delta }{x}^k=B_k\mathrm{\Delta }{x}^k-\frac{B_k\mathrm{\Delta }{x}^k{\left(\mathrm{\Delta }{x}^k\right)}^{\text{T}}{B}_k\mathrm{\Delta }{x}^k}{{\left(\mathrm{\Delta }{x}^k\right)}^{\text{T}}{B}_k\mathrm{\Delta }{x}^k}+\frac{\mathrm{\Delta }{y}^k{\left(\mathrm{\Delta }{y}^k\right)}^{\text{T}}\mathrm{\Delta }{x}^k}{{\left(\mathrm{\Delta }{x}^k\right)}^{\text{T}}\mathrm{\Delta }{y}^k}=\mathrm{\Delta }{y}^k.\end{array}

As we all know, if h_x(\cdot ) is a strongly convex function on {\mathbb{R}}^n, then H_{\cdot }^M(\cdot ) is also strongly convex on {\mathbb{R}}^n [10]. Now assume that H_{\cdot }^M(\cdot ) is strongly convex on {\left\{D_k\right\}}_{k⩾0}. Then h_x(\cdot ) has a unique minima \overline{x} on {\left\{D_k\right\}}_{k⩾0}. Introduce the label K:=\{0\}\cup \{i\mid (20) or (21) does not hold for k=i-1\}\equiv \left\{k_0,k_1,\dots ,k_j,\dots \right\}. The following theorem shows that the assumption of boundedness of \left\{B_k\right\} and \left\{B_k^{-1}\right\} in Theorem 4.2 holds under certain conditions.

Theorem 5.1　 Suppose that H_{\cdot }^M(\cdot ) is strongly convex on {\left\{D_k\right\}}_{k⩾0}. The sequences \left\{x^k\right\} and \left\{B_k\right\} are generated by Algorithm 5.1. Suppose that \mathrm{\forall }k⩾0,\mathrm{\Delta }{x}^k and \mathrm{\Delta }{y}^k satisfy

\begin{array}{r}\frac{\Vert \mathrm{\Delta }{\overline{y}}^k-M_{\ast }\mathrm{\Delta }{x}^k\Vert }{\Vert \mathrm{\Delta }{x}^k\Vert }⩽{\delta }_k^{\mathrm{\prime }},\end{array}

where \mathrm{\Delta }{\overline{y}}^k=G_{x^{k+1}}^M\left(x^{k+1}\right)-G_{x^k}^M\left(x^k\right),M_{\ast } is some symmetric positive definite matrix and {\left\{{\delta }_k^{\mathrm{\prime }}\right\}}_{k=0}^{\mathrm{\infty }} is a sequence of positive numbers satisfying . Then \Vert B_k\Vert and \Vert B_k^{-1}\Vert are bounded.

Proof　If (20) and (21) hold, we have

\begin{array}{r}\begin{array}{rl}{\left(\mathrm{\Delta }{x}^k\right)}^{\text{T}}\mathrm{\Delta }{y}^k⩾& {\left(\mathrm{\Delta }{x}^k\right)}^{\text{T}}\mathrm{\Delta }{\overline{y}}^k-{\Vert \mathrm{\Delta }{x}^k\Vert }_M{\Vert \mathrm{\Delta }{y}^k-\mathrm{\Delta }{\overline{y}}^k\Vert }_{M^{-1}}\\ ⩾& {\left(\mathrm{\Delta }{x}^k\right)}^{\text{T}}\mathrm{\Delta }{\overline{y}}^k-{\Vert \mathrm{\Delta }{x}^k\Vert }_M({\Vert {\tilde{G}}_{x^k}^M\left(x^k\right)-G_{x^k}^M\left(x^k\right)\Vert }_{M^{-1}}\\ & +{\Vert {\tilde{G}}_{x^{k+1}}^M\left(x^{k+1}\right)-G_{x^{k+1}}^M\left(x^{k+1}\right)\Vert }_{M^{-1}})\\ ⩾& {\left(\mathrm{\Delta }{x}^k\right)}^{\text{T}}\mathrm{\Delta }{\overline{y}}^k-{\Vert \mathrm{\Delta }{x}^k\Vert }_M\left(\sqrt{2{\epsilon }_k}+\sqrt{2{\epsilon }_{k+1}}\right).\end{array}\end{array}

After organizing, we get

(22)\begin{array}{r}{\left(\mathrm{\Delta }{x}^k\right)}^{\text{T}}\mathrm{\Delta }{y}^k⩾\frac{1}{1+a_3}{\left(\mathrm{\Delta }{x}^k\right)}^{\text{T}}\mathrm{\Delta }{\overline{y}}^k.\end{array}

In addition,

\begin{array}{r}\begin{array}{rl}{\Vert \mathrm{\Delta }{\overline{y}}^k\Vert }^2& ⩾{\Vert \mathrm{\Delta }{y}^k\Vert }^2-2{\Vert \mathrm{\Delta }{y}^k\Vert }_M{\Vert \mathrm{\Delta }{\overline{y}}^k-\mathrm{\Delta }{y}^k\Vert }_{M^{-1}}\\ & ⩾{\Vert \mathrm{\Delta }{y}^k\Vert }^2-2{\Vert \mathrm{\Delta }{y}^k\Vert }_M\left(\sqrt{2{\epsilon }_k}+\sqrt{2{\epsilon }_{k+1}}\right).\end{array}\end{array}

We reorganized and get

(23)\begin{array}{r}{\Vert \mathrm{\Delta }{\overline{y}}^k\Vert }^2⩾\left(1-a_4\right){\Vert \mathrm{\Delta }{y}^k\Vert }^2.\end{array}

For i\in K, we have B_i=M. The i\notin K,B_i is obtained by B_{i-1} using the BFGS correction formula. Assuming that K has an infinite number of elements, for , where for a certain j⩾1,k_{j-1},k_j\in K, then (20) and (21) hold, i.e.,

\begin{array}{r}{\left(\mathrm{\Delta }{x}^k\right)}^{\text{T}}\mathrm{\Delta }{y}^k⩾\frac{1}{1+a_3}{\left(\mathrm{\Delta }{x}^k\right)}^{\text{T}}\mathrm{\Delta }{\overline{y}}^k>0.\end{array}

Then, by Lemma 2.4, (20) and (21), for all , we have

\begin{array}{r}\begin{array}{rl}\frac{\Vert \mathrm{\Delta }{y}^k-M_{\ast }\mathrm{\Delta }{x}^k\Vert }{\Vert \mathrm{\Delta }{x}^k\Vert }& ⩽{\delta }_k^{\mathrm{\prime }}+\frac{\Vert \mathrm{\Delta }{y}^k-\mathrm{\Delta }{\overline{y}}^k\Vert }{\Vert \mathrm{\Delta }{x}^k\Vert }\\ & ⩽{\delta }_k^{\mathrm{\prime }}+\frac{\Vert M\Vert \left(\sqrt{2{\epsilon }_k}+\sqrt{2{\epsilon }_{k+1}}\right)}{\Vert \mathrm{\Delta }{x}^k\Vert }\\ & ⩽{\delta }_k^{\mathrm{\prime }}+\frac{\sqrt{\Vert M\Vert \Vert M^{-1}\Vert }}{2\sqrt{1-a_4}}\left({\delta }_k^{\frac{1}{3}}+{\delta }_{k+1}^{\frac{1}{3}}\right)\frac{\Vert \mathrm{\Delta }{\overline{y}}^k\Vert }{\Vert \mathrm{\Delta }{x}^k\Vert }\\ & ⩽{\delta }_k^{\mathrm{\prime }}+\frac{\sqrt{\Vert M{\Vert }^3\Vert M^{-1}\Vert }}{2\sqrt{1-a_4}}\left({\delta }_k^{\frac{1}{3}}+{\delta }_{k+1}^{\frac{5}{2}}\right)\\ & =:{\overline{\delta }}_k.\end{array}.\end{array}

Then for all k, where k satisfying , there are

\begin{array}{r}\frac{\Vert \mathrm{\Delta }{y}^k-M_{\ast }\mathrm{\Delta }{x}^k\Vert }{\Vert \mathrm{\Delta }{x}^k\Vert }⩽{\overline{\delta }}_k.\end{array}

By assumption , we have . Then according to [2, Theorem 3.2], it follows that for all k satisfying , \Vert B_k\Vert and \Vert B_k^{-1}\Vert are bounded depending on B_{k_{j-1}}. Since B_{k_j}=M for all j⩾0, then the sequences \Vert B_k\Vert and \Vert B_k^{-1}\Vert are bounded. When the number of elements in K is finite, a similar proof method can be used to obtain that \Vert B_k\Vert and \Vert B_k^{-1}\Vert are bounded.□

The following results show that without the assumption of boundedness of \Vert B_k\Vert and \Vert B_k^{-1}\Vert, it is still possible to prove that the iterative sequence generated by the algorithm still converges to the optimal solution of Problem (1) for some special problems.

Theorem 5.2　 Suppose H_{\cdot }^M(\cdot ) is strongly convex on {\left\{D_k\right\}}_{k⩾0}, and the sequence \left\{x^k\right\},\left\{B_k\right\} are generated by Algorithm 5.1 with x^k\ne \overline{x}, \mathrm{\forall }k⩾0. Then \left\{x^k\right\} R-linearly converges to the unique solution \overline{x} of Problem (1) further, with

Proof　The division K is both a finite set and an infinite set. In either case, we can use the strong convexity of H_{\cdot }^M(\cdot ), the global continuity of G_{\cdot }^M(\cdot ) modulo \Vert M\Vert, the positive definite correction criterion of B_k, Lemma 2.2, and to obtain that there exists the constant \theta \in (0,1),\overline{N}\in (0,\mathrm{\infty }) and a positive integer \overline{k} such that for all k⩾\overline{k}, there is

\begin{array}{r}{H}_{x^{k+1}}^M\left(x^{k+1}\right)-H_{\overline{x}}^M(\overline{x})⩽\overline{N}{\left({\theta }^{\frac{1}{2}}\right)}^{k-\overline{k}+1}\left(H_{x^k}^M\left(x^k\right)-H_{\overline{x}}^M(\overline{x})\right).\end{array}

By the strong convexity of H_{\cdot }^M(\cdot ) and [1, Lemma 4.3], for all k⩾0, there exists a constant \alpha >0 such that

\begin{array}{r}\frac{1}{2}\alpha {\Vert x^k-\overline{x}\Vert }^2⩽H_{x^k}^M\left(x^k\right)-H_{\overline{x}}^M(\overline{x})⩽\frac{\alpha }{2}{\Vert G_{x^k}^M\left(x^k\right)\Vert }^2.\end{array}

So, there is

\begin{array}{r}\begin{array}{rl}\sum _{k=\overline{k}}^{\mathrm{\infty }}{\Vert x^k-\overline{x}\Vert }^2& ⩽{\left(\frac{2}{\alpha }\right)}^{\frac{1}{2}}\sum _{k=\overline{k}}^{\mathrm{\infty }}{\left(H_{x^k}^M\left(x^k\right)-H_{\overline{x}}^M(\overline{x})\right)}^{\frac{1}{2}}\\ & ⩽{\left[\frac{2\overline{N}\left(H_{x^k}^M\left(x^k\right)-H_{\overline{x}}^M(\overline{x})\right)}{\alpha }\right]}^{\frac{1}{2}}\sum _{k=\overline{k}}^{\mathrm{\infty }}{\left({\theta }^{\frac{1}{4}}\right)}^{k-\overline{k}}\\ & ⩽\mathrm{\infty }.\end{array}\end{array}

The conclusion is proofed.□

6 Conclusions

In this paper, we combine the infeasible bundle method of Sagastizdbal with the Moreau-Yosida regularization idea to construct a quasi-Newton bundle method for solving convex constrained optimization problems. The global convergence result of the algorithm is derived without guaranteeing that every iteration point is feasible. Although the method is essentially similar to the unconstrained bundle method, it differs from the previous methods in that the objective function of the subproblem to be solved at each iteration changes with the generation of descent steps, which requires necessary adjustments to ensure the convergence of the algorithm. For the special BFGS bundle method, we briefly discuss the conditions that guarantee the boundedness of the quasi-Newton matrix.

References

[1] Bonnans J F, Gilbert J Ch, Lemaréchal C, Sagastizábal C A. A family of variable metric proximal methods. Math Program 1995; 68(1/2/3): 15–47

[2] Byrd R H, Nocedal J. A tool for the analysis of quasi-Newton methods with application to unconstrained minimization. SIAM J Numer Anal 1989; 26(3): 727–739

[3] FletcherRLeyfferS. A bundle filter method for nonsmooth nonlinear optimization. Numerical Analysis Report NA/195, Dundee: Department of Mathematics, The University of Dundee, 1999

[4] Frangioni A. Generalized bundle methods. SIAM J Optim 2002; 13(1): 117–156

[5] Hare W, Sagastizábal C A. A redistributed proximal bundle method for nonconvex optimization. SIAM J Optim 2010; 20(5): 2442–2473

[6] KiwielK C. Methods of Descent for Nondifferentiable Optimization. Lecture Notes in Math, Vol 1133. Berlin: Springer-Verlag, 1985

[7] Kiwiel K C. Exact penalty functions in proximal bundle methods for constrained convex nondifferentiable minimization. Math Program 1991; 52(1/2/3): 285–302

[8] Kiwiel K C. A proximal bundle method with approximate subgradient linearizations. SIAM J Optim 2006; 16(4): 1007–1023

[9] Kiwiel K C. A proximal-projection bundle method for Lagrangian relaxation, including semidefinite programming. SIAM J Optim 2006; 17(4): 1015–1034

[10] Lemaréchal C, Sagastizábal C A. Practical aspects of the Moreau-Yosida regularization: theoretical preliminaries. SIAM J Optim 1997; 7(2): 367–385

[11] Miffin R. An algorithm for constrained optimization with semismooth functions. Math Oper Res 1977; 2(2): 191–207

[12] MifflinR. A modification and an extension of Lemarechal’s algorithm for nonsmooth minimization. In: Nondifferential and Variational Techniques in Optimization, Math Program Stud, Vol 17. Berlin: Springer-Verlag, 1982, 77–90

[13] Miflin R, Sun D F, Qi L Q. Quasi-Newton bundle type methods for nondifferentiable convex optimization. SIAM J Optim 1998; 8(2): 583–603

[14] OliveiraWSolodovM. A doubly stabilized bundle method for nonsmooth convex optimization. 2013

[15] Pironneau O, Polak E. Rate of convergence of a class of methods of feasible directions. SIAM J Numer Anal 1973; 10(1): 161–174

[16] Sagastizabal C A, Solodov M. An infeasible bundle method for nonsmooth convex constrained optimization without a penalty function or a filter. SIAM J Optim 2005; 16(1): 146–169

RIGHTS & PERMISSIONS

Higher Education Press 2023