Accelerating two projection methods via perturbations with application to Intensity-Modulated Radiation Therapy

TL;DR

This paper introduces two novel perturbation techniques to accelerate projection methods in solving convex feasibility problems, demonstrating significant speed-ups and improved solutions in applications like intensity-modulated radiation therapy.

Contribution

The paper develops and applies two new perturbation strategies leveraging bounded perturbation resilience to enhance projection methods for convex optimization, especially in ill-conditioned problems.

Findings

Perturbed methods accelerate convergence up to 4 times faster.

Perturbed methods achieve lower objective function values by 0.5% to 5.1%.

Effective in applications like IMRT treatment planning.

Abstract

Constrained convex optimization problems arise naturally in many real-world applications. One strategy to solve them in an approximate way is to translate them into a sequence of convex feasibility problems via the recently developed level set scheme and then solve each feasibility problem using projection methods. However, if the problem is ill-conditioned, projection methods often show zigzagging behavior and therefore converge slowly. To address this issue, we exploit the bounded perturbation resilience of the projection methods and introduce two new perturbations which avoid zigzagging behavior. The first perturbation is in the spirit of $k$-step methods and uses gradient information from previous iterates. The second uses the approach of surrogate constraint methods combined with relaxed, averaged projections. We apply two different projection methods in the unperturbed version, as well as the two perturbed versions, to linear feasibility problems along with nonlinear optimization problems arising from intensity-modulated radiation therapy (IMRT) treatment planning. We demonstrate that for all the considered problems the perturbations can significantly accelerate the convergence of the projection methods and hence the overall procedure of the level set scheme. For the IMRT optimization problems the perturbed projection methods found an approximate solution up to 4 times faster than the unperturbed methods while at the same time achieving objective function values which were 0.5 to 5.1% lower.

Tables

Table 1. Table 1: Iterations needed to find a feasible solution with λ ​ ( x ) ≡ 1.9 𝜆 𝑥 1.9 \lambda(x)\equiv 1.9 . SP converges slowly. Using perturbations accelerates the convergence significantly.

	SP	SP + HB			SP + SC
		${\lambda }_k^{HB}=8$	${\lambda }_k^{HB}=80$	${\lambda }_k^{HB}=800$
K	449	58	17	4	4

Table 2. Table 2: Iterations needed to find a feasible solution with λ ​ ( x ) ≡ 1.9 𝜆 𝑥 1.9 \lambda(x)\equiv 1.9 and the control sequence { j 1 ​ ( ν ) } subscript 𝑗 1 𝜈 \{j_{1}(\nu)\} . Not all perturbed versions of CP converge faster than unperturbed CP.

	CP	CP + HB			CP + SC
		${\lambda }_k^{HB}=8$	${\lambda }_k^{HB}=80$	${\lambda }_k^{HB}=800$
K	20	34	26	9	4

Table 3. Table 3: Iterations needed to find a feasible solution with λ ​ ( x ) ≡ 1.9 𝜆 𝑥 1.9 \lambda(x)\equiv 1.9 and the control sequence { j 2 ​ ( ν ) } subscript 𝑗 2 𝜈 \{j_{2}(\nu)\} . Perturbed CP converges faster than unperturbed CP.

	CP	CP + HB			CP + SC
		${\lambda }_k^{HB}=8$	${\lambda }_k^{HB}=80$	${\lambda }_k^{HB}=800$
K	32	29	20	7	3

Table 4. Table 4: Dose evaluation functions of the considered biological structures

	structure	function name	parameters
$f_1$	Left parotis	EUD	$p=2$
$f_2$	Right parotis	EUD	$p=2$
$f_3$	Myelon	EUD	$p=2$
$f_4$	Unclassified healthy tissue	EUD	$p=2$
$f_5$	Tumor volume	Tumor conformity	$p=2,d^{ref}=60$
$g_1$	Tumor volume	lower tail penalty	$L=55$
$g_2$	Tumor volume	upper tail penalty	$U=66$
$g_3$	Myelon	upper tail penalty	$U=45$

Table 5. Table 5: Iteration numbers and lowest objective function values achieved using perturbed and unperturbed SP. Perturbed SP reaches solutions with lower objective function values within less iterations.

	$f_{\text{SP}}^{\ast }$	$f_{\text{SP, HB}}^{\ast }$	$f_{\text{SP, SC}}^{\ast }$	$K_{\text{SP}}$	$K_{\text{SP, HB}}$	$K_{\text{SP, SC}}$
Case 1	3480	3464 (99.54%)	3387 (97.33%)	7159	3437 (48.00%)	2155 (30.10%)
Case 2	2378	2356 (99.07%)	2317 (97.43%)	3523	2340 (66.42%)	1108 (31.45%)
Case 3	3129	3056 (97.67%)	3012 (96.26%)	4773	3458 (72.45%)	1171 (24.53%)
Case 4	3098	2980 (96.19%)	2941 (94.93%)	4496	3898 (86.70%)	1292 (28.74%)

Table 6. Table 6: Iteration numbers and lowest objective function values achieved using perturbed and unperturbed CP. In only two of four cases perturbations are used and achieve lower objective function values than unperturbed CP.

	$f_{\text{CP}}^{\ast }$	$f_{\text{CP, HB}}^{\ast }$	$f_{\text{CP, SC}}^{\ast }$	$K_{\text{CP}}$	$K_{\text{CP, HB}}$	$K_{\text{CP, SC}}$
Case 1	3563	3490 (97.95%)	3420 (95.99%)	6665	6103 (91.57%)	4044 (60.66%)
Case 2	2424	-	-	4484	-	-
Case 3	3178	3083 (97.01%)	3052 (96.04%)	4238	5929 (139.90%)	5035 (118.81%)
Case 4	3093	-	-	5280	-	-

Table 7. Table 7: The number of perturbations used by the perturbed projection methods. For all cases HB perturbation is used more often than SC perturbation.

	SP + HB	SP + SC	CP + HB	CP + SC
Case 1	803	283	789	374
Case 2	534	101	0	0
Case 3	922	166	451	349
Case 4	1074	180	0	0

Table 8. Table 8: The fraction of K SP subscript 𝐾 SP K_{\text{SP}} or K CP subscript 𝐾 CP K_{\text{CP}} needed by the perturbed methods to find a solution with objective function value ≤ f SP ∗ absent subscript superscript 𝑓 ∗ SP \leq f^{\ast}_{\text{SP}} or f CP ∗ subscript superscript 𝑓 ∗ CP f^{\ast}_{\text{CP}} . SC perturbation achieves lower values than HB perturbation, in particular if the perturbed method is SP.

	SP + HB	SP + SC	CP + HB	CP + SC
Case 1	0.4801	0.2122	0.7685	0.4308
Case 2	0.6642	0.2191	-	-
Case 3	0.4324	0.1402	0.9122	0.8554
Case 4	0.4121	0.1417	-	-

Equations

$$f(x)$$

$$g_j(x)\le 0j=1,...,m$$

$$\text{Minimize }t$$

$$\text{s.t. }f(x)-t$$

$$g_j(x)$$

$$\text{Find }x\in C:=\bigcap _{i\in I}{C}_i:=\left\{x\in {\mathbb{R}}^n\mid {\phi }_i(x)\le 0\right\}$$

$${\phi }_1^s$$

$$\text{Find }x\in C=\bigcap _{i\in I}{C}_i^s:=\left\{x\in {\mathbb{R}}^n\mid {\phi }_i^s(x)\le 0\right\}$$

$$x^{k+1}$$

$$p(x^k)=$$

$$x^{k+1}$$

$$p(x^k)$$

$$y^{k+1}$$

$${\tilde{b}}^k$$

$$y^{k+1}$$

