R-Linear and Q-Linear Convergence of a Double Golden-Ratio Tseng Method Beyond Monotonicity

1. Introduction

Let C be a nonempty, closed, and convex subset of a real Hilbert space H, and let $\mathcal{F}:H\to H$ be an operator. We consider the classical variational inequality problem (VIP) associated with $\mathcal{F}$ on C, which consists of finding a point $p^{*}\in C$ such that

$$\langle \mathcal{F}\left(p^{*}\right),p-p^{*}\rangle \ge 0,\forall p\in C.$$

(1)

The solution set of (1) is denoted by $\Omega$.

The corresponding dual variational inequality problem (DVIP) is to find a point $p^{*}\in C$ satisfying

$$\langle \mathcal{F}\left(p\right),p-p^{*}\rangle \ge 0,\forall p\in C.$$

(2)

The solution set of (2) is denoted by ${\Omega }_D$. It is well known that ${\Omega }_D$ is a closed and convex set, although it may be empty. Moreover, if $\mathcal{F}$ is continuous and C is convex, then ${\Omega }_D\subseteq \Omega .$ Furthermore, if $\mathcal{F}$ is continuous and pseudomonotone, then $\Omega ={\Omega }_D,$ see, for example, Lemma 2.1 in [14]. However, this identity generally fails for continuous quasimonotone operators; see, for instance, Example 4.2 in [53].

Variational inequalities provide a useful framework for problems in economics, engineering mechanics, mathematical programming, transportation science, image reconstruction, signal processing, machine learning, and network equilibrium; see, for example, [1,6,17,24,25,38,55]. This broad range of applications has motivated extensive work on efficient algorithms and their convergence properties. The resulting approaches include projection, extragradient, proximal, splitting, and inertial methods in finite- and infinite-dimensional Hilbert spaces; see, for example, [8,9,10,31,52] and the references therein.

Among the earliest and most fundamental iterative schemes for solving (1) is the gradient projection method, given by