$$\{\begin{array}{l}{z}^0=\mathbf{T}(x^0)\\ z^{k+1}=\mathbf{T}(z^k+{\tilde{b}}^k)\end{array}$$

$$z^0+{\tilde{b}}^0$$

$$z^k+{\tilde{b}}^k$$

$$=\mathbf{T}(y^k)+{\tilde{b}}^k=y^{k+1}$$

$$\mathbf{T}(x)=x+\lambda (x)p(x)$$

$$\tilde{c}(y^k,y^{k-1})$$

$$y^{k+1}$$

$$v^k$$

$$z^{k+1}$$

$$v^k$$

$${\beta }_k$$

$$z^{l+1}$$

$$y^{l+1}$$

$$=y^l+(1-{\beta }_l)\lambda (y^l)p(y^l)+{\beta }_l{b}^l=y^l+b^l.$$

$$y^{l+2}$$

$$z^{l+2}$$

$$=y^{l+2}+\lambda (y^{l+2})p(y^{l+2})=\mathbf{T}(y^{l+2})=y^{l+3}$$

Full text

To address this issue, we exploit the bounded perturbation resilience of the projection methods and introduce two new perturbations which avoid zigzagging behavior. The first perturbation is in the spirit of $k$-step methods and uses gradient information from previous iterates. The second uses the approach of surrogate constraint methods combined with relaxed, averaged projections.

We apply two different projection methods in the unperturbed version, as well as the two perturbed versions, to linear feasibility problems along with nonlinear optimization problems arising from intensity-modulated radiation therapy (IMRT) treatment planning. We demonstrate that for all the considered problems the perturbations can significantly accelerate the convergence of the projection methods and hence the overall procedure of the level set scheme. For the IMRT optimization problems the perturbed projection methods found an approximate solution up to 4 times faster than the unperturbed methods while at the same time achieving objective function values which were 0.5 to 5.1% lower.

††journal: Applied Mathematics and Optimization

1 Introduction

The problem we consider in this paper is the general constrained convex optimization problem

[TABLE]

where $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ and $g_{j}:\mathbb{R}^{n}\rightarrow\mathbb{R}$ for $j\in J=\{1,\ldots,m\}$ are convex functions.

In order to solve this problem we translate it to an equivalent epigraph form (2.1) and then use a procedure called the level set scheme introduced in gkrs18 to construct a sequence of convex feasibility problems from it. In this procedure, we define an upper bound on the function $f$ and thereby generate an additional constraint. This constraint together with the constraints given in (1.1) formulates a feasibility problem. Each time the feasibility problem can be solved by a suitable algorithmic operator $\boldsymbol{T}$, we reduce the upper bound on $f$ and use it to formulate the next feasibility problem.

The solution $x^{\ast}_{\varepsilon}$ produced by the level set scheme is an approximation of the solution $x^{\ast}$ of the original optimization problem (1.1) in the sense that $f(x^{\ast}_{\varepsilon})<f(x^{\ast})+\varepsilon_{s}$ for some $\varepsilon_{s}>0$, while at the same time both $x^{\ast}_{\varepsilon}$ and $x^{\ast}$ fulfill $g_{j}(x)\leq 0$ for all $j\in J$. $x^{\ast}_{\varepsilon}$ is therefore called an $\varepsilon_{s}$-optimal solution of (1.1).

In this work we choose the simultaneous subgradient projection method and the cyclic subgradient projection method as options for the operator $T$ used to solve each feasibility problem and compare their behavior. Both methods, as well as many others suitable to solve a convex feasibility problem, are bounded perturbation resilient according to Definition 2 special-issue . This means that certain changes can be made at each iteration step and $T$ still produces a solution of the feasibility problem. One favorable property of bounded perturbation resilient algorithms is the fact that the convergence rate of the unperturbed variant is preserved by the perturbed algorithm (see e.g. (brz17, , Theorem 11 and Corollary 14)).

There exist two kinds of perturbations, which differ in the way the changes are applied. They are called inner and outer perturbations. Combettes in combettes01 studied the concept of perturbation resilience of projection methods using inner perturbations, which they and others refer to as convergence under summable errors tw13 . Algorithms using outer perturbations are also referred to as inexact methods sv12 , because the perturbation vectors can be interpreted e.g. as calculation error if the evaluation of $\boldsymbol{T}(x)$ cannot be performed in an exact way. These kind of perturbations are used e.g. in dgj18 .

One way to exploit bounded perturbation resilience is implemented by the superiorization methodology, which uses perturbations pro-actively during the performance of the iteration scheme to steer the algorithm to an output that still solves the given problem, but is superior with respect to a given secondary criterion. See e.g. censor18 ; censor15 ; censor17 ; cdh10 ; cdhst14 ; cz15 . This approach has been used in the same application as the one we consider in this paper: Intensity-modulated radiation therapy (IMRT) treatment planning bgks17 ; gkrs18 . Gibali et al. in gkrs18 presented a detailed comparison of superiorized and unsuperiorized projection methods used in the level set scheme to solve optimization problems in the field of IMRT. In bgks17 these optimization problems are reformulated as lexicographic optimization problems. The perturbations used in a given level of the lexicographic optimization scheme are chosen to be descent directions with respect to the objective function of the subsequent optimization level. This approach exhibits faster convergence than the unperturbed algorithm.

Our motivation for using perturbations, however, is a different one. When we solve IMRT optimization problems as described previously, we observe zigzagging behavior of the projection methods. This behavior is not surprising as we attempt to satisfy two main goals, which conflict with each other: Irradiating the tumor and at the same time sparing surrounding healthy organs. The zigzagging behavior is similar to the one known to occur in projected gradient methods, which arises from the conflict between seeking feasibility and reduction of the objective function, see wang2008 and the references therein.

In this paper we introduce the heavy ball and the surrogate constraint perturbation, which address the zigzagging behavior. We offer formulations of them both as inner and outer perturbations. We show that all of these formulations describe bounded perturbations. Finally, we use the heavy ball and surrogate constraint perturbation together with both the cyclic and simultaneous subgradient projection method.

The heavy ball perturbation uses the approach of the heavy ball method, which belongs to the group of $k$-step or inertial methods. These algorithms use gradient information from $k$ previous iteration steps to calculate their next iterate. Incorporating previous gradient information alleviates zigzagging behavior compared to methods, which only use current gradient information. In recent years authors demonstrated that many optimization algorithms of this kind can be seen as a discretization of the trajectories of differential equations derived from the field of continuous dynamical systems, e.g., acpr:18 ; frv:18 ; Nesterov:83 ; Polyak:64 . It is shown that inertial approaches are very effective and demonstrate good convergence properties.

The second perturbation we discuss is the surrogate constraint perturbation. It is inspired by the algorithm presented in Dudek2007 . This algorithm combines the approach of classical surrogate constraint methods (see e.g. Yang1992 ) with relaxed, averaged projections. The author considers linear feasibility problems and demonstrates that especially for flat solution sets (which are associated with zigzagging behavior by methods using gradient information), their algorithm converges faster than a classical surrogate constraint method. In our work we illustrate that their idea can be successfully translated to the nonlinear context of IMRT optimization problems.

The paper is organized as follows. In Section 2 we present definitions of methods and concepts, which we will use in our algorithm. In Section 3 we present some results on the convergence of algorithmic operators using inner and outer perturbations. Then we introduce the perturbations developed by us and show that they are bounded. We end this section by summarizing our algorithm, which incorporates all of the concepts described before. Finally, in Section 4 we present numerical results from both linear feasibility problems and nonlinear optimization problems arising from IMRT treatment planning.

2 Preliminaries

We now present concepts and methods as well as mathematical background which is essential for the introduction of our results.

In this paper we consider the general constrained convex optimization problem (1.1). We reformulate it into an equivalent epigraph form

[TABLE]

We denote the optimal value of (2.1) by $t^{\ast}$. We assume that $t^{\ast}$ is finite and that there exists a point $x^{\ast}$ with $f(x^{\ast})=t^{\ast}$, which is feasible for (2.1).

Our approach to solving (2.1) is to transform it into a sequence of convex feasibility problems (CFP), see also (Bertsekas99, , Subsection 2.1.2). The general formulation of a CFP is the following:

[TABLE]

where the functions $\varphi_{i}$ are convex for all $i\in I$.

In the following we describe the procedure we use to solve (2.1), which is called the level set scheme gkrs18 . We use a decreasing sequence $\left\{t_{s}\right\}_{s=1}^{\infty}$. This sequence is not known a priori. The first element $t_{1}=\infty$ is fixed and every subsequent element of the sequence is determined later during the iteration process. We define the functions

[TABLE]

Then we try to solve the CFP $\boldsymbol{P}^{s}$:

[TABLE]

with $I=\{1,...,m+1\}$. If $C\neq\emptyset$, a solution of $\boldsymbol{P}^{s}$ found by a suitable algorithmic operator $\boldsymbol{T}$ is denoted by $x^{\ast}_{s}$. Now $t_{s+1}$ is calculated according to a user defined update rule like for example $t_{s+1}=f(x^{\ast}_{s})-\varepsilon_{s}$ or $t_{s+1}=f(x^{\ast}_{s})(1-\varepsilon_{s})$ where $\{\varepsilon_{s}\}$ is some user chosen sequence with $\varepsilon_{s}>0$ for all $s\geq 0$. Using $t_{s+1}$ the next CFP $\boldsymbol{P}^{s+1}$ can be formulated and the iteration proceeds.

If $C=\emptyset$ or the algorithmic operator $\boldsymbol{T}$ is unable to find a solution of $\boldsymbol{P}^{s}$ for other reasons, we consider $x^{\ast}_{s-1}$ to be the result of the level set scheme.

One successful class of algorithmic operators for solving CFPs is the class of projection methods. These are iterative algorithms that use projections onto sets, relying on the principle that when a family of sets is present, then projections onto the given individual sets are easier to perform than projections onto other sets (intersections, image sets under some transformation, etc.) that are derived from the given individual sets. Projection methods can have various algorithmic structures some of which are particularly suitable for parallel computing, and possess desirable convergence properties. See e.g. cz97 ; combettes97 ; znly18 . A special case of a CFP is the linear feasibility problem $Ax=b$. An illustration of different projection methods for this problem can be found in ceh2000 .

In this paper we present two different projection methods as options for the algorithmic operator $\boldsymbol{T}$. The first one is the cyclic subgradient projection method, which is also used in combination with the level set scheme in gkrs18 . Like some other projection methods, it utilizes the concept of control sequences. These are sequences $\{i(\nu)\}_{\nu=0}^{\infty}$, which determine an ordering of the indices of the sets $C_{i}$ mentioned in (2.2). The kind of control sequence used by the cyclic subgradient projection method is defined below.

Definition 1

The sequence $\{i(\nu)\}_{\nu=0}^{\infty}$ is called a cyclic control sequence, if $i(\nu)=\left(\nu\operatorname{mod}n\right)+1$, where $n$ is the number of sets in (2.2).

The cyclic subgradient projection method can be written as algorithmic operator $\boldsymbol{T}$ in the following way:

Let $x^{0}\in\mathbb{R}^{n}$ be an arbitrary starting point. Given the current iterate $x^{k}$, the next iterate $x^{k+1}$ can be calculated via

[TABLE]

with

[TABLE]

where $\xi^{k}\in\partial\varphi_{i(k)}(x^{k})$ (subgradient of $\varphi_{i(k)}$ at $x^{k}$ ) is arbitrary, $\lambda_{k}\in[\epsilon_{1},2-\epsilon_{2}]$ (relaxation parameters) for arbitrary $\epsilon_{1},\epsilon_{2}\in(0,1]$ and $\{i(k)\}$ is a cyclic control sequence.

The second option we present for the algorithmic operator $\boldsymbol{T}$ is the simultaneous subgradient projection method. We can write it in terms of $\boldsymbol{T}$ as follows:

Let $x^{0}\in\mathbb{R}^{n}$ be an arbitrary starting point. Given the current iterate $x^{k}$, the next iterate $x^{k+1}$ can be calculated via

[TABLE]

with

[TABLE]

where $\xi^{k}\in\partial\varphi_{i}(x^{k})$ (subgradient of $\varphi_{i}$ at $x^{k}$) is arbitrary, $w_{i}>0$ are weights with $\sum_{i\in I}w_{i}=1$ and $\lambda_{k}\in[\epsilon_{1},2-\epsilon_{2}]$ (relaxation parameters) for arbitrary $\epsilon_{1},\epsilon_{2}\in(0,1]$.

In the following, we use the notation $\bar{p}(x^{k})\colonequals p(x^{k})/\|p(x^{k})\|$, refer to the simultaneous subgradient projection method simply as simultaneous projection and to the cyclic subgradient projection method as cyclic projection.

The convergence proof of the level set scheme (gkrs18, , Theorem 3.6) relies on using a finite convergence method. We can transform both the cyclic and the simultaneous projection to fulfill this property using the approach of pi1988 . Note that (gkrs18, , Theorem 3.6), in contrast to convergence proofs for other finite convergent projection methods, e.g. ccp11 , is not based on the assumption that the Slater Condition holds. In our context this means that we do not rely on the existence of $x\in\mathbb{R}^{n}$ with $\varphi_{i}(x)<0$ for all $i\in I$ for the level set scheme to converge.

Perturbations can be applied to any algorithmic operator $\boldsymbol{T}$ which generates a sequence $\{x^{k}\}_{k=0}^{\infty}$ via $x^{k+1}=\boldsymbol{T}(x^{k})$. A central concept in the context of perturbations is bounded perturbation resilience. This concept is defined in bdhk07 in a general way, i.e. for a general mathematical problem $\boldsymbol{P}$ and an algorithmic operator $\boldsymbol{T}$, which is suited to solve $\boldsymbol{P}$.

In dgj18 ; dgjt17 the concept of bounded perturbation resilience is considered with the problem $\boldsymbol{P}$ being the variational inequality problem. The operators $\boldsymbol{T}$ used to solve $\boldsymbol{P}$ in dgj18 ; dgjt17 are the extragradient method, the subgradient extragradient method and the projection and contraction method. In gh14 ; jcj13 the problem is the bioluminescence imaging problem, which can be phrased as a constrained optimization problem, and the operator used to solve it is the expectation maximization method. For a detailed overview and more examples see censor18 .

In this paper, the problem $\boldsymbol{P}$ is a CFP and $\boldsymbol{T}$ can be read as one of the projection methods we presented before.

Definition 2

Given a problem $\boldsymbol{P}$, an algorithmic operator $\boldsymbol{T}$ and a starting point $x^{0}$ such that the sequence $\{x^{k}\}_{k=0}^{\infty},$ generated by $x^{k+1}=\boldsymbol{T}(x^{k})$ converges to a solution of $\boldsymbol{P}$. Then $\boldsymbol{T}$ is called bounded perturbation resilient if any sequence $\{y^{k}\}_{k=0}^{\infty}$ with $y^{0}=x^{0}$ generated using either inner perturbations via

[TABLE]

or using outer perturbations via

[TABLE]

where $\beta_{k}v^{k}$ are bounded perturbations (i.e. $\beta_{k}\in\mathbb{R}_{\geq 0}$ for all $k\geq 0$, ${\displaystyle\sum\limits_{k=0}^{\infty}}\beta_{k}\,<\infty$, $v^{k}\in\mathbb{R}^{n}$ and $\|v^{k}\|\leq M\in\mathbb{R}$ for all $k\geq 0$) also converges to a solution of $\boldsymbol{P}$.

Both the cyclic and the simultaneous projection method are known to be bounded perturbation resilient special-issue .

3 Perturbation of the simultaneous and cyclic subgradient projection method

In this section we first present results on the convergence of sequences, which were generated using inner or outer perturbations. In Subsections 3.2 and 3.3 we introduce two specific perturbations developed by us, show how they can be formulated both as inner and outer perturbation and that they are bounded in the sense of Definition 2. Finally, we summarize how we combined the tools presented in this section and Section 2 in our algorithm.

3.1 The relation of inner and outer perturbations

In this section we explore the relation of iteration sequences resulting from methods using inner and outer perturbations concerning their convergence behavior.

First, we are extending Proposition 5 from CensorReem2014 to include the simultaneous subgradient projection method (2.7) as algorithmic operator $\boldsymbol{T}$.

Let

[TABLE]

be perturbation vectors where $c:\mathbb{R}^{n}\rightarrow\{$false, true$\}$ is a function that gives a condition, which determines whether or not to perturb the current iterate.

Consider the sequences $\{y^{k}\}_{k=0}^{\infty}$ using outer perturbations and $\{z^{k}\}_{k=0}^{\infty}$ using inner perturbations, which are defined as follows:

[TABLE]

and

[TABLE]

Proposition 1

Suppose that $\{b^{k}\}_{k=0}^{\infty}$ is a sequence in $\mathbb{R}^{n}$ satisfying $\lim_{k\rightarrow\infty}b^{k}\rightarrow 0$. If $\{y^{k}\}_{k=0}^{\infty}$ as defined in (3.2) converges weakly to some $y^{\ast}$, then also $\{z^{k}\}_{k=0}^{\infty}$ as defined in (3.3) converges weakly to $y^{\ast}$ and vice versa. If $\{y^{k}\}_{k=0}^{\infty}$ converges strongly, then $\{z^{k}\}_{k=0}^{\infty}$ converges strongly to the same limit and vice versa.

Proof

We show that by induction $y^{k+1}=z^{k}+\tilde{b}^{k}$ for all $k\in\mathbb{N}\cup\{0\}$. The rest of the statement follows from the proof of Proposition 5 in CensorReem2014 .

k=0

Show that $y^{1}=z^{0}+\tilde{b}^{0}$.

[TABLE]

Induction step

Show that $y^{k}=z^{k-1}+\tilde{b}^{k-1}\Rightarrow y^{k+1}=z^{k}+\tilde{b}^{k}$.

[TABLE]

Proposition 1 implies that inner and outer perturbations of the simultaneous projection method using the same perturbation vectors $\tilde{b}^{k}$ as described before can be used interchangeably with regard to weak and strong convergence.

In Sections 3.2 and 3.3 we present two perturbations developed by us and formulate them both as inner and outer perturbation. In contrast to (3.2) and (3.3), we do not use the same perturbation vectors $\tilde{b}^{k}$ for both formulations there. Instead, we formulate them according to the following schemes. We assume specific choices of $\lambda_{k},\beta_{k}$, which will be explained in more detail, and present a result regarding the convergence of the sequence of iterates produced by methods using these perturbations.

Let again $\boldsymbol{T}$ be an algorithmic operator with

[TABLE]

where $p(x)$ is either the simultaneous (2.8) or the cyclic projection step (2.6). Note that $\lambda:\mathbb{R}^{n}\rightarrow[\epsilon_{1},2-\epsilon_{2})$ with arbitrary $\epsilon_{1},\epsilon_{2}\in(0,1]$ is a function here.

Let $c:\mathbb{R}^{n}\rightarrow\{$false, true$\}$ and $\tilde{c}:\mathbb{R}^{n}\times\mathbb{R}^{n}\rightarrow\{$false, true$\}$ with

[TABLE]

be functions, which determine, whether perturbation is applied in a certain iteration. $\tilde{c}(y^{k},y^{k-1})$ is by construction only true, if $c(y^{k-1})$ is false and $c(y^{k})$ is true.

Let $\{y^{k}\}_{k=0}^{\infty}$ be the sequence of iterates produced by the outer perturbation scheme, which is defined as follows.

[TABLE]

where $\{\beta_{k}\}_{k=0}^{\infty}$ is a sequence with $\beta_{k}\in\mathbb{R}_{\geq 0}$ for all $k\geq 0$ and ${\displaystyle\sum\limits_{k=0}^{\infty}}\beta_{k}\,<\infty$,

[TABLE]

and the sequence $\{b^{k}\}_{k=0}^{\infty}$ is bounded. Furthermore assume that $\|p(y^{k})\|<q\in\mathbb{R}$ for all $k\geq 0$.

Let $\{z^{k}\}_{k=0}^{\infty}$ be the sequence of iterates produced by the inner perturbation scheme, which is defined as follows.

[TABLE]

where $\{\beta_{k}\}_{k=0}^{\infty}$ is a sequence with $\beta_{k}\in\mathbb{R}_{\geq 0}$ for all $k\geq 0$ and ${\displaystyle\sum\limits_{k=0}^{\infty}}\beta_{k}\,<\infty$,

[TABLE]

and the sequence $\{b^{k}\}_{k=0}^{\infty}$ is bounded.

Now we choose

[TABLE]

for some $K\in\mathbb{N}$. Choosing $\beta_{k}$ like this means that for $k>K$ no more perturbations will be applied. In our computations we choose $K$ as the maximum number of iterations allowed for the projection method to solve the current CFP.

Lemma 1

If $y^{0}=z^{0}$ and $\{\beta_{k}\}_{k=0}^{\infty}$ is chosen as in (3.10), the following statements are true.

• (a)

If for some $l<K$ it holds that $c(z^{k})=$ false for all $k=0,...,l-1$ and $c(z^{l})=$ true, we have $z^{k}=y^{k}$ for all $k=0,...,l$.

• (b)

$\{z^{k}\}_{k=0}^{\infty}\subseteq\{y^{k}\}_{k=0}^{\infty}$.

• (c)

For all $z^{k}\in\{z^{k}\}_{k=0}^{\infty}$ exists $N\in\mathbb{N}_{0},0\leq N\leq K$ such that $z^{k}=y^{k+N}$.

Proof

Statement (a) follows directly from equations (3.9) and (3.7).

Now suppose, that $l$ with $0\leq l<K$ is the first iteration index where the condition for applying perturbations is fulfilled. This means that $c(z^{l})=\tilde{c}(y^{l},y^{l-1})=$ true and $c(z^{k})=\tilde{c}(y^{k},y^{k-1})=$ false for all $k<l$. We know from (a) that then $z^{k}=y^{k}$ for all $k=0,...,l$.

We have $\beta_{l}=1$ because $l<K$. The iteration schemes give us

[TABLE]

and

[TABLE]

$c(y^{l})=$ true, so $\tilde{c}(y^{l+1},y^{l})=$ false, no matter what $c(y^{l+1})$ is. The following iterate is therefore an unperturbed one.

[TABLE]

Because $\lambda$ depends on the iterate (and not on the iteration index) and $y^{l+2}=z^{l+1}$, we get

[TABLE]

assuming that $c(z^{l+1})=$ false. Otherwise, the same argument as before holds. Statement (b) follows from this. Before any perturbations are applied, the number $N$ from statement (c) is 0. Each time perturbations are applied in an iteration, $N$ increases by 1.

Proposition 2

Suppose that $\{\beta_{k}\}_{k=0}^{\infty}$ is chosen as in (3.10) and $\{b^{k}\}_{k=0}^{\infty}$ is a bounded sequence in $\mathbb{R}^{n}$. If $\{y^{k}\}_{k=0}^{\infty}$ converges weakly to some $y^{\ast}$, then also $\{z^{k}\}_{k=0}^{\infty}$ converges weakly to $y^{\ast}$ and vice versa. If $\{y^{k}\}_{k=0}^{\infty}$ converges strongly, then $\{z^{k}\}_{k=0}^{\infty}$ converges strongly to the same limit and vice versa.

Proof

It follows from Lemma 1(c) that there exists $0\leq N\leq K$ such that $z^{k}=y^{k+N}$ for all $k=l,...\infty$ for some $l\in\mathbb{N}_{0}$.

Proposition 2 implies that inner and outer perturbations of both the simultaneous and cyclic projection method using the perturbation vectors (3.9) or (3.7) can be used interchangeably with regard to weak and strong convergence.

3.2 The heavy ball perturbation

An early example of inertial-type methods (also known as $k$-step methods) is the heavy ball method of Polyak:64 . Inertial-type methods are a time-discretization of an ordinary differential equation defining a continuous-time dynamical system (in general are easier to understand than their discrete-time counterparts, see Borwein:2018aa ). These methods incorporate gradient information from the last $k$ iterates into the calculation of the iteration step towards the next iterate and can be used to avoid zigzagging behavior. It is shown that methods using such inertial terms progress converge faster than methods using only current gradient information. For a deeper discussion of this matter, see e.g. acpr:18 ; Nesterov:83 as well as dgj18 ; dgjt17 and the many references therein.

Our approach is to use such terms as perturbations for algorithmic operators $\boldsymbol{T}$ as defined in (3.4) and show that such interference can accelerate the convergence of the algorithm.

3.2.1 Formulation as inner perturbation

The iteration scheme using the heavy ball perturbation as inner perturbation is the following:

[TABLE]

where

[TABLE]

$\{\beta_{k}\}$ is a sequence with $\beta_{k}\in\mathbb{R}_{\geq 0}$ for all $k\geq 0$ and $\sum_{k=1}^{\infty}\beta_{k}<\infty$, $\{\lambda^{HB}_{k}\}$ is a bounded user-chosen sequence of step lengths and $d^{HB}=\bar{p}(x^{k-1})+\bar{p}(x^{k})$, where $p(x)$ is either (2.6) or (2.8). The perturbation direction $d^{HB}$ is chosen like this, similar to iteration rules of $k$-step methods, in order to make opposing parts of $\bar{p}(x^{k-1})$ and $\bar{p}(x^{k})$ cancel each other, because those parts are provoking the zigzagging behavior. Instead we get a direction which contains what both vectors have in common. The function $c$, which determines, whether perturbations are applied or not, is defined as follows.

[TABLE]

Choosing $c$ like this means, that perturbations are applied, when $\langle\bar{p}(x^{k-1}),\bar{p}(x^{k})\rangle$ is close to $-1$, which translates to $\bar{p}(x^{k-1})$ and $\bar{p}(x^{k})$ pointing into almost opposite directions and therefore zigzagging behavior. If that is the case, the convergence of simultaneous and cyclic projection is slow and our aim is to accelerate it by using perturbations.

Recalling Definition 2, we now show that the perturbations (3.12) are bounded.

Lemma 2

Let $c(x^{k})=$ true. Choose a bounded sequence $\{\lambda^{HB}_{k}\}$ with $\lambda^{HB}_{k}\geq 0$ for all $k\geq 0$. Then, the perturbations $\beta_{k}v^{k}$ defined in (3.12) are bounded.

Proof

With $\{\beta_{k}\},\{\lambda^{HB}_{k}\}$ chosen as described in the assumptions, it suffices to show that $\|d^{HB}\|$ is bounded.

[TABLE]

3.2.2 Formulation as outer perturbation

The iteration scheme using the heavy ball perturbation as outer perturbation is the following:

[TABLE]

where

[TABLE]

$\{\beta_{k}\}$ is a sequence with $\beta_{k}\in\mathbb{R}_{\geq 0}$ for all $k\geq 0$ and $\sum_{k=1}^{\infty}\beta_{k}<\infty$, $\{\lambda^{HB}_{k}\}$ is a user-chosen bounded sequence of step lengths and $d^{HB}=\bar{p}(x^{k-1})+\bar{p}(x^{k})$, where $p(x)$ is either (2.6) or (2.8). The function $\tilde{c}$ is defined as in (3.5) and (3.13).

Recalling Definition 2, we now show that the perturbations (3.15) are bounded.

Lemma 3

Assume, that $\|p(x^{k})\|\leq q\in\mathbb{R}$ for all $k\geq 0$ and $\{\lambda^{HB}_{k}\}$ is a bounded sequence which fulfills $\lambda^{HB}_{k}\geq 0$ for all $k\geq 0$.

Then, the perturbations $\beta_{k}v^{k}$ defined in (3.15) are bounded.

Proof

With $\{\beta_{k}\}$ chosen as described in the assumptions, it suffices to show that $\|v^{k}\|<\bar{M}\in\mathbb{R}$ for all $k\geq 0$. We have shown before, that $\|d^{HB}\|<M\in\mathbb{R}$.

[TABLE]

3.3 The surrogate constraint perturbation

The approach we used for the surrogate constraint perturbation originates from a special kind of surrogate constraint method presented in Dudek2007 . In this method, the author combines the traditional idea of surrogate constraint methods (see e.g. Yang1992 ) with relaxed, averaged projections. This algorithm has the ability to prevent zigzagging behavior in solving a linear feasibility problem.

Again, we consider an algorithmic operator $\boldsymbol{T}$ as defined in (3.4).

We modified the idea of Dudek2007 in such a manner that we did not project the individual gradients of the violated constraint functions onto the surrogate half-space $H^{SC}$. Instead, we calculated $p(x^{k})$ as defined in (2.7) and then projected it onto the half-space $H^{SC}$.

In Figure 1 the vectors $\bar{p}(x^{k-1})$ and $p(x^{k})$ are depicted in red. We define the half-spaces

[TABLE]

The vector

[TABLE]

is the vector pointing from $x^{k}$ to the metric projection of the point $x^{k}+p(x^{k})$ onto the surrogate half-space $H^{SC}$. It is depicted in cyan in Figure 1. Projecting onto $H^{SC}$ eliminates the part of $p(x^{k})$ which points into the direction of $-p(x^{k-1})$. This means that $d^{SC}\perp p(x^{k-1})$ and steering the iteration into the direction of $d^{SC}$ reduces zigzagging behavior.

By $A$ we denote the intersection point of $H^{k-1}$, $H^{k}$ and $span\{p(x^{k-1}),p(x^{k}))\}$. Let $\alpha:=\cos^{-1}(1-\epsilon)$ be the angle between $-p(x^{k-1})$ and $p(x^{k})$ and $\beta\colonequals\pi/2-\alpha$. Let $B\colonequals P_{H^{SC}}(x^{k}+p(x^{k}))$ and $C\colonequals x^{k}+p(x^{k})$ and . The triangles $\Delta CBx^{k}$ and $\Delta CAx^{k}$ are similar, so it is true that

[TABLE]

Choosing the step length

[TABLE]

we get

[TABLE]

These equations yield $\|A-x^{k}\|=\|\lambda^{SC}_{k}d^{SC}\|$. The construction of $d^{SC}$ implies that $d^{SC}/\|d^{SC}\|=(A-x^{k})/\|A-x^{k}\|$ and therefore $A-x^{k}=\lambda^{SC}_{k}d^{SC}$. This means that the point $x^{k}+\lambda^{SC}_{k}d^{SC}$ is actually the same as the point $A$ depicted in Figure 1. Note that

[TABLE]

As before using the heavy ball perturbation, we perturb the iteration sequence generated by $\boldsymbol{T}$ if $\langle\bar{p}(x^{k}),\bar{p}(x^{k-1})\rangle=-1+\epsilon$ for $\epsilon\in[\epsilon_{\min},\epsilon_{\max}]\subset(0,1)$.

3.3.1 Formulation as inner perturbation

The iteration scheme using the surrogate constraint perturbation as inner perturbation is

[TABLE]

where $\{\beta_{k}\}$ is a sequence with $\beta_{k}\in\mathbb{R}_{\geq 0}$ for all $k\geq 0$ and ${\displaystyle\sum\limits_{k=0}^{\infty}}\beta_{k}\,<\infty$ and