$$\left\{\begin{array}{c}{v}_0\in C,\hfill \\ v_{k+1}=P_C\left(v_k-\tau \mathcal{F}\left(v_k\right)\right),k\ge 0,\hfill \end{array}\right.$$

(3)

where $P_C$ denotes the metric projection of H onto C and $\tau >0$ is a stepsize parameter. The convergence of the classical gradient projection method typically requires the operator $\mathcal{F}$ to be Lipschitz continuous and strongly monotone, or to be inverse strongly monotone.

To weaken these requirements, Korpelevich [26] and independently Antipin [5] introduced the celebrated extragradient method:

$$\left\{\begin{array}{c}{v}_0\in C,\hfill \\ u_k=P_C\left(v_k-\tau \mathcal{F}\left(v_k\right)\right),\hfill \\ v_{k+1}=P_C\left(v_k-\tau \mathcal{F}\left(u_k\right)\right),k\ge 0,\hfill \end{array}\right.$$

(4)

where $\mathcal{F}:C\to H$ is monotone and L-Lipschitz continuous, and the stepsize satisfies $0<\tau <{\displaystyle \frac{1}{L}}.$

The extragradient method has inspired many extensions; see, for example, [16,30,43,50] and the references therein. Its main computational drawback is that each iteration requires two metric projections onto C. If projecting onto C is expensive, the second projection can significantly increase the overall cost.

Several methods have therefore been developed to reduce the projection cost, including the subgradient extragradient method [9], Tseng’s forward–backward–forward method [47], and the projection-and-contraction method [21]. These methods reduce the number of projections while retaining desirable convergence properties under suitable assumptions. Even so, most available convergence analyses still rely on monotonicity-type assumptions together with global Lipschitz continuity of the underlying operator.

To enlarge the class of admissible problems, several authors have recently proposed new solution conditions together with extragradient-type methods for non-monotone and non-Lipschitz variational inequalities in real Hilbert spaces; see, for example, [4,46,49]. Moreover, a number of existing convergence analyses impose the additional requirement that the set $\Omega \setminus {\Omega }_D$ be finite; see, for instance, Condition (A5’) in [27]. This requirement is restrictive because it excludes variational inequality problems for which $\Omega \setminus {\Omega }_D$ may contain infinitely many points. Here we avoid the finiteness assumption by using the solution-oriented condition (A2), which relates $\mathcal{F}$ on the feasible set only to points in the solution set. This condition is implied by monotonicity and pseudomonotonicity, but is less restrictive than either of these assumptions [46].

In a parallel direction aimed at reducing the projection cost per iteration, Malitsky [32] introduced an elegant iterative framework based on the golden ratio for solving the mixed variational inequality problem of finding $p^{*}\in H$ such that

$$\langle \mathcal{F}\left(p^{*}\right),p-p^{*}\rangle +g\left(p\right)-g\left(p^{*}\right)\ge 0,\forall p\in H,$$

(5)

where $\mathcal{F}:H\to H$ is a monotone operator and $g:H\to (-\infty ,+\infty ]$ is a proper, convex, and lower semicontinuous function. To solve (5), he proposed the following golden-ratio algorithm:

$$\left\{\begin{array}{c}{\displaystyle {\overline{z}}_k=\frac{(\phi -1)z_k+{\overline{z}}_{k-1}}{\phi },}\hfill \\ z_{k+1}={prox}_{\lambda g}\left({\overline{z}}_k-\lambda \mathcal{F}\left(z_k\right)\right),\hfill \end{array}\right.$$

(6)

where

$$\phi ={\displaystyle \frac{1+\sqrt{5}}{2}}$$

is the golden ratio, and ${prox}_{\lambda g}$ denotes the proximal operator associated with g.

When $g={\iota }_C$, the indicator function of a nonempty closed and convex set C, problem (5) reduces to the classical variational inequality problem (1), while the proximal mapping becomes the metric projection:

$${prox}_{\lambda {\iota }_C}=P_C.$$

A useful feature of the golden-ratio algorithm is that the auxiliary iterate ${\overline{z}}_k$ is generated through a golden-ratio convex combination of previously computed points. The method therefore uses information from earlier iterates while requiring only one proximal or projection evaluation per step. Since its introduction, golden-ratio techniques have inspired numerous extensions for variational inequalities, equilibrium problems, fixed-point problems, and monotone inclusions; see, for example, [2,3,11,12,13,23,39,40,54] and the references therein.

Most existing golden-ratio methods use a single extrapolation step. To the best of our knowledge, no double golden-ratio Tseng-type extragradient method has yet been developed. This leads us to ask whether two golden-ratio averaging mechanisms can make better use of previous iterates while preserving a low per-iteration cost.

Linear convergence is important both theoretically and computationally. Existing R-linear convergence results for variational inequality algorithms are mainly established under relatively strong assumptions such as strong monotonicity or strong pseudomonotonicity, typically together with global Lipschitz continuity of the underlying operator; see, for example, [2,3,13,20,34,40,45,54]. Although these assumptions yield linear convergence, they restrict the range of admissible problems and force the solution set to be a single point. Below we obtain R-linear convergence of $dist(p_k,\Omega )$ and a Q-linear contraction of an associated weighted paired distance under a local error-bound condition of Robinson/Luo–Tseng type, without imposing strong monotonicity or requiring $\Omega$ to be a singleton.

Motivated by these observations, we develop a double golden-ratio Tseng-type extragradient method for variational inequalities in real Hilbert spaces. Its two averaging steps retain information from previous iterates without adding another projection. The analysis uses a solution-oriented condition in place of the usual monotonicity assumptions. The method also updates its stepsize adaptively; unlike (4), where $\tau <1/L$ must be chosen in advance, it can be implemented without knowing the Lipschitz constant.

The main contributions of this paper are summarized as follows.

(i)

We introduce a double golden-ratio Tseng-type extragradient algorithm that uses two averaging steps to retain information from previous iterates while requiring only one projection per iteration.

(ii)

Unlike [2,3,11,12,13,23,39,40,54], convergence is established under the solution-oriented condition (A2) rather than monotonicity or pseudomonotonicity assumptions. The analysis also does not require $\mathcal{F}\left(p\right)\ne 0$ throughout the feasible set; this case is handled directly by the stopping rule of the algorithm.

(iii)

An adaptive stepsize removes the need to know the Lipschitz constant of $\mathcal{F}$ in advance. The operator is still assumed globally Lipschitz continuous, but this constant never has to be computed or estimated to run the algorithm.

(iv)

Unlike the R-linear results in [2,3,13,20,34,40,45,54], we prove weak convergence of $\left\{p_k\right\}$ to a point of $\Omega$, and R-linear convergence of $dist(p_k,\Omega )$ to 0, under a Robinson/Luo–Tseng-type local error bound without assuming strong monotonicity or strong pseudomonotonicity and without forcing $\Omega$ to be a singleton. The error bound is weaker than the strong monotonicity-type assumptions commonly used to obtain linear convergence and does not require the solution set to be a singleton.

(v)

Under the same error bound, we obtain a Q-linear contraction of a weighted paired distance for all sufficiently large iterations. When $\Omega$ is a singleton, the weighted paired error remains Q-linear, while the iterates satisfy an explicit R-linear estimate.

(vi)

Numerical experiments on sequence-space and function-space variational inequality problems, together with applications to sparse signal reconstruction and multi-OD urban traffic network equilibrium, demonstrate the computational effectiveness of DGR–Tseng. In comparison with several well-known single golden-ratio methods, the proposed method exhibits favorable performance in terms of iteration counts, CPU time, and solution accuracy, as well as reconstruction and traffic-equilibrium performance measures in the respective applications.

The rest of the paper is organized as follows. Section 2 collects the preliminary facts, identities, and lemmas used throughout. Section 3 presents the proposed algorithm, while Section 4 establishes its weak convergence. Section 5 develops R-linear convergence under a local error-bound condition, and Section 6 sharpens this to Q-linear convergence. Section 7 reports numerical experiments, Section 8 presents an application to sparse signal reconstruction, and Section 9 applies the method to a multi-OD urban traffic network equilibrium problem. Finally, Section 10 concludes the paper.

2. Preliminaries

Throughout, $\mathbb{H}$ is a real Hilbert space with inner product $\langle ·,·\rangle$ and induced norm $\parallel ·\parallel$, and $\mathcal{D}$ is a nonempty closed convex subset of $\mathbb{H}$. We write $u_k\rightharpoonup u$ for weak convergence and $u_k\to u$ for strong convergence of a sequence $\left\{u_k\right\}\subset \mathbb{H}$. Recall the identity

$${\parallel \xi -\eta \parallel }^2={\parallel \xi \parallel }^2-2\langle \xi ,\eta \rangle +{\parallel \eta \parallel }^2,\xi ,\eta \in \mathbb{H},$$

(7)

and that for every $\xi \in \mathbb{H}$ there is a unique nearest point $P_{\mathcal{D}}\xi \in \mathcal{D}$; the operator $P_{\mathcal{D}}$ is nonexpansive.

Lemma 1. 

For every $\xi ,\eta \in \mathbb{H}$ and every $\lambda \in \mathbb{R}$,

$${\parallel \lambda \xi +(1-\lambda )\eta \parallel }^2={\lambda \parallel \xi \parallel }^2+{(1-\lambda )\parallel \eta \parallel }^2-\lambda (1-\lambda ){\parallel \xi -\eta \parallel }^2.$$

Lemma 2 

([19]). Let $z\in \mathcal{D}$ and $\xi \in \mathbb{H}$. Then

$$z=P_{\mathcal{D}}\xi \iff \langle \xi -z,z-\eta \rangle \ge 0\forall \eta \in \mathcal{D}.$$

Moreover, for every $\xi \in \mathbb{H}$ and $\eta \in \mathcal{D}$,

$$\parallel P_{\mathcal{D}}{\xi -\eta \parallel }^2\le {\parallel \xi -\eta \parallel }^2-{\parallel \xi -P_{\mathcal{D}}\xi \parallel }^2.$$

We omit the proof, which is standard.

Lemma 3 

([35]). Let $\left\{u_k\right\}\subset \mathbb{H}$ and let $\Omega \subset \mathbb{H}$ be nonempty. If

(a) 

${lim}_{k\to \infty }\parallel u_k-p\parallel$ exists for every $p\in \Omega$, and

(b) 

every weak sequential cluster point of $\left\{u_k\right\}$ lies in Ω,

then $\left\{u_k\right\}$ converges weakly to a point of Ω.

Lemma 4 

([44]). Let $\left\{a_k\right\}$, $\left\{b_k\right\}$ be nonnegative real sequences with $a_{k+1}\le a_k+b_k$ for all k and ${\sum }_{k=1}^{\infty }{b}_k<\infty$. Then ${lim}_{k\to \infty }{a}_k$ exists.

Lemma 5 

(Jensen’s inequality, [42] (Theorem 4.3)). Let $g:\mathbb{H}\to \mathbb{R}\cup \{+\infty \}$ be a convex function. If ${\xi }_1,\dots ,{\xi }_k\in domg$ and $t_1,\dots ,t_k\ge 0$ satisfy ${\sum }_{i=1}^k{t}_i=1$, then

$$g\left({\displaystyle \sum _{i=1}^k}{t}_i{\xi }_i\right)\le {\displaystyle \sum _{i=1}^k}{t}_{i}g\left({\xi }_i\right).$$

In particular, if $\Omega \subset \mathbb{H}$ is nonempty, closed and convex, the distance function $dist(·,\Omega ):\mathbb{H}\to [0,\infty )$ is convex (indeed 1-Lipschitz), so for every $\xi ,\eta \in \mathbb{H}$ and $t\in [0,1]$,

$$dist\left(t\xi +(1-t)\eta ,\Omega \right)\le tdist(\xi ,\Omega )+(1-t)dist(\eta ,\Omega ).$$

(8)

Lemma 6 

(Perron–Frobenius theorem, [22] (Theorem 8.4.4)). Let $N\in {\mathbb{R}}^{n\times n}$ have nonnegative entries and be irreducible. Then there is a real number $\rho \left(N\right)>0$, the Perron root of N, such that:

(i) 

$\rho \left(N\right)\ge \left|\xi \right|$ for every eigenvalue ξ of N (i.e. $\rho \left(N\right)$ equals the spectral radius of N);

(ii) 

$\rho \left(N\right)$ is a simple eigenvalue of N;

(iii) 

N has an eigenvector associated with $\rho \left(N\right)$ all of whose entries are strictly positive.

Lemma 7 

(Gelfand’s spectral radius formula, [18,22] (Corollary 5.6.14)). Let $N\in {\mathbb{R}}^{n\times n}$ and let $\parallel ·\parallel$ be any submultiplicative matrix norm. Then

$$\rho \left(N\right)=\underset{j\to \infty }{lim}{\parallel N^j\parallel }^{1/j}.$$

Consequently, for every $\sigma \in \left(\rho \right(N),1)$ there exists a constant $C_{\sigma }\ge 1$, depending only on N, σ and the chosen norm, such that

$$\parallel N^j\parallel \le C_{\sigma }{\sigma }^j,\forall j\ge 0.$$

Definition 1 

(Robinson/Luo–Tseng-type error bound, [28,29,41]). Let $\mathcal{D}\subset \mathbb{H}$ be nonempty, closed and convex, let $\mathcal{F}:\mathcal{D}\to \mathbb{H}$, and let $\Omega =Sol(\mathcal{D},\mathcal{F})\ne \varnothing$. The pair $(\mathcal{D},\mathcal{F})$ is said to admit aRobinson/Luo–Tseng-type local error boundat Ω if there exist $\kappa >0$, $\delta \in (0,\infty ]$ and a compact interval $[m,\Phi ]\subset (0,\infty )$ such that

$$dist(v,\Omega )\le \kappa ∥v-P_{\mathcal{D}}(v-\varphi \mathcal{F}v)∥$$

for every $\varphi \in [m,\Phi ]$ and every $v\in \mathcal{D}$ satisfying $\parallel v-P_{\mathcal{D}}(v-\varphi \mathcal{F}v)\parallel \le \delta$. The bound is calledglobalwhen $\delta =\infty$. Robinson’s polyhedral error-bound theorem [41] shows that a global bound of this type holds whenever $\mathcal{D}$ is polyhedral and $\mathcal{F}$ is affine; Luo and Tseng [28,29] extended this to piecewise-affine and, more generally, essentially smooth $\mathcal{F}$ on polyhedral $\mathcal{D}$, and to the associated variational-inequality residuals $v-P_{\mathcal{D}}(v-\varphi \mathcal{F}v)$ used above.

3. Proposed Method

We work under the following assumptions.

Assumption 1.

(A1) 

$\mathcal{F}$ is L-Lipschitz continuous for some $L>0$:

$$\parallel \mathcal{F}u-\mathcal{F}v\parallel \le L\parallel u-v\parallel ,\forall u,v\in \mathbb{H};$$

(9)

(A2) 

$\langle \mathcal{F}v,v-p\rangle \ge 0$ for every $p\in \Omega$ and $v\in \mathcal{D}$;

(A3) 

$u_k\rightharpoonup u$ implies ${lim\; sup}_{k\to \infty }\langle \mathcal{F}{u}_k,u_k-p\rangle \ge \langle \mathcal{F}u,u-p\rangle$ for every $p\in \mathcal{D}$;

(A4) 

$\Omega =Sol(\mathcal{D},\mathcal{F})$ is nonempty, closed and convex.

The proposed method is given below.

Algorithm 1 Adaptive Double Golden-Ratio Tseng Method (DGR-Tseng)

Initialization: Choose $p_0,p_1,w_0\in \mathbb{H}$, ${\varphi }_1>0$, $1<{\phi }_1,{\phi }_2\le {\displaystyle \frac{1+\sqrt{5}}{2}}$, $\mu \in (0,1)$, $\left\{{\lambda }_k\right\}\subset [0,\infty )$, $\left\{{\theta }_k\right\}\subset [1,\infty )$, $\left\{{\gamma }_k\right\}\subset [0,\infty )$ such that ${\lambda }_k\to 0$, ${\theta }_k\to 1$, ${\displaystyle {\displaystyle \sum _{k=1}^{\infty }}{\gamma }_k<\infty }$. Set $k=1$. Iterative steps: Given the current iterates $p_k$, $w_{k-1}$, and ${\varphi }_k$, calculate $w_k,q_k,r_k,p_{k+1},{\varphi }_{k+1}$ as follows: Step 1. Compute $$w_k={\displaystyle \frac{{\phi }_1-1}{{\phi }_1}}{p}_k+{\displaystyle \frac{1}{{\phi }_1}}{w}_{k-1}.$$ (10) Step 2. Compute $$q_k=P_{\mathcal{D}}\left(w_k-{\varphi }_k\mathcal{F}{w}_k\right).$$ (11) If $q_k=w_k$, then stop. Otherwise, go to Step 3. Step 3. Compute $$r_k=q_k-{\varphi }_k\left(\mathcal{F}{q}_k-\mathcal{F}{w}_k\right).$$ (12) Step 4. Set $$p_{k+1}={\displaystyle \frac{{\phi }_2-1}{{\phi }_2}}{r}_k+{\displaystyle \frac{1}{{\phi }_2}}{p}_k.$$ (13) Step 5. Update $${\varphi }_{k+1}=\left\{\begin{array}{cc}min\left\{({\lambda }_k+{\theta }_k\mu ){\displaystyle \frac{\parallel w_k-q_k\parallel }{\parallel \mathcal{F}{w}_k-\mathcal{F}{q}_k\parallel }},(1+{\gamma }_k){\varphi }_k\right\},\hfill & \mathcal{F}{w}_k\ne \mathcal{F}{q}_k,\hfill \\ (1+{\gamma }_k){\varphi }_k,\hfill & \mathcal{F}{w}_k=\mathcal{F}{q}_k.\hfill \end{array}\right.$$ (14) Set $k\leftarrow k+1$ and return to Step 1.

Remark 1. 

We fix throughout the constants attached to ${\phi }_1,{\phi }_2\in (1,\phi ]$, $\phi =(1+\sqrt{5})/2$:

$${\lambda }_1:={\displaystyle \frac{{\phi }_1-1}{{\phi }_1}},{\lambda }_2:={\displaystyle \frac{{\phi }_2-1}{{\phi }_2}},{\nu }_1:={\lambda }_1(1-{\lambda }_1),{\nu }_2:={\lambda }_2(1-{\lambda }_2),c:={\displaystyle \frac{{\lambda }_2(1-{\lambda }_1)}{{\lambda }_1}}.$$

As ${\phi }_1,{\phi }_2>1$, all five are strictly positive and independent of k, and (10), (13) take the form

$$w_k={\lambda }_1{p}_k+(1-{\lambda }_1)w_{k-1},p_{k+1}={\lambda }_2{r}_k+(1-{\lambda }_2)p_k.$$

Remark 2. 

• If $q_k=w_k$, then $w_k=P_{\mathcal{D}}(w_k-{\varphi }_k\mathcal{F}{w}_k)$, and Lemma 2 gives $\langle \mathcal{F}{w}_k,v-w_k\rangle \ge 0$ for all $v\in \mathcal{D}$, so $w_k\in \Omega$ and the stopping rule is well defined.
• Condition (A2) is less restrictive than monotonicity or pseudomonotonicity, since it only relates points of $\mathcal{D}$ to points of the solution set.
• When ${\phi }_1={\phi }_2=\phi$, the two averaging steps use the same classical golden-ratio weighting, ${\lambda }_1={\lambda }_2=1/{\phi }^2$, and $c=1-{\lambda }_1=1/\phi$. Accordingly, the Lyapunov function of Lemma 10 becomes $\parallel p_k{-u\parallel }^2+{\displaystyle \frac{1}{\phi }}{\parallel w_{k-1}-u\parallel }^2$.

4. Weak Convergence Analysis

We begin with the behaviour of the step-size sequence, which under the global Lipschitz assumption (A1) can be bounded away from zero. The update rule (14) does not multiply the ratio $\parallel w_k-q_k\parallel /\parallel \mathcal{F}{w}_k-\mathcal{F}{q}_k\parallel$ by the bare constant $\mu$, but by the k-dependent factor ${\lambda }_k+{\theta }_k\mu$, so it is this factor, not $\mu$ alone, that has to be tracked. To keep the notation light we write ${\mu }_k:={\lambda }_k+{\theta }_k\mu ,k\ge 1.$ Because ${\lambda }_k\ge 0$ and ${\theta }_k\ge 1$, we always have ${\mu }_k\ge \mu$; and because ${\lambda }_k\to 0$, ${\theta }_k\to 1$, we also have ${\mu }_k\to \mu$. Both facts are used repeatedly below.

Lemma 8. 

Let $\left\{{\varphi }_k\right\}$ be generated by (14) and set $m:=min\{{\varphi }_1,\mu /L\}$ and $\Phi :={\varphi }_{1}exp\left({\sum }_{j=1}^{\infty }{\gamma }_j\right)$. Then $m\le {\varphi }_k\le \Phi$ for every k, and $\left\{{\varphi }_k\right\}$ converges to some $\overline{\varphi }\in [m,\Phi ]$. Furthermore

$$\parallel \mathcal{F}{w}_k-\mathcal{F}{q}_k\parallel \le {\displaystyle \frac{{\mu }_k}{{\varphi }_{k+1}}}\parallel w_k-q_k\parallel ,\forall k\ge 1.$$

(15)

Proof. 

We induct on the floor bound. At $k=1$ it holds because $m\le {\varphi }_1$. Assume ${\varphi }_k\ge m$. If $\mathcal{F}{w}_k=\mathcal{F}{q}_k$, the second branch of (14) gives ${\varphi }_{k+1}=(1+{\gamma }_k){\varphi }_k\ge {\varphi }_k\ge m$ outright. If instead $\mathcal{F}{w}_k\ne \mathcal{F}{q}_k$, then $w_k\ne q_k$, and (A1) bounds the numerator: $\parallel \mathcal{F}{w}_k-\mathcal{F}{q}_k\parallel \le L\parallel w_k-q_k\parallel$, so

$${\mu }_k{\displaystyle \frac{\parallel w_k-q_k\parallel }{\parallel \mathcal{F}{w}_k-\mathcal{F}{q}_k\parallel }}\ge {\displaystyle \frac{{\mu }_k}{L}}\ge {\displaystyle \frac{\mu }{L}},$$

the last step using ${\mu }_k\ge \mu$ noted above. Since (14) takes the minimum of this quantity with $(1+{\gamma }_k){\varphi }_k\ge {\varphi }_k\ge m$, we get ${\varphi }_{k+1}\ge min\{\mu /L,m\}=m$, as $m\le \mu /L$ by the choice of m. This closes the induction, and in particular the correction term ${\lambda }_k$ and the slack factor ${\theta }_k$ never push the stepsize below the same floor m that a fixed constant $\mu$ would have given.

For the upper bound, (14) gives ${\varphi }_{k+1}\le (1+{\gamma }_k){\varphi }_k$ in either branch. Hence ${\varphi }_k\le {\varphi }_1{\prod }_{j=1}^{k-1}(1+{\gamma }_j)\le {\varphi }_{1}exp\left({\sum }_{j=1}^{\infty }{\gamma }_j\right)=\Phi$. Therefore ${\varphi }_{k+1}\le {\varphi }_k+\Phi {\gamma }_k$. Since ${\sum }_{k=1}^{\infty }\Phi {\gamma }_k<\infty$, Lemma 4 yields the convergence of $\left\{{\varphi }_k\right\}$ to some $\overline{\varphi }\in [m,\Phi ]$.

Finally, (15) is immediate when $\mathcal{F}{w}_k=\mathcal{F}{q}_k$, and otherwise follows directly from the definition ${\varphi }_{k+1}\le {\mu }_k\parallel w_k-q_k\parallel /\parallel \mathcal{F}{w}_k-\mathcal{F}{q}_k\parallel$ by rearranging. □

Remark 3. 

The global Lipschitz condition (A1) is essential for this lower bound. Under the weaker uniform continuity assumption used in [4,46], one can only show that ${lim}_k{\varphi }_k\ge 0$, leaving open the possibility that the step size decays to zero. That possibility must then be addressed separately in the convergence proof. Here, the explicit bound ${\varphi }_k\ge m>0$ rules it out.

Lemma 9. 

Suppose (A2) holds. For every $u\in \Omega$ and $k\ge 1$,

$$\parallel r_k{-u\parallel }^2\le \parallel w_k{-u\parallel }^2-\left(1-{\left({\mu }_k\right)}^2{\left({\displaystyle \frac{{\varphi }_k}{{\varphi }_{k+1}}}\right)}^2\right){\parallel w_k-q_k\parallel }^2.$$

(16)

Proof. 

Since $q_k=P_{\mathcal{D}}(w_k-{\varphi }_k\mathcal{F}{w}_k)$ and $u\in \Omega \subset \mathcal{D}$, Lemma 2 gives

$$\begin{array}{cc}\hfill \parallel q_k{-u\parallel }^2& \le \parallel w_k-{\varphi }_k\mathcal{F}{w}_k{-u\parallel }^2-{\parallel w_k-{\varphi }_k\mathcal{F}{w}_k-q_k\parallel }^2\hfill \\ \hfill & =\parallel w_k{-u\parallel }^2-2{\varphi }_k\langle \mathcal{F}{w}_k,w_k-u\rangle +{\varphi }_k^2{\parallel \mathcal{F}{w}_k\parallel }^2\hfill \\ \hfill & -\parallel w_k-q_k{\parallel }^2+2{\varphi }_k\langle \mathcal{F}{w}_k,w_k-q_k\rangle -{\varphi }_k^2{\parallel \mathcal{F}{w}_k\parallel }^2\hfill \\ \hfill & =\parallel w_k{-u\parallel }^2-{\parallel w_k-q_k\parallel }^2-2{\varphi }_k\langle \mathcal{F}{w}_k,q_k-u\rangle .\hfill \end{array}$$

(17)

By (12), we have

$$\begin{array}{cc}\hfill \parallel r_k{-u\parallel }^2& =\parallel (q_k-u)-{\varphi }_k(\mathcal{F}{q}_k-\mathcal{F}{w}_k){\parallel }^2\hfill \\ \hfill & =\parallel q_k{-u\parallel }^2-2{\varphi }_k\langle \mathcal{F}{q}_k-\mathcal{F}{w}_k,q_k-u\rangle +{\varphi }_k^2{\parallel \mathcal{F}{q}_k-\mathcal{F}{w}_k\parallel }^2.\hfill \end{array}$$

Inserting (17) and noting that $-2{\varphi }_k\langle \mathcal{F}{w}_k,q_k-u\rangle -2{\varphi }_k\langle \mathcal{F}{q}_k-\mathcal{F}{w}_k,q_k-u\rangle =-2{\varphi }_k\langle \mathcal{F}{q}_k,q_k-u\rangle$, we get

$$\parallel r_k{-u\parallel }^2\le \parallel w_k{-u\parallel }^2-\parallel w_k-q_k{\parallel }^2-2{\varphi }_k\langle \mathcal{F}{q}_k,q_k-u\rangle +{\varphi }_k^2{\parallel \mathcal{F}{q}_k-\mathcal{F}{w}_k\parallel }^2.$$

By (A2), $\langle \mathcal{F}{q}_k,q_k-u\rangle \ge 0$ because $q_k\in \mathcal{D}$ and $u\in \Omega$, so the corresponding term may be discarded:

$$\parallel r_k{-u\parallel }^2\le \parallel w_k{-u\parallel }^2-\parallel w_k-q_k{\parallel }^2+{\varphi }_k^2{\parallel \mathcal{F}{q}_k-\mathcal{F}{w}_k\parallel }^2.$$

(18)

Finally, (15) bounds the last term by ${\left({\mu }_k\right)}^2{({\varphi }_k/{\varphi }_{k+1})}^2{\parallel w_k-q_k\parallel }^2$, and substituting this into (18) gives (16). □

It is at this point that the correction terms ${\lambda }_k,{\theta }_k$ in the stepsize rule matter, and it is also here that they stop mattering: by Lemma 8, ${\varphi }_k\to \overline{\varphi }>0$, so ${\varphi }_k/{\varphi }_{k+1}\to 1$; and ${\mu }_k\to \mu$ by construction. Hence ${\left({\mu }_k\right)}^2{({\varphi }_k/{\varphi }_{k+1})}^2\to {\mu }^2$, and since $\mu \in (0,1)$,

$$1-{\left({\mu }_k\right)}^2{\left({\displaystyle \frac{{\varphi }_k}{{\varphi }_{k+1}}}\right)}^2⟶1-{\mu }^2>0.$$

Fixing any $\theta \in (0,1-{\mu }^2)$, this convergence gives some $k_0\in \mathbb{N}$ with

$$1-{\left({\mu }_k\right)}^2{\left({\displaystyle \frac{{\varphi }_k}{{\varphi }_{k+1}}}\right)}^2\ge \theta ,\forall k\ge k_0,$$

(19)

so that (16) yields, for all $k\ge k_0$ and $u\in \Omega$,

$$\parallel r_k{-u\parallel }^2\le \parallel w_k{-u\parallel }^2-\theta {\parallel w_k-q_k\parallel }^2.$$

(20)

Lemma 10. 

Suppose (A2) holds and let c be as in Remark 1. For $u\in \Omega$ and $k\ge 1$ set $L_k\left(u\right):=\parallel p_k{-u\parallel }^2+c{\parallel w_{k-1}-u\parallel }^2$. Then there are constants ${\sigma }_1,{\sigma }_2,{\sigma }_3>0$, not depending on k or u, such that

$$L_{k+1}\left(u\right)\le L_k\left(u\right)-{\sigma }_1\parallel p_k-w_{k-1}{\parallel }^2-{\sigma }_2\parallel w_k-q_k{\parallel }^2-{\sigma }_3{\parallel r_k-p_k\parallel }^2,\forall k\ge k_0,u\in \Omega .$$

(21)

Consequently $\left\{L_k\left(u\right)\right\}$ converges for every $u\in \Omega$,

$$\underset{k\to \infty }{lim}\parallel p_k-w_{k-1}\parallel =\underset{k\to \infty }{lim}\parallel w_k-q_k\parallel =\underset{k\to \infty }{lim}\parallel r_k-p_k\parallel =0,$$

(22)

the sequences $\left\{p_k\right\}$, $\left\{w_k\right\}$, $\left\{q_k\right\}$, $\left\{r_k\right\}$ are bounded, and ${lim}_{k\to \infty }{\parallel p_k-u\parallel }^2$ exists for every $u\in \Omega$.

Proof. 

Fix $u\in \Omega$ and put ${\Delta }_k:={\parallel p_k-u\parallel }^2$, ${\Theta }_k:={\parallel w_{k-1}-u\parallel }^2$, so $L_k\left(u\right)={\Delta }_k+c{\Theta }_k$. By Lemma 1, applied first to $w_k={\lambda }_1{p}_k+(1-{\lambda }_1)w_{k-1}$,

$$\parallel w_k{-u\parallel }^2={\lambda }_1{\Delta }_k+(1-{\lambda }_1){\Theta }_k-{\nu }_1{\parallel p_k-w_{k-1}\parallel }^2,$$

(23)

and then to $p_{k+1}={\lambda }_2{r}_k+(1-{\lambda }_2)p_k$, together with (20) for $k\ge k_0$,

$$\begin{array}{cc}\hfill {\Delta }_{k+1}& ={\lambda }_2\parallel r_k{-u\parallel }^2+(1-{\lambda }_2){\Delta }_k-{\nu }_2{\parallel r_k-p_k\parallel }^2\hfill \\ \hfill & \le {\lambda }_2\parallel w_k{-u\parallel }^2+(1-{\lambda }_2){\Delta }_k-{\lambda }_2\theta \parallel w_k-q_k{\parallel }^2-{\nu }_2{\parallel r_k-p_k\parallel }^2.\hfill \end{array}$$

Inserting (23) for $\parallel w_k{-u\parallel }^2$,

$${\Delta }_{k+1}\le \left[{\lambda }_2{\lambda }_1+1-{\lambda }_2\right]{\Delta }_k+{\lambda }_2(1-{\lambda }_1){\Theta }_k-{\lambda }_2{\nu }_1\parallel p_k-w_{k-1}{\parallel }^2-{\lambda }_2\theta \parallel w_k-q_k{\parallel }^2-{\nu }_2{\parallel r_k-p_k\parallel }^2,$$

and adding $c{\Theta }_{k+1}=c{\lambda }_1{\Delta }_k+c(1-{\lambda }_1){\Theta }_k-c{\nu }_1{\parallel p_k-w_{k-1}\parallel }^2$, which is again (23) shifted by one index, gives

$$\begin{array}{cc}\hfill L_{k+1}\left(u\right)& \le \left[{\lambda }_2{\lambda }_1+1-{\lambda }_2+c{\lambda }_1\right]{\Delta }_k+\left[{\lambda }_2(1-{\lambda }_1)+c(1-{\lambda }_1)\right]{\Theta }_k\hfill \\ \hfill & -({\lambda }_2+c){\nu }_1\parallel p_k-w_{k-1}{\parallel }^2-{\lambda }_2\theta \parallel w_k-q_k{\parallel }^2-{\nu }_2{\parallel r_k-p_k\parallel }^2.\hfill \end{array}$$

The two bracketed coefficients simplify considerably once we recall that c was chosen so that $c{\lambda }_1={\lambda }_2(1-{\lambda }_1)$. Indeed

$${\lambda }_2{\lambda }_1+1-{\lambda }_2+c{\lambda }_1=1-{\lambda }_2(1-{\lambda }_1)+c{\lambda }_1=1,$$

and

$${\lambda }_2(1-{\lambda }_1)+c(1-{\lambda }_1)=(1-{\lambda }_1)({\lambda }_2+c)=(1-{\lambda }_1){\lambda }_2+(1-{\lambda }_1)c=c{\lambda }_1+(1-{\lambda }_1)c=c,$$

where the last step uses $(1-{\lambda }_1){\lambda }_2=c{\lambda }_1$ once more. So the coefficients of ${\Delta }_k,{\Theta }_k$ collapse to $1,c$ exactly, and we arrive at (21) with ${\sigma }_1:=({\lambda }_2+c){\nu }_1$, ${\sigma }_2:={\lambda }_2\theta$, ${\sigma }_3:={\nu }_2$, all strictly positive by Remark 1.

Because ${\sigma }_1,{\sigma }_2,{\sigma }_3>0$, the sequence ${\left\{L_k\left(u\right)\right\}}_{k\ge k_0}$ is nonincreasing and nonnegative, hence convergent, and telescoping (21) from $k_0$ to ∞ gives

$${\sigma }_1\sum _{k\ge k_0}\parallel p_k-w_{k-1}{\parallel }^2+{\sigma }_2\sum _{k\ge k_0}\parallel w_k-q_k{\parallel }^2+{\sigma }_3\sum _{k\ge k_0}{\parallel r_k-p_k\parallel }^2\le L_{k_0}\left(u\right)<\infty ,$$

which yields (22).

Boundedness of $\left\{p_k\right\}$ follows from ${\Delta }_k\le L_k\left(u\right)\le L_{k_0}\left(u\right)$ for $k\ge k_0$, and boundedness of $\left\{w_k\right\}$ from (10), since $w_k$ is a convex combination of $p_k$ and $w_{k-1}$. Boundedness of $\left\{q_k\right\}$ and $\left\{r_k\right\}$ is then immediate from (22).

For the last claim, note $\left\{p_k\right\},\left\{w_{k-1}\right\}$ bounded gives some $M>0$ with $\parallel w_{k-1}-u\parallel +\parallel p_k-u\parallel \le M$, and hence

$$|{\Theta }_k-{\Delta }_k|=\left|\parallel w_{k-1}-u\parallel -\parallel p_k-u\parallel \right|·\left(\parallel w_{k-1}-u\parallel +\parallel p_k-u\parallel \right)\le M\parallel p_k-w_{k-1}\parallel \to 0.$$

Writing $L_k\left(u\right)=(1+c){\Delta }_k+c({\Theta }_k-{\Delta }_k)$, convergence of $\left\{L_k\left(u\right)\right\}$ together with ${\Theta }_k-{\Delta }_k\to 0$ forces $\left\{{\Delta }_k\right\}$ to converge as well. □

Lemma 11. 

Suppose (A1)–(A4) hold, and let $\left\{w_{k_j}\right\}$ be a subsequence of $\left\{w_k\right\}$ with $w_{k_j}\rightharpoonup z$ and $\parallel w_{k_j}-q_{k_j}\parallel \to 0$. Then $z\in \Omega$.

Proof. 

From $\parallel w_{k_j}-q_{k_j}\parallel \to 0$ and $w_{k_j}\rightharpoonup z$ we get $q_{k_j}\rightharpoonup z$, and since $\left\{q_k\right\}\subset \mathcal{D}$ with $\mathcal{D}$ weakly closed, $z\in \mathcal{D}$.

Let $\xi \in \mathcal{D}$ be arbitrary. Since $q_{k_j}=P_{\mathcal{D}}(w_{k_j}-{\varphi }_{k_j}\mathcal{F}{w}_{k_j})$, Lemma 2 gives

$$〈w_{k_j}-{\varphi }_{k_j}\mathcal{F}{w}_{k_j}-q_{k_j},\xi -q_{k_j}〉\le 0,$$

i.e. ${\varphi }_{k_j}\langle \mathcal{F}{w}_{k_j},\xi -q_{k_j}\rangle \ge \langle w_{k_j}-q_{k_j},\xi -q_{k_j}\rangle$. Dividing by ${\varphi }_{k_j}$ and splitting $\xi -q_{k_j}=(\xi -w_{k_j})+(w_{k_j}-q_{k_j})$ inside the left-hand inner product gives

$${\displaystyle \frac{1}{{\varphi }_{k_j}}}\langle q_{k_j}-w_{k_j},\xi -q_{k_j}\rangle +\langle \mathcal{F}{w}_{k_j},w_{k_j}-q_{k_j}\rangle \ge \langle \mathcal{F}{w}_{k_j},w_{k_j}-\xi \rangle .$$

(24)

The sequence $\left\{w_{k_j}\right\}$ is bounded, so by (A1) so is $\left\{\mathcal{F}{w}_{k_j}\right\}$; combined with $\parallel w_{k_j}-q_{k_j}\parallel \to 0$ this gives a bounded $\left\{q_{k_j}\right\}$ and $\parallel \mathcal{F}{w}_{k_j}-\mathcal{F}{q}_{k_j}\parallel \to 0$. By Lemma 8, ${\varphi }_{k_j}\ge m>0$ for every j, so $1/{\varphi }_{k_j}$ is bounded, and Cauchy–Schwarz applied to both terms on the left of (24) shows that each vanishes as $j\to \infty$:

$${\displaystyle \frac{1}{{\varphi }_{k_j}}}\left|\langle q_{k_j}-w_{k_j},\xi -q_{k_j}\rangle \right|\le {\displaystyle \frac{1}{m}}\parallel q_{k_j}-w_{k_j}\parallel \parallel \xi -q_{k_j}\parallel \to 0,$$

$$\left|\langle \mathcal{F}{w}_{k_j},w_{k_j}-q_{k_j}\rangle \right|\le \parallel \mathcal{F}{w}_{k_j}\parallel \parallel w_{k_j}-q_{k_j}\parallel \to 0.$$

Hence

$$\underset{j\to \infty }{lim\; sup}\langle \mathcal{F}{w}_{k_j},w_{k_j}-\xi \rangle \le 0.$$

(25)

Applying (A3) to $u^k=w_{k_j}$, $u=z$, $p=\xi$, and using (25),

$$\langle \mathcal{F}z,z-\xi \rangle \le \underset{j\to \infty }{lim\; sup}\langle \mathcal{F}{w}_{k_j},w_{k_j}-\xi \rangle \le 0.$$

As $\xi \in \mathcal{D}$ was arbitrary, $z\in Sol(\mathcal{D},\mathcal{F})=\Omega$. □

Theorem 1. 

Under Assumption 1, if $\mu \in (0,1)$, ${\phi }_1,{\phi }_2\in (1,\phi ]$ and ${\sum }_{k=1}^{\infty }{\gamma }_k<\infty$, the sequence $\left\{p_k\right\}$ generated by Algorithm 1 converges weakly to a point of Ω.

Proof. 

Condition (a) of Lemma 3 is exactly what Lemma 10 gives us: ${lim}_{k\to \infty }{\parallel p_k-u\parallel }^2$ exists for every $u\in \Omega$.

For condition (b), take any weakly convergent subsequence $p_{k_j}\rightharpoonup \zeta$. We would like to place $\zeta$ in $\Omega$ via Lemma 11, but that lemma speaks about $w_k$, not $p_k$ – so we first move the subsequence over. Since $\parallel p_k-w_{k-1}\parallel \to 0$ by (22), $w_{k_j-1}$ inherits the same weak limit $\zeta$; and since $\parallel w_k-q_k\parallel \to 0$ as well, the residual $\parallel w_{k_j-1}-q_{k_j-1}\parallel \to 0$ too. Both hypotheses of Lemma 11 are met by the subsequence indexed by $k_j-1$, so $\zeta \in \Omega$.

With both conditions of Lemma 3 in hand, $\left\{p_k\right\}$ converges weakly to a point of $\Omega$. □

5. R-Linear Convergence

By Lemma 8, every stepsize generated by the algorithm satisfies ${\varphi }_k\in [m,\Phi ]$, where $m:=min\{{\varphi }_1,\mu /L\}$ and $\Phi :={\varphi }_{1}exp\left({\sum }_{j=1}^{\infty }{\gamma }_j\right)$. We impose a local Robinson/Luo–Tseng-type error bound uniformly over this fixed interval.

Assumption A2. 

There exist $\kappa >0$ and $\delta >0$ such that, for every $\varphi \in [m,\Phi ]$ and every $v\in \mathcal{D}$ satisfying $\parallel v-P_{\mathcal{D}}(v-\varphi \mathcal{F}v)\parallel \le \delta$,

$$dist(v,\Omega )\le \kappa ∥v-P_{\mathcal{D}}(v-\varphi \mathcal{F}v)∥.$$

(26)

Example 1. 

Let $\mathbb{H}={\mathbb{R}}^2$, $\mathcal{D}={\mathbb{R}}^2$, and define $\mathcal{F}:{\mathbb{R}}^2\to {\mathbb{R}}^2$ by

$$\mathcal{F}(p_1,p_2)=(p_1,0).$$

Since $\mathcal{D}=\mathbb{H}$, we have $P_{\mathcal{D}}=I$, and the variational inequality is equivalent to $\mathcal{F}p=0$. Hence, $\Omega =\left\{\right(0,t):t\in \mathbb{R}\},$ which is not a singleton. For every $p=(p_1,p_2)\in {\mathbb{R}}^2$, $dist(p,\Omega )=|p_1|.$ Moreover, for every $\varphi >0$,

$$\begin{array}{cc}\hfill p-P_{\mathcal{D}}(p-\varphi \mathcal{F}p)& =p-(p-\varphi \mathcal{F}p)\hfill \\ & =\varphi \mathcal{F}p\hfill \\ & =(\varphi p_1,0),\hfill \end{array}$$

and therefore

$$∥p-P_{\mathcal{D}}(p-\varphi \mathcal{F}p)∥=\varphi \left|p_1\right|.$$

Thus, for every $\varphi \in [m,\Phi ]$,

$$dist(p,\Omega )\le {\displaystyle \frac{1}{m}}∥p-P_{\mathcal{D}}(p-\varphi \mathcal{F}p)∥,$$

so Assumption 2 holds globally with $\kappa =1/m$. On the other hand, $\mathcal{F}$ is not strongly monotone. Indeed, for $p=(0,1)$ and $q=(0,0)$,

$$\langle \mathcal{F}p-\mathcal{F}q,p-q\rangle ={0,\parallel p-q\parallel }^2=1.$$

Hence no $\rho >0$ can satisfy

$$\langle \mathcal{F}p-\mathcal{F}q,p-q\rangle \ge {\rho \parallel p-q\parallel }^2$$

for all $p,q\in {\mathbb{R}}^2$. Therefore, the error-bound condition may hold even when strong monotonicity fails and the solution set is not a singleton.

Remark 4. 

Example 1 shows that Assumption 2 may hold even when strong monotonicity fails and the solution set is not a singleton. Thus the linear-rate analysis below is based on an error-bound property rather than on the strong monotonicity-type assumptions commonly used to obtain linear convergence.

Recall $dist(p,\Omega )={inf}_{u\in \Omega }\parallel p-u\parallel =\parallel p-P_{\Omega }p\parallel$, and that $dist(·,\Omega )$ is convex and 1-Lipschitz because $\Omega$ is closed and convex (Lemma 5).

Lemma 12. 

Suppose (A2) and Assumption 2 hold. Let

$$\theta :=\frac{1}{2}(1-{\mu }^2)\in (0,1-{\mu }^2),$$

and let $k_0$ be as in (19) for this choice of θ. Then there exists $k_1\ge k_0$ such that, for every $k\ge k_1$,

$$dist(r_k,\Omega )\le \gamma dist(w_k,\Omega ),\gamma :=\sqrt{1-\theta {\widehat{\kappa }}^{-2}}\in (0,1),$$

(27)

where $\widehat{\kappa }:=1+\kappa (1+\Phi L)$.

Proof. 

For $v\in \mathcal{D}$ and $\varphi >0$, write $R_{\varphi }\left(v\right):=\parallel v-P_{\mathcal{D}}(v-\varphi \mathcal{F}v)\parallel$. Since $q_k=P_{\mathcal{D}}(w_k-{\varphi }_k\mathcal{F}{w}_k)$, the nonexpansiveness of $P_{\mathcal{D}}$ and (A1) give

$$R_{{\varphi }_k}\left(q_k\right)=\parallel q_k-P_{\mathcal{D}}(q_k-{\varphi }_k\mathcal{F}{q}_k)\parallel \le (1+{\varphi }_{k}L)\parallel w_k-q_k\parallel \le (1+\Phi L)\parallel w_k-q_k\parallel .$$

By (22), $\parallel w_k-q_k\parallel \to 0$. Hence there exists $k_1\ge k_0$ such that $R_{{\varphi }_k}\left(q_k\right)\le \delta$ for every $k\ge k_1$. Because $q_k\in \mathcal{D}$ and ${\varphi }_k\in [m,\Phi ]$, Assumption 2 yields

$$dist(q_k,\Omega )\le \kappa R_{{\varphi }_k}\left(q_k\right)\le \kappa (1+\Phi L)\parallel w_k-q_k\parallel ,k\ge k_1.$$

Therefore

$$dist(w_k,\Omega )\le \parallel w_k-q_k\parallel +dist(q_k,\Omega )\le \widehat{\kappa }\parallel w_k-q_k\parallel ,\widehat{\kappa }:=1+\kappa (1+\Phi L).$$

Write ${\overline{u}}_k:=P_{\Omega }\left(w_k\right)$, so $\parallel w_k-{\overline{u}}_k\parallel =dist(w_k,\Omega )=:D_k^w$. Since ${\overline{u}}_k\in \Omega$, inequality (20) gives, for $k\ge k_1$,

$$\parallel r_k-{\overline{u}}_k{\parallel }^2\le {\left(D_k^w\right)}^2-\theta {\parallel w_k-q_k\parallel }^2\le \left(1-\theta {\widehat{\kappa }}^{-2}\right){\left(D_k^w\right)}^2.$$

Since $\widehat{\kappa }\ge 1$ and $0<\theta <1$, the number $\gamma =\sqrt{1-\theta {\widehat{\kappa }}^{-2}}$ belongs to $(0,1)$. Finally, $dist(r_k,\Omega )\le \parallel r_k-{\overline{u}}_k\parallel \le \gamma D_k^w$, which proves (27). □

Lemma 13. 

Set $D_k^p:=dist(p_k,\Omega )$, $D_k^w:=dist(w_k,\Omega )$, and $y_k:={(D_k^p,D_{k-1}^w)}^{\top }$. Then, for all $k\ge k_1$,

$$y_{k+1}\le N{y}_k\left(\mathrm{entrywise}\right),N:=\left(\begin{array}{cc}(1-{\lambda }_2)+{\lambda }_2{\lambda }_1\gamma & {\lambda }_2(1-{\lambda }_1)\gamma \\ {\lambda }_1& 1-{\lambda }_1\end{array}\right),$$

(28)

with ${\lambda }_1,{\lambda }_2$ as in Remark 1 and γ as in Lemma 12.

Proof. 

Since $\Omega$ is closed and convex, $dist(·,\Omega )$ is convex, so Jensen’s inequality (Lemma 5, specifically (8)) applied to $w_k={\lambda }_1{p}_k+(1-{\lambda }_1)w_{k-1}$ gives

$$D_k^w\le {\lambda }_1{D}_k^p+(1-{\lambda }_1)D_{k-1}^w.$$

(29)

Applying Jensen’s inequality (Lemma 5) to $p_{k+1}={\lambda }_2{r}_k+(1-{\lambda }_2)p_k$ and then Lemma 12,

$$D_{k+1}^p\le {\lambda }_{2}dist(r_k,\Omega )+(1-{\lambda }_2)D_k^p\le {\lambda }_2\gamma D_k^w+(1-{\lambda }_2)D_k^p.$$

Since ${\lambda }_2\gamma \ge 0$, inserting (29) preserves the inequality:

$$D_{k+1}^p\le \left[(1-{\lambda }_2)+{\lambda }_2{\lambda }_1\gamma \right]D_k^p+{\lambda }_2(1-{\lambda }_1)\gamma D_{k-1}^w.$$

Together with (29) shifted by one index, this is exactly (28). □

Lemma 14. 

The matrix N of Lemma 13 has nonnegative, irreducible entries, and $\rho \left(N\right)<1$.

Proof. 

Nonnegativity is clear since ${\lambda }_1,{\lambda }_2\in (0,1)$ and $\gamma \in (0,1)$; irreducibility holds because both off-diagonal entries ${\lambda }_2(1-{\lambda }_1)\gamma$ and ${\lambda }_1$ are strictly positive. A direct computation gives

$$detN=(1-{\lambda }_1)(1-{\lambda }_2)=:D_0\in (0,1),trN=T\left(\gamma \right):=2-{\lambda }_1-{\lambda }_2+{\lambda }_1{\lambda }_2\gamma ,$$

where T is strictly increasing in $\gamma$. By the Perron–Frobenius theorem (Lemma 6), N has a real, positive, simple dominant eigenvalue ${\xi }_1\left(\gamma \right)$ at least as large in modulus as any other eigenvalue; for a real $2\times 2$ matrix the remaining eigenvalue ${\xi }_2\left(\gamma \right)=T\left(\gamma \right)-{\xi }_1\left(\gamma \right)$ is then automatically real, so ${\xi }_1,{\xi }_2$ are the two roots of ${\xi }^2-T\left(\gamma \right)\xi +D_0=0$ and, by part (i) of Lemma 6,

$$\rho \left(N\left(\gamma \right)\right)={\xi }_1\left(\gamma \right)={\displaystyle \frac{T\left(\gamma \right)+\sqrt{T{\left(\gamma \right)}^2-4D_0}}{2}},$$

which is strictly increasing in T (its T-derivative is $\frac{1}{2}+\frac{T}{2\sqrt{T^2-4D_0}}>0$). At $\gamma =1$, one checks directly (as in the coefficient identity of Lemma 10) that $\xi =1$ solves ${\xi }^2-T\left(1\right)\xi +D_0=0$, so $\{{\xi }_1\left(1\right),{\xi }_2\left(1\right)\}=\{1,D_0\}$ and $\rho \left(N\right(1\left)\right)=1$. Since $\gamma <1$ implies $T\left(\gamma \right)<T\left(1\right)$, monotonicity gives $\rho \left(N\left(\gamma \right)\right)={\xi }_1\left(\gamma \right)<{\xi }_1\left(1\right)=1$. □

Theorem 2. 

Suppose (A1)–(A4) and Assumption 2 hold, and let $\left\{p_k\right\}$ be generated by Algorithm 1 with $\mu \in (0,1)$, ${\phi }_1,{\phi }_2\in (1,\phi ]$ and ${\sum }_{k=1}^{\infty }{\gamma }_k<\infty$. Then there exist $C>0$ and $\sigma \in \left(\rho \right(N),1)$ such that

$$dist(p_k,\Omega )\le C{\sigma }^{k-k_1},\forall k\ge k_1.$$

(30)

That is, $\left\{p_k\right\}$ converges to Ω R-linearly. Consequently, together with Theorem 1, $\left\{p_k\right\}$ converges weakly to some $p^{*}\in \Omega$ with $dist(p_k,\Omega )\to 0$ at the geometric rate (30); if in addition $\Omega =\left\{u^{*}\right\}$ is a singleton, this reads $\parallel p_k-u^{*}\parallel \le C{\sigma }^{k-k_1}\to 0$, i.e. strong convergence at a linear rate.

Proof. 

Since N has nonnegative entries, entrywise inequalities propagate under repeated multiplication: from Lemma 13, $y_{k_1+j}\le N^j{y}_{k_1}$ entrywise for every $j\ge 0$. By Lemma 14, $\rho \left(N\right)<1$, so by Gelfand’s spectral radius formula (Lemma 7), for any $\sigma \in \left(\rho \right(N),1)$ there is $C_{\sigma }\ge 1$ with $\parallel N^j\parallel \le C_{\sigma }{\sigma }^j$ for all $j\ge 0$ (in any fixed matrix norm). The vector $y_{k_1}$ is finite and nonnegative, since $\left\{p_k\right\}$ and $\left\{w_k\right\}$ are bounded by Lemma 10 (which uses only (A2)). Hence each entry of $N^j{y}_{k_1}$ is bounded by $C{\sigma }^j$ for a suitable constant $C=C(y_{k_1},C_{\sigma })$; in particular $D_{k_1+j}^p\le C{\sigma }^j$ for all $j\ge 0$. Writing $k=k_1+j$ gives (30). □

Remark 5. 

Assumption 2 is imposed only on points of $\mathcal{D}$. Lemma 12 applies it at $q_k\in \mathcal{D}$ and then transfers the resulting estimate to $w_k$ by nonexpansiveness of the projection and Lipschitz continuity of $\mathcal{F}$. Thus the error bound is never applied to an iterate outside its stated domain. Moreover, Lemma 13 uses Jensen’s inequality (Lemma 5) for the convex function $dist(·,\Omega )$, so the matrix recursion propagates distances rather than squared distances.

6. Q-Linear Convergence

Theorem 2 of Section 5 is an R-linear statement: it only bounds $dist(p_k,\Omega )$ by a fixed geometric envelope, for some $\sigma$ strictly larger than $\rho \left(N\right)$, via Gelfand’s formula (Lemma 7), which is inherently asymptotic. We now show that the same matrix inequality (28) in fact yields a genuine, non-asymptotic Q-linear contraction with the explicit Perron factor $\rho \left(N\right)$ once $D_k^p$ and $D_{k-1}^w$ are combined through the left Perron eigenvector of N. Throughout this section we retain Assumption 2 and all the notation of Section 5 ($m,\Phi ,N,\rho \left(N\right),k_1,D_k^p,D_k^w$, etc.).

Lemma 15. 

Let N be the matrix of Lemma 13. The transpose $N^{\top }$ is again nonnegative and irreducible, since transposition preserves both properties. Hence Lemma 6, applied to $N^{\top }$, yields $\rho \left(N^{\top }\right)=\rho \left(N\right)$ together with an eigenvector $v={(v_1,v_2)}^{\top }$ of $N^{\top }$ associated with $\rho \left(N\right)$ all of whose entries are strictly positive, i.e.

$$v^{\top }N=\rho \left(N\right)v^{\top },v_1,v_2>0.$$

(31)

Theorem 3. 

Suppose (A1)–(A4) and Assumption 2 hold, let $\left\{p_k\right\}$ be generated by Algorithm 1 as in Theorem 2, and let $v={(v_1,v_2)}^{\top }>0$ be the left Perron eigenvector of Lemma 15. Define, for $k\ge k_1$,

$$Y_k:=v_1{D}_k^p+v_2{D}_{k-1}^w=v_{1}dist(p_k,\Omega )+v_{2}dist(w_{k-1},\Omega )\ge 0.$$

(32)

Then

$$Y_{k+1}\le \rho \left(N\right)Y_k,\forall k\ge k_1,$$

(33)

that is, ${\left\{Y_k\right\}}_{k\ge k_1}$ converges Q-linearly to 0 with ratio $\rho \left(N\right)\in (0,1)$. Consequently

$$dist(p_k,\Omega )\le {\displaystyle \frac{1}{v_1}}{Y}_{k_1}\rho {\left(N\right)}^{k-k_1},\forall k\ge k_1,$$

(34)

which provides an explicit R-linear estimate for $dist(p_k,\Omega )$ with geometric factor $\rho \left(N\right)$ and constant $Y_{k_1}/v_1$, instead of an arbitrary $\sigma \in \left(\rho \right(N),1)$ and the Gelfand constant $C_{\sigma }$.

Proof. 

Fix $k\ge k_1$. By Lemma 13, $y_{k+1}\le N{y}_k$ entrywise, where $y_k={(D_k^p,D_{k-1}^w)}^{\top }$. Since $v_1,v_2\ge 0$, left-multiplying an entrywise inequality between nonnegative vectors by $v^{\top }$ preserves it:

$$v^{\top }{y}_{k+1}\le v^{\top }N{y}_k.$$

By (31), $v^{\top }N=\rho \left(N\right)v^{\top }$, so the right-hand side equals $\rho \left(N\right)v^{\top }{y}_k$. Since $v^{\top }{y}_k=Y_k$ by (32), this is precisely (33).

Iterating (33) from $k_1$ gives $Y_k\le \rho {\left(N\right)}^{k-k_1}{Y}_{k_1}$ for all $k\ge k_1$. Since $v_2{D}_{k-1}^w\ge 0$, we have $v_1{D}_k^p\le Y_k$, and dividing by $v_1>0$ gives (34). □

Remark 6. 

Neither Theorem 2 nor Theorem 3 assumes strong monotonicity or strong pseudomonotonicity of $\mathcal{F}$. Theorem 2 gives R-linear convergence of $dist(p_k,\Omega )$, whereas Theorem 3 gives Q-linear convergence of the weighted paired distance $Y_k$. Assumption 2 does not force Ω to be a singleton, so both conclusions remain meaningful for a non-singleton solution set.

Corollary 1. 

Suppose, in addition to the hypotheses of Theorem 3, that $\Omega =\left\{u^{★}\right\}$ is a singleton. Then $D_k^p=\parallel p_k-u^{★}\parallel$ and $D_{k-1}^w=\parallel w_{k-1}-u^{★}\parallel$, so

$$Y_k=v_1\parallel p_k-u^{★}\parallel +v_2\parallel w_{k-1}-u^{★}\parallel$$

satisfies the Q-linear contraction (33) for every $k\ge k_1$. Thus the weighted paired error contracts by a factor of at most $\rho \left(N\right)$ at each such iteration, while (34) gives the explicit R-linear estimate $\parallel p_k-u^{★}\parallel \le (Y_{k_1}/v_1)\rho {\left(N\right)}^{k-k_1}$ for all $k\ge k_1$.

7. Numerical Examples

This section illustrates the convergence results established in Section 5Section 6 and examines the numerical performance of the proposed Adaptive Double Golden-Ratio Tseng method (DGR-Tseng), given in Algorithm 1. We compare DGR-Tseng with three golden-ratio-type methods from the literature, namely:

• the Golden Ratio Projection Algorithm (GRPA) of Chu and Zhang [13];
• the Golden-Ratio Tseng Method (GRT-Tseng) of Abiodun et al. [3];
• the Golden-Ratio Subgradient Extragradient Method (SEM-GRT) of Oyewole and Reich [36].

Unlike DGR-Tseng, which uses the double golden-ratio mechanism (10)–(13), the competing methods employ a single golden-ratio extrapolation step.

All computations were carried out in MATLAB R2018a. In Example 2, the space ${\ell }^2$ was approximated by its first $N=500$ coordinates. In Example 3, the interval $[0,1]$ was discretized using 501 equally spaced grid points, and the $L^2$-norm and the required integrals were approximated by the composite trapezoidal rule.

For each test case, the same initial data were used for all four algorithms. The common stopping criterion was $\parallel p_{k+1}-p_k\parallel <{10}^{-8}$, with a maximum of 1500 iterations.

Let $\phi ={\displaystyle \frac{1+\sqrt{5}}{2}}.$ The control parameters used in the experiments are listed in Table 1.

The selected sequences satisfy ${\lambda }_k\to 0,{\theta }_k\to 1,{\sum }_{k=1}^{\infty }{\gamma }_k<\infty ,$ for DGR-Tseng, and ${\sum }_{k=1}^{\infty }{\beta }_k<\infty ,{\sum }_{k=1}^{\infty }({\delta }_k-1)<\infty$ for GRT-Tseng and SEM-GRT.

We consider two infinite-dimensional test problems. The first is posed in ${\ell }^2$, while the second is formulated in $L^2(0,1)$. In both examples, the operator is globally Lipschitz continuous but nonmonotone, and the solution set is $\Omega =\left\{0\right\}$.

Example 2. 

Let $\mathbb{H}=\mathcal{D}={\ell }^2=\left\{p=(p_1,p_2,\dots ):{\sum }_{j=1}^{\infty }{\left|p_j\right|}^2<\infty \right\},$ and define $h\left(t\right)=t+2sint,t\in \mathbb{R}.$ Consider $\mathcal{F}\left(p\right)=\left(h\left(p_1\right),p_2,p_3,\dots \right).$ Set $c_0:=1-{\displaystyle \frac{2}{\pi }}>0.$

We first note that

$$th\left(t\right)\ge c_0{t}^2,\left|h\left(t\right)\right|\ge c_0\left|t\right|,t\in \mathbb{R}.$$

(35)

If $\left|t\right|\le \pi$, then $tsint\ge 0$, so $th\left(t\right)=t^2+2tsint\ge t^2\ge c_0{t}^2.$ If $\left|t\right|\ge \pi$, then $2tsint\ge -2\left|t\right|\ge -{\displaystyle \frac{2}{\pi }}{t}^2,$ and again $th\left(t\right)\ge c_0{t}^2.$ Since $h^{\prime }\left(t\right)=1+2cost,\left|h^{\prime }\left(t\right)\right|\le 3,$ we have

$$\begin{array}{cc}\hfill {\parallel \mathcal{F}p-\mathcal{F}y\parallel }^2& =|h\left(p_1\right)-h\left(y_1\right){|}^2+{\displaystyle \sum _{j=2}^{\infty }}{|p_j-y_j|}^2\hfill \\ \hfill & \le 9|p_1-y_1{|}^2+{\displaystyle \sum _{j=2}^{\infty }}{|p_j-y_j|}^2\hfill \\ \hfill & \le {9\parallel p-y\parallel }^2.\hfill \end{array}$$

Thus $\mathcal{F}$ is 3-Lipschitz continuous.

Because $\mathcal{D}=\mathbb{H}$, the variational inequality reduces to $\mathcal{F}p=0$. From (35), $h\left(p_1\right)=0⟺p_1=0,$ and hence $p_j=0,j\ge 2.$ Therefore, $\Omega =\left\{0\right\}.$ Furthermore,

$$\begin{array}{cc}\hfill \langle \mathcal{F}p,p\rangle & =p_{1}h\left(p_1\right)+{\displaystyle \sum _{j=2}^{\infty }}{\left|p_j\right|}^2\hfill \\ \hfill & \ge c_0|p_1{|}^2+{\displaystyle \sum _{j=2}^{\infty }}{\left|p_j\right|}^2\hfill \\ \hfill & \ge c_0{\parallel p\parallel }^2.\hfill \end{array}$$

Hence Assumption(A2)holds.

Let $p^k\rightharpoonup p$ in ${\ell }^2$. Since $p_1^k\to p_1$, and with $Qp=(0,p_2,p_3,\dots ),$ we have $Q{p}^k\rightharpoonup Qp$. For $u=(u_1,u_2,\dots )\in {\ell }^2$,

$$\langle \mathcal{F}{p}^k,p^k-u\rangle =h\left(p_1^k\right)(p_1^k-u_1)+{\parallel Q{p}^k\parallel }^2-\langle Q{p}^k,Qu\rangle .$$

Using the weak lower semicontinuity of the norm gives

$$\underset{k\to \infty }{lim\; sup}\langle \mathcal{F}{p}^k,p^k-u\rangle \ge \langle \mathcal{F}p,p-u\rangle .$$

Thus Assumption(A3)holds, while Assumption(A4)follows from $\Omega =\left\{0\right\}$.

The operator is nonmonotone since $h^{\prime }\left(\pi \right)=-1<0.$ Thus there exist $a<b$ such that $\left(h\right(a)-h(b\left)\right)(a-b)<0.$ Setting $p=(a,0,0,\dots ),y=(b,0,0,\dots ),$ gives $\langle \mathcal{F}p-\mathcal{F}y,p-y\rangle <0.$

Finally, $\parallel \mathcal{F}p\parallel \ge c_0\parallel p\parallel .$ Since $P_{\mathcal{D}}=I$, for every $\varphi \in [m,\Phi ]$,

$$\parallel p-P_{\mathcal{D}}(p-\phi \mathcal{F}p)\parallel \ge m{c}_{0}dist(p,\Omega ).$$

Hence Assumption 2 holds globally with $\kappa ={\displaystyle \frac{1}{m{c}_0}}.$

For the numerical test, ${\ell }^2$ was truncated to its first $N=500$ coordinates. For $j=1,\dots ,N$, we used

$${\displaystyle \begin{array}{cc}\mathbf{Case}\mathbf{1}:\hfill & {\displaystyle {\left(p_0\right)}_j=\frac{3}{j},{\left(p_1\right)}_j=\frac{2}{j},{\left(w_0\right)}_j=\frac{1}{j},}\hfill \end{array}\mathbf{Case}\mathbf{2}:{\left(p_0\right)}_j=\frac{4{(-1)}^{j+1}}{j},{\left(p_1\right)}_j=\frac{3{(-1)}^{j+1}}{j},{\left(w_0\right)}_j=\frac{1.5{(-1)}^{j+1}}{j},\left[3mm\right]}$$

Case 3:(p0)j=3.5(0.35j)j, (p1)j=2.5(0.25j)j, (w0)j=1.5(0.50j)j,

Case 4:(p0)j= 3(0.20j)+(0.60j)j,

(p1)j= 2.5(0.30j)-0.8(0.70j)j, (w0)j= 1.5(0.40j)+0.5(0.90j)j.

Table 2 and Figure 1, Figure 2 and Figure 3 show that DGR-Tseng reaches the stopping tolerance with fewer iterations in all four cases. The CPU times and residual curves display the same overall trend.

Example 3. 

Let $\mathbb{H}=\mathcal{D}=L^2(0,1),$ e(t)=1. Then ${\parallel e\parallel }_{L^2}=1$. Define

$$M=span\left\{e\right\},M^{\perp }=\left\{p\in L^2(0,1):{\int }_0^{1}p\left(t\right)dt=0\right\}.$$

Every $p\in L^2(0,1)$ can be written uniquely as $p=a_{p}e+z_p,a_p={\int }_0^{1}p\left(t\right)dt,z_p\in M^{\perp }.$ Let $h\left(s\right)=s+2sins,$ and define

$$\mathcal{F}p=h\left(a_p\right)e+z_p.$$

(36)

For $p=a_{p}e+z_p$ and $y=a_{y}e+z_y$,

$$\begin{array}{cc}\hfill {\parallel \mathcal{F}p-\mathcal{F}y\parallel }_{L^2}^2& =|h\left(a_p\right)-h\left(a_y\right){|}^2+{\parallel z_p-z_y\parallel }_{L^2}^2\hfill \\ \hfill & \le {9\parallel p-y\parallel }_{L^2}^2.\hfill \end{array}$$

Hence $\mathcal{F}$ is 3-Lipschitz continuous.

Since $\mathcal{D}=\mathbb{H}$, the variational inequality reduces to $\mathcal{F}p=0$. Using $\left|h\left(s\right)\right|\ge c_0\left|s\right|,c_0=1-{\displaystyle \frac{2}{\pi }},$ we obtain $\Omega =\left\{0\right\}.$

Moreover,

$$\langle \mathcal{F}p,p\rangle \ge c_0{\parallel p\parallel }_{L^2}^2,$$

(37)

so Assumption(A2)holds.

If $p^k\rightharpoonup p$, write $p^k=a_{k}e+z_k,p=ae+z.$ Then $a_k\to a$ and $z_k\rightharpoonup z$. For $u=be+v$, $v\in M^{\perp }$,

$$\langle \mathcal{F}{p}^k,p^k-u\rangle =h\left(a_k\right)(a_k-b)+{\parallel z_k\parallel }_{L^2}^2-\langle z_k,v\rangle .$$

Weak lower semicontinuity yields the required limit inequality, so Assumption(A3)follows. Assumption (A4)follows from $\Omega =\left\{0\right\}$.

The operator is nonmonotone. Indeed, for sufficiently small $\epsilon >0$, take $p=(\pi -\epsilon )e,y=(\pi +\epsilon )e.$ Since $h^{\prime }\left(\pi \right)=-1$, $\langle \mathcal{F}p-\mathcal{F}y,p-y\rangle <0.$ Furthermore, ${\parallel \mathcal{F}p\parallel }_{L^2}\ge c_0{\parallel p\parallel }_{L^2},$ and therefore

$$\parallel p-P_{\mathcal{D}}{(p-\phi \mathcal{F}p)\parallel }_{L^2}\ge m{c}_{0}dist(p,\Omega ),\varphi \in [m,\Phi ].$$

Hence Assumption 2 holds globally with $\kappa ={\displaystyle \frac{1}{m{c}_0}}.$

For the numerical experiment, four starting-point cases were used:

$$\begin{array}{cc}\hfill \mathrm{Case}1\text{:}& \parallel p_0{\parallel }_2=3.000000,\parallel p_1{\parallel }_2=2.000000,{\parallel w_0\parallel }_2=1.000000,\hfill \\ \hfill \mathrm{Case}2\text{:}& \parallel p_0{\parallel }_2=2.912912,\parallel p_1{\parallel }_2=1.749286,{\parallel w_0\parallel }_2=1.236932,\hfill \\ \hfill \mathrm{Case}3\text{:}& \parallel p_0{\parallel }_2=2.478576,\parallel p_1{\parallel }_2=1.950215,{\parallel w_0\parallel }_2=1.134903,\hfill \\ \hfill \mathrm{Case}4\text{:}& \parallel p_0{\parallel }_2=2.598582,\parallel p_1{\parallel }_2=1.834649,{\parallel w_0\parallel }_2=1.633381.\hfill \end{array}$$

Table 3 and Figure 4, Figure 5 and Figure 6 show the same pattern as in Example 2. DGR-Tseng requires fewer iterations in all four cases and records the smallest CPU times in the reported experiments.

8. Application to Sparse Signal Reconstruction

We next consider an application to sparse signal reconstruction. Let $p^{†}\in {\mathbb{R}}^N$ denote the unknown sparse signal and suppose that

$$b=A{p}^{†}+\epsilon ,$$

(38)

where $A\in {\mathbb{R}}^{M\times N}$, $M<N$, is the sensing matrix and $\epsilon$ represents measurement noise.

We consider the regularized problem

$$\underset{p\in \mathcal{D}}{min}\left\{{\displaystyle \frac{1}{2}}{\parallel Ap-b\parallel }^2+\lambda {\mathcal{R}}_{ϵ}\left(p\right)\right\},$$

(39)

where $\mathcal{D}\subset {\mathbb{R}}^N$ is nonempty, closed, and convex, and ${\mathcal{R}}_{ϵ}\left(p\right)={\sum }_{i=1}^N{\displaystyle \frac{p_i^2}{p_i^2+{ϵ}^2}},ϵ>0$. Its gradient is given by ${[\nabla {\mathcal{R}}_{ϵ}\left(p\right)]}_i={\displaystyle \frac{2{ϵ}^2{p}_i}{{(p_i^2+{ϵ}^2)}^2}}.$ Define $\mathcal{F}\left(p\right)=A^{\top }(Ap-b)+\lambda \nabla {\mathcal{R}}_{ϵ}\left(p\right).$ The corresponding first-order condition can be written as

$$\mathrm{find}{p}^{*}\in \mathcal{D}\mathrm{such}\mathrm{that}\langle \mathcal{F}\left(p^{*}\right),p-p^{*}\rangle \ge 0,p\in \mathcal{D}.$$

(40)

The gradient of the least-squares term, $p\mapsto A^{\top }(Ap-b)$, is Lipschitz continuous with constant ${\parallel A\parallel }^2$. Moreover, $\nabla {\mathcal{R}}_{ϵ}$ is globally Lipschitz continuous because the derivative of each scalar component is bounded in modulus by $2/{ϵ}^2$. Hence $\mathcal{F}$ is globally Lipschitz continuous with Lipschitz constant at most ${\parallel A\parallel }^2+2\lambda /{ϵ}^2$. Since ${\mathcal{R}}_{ϵ}$ is nonconvex, the resulting operator need not be monotone.

We consider $N\in \{256,512,1024,2048\}$, and $M=\lfloor 0.5N\rfloor .$ The sensing matrix is generated according to $A_{ij}\sim {\displaystyle \frac{1}{\sqrt{M}}}\mathcal{N}(0,1),$ and additive Gaussian noise is included in the measurements. For visual comparison, the back-projected signal $p^{\mathrm{meas}}=A^{\top }b$ is also displayed.

The reconstruction quality is measured using

$$MSE\left(\widehat{p}\right)={\displaystyle \frac{1}{N}}{\parallel \widehat{p}-p^{†}\parallel }^2$$

(41)

and

$$SNR\left(\widehat{p}\right)=20{log}_{10}\left({\displaystyle \frac{\parallel p^{†}\parallel }{\parallel \widehat{p}-p^{†}\parallel }}\right)\mathrm{dB}.$$

(42)

Thus, smaller MSE and larger SNR indicate better reconstruction.

For each dimension, Figure 7, Figure 8, Figure 9 and Figure 10 show the original signal, the back-projected noisy measurement, and the reconstructions produced by the four methods.

The MSE and SNR comparisons are shown in Figure 11 and Figure 12, respectively.

Table 4 shows that DGR-Tseng gives the lowest MSE and the highest SNR in all four cases. For $N=256$, DGR-Tseng attains an MSE of $3.5088\times {10}^{-5}$ and an SNR of $39.300$ dB, while GRT-Tseng, the next-best method in this case, records $2.0888\times {10}^{-4}$ and $31.553$ dB.

The same pattern persists as the dimension increases. For $N=2048$, DGR-Tseng achieves an MSE of $2.9318\times {10}^{-5}$ and an SNR of $37.986$ dB, compared with SNR values of $29.244$, $24.578$, and $22.342$ dB for GRT-Tseng, SEM-GRT, and GRPA, respectively. Together with the reconstruction plots, these results show that DGR-Tseng more closely reproduces the significant coefficients of the original sparse signal in the reported experiment.

9. Application to a Multi-OD Urban Traffic Network Equilibrium Problem

Traffic network equilibrium is one of the classical applications of variational inequality theory. According to Wardrop’s user-equilibrium principle [51], a traffic flow is at equilibrium if no traveller can reduce the experienced travel cost by unilaterally switching to another admissible route connecting the same origin and destination. The classical optimization formulation of traffic assignment was developed by Beckmann et al. [7], while its more general variational inequality formulation was studied by Dafermos [15]; see also [33,37].

In this section, we apply the proposed DGR–Tseng method to a multi-origin–destination urban traffic network involving nonlinear congestion, asymmetric route interaction, and incident-induced capacity reductions. Besides providing a practical illustration of the variational inequality framework, the example permits a route-level and network-level comparison of DGR–Tseng with GRPA, GRT–Tseng, and SEM–GRT.

Let $\mathcal{G}=(\mathcal{N},\mathcal{A})$ be the directed traffic network with node set $\mathcal{N}=\{W,S,A,B,C,J,\mathrm{CBD}\}.$ The nodes W and S are the two origins, while $\mathrm{CBD}$ is the common destination. We consider the two origin–destination pairs $W⟶\mathrm{CBD},S⟶\mathrm{CBD},$ with traffic demands $d_1=1200\mathrm{veh}/\mathrm{h},d_2=900\mathrm{veh}/\mathrm{h}.$

The complete network consists of fourteen directed links and is displayed in Figure 13.

The free-flow travel times and nominal capacities are given in Table 5.

For the OD pair $W\to \mathrm{CBD}$, the admissible routes are

$$\begin{array}{cccc}\hfill P_1& :W-A-\mathrm{CBD},\hfill & \hfill P_2& :W-B-\mathrm{CBD},\hfill \\ \hfill P_3& :W-B-A-\mathrm{CBD},\hfill & \hfill P_4& :W-C-\mathrm{CBD},\hfill \end{array}$$

while for $S\to \mathrm{CBD}$ we consider

$$\begin{array}{cccc}\hfill P_5& :S-B-\mathrm{CBD},\hfill & \hfill P_6& :S-J-\mathrm{CBD},\hfill \\ \hfill P_7& :S-C-\mathrm{CBD},\hfill & \hfill P_8& :S-J-B-\mathrm{CBD}.\hfill \end{array}$$

These routes are displayed in Figure 14.

Let $p={(p_1,p_2,\dots ,p_8)}^{\top }\in {\mathbb{R}}^8$ denote the path-flow vector, where $p_i$ is the flow assigned to route $P_i$. The feasible path-flow set is

$$\mathcal{D}=\left\{p\in {\mathbb{R}}_{+}^8:{\displaystyle \sum _{i=1}^4}{p}_i=d_1,{\displaystyle \sum _{i=5}^8}{p}_i=d_2\right\}.$$

(43)

Thus $\mathcal{D}=D_1\times D_2,$ where $D_1$ and $D_2$ are simplices. In particular, $\mathcal{D}$ is nonempty, closed, bounded, convex, and polyhedral.

Let $\Delta \in {\mathbb{R}}^{14\times 8}$ be the link–path incidence matrix defined by

$${\delta }_{ai}=\left\{\begin{array}{cc}1,\hfill & \mathrm{if}\mathrm{link}a\mathrm{belongs}\mathrm{to}\mathrm{route}{P}_i,\hfill \\ 0,\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

For the eight routes under consideration,

$$\Delta =\left(\begin{array}{cccccccc}1& 0& 0& 0& 0& 0& 0& 0\\ 0& 1& 1& 0& 0& 0& 0& 0\\ 1& 0& 1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 1& 0& 0& 1\\ 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0& 1& 0\\ 0& 0& 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 1\\ 0& 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1\end{array}\right).$$

Hence, the link-flow vector corresponding to p is $v\left(p\right)=\Delta p.$

For every link $a\in \mathcal{A}$, its travel time is modelled by the BPR-type volume-delay function [48]

$$t_a\left(v_a\right)=t_a^0\left[1+0.15{\left({\displaystyle \frac{v_a}{c_a^{\mathrm{eff}}}}\right)}^4\right].$$

(44)

To represent an incident-sensitive bottleneck, the effective capacity of $B\to \mathrm{CBD}$ is reduced to $c_4^{\mathrm{eff}}=0.72c_4,$ whereas the effective capacity of the merging link $J\to B$ is reduced to $c_{14}^{\mathrm{eff}}=0.85c_{14}.$ For all other links, $c_a^{\mathrm{eff}}=c_a.$ The route-cost component generated by the link travel times is $C\left(p\right)={\Delta }^{\top }t(\Delta p).$

To incorporate interaction among different traffic streams, define the operator $\mathcal{F}:{\mathbb{R}}^8\to {\mathbb{R}}^8$ by

$$\mathcal{F}\left(p\right)={\Delta }^{\top }t(\Delta p)+Hp,$$

(45)

where $H\in {\mathbb{R}}^{8\times 8}$ has the nonzero entries

$$\begin{array}{cc}& H_{2,5}=0.0020,H_{5,2}=0.0010,\hfill \\ & H_{3,5}=0.0015,H_{5,3}=0.0008,\hfill \\ & H_{6,8}=0.0018,H_{8,6}=0.0010,\hfill \\ & H_{4,7}=0.0012,H_{7,4}=0.0006.\hfill \end{array}$$

The asymmetric term $Hp$ provides a simplified representation of cross-route interaction, merging interference, and spillback.

The traffic equilibrium problem is therefore to find $p^{*}\in \mathcal{D}$ such that

$$〈\mathcal{F}\left(p^{*}\right),p-p^{*}〉\ge 0,\forall p\in \mathcal{D}.$$

(46)

Thus, (46) is precisely a finite-dimensional instance of the variational inequality problem studied throughout this paper, with solution set $\Omega =Sol(\mathcal{D},\mathcal{F}).$ Since $\mathcal{D}=D_1\times D_2$, its metric projection satisfies $P_{\mathcal{D}}\left(z\right)=\left(\begin{array}{c}{P}_{D_1}(z_1,\dots ,z_4)\\ P_{D_2}(z_5,\dots ,z_8)\end{array}\right).$ Consequently, the single projection required at each DGR–Tseng iteration reduces to two standard simplex projections.

Since $\mathcal{D}$ is compact and the BPR travel-time functions are continuously differentiable on the relevant bounded traffic range, $\mathcal{F}$ is Lipschitz continuous on $\mathcal{D}$. However, the quartic BPR extension is not globally Lipschitz continuous on all of ${\mathbb{R}}^8$, and the iterates $p^k$ need not remain in $\mathcal{D}$. Accordingly, this section is intended as a computational application of the proposed method; we do not claim that this particular traffic model satisfies the global Lipschitz assumption (A1) or the error-bound hypothesis of Section 5Section 6 on the whole ambient space.

Algorithm 1 generates the sequence $\left\{p^k\right\}$ together with $\left\{w^k\right\}$, $\left\{q^k\right\}$, and $\left\{r^k\right\}$. Since $q^k\in \mathcal{D}$, whereas $p^k$ need not necessarily belong to $\mathcal{D}$, we use the feasible diagnostic path-flow vector ${\widehat{p}}^k:=P_{\mathcal{D}}\left(p^k\right)$ when reporting traffic-specific quantities. The normalized natural residual is defined by $\mathcal{R}\left(p\right)={\displaystyle \frac{∥p-P_{\mathcal{D}}\left(p-\tau \mathcal{F}\left(p\right)\right)∥}{1+\parallel p\parallel }},\tau =1.$ Accordingly, the residual reported at iteration k is $\mathcal{R}\left({\widehat{p}}^k\right)$.

For each OD pair, define $C_{min}^{\left(1\right)}\left(p\right)={min}_{1\le i\le 4}{\mathcal{F}}_i\left(p\right),C_{min}^{\left(2\right)}\left(p\right)={min}_{5\le i\le 8}{\mathcal{F}}_i\left(p\right).$ The associated Wardrop relative gaps are

$$R{G}_1\left(p\right)={\displaystyle \frac{{\displaystyle {\displaystyle \sum _{i=1}^4}{p}_i{\mathcal{F}}_i\left(p\right)-d_1{C}_{min}^{\left(1\right)}\left(p\right)}}{{\displaystyle {\displaystyle \sum _{i=1}^4}{p}_i{\mathcal{F}}_i\left(p\right)}}}$$

and

$$R{G}_2\left(p\right)={\displaystyle \frac{{\displaystyle {\displaystyle \sum _{i=5}^8}{p}_i{\mathcal{F}}_i\left(p\right)-d_2{C}_{min}^{\left(2\right)}\left(p\right)}}{{\displaystyle {\displaystyle \sum _{i=5}^8}{p}_i{\mathcal{F}}_i\left(p\right)}}}.$$

We set $RG\left(p\right)=max\{R{G}_1\left(p\right),R{G}_2\left(p\right)\}.$

The two network-level traffic indicators are $TSTT\left(p\right)={\sum }_{a\in \mathcal{A}}{v}_a\left(p\right)t_a\left(v_a\left(p\right)\right),$ and $MaxVC\left(p\right)={max}_{a\in \mathcal{A}}{\displaystyle \frac{v_a\left(p\right)}{c_a^{\mathrm{eff}}}}.$ Thus, the reported performance quantities are

$$\mathrm{Iter}.,\mathrm{CPU},\mathcal{R}\left({\widehat{p}}^k\right),RG\left({\widehat{p}}^k\right),TSTT\left({\widehat{p}}^k\right),MaxVC\left({\widehat{p}}^k\right).$$

The common tolerance is $\epsilon ={10}^{-7},$ and the traffic stopping test is $max\left\{\mathcal{R}\left({\widehat{p}}^k\right),RG\left({\widehat{p}}^k\right)\right\}<\epsilon .$ To examine sensitivity to initialization, we define the four feasible base path-flow vectors

$$\begin{array}{cc}\hfill {\overline{p}}^{\left(1\right)}& ={(300,300,300,300,225,225,225,225)}^{\top },\hfill \\ \hfill {\overline{p}}^{\left(2\right)}& ={(660,180,240,120,90,450,225,135)}^{\top },\hfill \\ \hfill {\overline{p}}^{\left(3\right)}& ={(120,240,660,180,315,90,135,360)}^{\top },\hfill \\ \hfill {\overline{p}}^{\left(4\right)}& ={(180,132,408,480,180,270,360,90)}^{\top }.\hfill \end{array}$$

For each case $j=1,2,3,4$, Algorithm 1 is initialized by $p_0=p_1=w_0={\overline{p}}^{\left(j\right)}.$ This choice is admissible under the initialization of Algorithm 1 and provides a uniform starting state for the main and auxiliary sequences. Moreover, ${\overline{p}}^{\left(j\right)}\in \mathcal{D},j=1,2,3,4,$ since ${\overline{p}}_i^{\left(j\right)}\ge 0,{\sum }_{i=1}^4{\overline{p}}_i^{\left(j\right)}=1200,{\sum }_{i=5}^8{\overline{p}}_i^{\left(j\right)}=900.$ For fairness, the same base path-flow vector ${\overline{p}}^{\left(j\right)}$ is used to initialize each competing method in Case j, according to its respective initialization rule.

For the experiment in this section, we maintain the same control parameters as in Section 7. The high-accuracy reference equilibrium is computed using tolerance ${10}^{-10}$ and a maximum of 20000 iterations.

For a fair numerical comparison, DGR–Tseng, GRPA, GRT–Tseng, and SEM–GRT are considered on the same traffic equilibrium problem using the same starting-point cases and stopping criterion. The reported iteration counts, CPU times, terminal natural residuals, Wardrop relative gaps, total system travel times (TSTT), and maximum volume-to-capacity ratios (MaxVC) are used to compare the methods in terms of computational efficiency, equilibrium accuracy, and the quality of the resulting network traffic states.

The numerical results for the four initial path-flow cases are reported in Table 6.

The results in Table 6 show that DGR–Tseng attains the smallest iteration count, CPU time, natural residual, and Wardrop relative gap in all four cases. In addition, the traffic states obtained by DGR–Tseng yield the smallest TSTT and MaxVC values among the four methods.

Figure 15 displays the natural residual and Wardrop relative gap for Case 1.

Both quantities are displayed on logarithmic scales. The trajectories exhibit method-dependent curvature, transient behaviour, changes of slope, and different termination points. In particular, the DGR–Tseng profile reaches the terminal accuracy in substantially fewer iterations.

Figure 16 compares the iteration counts and CPU times for the four initial path-flow cases.

Figure 17 shows the corresponding terminal accuracy.

Since $-{log}_{10}\left(\epsilon \right)$ increases as $\epsilon$ decreases, a taller bar represents a smaller error. Thus, DGR–Tseng attains the highest terminal accuracy in all four cases.

Let $p^{\mathrm{ref}}\in \mathcal{D}$ denote the high-accuracy reference equilibrium computed by DGR–Tseng. To supplement the direct solver comparison with a route-level sensitivity analysis around this reference equilibrium, define

$$p_{\mathrm{bud}}^{\left(m\right)}=P_{\mathcal{D}}\left(p^{\mathrm{ref}}+{\alpha }_m{d}^{\left(m\right)}\right),$$

(47)

where $d^{\left(m\right)}$ is a demand-preserving perturbation direction and $({\alpha }_1,{\alpha }_2,{\alpha }_3,{\alpha }_4)=(3,60,34,46).$ The four components are associated, respectively, with DGR–Tseng, GRPA, GRT–Tseng, and SEM–GRT. These reference-based states are used only for the supplementary route- and link-level sensitivity comparisons below; the direct algorithmic results are those reported in Table 6.

The associated generalized route-cost vectors are $C_{\mathrm{bud}}^{\left(m\right)}=\mathcal{F}\left(p_{\mathrm{bud}}^{\left(m\right)}\right).$

Figure 18 compares the resulting route flows and generalized route costs.

Neither a smaller nor a larger individual route flow is intrinsically better, since the route flows jointly satisfy the fixed OD demands. For this supplementary sensitivity comparison, the relevant quantity is proximity to the reference equilibrium. In particular, the selected reference-based states satisfy

$${\displaystyle \frac{\parallel p_{\mathrm{bud}}^{\mathrm{DGR}}-p^{\mathrm{ref}}\parallel }{\parallel p^{\mathrm{ref}}\parallel }}<{\displaystyle \frac{\parallel p_{\mathrm{bud}}^{\left(m\right)}-p^{\mathrm{ref}}\parallel }{\parallel p^{\mathrm{ref}}\parallel }}$$

for each competing method m, together with the analogous route-cost comparison.

The corresponding link-flow vectors are $v_{\mathrm{bud}}^{\left(m\right)}=\Delta p_{\mathrm{bud}}^{\left(m\right)},$ and their link volume-to-capacity ratios are $V{C}_a^{\left(m\right)}={\displaystyle \frac{v_{a,\mathrm{bud}}^{\left(m\right)}}{c_a^{\mathrm{eff}}}}.$

Figure 19 displays the corresponding network-loading patterns.

The threshold ${\displaystyle \frac{v_a}{c_a^{\mathrm{eff}}}}=1$ corresponds to effective link capacity; values exceeding one therefore indicate overloaded links. As in the route-flow comparison, the magnitude of an individual link flow is not itself an algorithmic accuracy measure. Rather, proximity to the reference loading and the resulting level of congestion provide the appropriate interpretation.

The aggregate traffic-state indicators are reported in Figure 20.

For the computed traffic states, a smaller TSTT represents a smaller aggregate network travel time, whereas a smaller MaxVC indicates a less severe peak congestion level. These are descriptive network indicators rather than variational-inequality residuals, and TSTT is not a system-optimality certificate for a Wardrop equilibrium. The DGR–Tseng traffic state gives the smallest values of both quantities in each of the four cases.

Remark 7. 

The traffic experiment provides complementary algorithmic and transportation-level information. Figure 13 and Figure 14 describe the network structure and the eight admissible routes. Figure 15 illustrates the natural-residual and Wardrop-gap trajectories, while Figure 16 and Figure 17 compare the computational effort and terminal accuracy across the four initial path-flow cases.

At the transportation level, Figure 18 compares route flows and generalized route costs with the high-accuracy reference equilibrium, whereas Figure 19 shows the induced link loading and congestion. Figure 20 summarizes the corresponding total system travel time and maximum volume-to-capacity ratio.

The numerical results show that DGR–Tseng exhibits the most favorable iteration count, CPU time, natural residual, and Wardrop relative gap in all four cases. The supplementary reference-based route and link states associated with DGR–Tseng are closest to the high-accuracy reference equilibrium, and the corresponding traffic states have the smallest TSTT and MaxVC values.

The experiment therefore illustrates the potential usefulness of the double golden-ratio mechanism for a structured multi-OD traffic equilibrium problem involving nonlinear congestion, asymmetric route interaction, and incident-sensitive capacities.

10. Conclusions

In this paper, we introduced an adaptive double golden-ratio Tseng-type extragradient method, called DGR–Tseng, for solving variational inequality problems in real Hilbert spaces. The method employs two golden-ratio averaging steps while requiring only one metric projection per iteration. Its adaptive stepsize also removes the need for prior knowledge of the Lipschitz constant.

Under the solution-oriented assumptions adopted in this work, we established weak convergence of the generated sequence to a point of $\Omega$. Under an additional Robinson/Luo–Tseng-type local error-bound condition, we proved R-linear convergence of $dist(p_k,\Omega )$ and a Q-linear contraction for the associated weighted paired distance. These linear convergence results do not require strong monotonicity, strong pseudomonotonicity, or a singleton solution set.

The theoretical results were supported by numerical experiments in ${\ell }^2$ and $L^2(0,1)$, together with applications to sparse signal reconstruction and multi-OD urban traffic network equilibrium. The reported results show that DGR–Tseng compares favorably with GRPA, GRT–Tseng, and SEM–GRT in terms of iteration counts, CPU time, and convergence accuracy. It also produced favorable reconstruction and traffic-equilibrium performance measures in the respective applications.

These results demonstrate the effectiveness of the double golden-ratio mechanism for projection-efficient variational inequality methods beyond the standard monotonicity setting. Future work may consider extensions to non-Lipschitz operators, split and bilevel variational inequalities, and related equilibrium and inclusion problems.