[TABLE]

Again, $c$ is defined as in (3.13).

As before in Section 3.2, we now show that the perturbations (3.23) are bounded in the sense of the definition 2 of bounded perturbation resilience.

Lemma 4

Let $\langle\bar{p}(x^{k}),\bar{p}(x^{k-1})\rangle=-1+\epsilon$ for $\epsilon\in[\epsilon_{\min},\epsilon_{\max}]\subset(0,1)$. Further assume that $\|p(x^{k})\|\leq q\in\mathbb{R}$ holds for all $k\geq 0$.

Then, the perturbations $\beta_{k}v^{k}$ as defined in (3.23) are bounded.

Proof

With $\{\beta_{k}\}$ chosen as described in the assumptions, it suffices to show that $\|v^{k}\|\leq M\in\mathbb{R}$ for all $k\geq 0$. Recall that $\alpha=\cos^{-1}(1-\epsilon)$ as in Figure 1. We know from (3.21)

[TABLE]

3.3.2 Formulation as outer perturbation

The iteration scheme using the surrogate constraint perturbation as outer perturbation is

[TABLE]

where $\{\beta_{k}\}$ is a sequence with $\beta_{k}\in\mathbb{R}_{\geq 0}$ for all $k\geq 0$ and ${\displaystyle\sum\limits_{k=0}^{\infty}}\beta_{k}\,<\infty$ and

[TABLE]

The function $\tilde{c}$ is the same as in (3.5) with $c$ as in (3.13).

As before in Section 3.2, we now show that the perturbations (3.25) are bounded in the sense of the Definition 2.

Lemma 5

Let $\langle\bar{p}(x^{k-1}),\bar{p}(x^{k})\rangle=-1+\epsilon$ for $\epsilon\in[\epsilon_{\min},\epsilon_{\max}]\subset(0,1)$. Further assume that $\|p(x^{k})\|\leq q\in\mathbb{R}$ holds for all $k\geq 0$.

*Then, the perturbations $\beta_{k}v^{k}$ as defined in (3.25) are bounded.

Proof

With $\{\beta_{k}\}$ chosen as described in the assumptions, it suffices to show that $\|v^{k}\|\leq\bar{M}\in\mathbb{R}$ for all $k\geq 0$. We have shown before that $\|\lambda^{SC}d^{SC}\|\leq M\in\mathbb{R}$.

Let $\beta=\cos^{-1}(\langle\bar{p}(x^{k-1}),\bar{p}(x^{k})\rangle)$. We have $d^{SC}\in span\{p(x^{k-1}),p(x^{k})\},d^{SC}\perp p(x^{k-1})$ and $1/\|d^{SC}\|\langle d^{SC},\bar{p}(x^{k})\rangle>0$. Therefore, $\angle(d^{SC},p(x^{k}))=\beta-\pi/2$, which lies in the open interval $(0,\pi/2)$.

We also know that

[TABLE]

[TABLE]

3.4 Convergence speed of the perturbed projection methods

We have shown that the property of boundedness in the sense of Definition 2 is fulfilled by the heavy ball and the surrogate constraint perturbation, no matter whether they are formulated as inner or outer perturbation. Because both simultaneous and cyclic projection are bounded perturbation resilient, we can use the heavy ball and the surrogate constraint perturbation together with these methods and retain convergence to a solution of $\boldsymbol{P}^{s}$, if such a solution exists. (brz17, , Theorem 11 and Corollary 14) states that we retain the convergence rate of the simultaneous or cyclic projection when we use them with bounded perturbations. In the Section 4 we show numerical results, which indicate that in some cases using the heavy ball or the surrogate constraint perturbation can even speed up the convergence. In this section we will elaborate on theoretical arguments, which illustrate in which cases the perturbed simultaneous projection method converges faster than its unperturbed counterpart.

Statements about the convergence speed of both the simultaneous projection method and the method introduced by Dudek Dudek2007 exist in the literature. Cegielski (cegielski13, , Theorem 4.4.5) states the following result about the convergence speed of the simultaneous projection method. It is phrased in terms of metric projections $P_{C_{i}}$ onto closed convex sets $C_{i}$.

Theorem 3.1

If $z\in C\colonequals\bigcap_{i\in I}C_{i},C\neq\emptyset,\lambda\in(0,2]$ and

[TABLE]

it holds that

[TABLE]

The essential property of the operators $P_{C_{i}}$ needed for the proof of this result is that they are cutters, i.e.

[TABLE]

(cegielski13, , Corrolary 4.2.6) states that not only metric projections but also subgradient projections are cutters. Therefore (3.27) also holds for the subgradient projections we introduced in (2.8). To distinguish between metric and subgradient projections, we use the following notation:

[TABLE]

where $\xi$ is an arbitrary element of the subdifferential $\partial\varphi_{i}(x)$ of $\varphi_{i}$ at $x$.

Dudek in Dudek2007 states a result about the convergence of their method with a very similar structure to (3.27). Their method is phrased as an algorithm to solve a system of linear inequalities and the projections considered are metric projections. Nevertheless, the result also applies to subgradient projections, because the proof relies on the fact that the projections are cutters.

Dudeks result rephrased in terms of subgradient projection operators $\tilde{P}_{C_{i}}$ is the following:

Theorem 3.2

If $z\in C\colonequals\bigcap_{i\in I}C_{i},,C\neq\emptyset,\lambda\in[0,2),H^{SC}$ defined as in (3.16) and

[TABLE]

with

[TABLE]

it holds that

[TABLE]

Comparing (3.27) and (3.34) we observe that the bound on the reduction $\|x^{k+1}-z\|^{2}-\|x^{k}-z\|^{2}$ provided by the simultaneous projection and the bound on the reduction $\|x^{SC}-z\|^{2}-\|x^{k}-z\|^{2}$ provided by Dudek’s method differs only by the factor

[TABLE]

Because $P_{H^{SC}}$ is a metric projection onto a half-space we know that for all $k\geq 0$

[TABLE]

It follows that

[TABLE]

and thus $\delta(x^{k})\geq 1$ for all $k\geq 0$. The immediate consequence is that Dudek’s method reduces the distance to any $z\in C$ at least as much as the simultaneous projection.

3.4.1 Acceleration of the convergence by perturbations

We consider the iteration step $k\rightarrow k+1$, which is perturbed using outer perturbations, i.e. we assume that $\tilde{c}(x^{k},x^{k-1})=$ true with the function $\tilde{c}$ as defined in (3.5).

Surrogate constraint perturbation

Lemma 6

If $\langle p(x^{k-1}),\tilde{P}_{C_{i}}(x^{k})-x^{k}\rangle\leq 0$ for all $i\in I$ then

[TABLE]

Proof

$-\bar{p}(x^{k-1})$ is the normal vector $\nu$ to the surrogate half-space $H^{SC}$ (with $\nu$ pointing away from $H^{SC}$).

Let $v=\nu\langle\nu,v\rangle+u$ and $\langle u,\nu\rangle=0$. From the definition of the metric projection it follows that

[TABLE]

Now let

[TABLE]

The comparison of

[TABLE]

and

[TABLE]

proves Lemma 6.

Lemma 6 implies that if $\langle p(x^{k-1}),\tilde{P}_{C_{i}}(x^{k})-x^{k}\rangle\leq 0$ is fulfilled for all $i\in I$, then the surrogate constraint perturbation (with $\beta_{k}=1$) and Dudek’s method (with $\lambda=1$) actually result in the same iterate $x^{SC}\colonequals x^{k}+\lambda^{SC}_{k}d^{SC}$. This means that the statement of Theorem 3.2 also applies for the next iterate resulting from the surrogate constraint perturbation. Note that the assumption $\langle p(x^{k-1}),\tilde{P}_{C_{i}}(x^{k})-x^{k}\rangle\leq 0$ for all $i\in I$ is a reasonable assumption if there are two main groups of conflicting goals, each iteration step fulfills one group of goals, and the groups alternate in each iteration step.

Furthermore we can quantify the factor $\delta(x^{k})$ in dependence of $\langle\bar{p}(x^{k}),\bar{p}(x^{k-1})\rangle$ under this assumption, as we will demonstrate in the following theorem. Let $\alpha\colonequals\pi-\cos^{-1}(\langle\bar{p}(x^{k}),\bar{p}(x^{k-1})\rangle)$.

Theorem 3.3

If $\langle p(x^{k-1}),\tilde{P}_{C_{i}}(x^{k})-x^{k}\rangle\leq 0$ for all $i\in I$ and

[TABLE]

with $\lambda^{SC}_{k}$ and $d^{SC}$ as defined in (3.18) and (3.19) then

[TABLE]

Proof

From Theorem 3.2 we know

[TABLE]

The inequality

[TABLE]

yields

[TABLE]

The result now follows from (3.20).

Heavy ball perturbation

By $x^{HB}\colonequals x^{k}+0.5\lambda^{HB}(\bar{p}(x^{k})+\bar{p}(x^{k-1}))$ we denote the next iterate generated by the perturbed iteration scheme using heavy ball perturbation. Let $\langle\bar{p}(x^{k}),\bar{p}(x^{k-1})\rangle=-1+\epsilon$.

The triangle in Figure 2 is the same as in Figure 1. The Figure illustrates the progress towards the triangle tip $A$ made by the perturbed iteration scheme using heavy ball perturbation. The point $x^{HB}$ lies on the dashed cyan line, its exact position depending on the step length $\lambda^{HB}$. Let $\gamma\colonequals\angle(A-x^{k},x^{HB}-x^{k})$. The triangle $x^{k},x^{k}+0.5\bar{p}(x^{k-1}),x^{HB}$ is isosceles, because $\|0.5\bar{p}(x^{k-1})\|=\|0.5\bar{p}(x^{k})\|=0.5$. Therefore $\gamma=\pi/2-(\pi-\alpha)/2=\alpha/2$.

We have

[TABLE]

The distance $\|x^{HB}-A\|$ of course depends on the step size $\lambda^{HB}$. It attains its minimum if $\langle x^{HB}-A,x^{HB}-x^{k}\rangle=0$, which occurs if

[TABLE]

which means that

[TABLE]

We get

[TABLE]

In Section 3.3 we have illustrated that $x^{SC}$ coincides with the triangle tip $A$. Recall that $\alpha$ is tied to the parameters $\epsilon_{\min},\epsilon_{\max}$ in (3.13), which are used to decide whether the perturbation is applied in a certain iteration $k$.

Theorem 3.4

If $\langle p(x^{k-1}),\tilde{P}_{C_{i}}(x^{k})-x^{k}\rangle\leq 0$ for all $i\in I$ then there exists $\tilde{\alpha}\in(0,\pi/3]$ such that if $\alpha=\pi-\cos^{-1}(\langle\bar{p}(x^{k}),\bar{p}(x^{k-1})\rangle)<\tilde{\alpha}$ is fulfilled and $\lambda^{HB}$ is chosen as in (3.42) the following statement holds:

[TABLE]

Proof

From (3.43) it follows that for $x^{SC}=x^{k}+\lambda^{SC}_{k}d^{SC}$ and $x^{HB}_{opt}=x^{k}+0.5\lambda^{HB}_{opt}(\bar{p}(x^{k})+\bar{p}(x^{k-1}))$ it is true that

[TABLE]

Theorem 3.3 yields

[TABLE]

To show that the heavy ball perturbation reduces the distance to $z$ even more than the unperturbed simultaneous projection, we examine

[TABLE]

[TABLE]

For $\alpha<\pi/3$ we have

[TABLE]

For all (finite) $\|p(x^{k})\|,\|x^{SC}-z\|$ there exists $\tilde{\alpha}\in(0,\pi/3]$ such that for all $\alpha\in[0,\tilde{\alpha})$:

[TABLE]

and thus $\|x^{HB}_{opt}-z\|^{2}<\|x^{k}-z\|^{2}-\sum w_{i}(x^{k})\|\tilde{P}_{C_{i}}(x^{k})-x^{k}\|^{2}$.

Theorem 3.4 implies that under the mentioned assumptions there always exists a value for $\epsilon_{\max}$ with $\epsilon_{\max}>0$ as parameter for the function $\tilde{c}$, which guarantees that if $\langle\bar{p}(x^{k}),\bar{p}(x^{k-1})\rangle<-1+\epsilon_{\max}$ the reduction achieved by the heavy ball perturbation is larger than the reduction achieved by the unperturbed simultaneous projection.

3.5 The algorithm

Now we combine the tools we presented in the previous sections in our algorithm.

To solve a constrained convex optimization problem (1.1) we translate it into its epigraph representation (2.1). Then we use the level set scheme to transform it into a sequence $\{\boldsymbol{P}^{s}\}$ of CFPs . Next, we use simultaneous (2.7) or cyclic projection (2.5) to solve each CFP $\boldsymbol{P}^{s}$. If the use of perturbations is specified, we alter the sequence of iterates generated by the projection methods using the heavy ball (3.15) or the surrogate constraint perturbation (3.25).

4 Numerical demonstrations

In this section we present our numerical results. These were achieved with an algorithm implemented according to the description in Section 3.5. It uses the heavy ball and the surrogate constraint perturbation with the outer perturbation scheme (3.6). As the results in this section will show, both perturbations are eminently useful to speed up the convergence of the algorithm towards a solution compared to the unperturbed methods.

4.1 Linear feasibility problem

We present a linear feasibility problem and demonstrate the behavior of both the simultaneous and cyclic projection method with and without perturbations when we use them to solve the problem stated below. We show that in this example, the unperturbed versions of the methods converge slower than the perturbed versions. Furthermore, we demonstrate that the control sequence for the cyclic projection method has an influence on whether the condition triggering the perturbations is fulfilled.

4.1.1 Problem formulation

We consider the system of linear inequalities

[TABLE]

where

[TABLE]

and $x=(x_{1},x_{2},x_{3})^{T}\in\mathbb{R}^{3}$. In the following, by inequality $i$ we refer to the inequality $\langle a_{i},x\rangle\leq b_{i}$, where $a_{i}$ denotes the $i$-th row of the matrix $A$ and $b_{i}$ is the $i$-th coordinate of $b$ for $i\in\{1,2,3,4\}$.

These linear inequalities define half-spaces. The separating hyperplanes $H_{i}$ of these half-spaces intersect the $x_{1}$-axis at $\pm\delta_{x_{1}}$, the $x_{2}$-axis at $\pm\delta_{x_{2}}$ and the $x_{3}$-axis at $\delta_{x_{3}}$.

In this example, we choose $\delta_{x_{3}}=100$, $\delta_{x_{1}}=\tan(\beta)\delta_{x_{3}}/\sin(\alpha)$ and $\delta_{x_{2}}=\tan(\beta)\delta_{x_{3}}/\cos(\alpha)$ with $\alpha=30^{\circ}$ and $\beta=5^{\circ}$.

Next we demonstrate the behavior of simultaneous and cyclic projection with and without perturbations when we use them to solve (4.1).

4.1.2 Results

In what follows, we occasionally (especially in the descriptions and legends of the figures and tables) use abbreviations for the simultaneous projection (SP) and the cyclic projection method (CP) as well as for the heavy ball perturbation (HB) and the surrogate constraint perturbation (SC). CP+HB means for example that cyclic projection was used together with heavy ball perturbation.

For all methods, we choose the starting point $x^{0}=(15,0,0)^{T}$ and the parameters $\epsilon_{max}=6\cdot 10^{-2},\epsilon_{min}=10^{-6}$ and $\lambda^{SC}_{k}=\|p(x^{k})\|^{2}/\|d^{SC}\|^{2}$ as described in Section 3.3. We consider an iterate $x^{\ast}$ to be a solution of (4.1), if $\|Ax^{\ast}-b\|_{\infty}\leq 10^{-10}$ and $K$ denotes the number of iterations needed to find $x^{\ast}$.

In contrast to the step size $\lambda^{SC}_{k}$ of the surrogate constraint perturbation, the choice of the step size $\lambda^{HB}_{k}$ of the heavy ball perturbation is not motivated in a geometrical way. It can only be chosen in an optimal way, if the solution to the problem is already known. We present results for three different values of $\lambda^{HB}_{k}$ to demonstrate the effects of choosing it in an empirical way.

First, we use the simultaneous projection to solve (4.1). We choose $\lambda(x)\equiv 1.9$ to generate a sequence of iterates, which alternates between fulfilling inequalities $\{1,4\}$ and $\{2,3\}$. In this way, we obtain a sequence of projection steps $\{p(x^{k})\}$ with $\langle\bar{p}(x^{k}),\bar{p}(x^{k-1})\rangle\in[-1+\epsilon_{min},-1+\epsilon_{max}]$. If we use perturbations in this setting, they are therefore triggered after every second unperturbed iteration.

The results of our calculations are given in Table 1. The unperturbed method exhibits comparably slow convergence speed due to the opposing projection steps, which offer little progress in the direction of the $x_{3}$-axis. Both the heavy ball and the surrogate constraint perturbation are able to speed up the iteration process significantly.

Next, we solve (4.1) using the cyclic projection method with $\lambda(x)\equiv 1$ and the control sequence $\{j_{1}(\nu)\}_{\nu=0}^{\infty}=\{1,2,3,4,1,2,3,4,...\}$. This results in a sequence $\{p(x^{k})\}$ of projection steps, which does not fulfill $\langle\bar{p}(x^{k}),\bar{p}(x^{k-1})\rangle\in[-1+\epsilon_{min},-1+\epsilon_{max}]$ until the constraint violations fall below the tolerance of $10^{-10}$. Perturbations as defined in (3.12) or (3.23) can therefore not be used in a meaningful way with this method and choice of parameters. The values of $\langle\bar{p}(x^{k}),\bar{p}(x^{k-1})\rangle$ do not reflect the conflict inherent in the system of linear inequalities due to the choice of the control sequence. It takes the unperturbed algorithm 1917 iterations to find a feasible solution $x^{\ast}$.

We now use the same control sequence in combination with $\lambda(x)\equiv 1.9$. This results in a sequence of iterates $\{x^{k}\}$ and projection steps $\{p(x^{k})\}$, which fulfill $\langle\bar{p}(x^{k}),\bar{p}(x^{k-1})\rangle\in[-1+\epsilon_{min},-1+\epsilon_{max}]$ for some $k$. Therefore, we are able to use perturbations with this method and choice of parameters. The results of our calculations are presented in Table 2. This method with this particular choice of parameters is the quickest among our experiments to converge. We observe that for this choice of parameters, not all step sizes $\lambda_{k}^{HB}$ result in accelerated convergence behavior. Cyclic projection using the surrogate constraint perturbation, however, converges faster than any of the other variants of cyclic projection.

Finally, we alter the feasibility problem (4.1) by enlarging the linear inequality system with duplicates of the matrix rows $a_{i}$ and the values $b_{i}$ of the right hand side vector. We make this alteration in order to be able to give a control sequence, which both fits Definition 1 and results in the phenomena we will describe in the following. The new matrix $\tilde{A}\in\mathbb{R}^{8\times 3}$ and vector $\tilde{b}\in\mathbb{R}^{8}$ are given by

[TABLE]

The system of linear inequalities $\tilde{A}x\leq\tilde{b}$ of course has the same set of solutions as $Ax\leq b$. We solve $\tilde{A}x\leq\tilde{b}$ using the cyclic projection method with $\lambda=1.9$ and the control sequence $\{j_{2}(\nu)\}_{\nu=0}^{\infty}=\{1,2,3,4,5,6,7,8,1,2,3,4,5,6,7,8,...\}$. By repeating pairs of almost opposing inequalities, we generate a sequence of projection steps $\{p(x^{k})\}$, which fulfills $\langle\bar{p}(x^{k}),\bar{p}(x^{k-1})\rangle\in[-1+\epsilon_{min},-1+\epsilon_{max}]$ for some $k$. The number of iterations needed by the variants of the cyclic projection to find a feasible solution are given in Table 3. Here, all perturbations are able to accelerate the convergence. Again, cyclic projection using the surrogate constraint perturbation was the fastest to converge. In Figure 3 we present the norm of the constraint violation $\|\max\{0,Ax^{k}-b\}\|$ plotted against the iteration index $k$. HB8, HB80 and HB800 denote the three versions of the heavy ball perturbation using $\lambda_{k}^{HB}=8,80$ or $800$.

4.2 Intensity-modulated radiation therapy

Intensity-modulated radiation therapy (IMRT) is used as one treatment option in clinical oncology. With it, physicians aim to destroy the tumor tissue with irradiation while sparing surrounding healthy organs as much as possible. These goals are conflicting because of physical limitations like the proximity of the tumor to healthy organs and the attainable decrease of the dose values over a given distance. Therefore, the treatment planner usually has to decide for a compromise between several goals.

The IMRT planning problem, as we consider it in this paper, is to determine optimal intensity maps, i.e. a set of fluence intensity values, which causes a dose distribution in the patient’s body which will best fulfill the clinical goals formulated by the physician. In contrast to the previous examples, the IMRT planning problem is a nonlinear one.

4.2.1 Tools for modeling

The vector $d$ of dose values received by each voxel in the patient’s body when the fluence intensity given by the vector $x$ is applied can be calculated using the so-called dose matrix $P$. Due to physical reasons, the vector $x$ is restricted to be non-negative, which results in dose values $d=P\cdot x$, which are also non-negative.

We use the following set of functions to represent the dose prescriptions given by the clinical goals. All functions refer to a biological structure and evaluate the dose received by the voxels contained in the structure. From a mathematical point of view that structure is a set $\mathcal{O}$ of indices which correspond to the voxels contained in the biological structure. By $d_{i}=\langle p_{i},x\rangle$ we denote the dose value received by the voxel with index $i$, where $p_{i}$ is the $i$-th row of $P$.

Quadratic upper tail penalty function

The upper tail penalty function penalizes dose values of $d$ corresponding to the structure $\mathcal{O}$ which exceed a given threshold $U\in\mathbb{R}$.

[TABLE]

Quadratic lower tail penalty function

The lower tail penalty function penalizes dose values of $d$ corresponding to the structure $\mathcal{O}$ which fall below a given threshold $L\in\mathbb{R}$.

[TABLE]

Equivalent uniform dose (EUD)

The EUD function is related to the well-known EUD concept of Niemierko niemierko1997 . It penalizes dose values of $d$ corresponding to the structure $\mathcal{O}$ which deviate from 0.

[TABLE]

Tumor conformity

The tumor conformity function is used to ensure an even dose distribution within the tumor volume. It penalizes dose values of $d$ corresponding to the structure $\mathcal{O}$ which deviate from a given reference value.

[TABLE]

All of these functions take into account the system of linear (in)equalities resulting from $d_{i}=\langle p_{i},x\rangle$ being less/ greater or equal than a right hand side value determined by the function parameters, where $i\in\mathcal{O}$. The functions measure the distance of the vector $x$ to the (separating) hyperplanes corresponding to the linear (in)equalities in a nonlinear way. Figure 4 illustrates this concept with a simple 2D example.

The given functions do not represent a complete list of dose evaluation functions used in clinical applications. A more complete survey of such functions can be found in rdl04 and sfom99 .

4.2.2 The IMRT optimization problem

In this work we choose four head neck cases for our numerical experiments. We consider a reduced set of biological structures and their dose evaluation functions, which focuses on the main conflict between irradiating the tumor volume on the one hand and sparing the myelon and parotids on the other hand. Additionally, we include the healthy tissue not associated with any of the structures mentioned above with the goal to keep the dose in this tissue as low as possible. The dose evaluation functions we use are given in Table 4 in detail.

We use a reduced set of biological structures for our calculations, because our focus is to demonstrate the effects of the mathematical methods. Therefore, the treatment plans resulting from our calculations do not fulfill all of the clinical goals a treatment planner would formulate for a full head neck case.

The IMRT optimization problem as we formulate it in this work is

[TABLE]

where $f=\sum_{i\in I}f_{i}$, $I=\{1,2,3,4,5\},J=\{1,2,3\}$ and $f_{i},g_{j}$ as listed in Table 4.

4.2.3 Results

We use the combination of methods described in Section 3.5 to solve the problem (4.3) for four different head neck cases. Results of these calculations are presented in Tables 5, 6, 7 and 8 and, in an exemplary manner for one of the cases, in Figures 5, 6 and 7. For our calculations we choose $\lambda(x)\equiv 1.9,\lambda^{SC}_{k}=1,\lambda^{HB}_{k}=1$ and $\epsilon_{max}=0.034,\epsilon_{min}=10^{-8}$. We stop the algorithm if the projection method we use is unable to find a feasible solution of the current CFP after $10^{3}$ iterations. Then, the current CFP is assumed to be infeasible and we consider the solution of the previous CFP to be the result of our algorithm.

By $f^{\ast}_{\text{method, perturbation}}$ we denote the lowest objective function value, for which the algorithm using this method and perturbation is able to find a solution within the given maximum number of iterations per CFP $\boldsymbol{P}^{s}$ as described in (2.4). In the same way we denote by $K_{\text{method, perturbation}}$ the total number of iterations it takes the specified method using the specified perturbation to find the solution $x^{\ast}_{\text{method, perturbation}}$ with the optimal objective function value $f^{\ast}_{\text{method, perturbation}}$.

Tables 5 and 6 present the lowest objective function values achieved by the simultaneous or the cyclic projection method in an unperturbed manner or using either the heavy ball perturbation or the surrogate constraint perturbation. The values are given with a precision of $10^{0}$. Furthermore, the number of iterations needed to find the solutions are presented. For easier comparison, the percentages with respect to the values corresponding to the unperturbed methods are given in brackets.

Note that the optimal objective function values of the perturbed simultaneous projection method $f^{\ast}_{\text{SP, HB}}$ and $f^{\ast}_{\text{SP, SC}}$ (given in Table 5) are smaller than the optimal objective function value $f^{\ast}_{\text{SP}}$ of the unperturbed method. The same is true for the cyclic projection method.

When we use cyclic projection, the condition required to apply the perturbations given in equations (3.12) and (3.23) is never met for two of the four head neck cases. We denote the values for $f^{\ast}_{\text{method, perturbation}}$ and $K_{\text{method, perturbation}}$ by ”$-$”.

We observe that for all cases $K_{\text{SP, HB}},K_{\text{SP, SC}}<K_{\text{SP}}$, but for some cases $K_{\text{CP, HB}}$, $K_{\text{CP, SC}}>K_{\text{CP}}$. This means that in these cases, the perturbed method continues beyond solutions with the objective value $f^{\ast}_{\text{CP}}$, but in total takes more iterations than $K_{\text{CP}}$ to get there.

To answer the question, whether perturbed methods converge faster, measured at the same objective function value, we present Table 8. There we give the number of iterations it takes the perturbed methods to find a solution with an objective function value less or equal than $f^{\ast}_{\text{SP}}$ or $f^{\ast}_{\text{CP}}$, normalized by $K_{\text{SP}}$ and $K_{\text{CP}}$. These values are given with a precision of $10^{-4}$ and are an indicator for the acceleration of the iteration process caused by the perturbations we used. Note that in all cases the surrogate constraint perturbation is able to speed up the iteration process even more than the heavy ball perturbation.

Figure 5 illustrates the progress of the different methods in an exemplary way for Case 1. It is notable that perturbations are used rather early in the iteration process when we use simultaneous projection and comparably late when we use cyclic projection. The indices of iterations in which perturbations are used by the different methods are illustrated in more detail by Figure 6 for the same case as in Figure 5.

This phenomenon occurs due to the fact that simultaneous projection uses a weighted sum of all function gradients corresponding to violated constraints, whereas cyclic projection uses only gradient information of the next (with respect to the control sequence) violated constraint function.

In our model, there are groups of functions which correspond to conflicting goals. By summing the function gradients, simultaneous projection incorporates the information about the conflict between these groups of functions from a very early stage of the iteration process. While this results in almost opposing subsequent projection steps $p(x^{k-1}),p(x^{k})$, which leads to slow convergence, it also triggers the perturbation of the iteration process. When we use cyclic projection, the conflict between the groups of functions only becomes obvious when the subsequently (with respect to the control sequence) violated constraints have opposing function gradients. Figure 6 shows that this happens for IMRT cases rather late in the iteration process. The values in Table 8 indicate that methods using perturbations early converge faster than methods using them later in the iteration process.

Table 7 presents the number of iterations, in which perturbations are used. We observe that in all cases the heavy ball perturbation is applied more often than the surrogate constraint perturbation. Together with $K_{\text{SP, SC}}<K_{\text{SP, HB}}$ and $K_{\text{CP, SC}}<K_{\text{CP, HB}}$, this indicates that in our computations the surrogate constraint perturbation is more effective than the heavy ball perturbation.

Finally, in Figure 7, we present a cumulative dose volume histogram (DVH) in which we compare the dose distributions calculated for case 1 by simultaneous projection without perturbation and with the surrogate constraint perturbation. DVHs are a tool used by treatment planners to evaluate the quality of a fluence map and the resulting dose distribution in the patient’s body. DVHs show, which percentage of the volume of a certain structure receives a dose greater or equal than the dose value on the horizontal axis. For the tumor volume, a treatment planner might want that 95% of the volume receives at least a dose of 55 Gy and for the myelon, they might want that at most 5% of the volume receives a dose greater than 45 Gy.

The solid lines in Figure 7 represent the solution resulting from simultaneous projection without perturbation and the dashed lines correspond to the solution resulting from simultaneous projection with the surrogate constraint perturbation. We observe that the curves for the tumor volume and both parotids do not differ much, but the curves for the myelon are significantly lower for the perturbed method and therefore represent a more desirable dose distribution than the distribution resulting from the unperturbed method.

5 Conclusions

In this paper we transformed a general constrained convex optimization problem into a sequence of feasibility problems via the level set scheme. We solved each feasibility problem using either the simultaneous or cyclic subgradient projection method. We exploited the fact that both projection methods are bounded perturbation resilient and modified the iteration scheme of the projection methods with the heavy ball and the surrogate constraint perturbation.

Our numerical results demonstrate that the perturbed projection methods converge faster than their unperturbed counterparts, both for the linear feasibility problems we discussed and those arising from IMRT treatment planning. Applying perturbations early in the iteration process, as it was done for the simultaneous projection, yields better results than applying them later, like for the cyclic projection.

In our computations the perturbed versions of the simultaneous projection method offer the biggest improvement in both iteration numbers and objective function values of the solutions. The simultaneous projection method using the surrogate constraint perturbation surpasses the objective function value achieved by the unperturbed simultaneous projection method after only 14-22% of the number of iterations needed by the unperturbed method. The least improvement is achieved by the cyclic projection using the heavy ball perturbation. This method surpasses the objective value achieved by the unperturbed cyclic projection after 43-91% of the iterations needed by the unperturbed method.

In our calculations the surrogate constraint perturbation outperforms the heavy ball perturbation every time. Initial experiments indicate that this discrepancy results from the choice of the step length $\lambda^{HB}$. More detailed results are left for future publications. In our problem setting the computational cost for the surrogate constraint perturbation is slightly higher than for the heavy ball perturbation (12$n$+1 FLOPS vs. 11$n$ FLOPS, where $n$ is the length of the decision variable vector $x\in\mathbb{R}^{n}$). In problem settings where this difference is critical the heavy ball perturbation might be a more attractive approach.

Our observations suggest that both perturbations introduced by us can be used as an acceleration technique when zigzagging behavior occurs. Hence these and other perturbations should be further investigated.

Acknowledgment

We wish to thank the editor and the anonymous referees for the thorough analysis and review, their comments and suggestions helped tremendously in improving the quality of this paper and made it suitable for publication.

Conflict of Interest: The authors declare that they have no conflict of interest.