**affiliationtext: Department of Mathematics, University of Tübingen,
Germany****affiliationtext: Normandie Université, ENSICAEN, UNICAEN, CNRS,
GREYC, France.

# Inertial Quasi-Newton Methods for Monotone Inclusion: Efficient Resolvent
Calculus and Primal-Dual Methods

Shida Wang  Jalal Fadili  Peter Ochs

###### Abstract

We introduce an inertial quasi-Newton Forward-Backward Splitting Algorithm to
solve a class of monotone inclusion problems. While the inertial step is
computationally cheap, in general, the bottleneck is the evaluation of the
resolvent operator. A change of the metric makes its computation hard even for
(otherwise in the standard metric) simple operators. In order to fully exploit
the advantage of adapting the metric, we develop a new efficient resolvent
calculus for a low-rank perturbed standard metric, which accounts exactly for
quasi-Newton metrics. Moreover, we prove the convergence of our algorithms,
including linear convergence rates in case one of the two considered operators
is strongly monotone. Beyond the general monotone inclusion setup, we
instantiate a novel inertial quasi-Newton Primal-Dual Hybrid Gradient Method
for solving saddle point problems. The favourable performance of our inertial
quasi-Newton PDHG method is demonstrated on several numerical experiments in
image processing.

##  1 Introduction

Nowadays, convex optimization appears in many modern disciplines, especially
when dealing with datasets of large scale. There is a strong need of efficient
optimization schemes. Unfortunately, the high dimension of the problems at
hand makes the use of second order methods intractable. A promising
alternative are quasi-Newton type methods, which aim for exploiting cheap and
accurate first order approximations of the second order information. In
particular, the so-called limited memory quasi-Newton method has proved very
successful for solving unconstrained large scale problems. However, many
practical problems in machine learning, image processing or statistics more
naturally have constraints or are non-smooth by construction.

A problem structure that can cover a broad class of non-smooth problems in
these applications is the monotone inclusion problem on a real Hilbert space
$\mathcal{H}$ as the following:

$\mathrm{find}\,\,x\in\mathcal{H}\ \ \mathrm{such\ that}\ \ (A+B)x\ni 0\,,$
(1)

where
$A\colon\mathcal{H}\mathrel{\vbox{\halign{#\cr\nointerlineskip\cr$\rightarrow$\crcr\nointerlineskip\cr$\rightarrow$\crcr}}}{\mathcal{H}}$
is a maximally monotone operator and $B\colon\mathcal{H}\to\mathcal{H}$ is a
single-valued $\beta$ cocoercive operator. As a special case (1) comprises the
setting of minimization problems of the form $\min_{x\in\mathcal{H}}f(x)+g(x)$
with a proper lower semi-continuous convex function $f$ and convex function
$g$ with Lipschitz continuous gradient by setting $A=\partial f$ and $B=\nabla
g$.

A fundamental algorithmic scheme to tackle the problem class (1) is Forward-
Backward Splitting (FBS). However, this algorithm can be quite slow for ill-
conditioned problems, for which one would like to exploit the second order
information to adapt to the local geometry of the problem. As a
computationally affordable approximation in this paper, we propose a quasi-
Newton variant that takes advantage from a variable metric that is computed
using first order information only. We manage to remedy the main computational
bottlenecks for these type of approach by developing an efficient low-rank
variable metric resolvent calculus.

Our approach is inspired by the proximal quasi-Newton method in [7]. We extend
the framework proposed in [7] to the resolvent setting with a $M$ $\pm$
rank-$r$ symmetric positive definite variable metric. We first study the
convergence of two variants of FBS algorithm with respect to this type of
metrics. One variant uses inertial step to accelerate its convergence and the
other uses relaxation to guarantee convergence with mild assumption. In
analogy to the proximal calculus proposed in [7], we develop a resolvent
calculus that allows for a efficient evaluation of resolvent operators with
respect to this type of metrics by splitting the evaluation into two
computationally simple steps: calculating a resolvent operator with respect to
a simple metric $M$ and solving a low dimensional root finding problem. This
allows for using popular quasi-Newton methods, such as the limited memory SR1
or BFGS method.

In order to exploit the power and to illustrate the variety of problems that
can be solved via the framework as (1), the developed algorithms are
instantiated for the following saddle point problem, which include quasi-
Newton Primal-Dual methods and enjoy an enormous number of potential
applications:

$\min_{x\in\mathcal{H}_{1}}\max_{y\in\mathcal{H}_{2}}\left\langle
Kx,y\right\rangle+g(x)+G(x)-f(y)-F(y)\,,$ (2)

where $g$ and $f$ are lower semi-continuous convex functions, $G$ and $F$ are
convex, differentiable with Lipschitz gradients and $K$ is a bounded linear
operator. The numerical performance of our algorithms is tested on several
numerical experiments and demonstrates a clear improvement when our quasi-
Newton methods are used.

###  1.1 Related Works

Smooth quasi-Newton. Quasi-Newton methods are widely studied and used for
optimization with objective functions sufficiently smooth [35]. Their
motivation is to build a cheap approximation to the Newton’s Method. If the
approximation of the second order information is given by a positive definite
matrix, quasi-Newton methods can be interpreted as adapting the metric of the
space locally to the objective function. The success of these methods requires
a good approximation of the second order information (Hessian) of the
objective by using the first order information, which is still actively
researched. Recently, [30] proposed to select greedily some basis vectors
instead of using the difference of successive iterates for updating the
Hessian approximations. Based on their works, [24] constructed an
approximation of the indefinite Hessian of a twice differentiable function.
All methods mentioned above require sufficient smoothness of the objective
functions.
Non-smooth quasi-Newton. A broad class of optimization problems is interpreted
as a composition of a smooth function $f$ and a non-smooth function $h$. To
deal with the non-smoothness of $h$ efficiently, many authors consider a
combination of FBS with the quasi-Newton methods. By using the forward-
backward envelope, [27, 33] reinterpreted the FBS algorithm as a variable
metric gradient method for a smooth optimization problem in order to apply the
classic Newton or quasi-Newton method. For non-smooth function $g$ as simple
as an indicator function of a convex set, [32, 31] proposed an elegant method
named projected quasi-Newton algorithm (PQN) which, however, requires either
solving a subproblem or a diagonal metric. [23] extended PQN to a more general
setting as long as the proximal operator with respect to $h$ is simple to
compute. For a class of low-rank perturbed metrics, [6, 7] developed a
proximal quasi-Newton method with a root finding problem as the subproblem
which can be solved easily and efficiently. This method can be extended to the
nonconvex setting [19]. Based on [6, 7], [20] incorporated a limited-memory
quasi-Newton update. Meanwhile, [21] developed a different algorithm to
evaluate the proximal operator of the separable $l_{1}$ norm with respect to a
low-rank metric $V=M-UU^{*}$. We extend the quasi-Newton approach of [6, 7]
from the nonsmooth convex minimization setting to monotone inclusion problems
of type (1). More related works for the minimization setting are intensively
discussed in [7]. Therefore, we focus here on approaches to solve the monotone
inclusion problem. It is generic to use a variable metric for solving a
monotone inclusion problem, which is discussed in [12, 13]. Its convergence
relies on quasi-Fejér monotonicity [12]. However the efficient calculation of
the resolvent is left as an open problem. Our approaches use a variable metric
to obtain a quasi-Newton method with an efficient resolvent calculus. From a
different perspective, [2] considered incorporating Newton methods by
introducing dynamic system with Hessian damping. Later, in [1, 3], time
discretization of some dynamics gives new insight into the Newton’s method for
solving monotone inclusions which is different from variable metric
approaches. Recently, there is a new approach that [22] proposed a generalized
damped Newton type which is based on second-order subdifferentials.
PDHG. Primal-Dual Hybrid Gradient (PDHG) is widely used for solving saddle
point problems of the form (2). PDHG can be interpreted as a proximal point
algorithm [17] applied on a monotone inclusion problem with a fixed metric.
Based on this idea, [26] proposed an inertial FBS method applied to the sum of
set-valued operators from which a generalization of PDHG method is derived.
Also based on similar ideas, [28] considers diagonal preconditioning to
accelerate PDHG. Their method can be regarded as using a fixed blocked matrix
as a metric. Later, [25] considers non-diagonal preconditioning and pointed
out if a special preconditioner is chosen, that kind of preconditioned PDHG
method will be a special form of the linearized ADMM. Their method requires an
inner loop due to the non-diagonal preconditioning. [16] introduces an
adaptive PDHG scheme which can also be understood as using a variable metric
with step size tuned automatically. However, [16] considered changes on the
diagonal of the metric. Our resolvent calculus in Section 4 provides another
possibility to change metrics at elements off the diagonal to deal with
resolvent operators.

##  2 Preliminaries

Let us recall some essential notations and definitions. Let $\mathcal{H}$ be a
Hilbert space equipped with inner product
$\left\langle\cdot,\cdot\right\rangle$ and induced norm
$\|\cdot\|=\sqrt{\left\langle\cdot,\cdot\right\rangle}$. An operator
$K\in\mathcal{B}(\mathcal{D},\mathcal{H})$ is a linear bounded mapping from a
Hilbert space $\mathcal{D}$ to $\mathcal{H}$. We abbreviate
$\mathcal{B}(\mathcal{H},\mathcal{H})$ to $\mathcal{B}(\mathcal{H})$. The
adjoint of $M\in\mathcal{B}(\mathcal{H})$ is denoted by $M^{*}$. We define
$\mathcal{S}(\mathcal{H})\coloneqq\\{M\in\mathcal{B}(\mathcal{H})|M=M^{*}\\}$
and the identity operator by $I\in\mathcal{S}(\mathcal{H})$. Without
ambiguity, we also use the notation $\|M\|$ for the operator norm of
$M\in\mathcal{S}(\mathcal{H})$ with respect to $\|\cdot\|$. The partial
ordering on $\mathcal{S}(\mathcal{H})$ is given by

$(\forall U\in\mathcal{S}(\mathcal{H}))(\forall
V\in\mathcal{S}(\mathcal{H}))\colon\quad U\succeq V\iff(\forall
x\in\mathcal{H})\left\langle Ux,x\right\rangle\geq\left\langle
Vx,x\right\rangle\,.$ (3)

For $\sigma\in[0,+\infty)$, we introduce
$\mathcal{S}_{\sigma}(\mathcal{H})\coloneqq\\{U\in\mathcal{S}(\mathcal{H})|U\succeq\sigma
I)\\}$. Similarly, we introduce
$\mathcal{S}_{++}(\mathcal{H})\coloneqq\\{U\in\mathcal{S}(\mathcal{H})|U\succ
0\\}$. The norm $\|\cdot\|_{M}$ is defined by $\sqrt{\left\langle
M\cdot,\cdot\right\rangle}$ for some $M\in\mathcal{S}_{++}(\mathcal{H})$.

A set valued operator
$A\colon\mathcal{H}\mathrel{\vbox{\halign{#\cr\nointerlineskip\cr$\rightarrow$\crcr\nointerlineskip\cr$\rightarrow$\crcr}}}{\mathcal{H}}$
is defined by its graph

$\mathrm{Graph}\,A\coloneqq\\{(x,y)\in\mathcal{H}|x\in\mathrm{Dom}(A),y\in
Ax\\}\,,$

and has a domain given by

$\mathrm{Dom}(A)\coloneqq\\{x\in\mathcal{H}|Ax\neq\emptyset\\}.$

Given two set-valued operators
$A,B\colon\mathcal{H}\mathrel{\vbox{\halign{#\cr\nointerlineskip\cr$\rightarrow$\crcr\nointerlineskip\cr$\rightarrow$\crcr}}}{\mathcal{H}}$,
we define
$A+B\colon\mathcal{H}\mathrel{\vbox{\halign{#\cr\nointerlineskip\cr$\rightarrow$\crcr\nointerlineskip\cr$\rightarrow$\crcr}}}{\mathcal{H}}$
as follows:

$\begin{split}\mathrm{Dom}(A+B)&=\mathrm{Dom}(A)\cap\mathrm{Dom}(B)\,,\\\
(A+B)x&=Ax+Bx\coloneqq\\{y\in\mathcal{H}|\exists y_{1}\in\mathcal{H},\exists
y_{2}\in\mathcal{H}\,\ \mathrm{such\ that}\,\ y=y_{1}+y_{2}\\}\,.\\\
\end{split}$

The inverse of $A$ is denoted by $A^{-1}$ given by
$A^{-1}(y)\coloneqq\\{x\in\mathcal{H}|y\in Ax\\}$ and the zero set of $A$ is
denoted by $zer(A+B)\coloneqq\\{x\in\mathcal{H}|(A+B)x\ni 0\\}$. We call $A$
is $\gamma_{A}$-strongly monotone with respect to norm $\|\cdot\|$ and
$\gamma_{A}\geq 0$ if $\left\langle
x-y,u-v\right\rangle\geq\gamma_{A}\|x-y\|^{2}$ for any pair
$(x,u)\,,(y,v)\in\mathrm{Graph}\,A$ and $\gamma_{A}=0$ if just monotone. The
resolvent of
$A\colon\mathcal{H}\mathrel{\vbox{\halign{#\cr\nointerlineskip\cr$\rightarrow$\crcr\nointerlineskip\cr$\rightarrow$\crcr}}}{\mathcal{H}}$
with respect to metric $M\in\mathcal{S}_{++}(\mathcal{H})$ is defined as

$J^{M}_{A}\coloneqq(I+M^{-1}A)^{-1}\quad\text{and we set}\ J_{A}\coloneqq
J_{A}^{I}\ \text{for the identity metric}\,I\,,$ (4)

which, as shown for example in [5], enjoys the following properties.

######  Proposition 2.1.

Let
$A\colon\mathcal{H}\mathrel{\vbox{\halign{#\cr\nointerlineskip\cr$\rightarrow$\crcr\nointerlineskip\cr$\rightarrow$\crcr}}}{\mathcal{H}}$
be maximally monotone, $M\in\mathcal{S}_{++}(\mathcal{H})$ and
$y\in\mathcal{H}$. Then, the following holds

$y=J^{M}_{\gamma A}(x)\iff x\in y+\gamma M^{-1}Ay\iff x-y\in\gamma
M^{-1}Ay\iff(y,\gamma^{-1}M(x-y))\in\mathrm{Graph}\,A\,.$ (5)

######  Proposition 2.2.

Let
$A:\mathcal{H}\mathrel{\vbox{\halign{#\cr\nointerlineskip\cr$\rightarrow$\crcr\nointerlineskip\cr$\rightarrow$\crcr}}}{\mathcal{H}}$
be maximally monotone. Then for every sequence
$(x_{k},u_{k})_{k\in\mathbb{N}}$ in $\mathrm{Graph}\,A$ and every
$(x,u)\in\mathcal{H}\times\mathcal{H}$, if $x_{k}\rightharpoonup x$ and
$u_{k}\to u$, then $(x,u)\in\mathrm{Graph}\,A$.

######  Lemma 2.3.

Let
$T\colon\mathcal{H}\mathrel{\vbox{\halign{#\cr\nointerlineskip\cr$\rightarrow$\crcr\nointerlineskip\cr$\rightarrow$\crcr}}}{\mathcal{H}}$
be a maximally monotone operator and let $M\in\mathcal{S}_{++}(\mathcal{H})$.
Then, for any $z\in\mathcal{H}$, we have $J^{M}_{T}(z)=M^{-1/2}\circ
J_{M^{-1/2}TM^{-1/2}}\circ M^{1/2}(z)$

###### Proof.

See Appendix A.1. ∎

######  Lemma 2.4.

If $A$ is $\gamma_{A}$-strongly monotone and $\gamma_{A}\geq 0$, then
$J_{A}^{M}$ is Lipschitz continuous with respect to $\|\cdot\|_{M}$ with
constant ${1}/(1+\tfrac{\gamma_{A}}{C})\in(0,1]$ for any $C$ satisfying
$\|M\|\leq C<\infty$.

###### Proof.

See Appendix A.2. ∎

######  Proposition 2.5 (Variable Metric quasi-Fejér monotone sequence [12]).

Let $\sigma\in(0,+\infty)$, let $\varphi\colon[0,+\infty)\to[0,+\infty)$ be
strictly increasing and such that $\lim_{t\to+\infty}\varphi(t)=+\infty$, let
$(M_{k})_{k\in\mathbb{N}}$ be in $\mathcal{S}_{\sigma}(\mathcal{H})$, let $C$
be a nonempty subset of $\mathcal{H}$, and let $(x_{k})_{k\in\mathbb{N}}$ be a
sequence in $\mathcal{H}$ such that

$\begin{split}(\exists(\eta_{k})_{k\in\mathbb{N}}\in\ell_{+}^{1}(\mathbb{N}))(\forall
z\in
C)(\exists(\epsilon_{k})_{k\in\mathbb{N}}\in\ell_{+}^{1}(\mathbb{N}))&(\forall
k\in\mathbb{N})\colon\\\
&\varphi(\|x_{k+1}-z\|_{M_{k+1}})\leq(1+\eta_{k})\varphi(\|x_{k}-z\|_{M_{k}})+\epsilon_{k}\,.\end{split}$
(6)

  1. (a)

Then $(x_{k})_{k\in\mathbb{N}}$ is bounded and, for every $z\in C$,
$(\|z_{k}-z\|_{M_{k}})_{k\in\mathbb{N}}$ converges.

  2. (b)

If additionally, there exists $M\in\mathcal{S}_{\sigma}(\mathcal{H})$ such
that $M_{k}\to M$ pointwise, as is the case when

$\sup_{k\in\mathbb{N}}\|M_{k}\|<+\infty\quad\mathrm{and}\quad(\exists(\eta_{k})_{k\in\mathbb{N}}\in\ell_{+}^{1}(\mathbb{N}))(\forall
k\in N)\colon\,(1+\eta_{k})M_{k}\succeq M_{k+1}\,,$ (7)

then $(x_{k})_{k\in\mathbb{N}}$ converges weakly to a point in $C$ if and only
if every weak sequential cluster point of $(x_{k})_{k\in\mathbb{N}}$ lies in
$C$.

A key result for our resolvent calculus in Section 4 is the following abstract
duality principle.

######  Lemma 2.6 (A duality result for operators [4]).

Let
$A\colon\mathcal{H}\mathrel{\vbox{\halign{#\cr\nointerlineskip\cr$\rightarrow$\crcr\nointerlineskip\cr$\rightarrow$\crcr}}}{\mathcal{H}}$
be an operator such that $A^{-1}$ is single-valued and let
$B\colon\mathcal{H}\to\mathcal{H}$ be a single-valued operator. Then, the
following holds for $x,u\in\mathcal{H}$:

$\begin{cases}0\in Ax+Bx\\\ 0\in
B^{-1}u-A^{-1}(-u)\end{cases}\iff\begin{cases}x\in B^{-1}u\\\ -u\in
Ax\end{cases}\iff\begin{cases}Bx=u\\\ x=A^{-1}(-u)\end{cases}\,.$ (8)

Moreover, if there exists $x\in\mathcal{H}$ such that $0\in Ax+Bx$ or there
exists $u\in\mathcal{H}$ such that $0\in B^{-1}u-A^{-1}(-u)$, then there
exists a unique primal-dual pair $(x,u)$ that satisfies the equivalent
conditions above.

We also need the following lemma which was stated as [29, Lemma 2.2.2].

######  Lemma 2.7.

Let $C_{k}\geq 0$ and let

$\begin{split}&C_{k+1}\leq(1+\nu_{k})C_{k}+\zeta_{k},\quad\nu_{k}\geq 0,\
\zeta_{k}\geq 0,\\\
&\sum_{k\in\mathbb{N}}\nu_{k}<\infty,\quad\sum_{k\in\mathbb{N}}\zeta_{k}<\infty.\\\
\end{split}$ (9)

Then, $C_{k}\to C\geq 0$ and $C<+\infty$.

##  3 Inertial Quasi-Newton Forward-Backward Splitting for Monotone Inclusion

###  3.1 Problem setting

Let
$A\colon\mathcal{H}\mathrel{\vbox{\halign{#\cr\nointerlineskip\cr$\rightarrow$\crcr\nointerlineskip\cr$\rightarrow$\crcr}}}{\mathcal{H}}$
be a maximally and $\gamma_{A}$-strongly monotone operator with
$\gamma_{A}\geq 0$ with respect to norm $\|\cdot\|$ and let
$B\colon\mathcal{H}\to\mathcal{H}$ be a single-valued, $\beta$ cocoercive and
$\gamma_{B}$-strongly monotone operator with $\gamma_{B}\geq 0$ with respect
to norm $\|\cdot\|$. We consider the following monotone inclusion problem:

$\mathrm{find}\,\,x\in\mathcal{H}\ \ \mathrm{such\ that}\ \ Ax+Bx\ni 0\,.$
(10)

Note that by setting $\gamma_{A}=0$ or $\gamma_{B}=0$, we include the general
case of monotone operators that are not necessarily strongly monotone. In this
paper, we propose two variants of an efficiently implementable quasi-Newton
Forward-Backward Splitting (FBS) algorithm in Algorithm 1 and Algorithm 2 to
solve (10). The main update step is a variable metric FBS step in both
algorithms as shown in (12) and (13), i.e., a forward step with respect to the
cocoercive operator $B$, followed by a proximal point step (computation of the
resolvent) with respect to the maximally monotone operator $A$, both evaluated
in an iteration dependent metric $M_{k}$ that is inspired by quasi-Newton
methods and therefore adapts to the local geometry of the problem. In contrast
to the related works, as discussed in Section 1.1, we emphasize the importance
of efficiently implementable resolvent operators (see Section 4). The two
variants allow for more or less flexibility for the choice of the metric (see
Section 5). Both variants account for potential (numerical) errors in the
evaluation of the forward-backward step. Algorithm 1 combines a FBS step with
an additional inertial step which has the potential of accelerating the
convergence, as we illuminate in our numerical experiments in Section 7. It is
a generalization of the algorithms in [10, 26] to a quasi-Newton variant. In
[26], Lorenz and Pock proposed an inertial Forward-Backward Splitting
algorithm but with a fixed metric which is different from our Algorithm 1
which uses a variable metric. Algorithm 2 combines FBS with a relaxation step
in (iii) which yields convergence under mild assumptions on the metric. It
generalizes the correction step introduced in [17] which can be retrieved by
metric $M_{k}\equiv M$ and $B\equiv 0$.

Algorithm 1 Inertial quasi-Newton Forward-Backward Splitting

$N\geq 0$, $(\|\epsilon_{k}\|)_{k\in\mathbb{N}}\in\ell_{+}^{1}(\mathbb{N})$

$x_{0}\in\mathcal{H}$

Update for $k=0,\cdots,N$:

  1. (i)

Compute $M_{k}$ according to a quasi-Newton framework.

  2. (ii)

Compute the inertial step (extrapolation step):

$\bar{z}_{k}\leftarrow{z}_{k}+\alpha_{k}({z}_{k}-z_{k-1})\,,$ (11)

  3. (iii)

and the forward-backward step:

$z_{k+1}\leftarrow
J_{A}^{M_{k}}(\bar{z}_{k}-M_{k}^{-1}B\bar{z}_{k})+\epsilon_{k}\,.$ (12)

End

Algorithm 2 Quasi-Newton Forward-Backward Splitting with Relaxation

$N\geq 0$, $(\|\epsilon_{k}\|)_{k\in\mathbb{N}}\in\ell_{+}^{1}(\mathbb{N})$

$x_{0}\in\mathcal{H}$

Update for $k=0,\cdots,N$

  1. (i)

Compute $M_{k}$ according to a quasi-Newton framework.

  2. (ii)

Compute the forward-backward step:

$\tilde{z}_{k}\leftarrow
J_{A}^{M_{k}}(z_{k}-M_{k}^{-1}Bz_{k})+\epsilon_{k}\,,$ (13)

  3. (iii)

and relaxation step:

$t_{k}=\frac{\left\langle
z_{k}-\tilde{z}_{k},(M_{k}-B)(z_{k}-\tilde{z}_{k})\right\rangle}{2\|(M_{k}-B)(z_{k}-\tilde{z}_{k})\|^{2}}\,,$
(14) $z_{k+1}\leftarrow z_{k}-t_{k}[(M_{k}-B)(z_{k}-\tilde{z}_{k})]\,.$ (15)

End

###  3.2 Convergence

In this subsection, we prove the convergence of Algorithm 1 and Algorithm 2.
The implementation details for the specific quasi-Newton features are deferred
to the upcoming section.

####  3.2.1 Algorithm 1: Inertial Quasi-Newton Forward-Backward Splitting

The following convergence result is a generalization of [11, Theorem 3.1] to
an inertial version of variable metric Forward-Backward Splitting.

######  Assumption 1.

Let $\sigma\in(0,+\infty)$. $(M_{k})_{k\in\mathbb{N}}$ is a sequence in
$\mathcal{S}_{\sigma}(\mathcal{H})$ such that

$\begin{cases}C\coloneqq\mathrm{sup}_{k\in\mathbb{N}}||M_{k}||<\infty\,,\\\
(\exists(\eta_{k})_{k\in\mathbb{N}}\in\ell_{+}^{1}(\mathbb{N}))(\forall
k\in\mathbb{N})\colon\quad(1+\eta_{k})M_{k}\succeq M_{k+1}\,,\\\ \end{cases}$
(16)

and $M_{k}-\tfrac{1}{2\beta}I\in\mathcal{S}_{\epsilon}(\mathcal{H})$ for all
$k\in\mathbb{N}$ and some $\epsilon>0$. Moreover,
$(\epsilon_{k})_{k\in\mathbb{N}}$ is a sequence in $\mathcal{H}$ such that
$\sum_{k\in\mathbb{N}}\|\epsilon_{k}\|<+\infty$.

######  Theorem 3.1.

Consider Problem (10) and let the sequence $(z_{k})_{k\in\mathbb{N}}$ be
generated by Algorithm 1 where Assumption 1 holds and
$(\alpha_{k})_{k\in\mathbb{N}}$ are selected such that
$\alpha_{k}\in(0,\Lambda]$ with $\Lambda<\infty$ and

$\sum_{k\in\mathbb{N}}\alpha_{k}\max\\{\|z_{k}-z_{k-1}\|_{M_{k}},\|z_{k}-z_{k-1}\|_{M_{k}}^{2}\\}<+\infty\,.$

Suppose that $\mathrm{zer}(A+B)\neq\emptyset$. Then $(z_{k})_{k\in\mathbb{N}}$
are bounded and weakly converge to a point $z^{*}\in zer(A+B)$, i.e.
$z_{k}\rightharpoonup z^{*}$ as $k\to\infty$.

Furthermore, if additionally we assume $\epsilon_{k}\equiv 0$ for any
$k\in\mathbb{N}$, $\gamma_{A}>0$ or $\gamma_{B}>0$ and
$M_{k}-\tfrac{1}{\beta}I\in\mathcal{S}_{\epsilon}(\mathcal{H})$, then there
exist some $\xi\in(0,1)$, some $\Theta>0$ and some $K_{0}\in\mathbb{N}$ such
that for any $k\geq K_{0}$,

$\|z_{k}-z^{*}\|_{M_{k}}^{2}\leq(1-\xi)^{k-K_{0}}\|z_{K_{0}}-z^{*}\|_{M_{K_{0}}}^{2}+\sum_{i=K_{0}}^{k-1}\Theta(1-\xi)^{k-i}\alpha_{i}\max\\{\|z_{i}-z_{i-1}\|_{M_{i}},\|z_{i}-z_{i-1}\|_{M_{i}}^{2}\\}\,.$
(17)

###### Proof.

See Appendix A.3. ∎

######  Remark 1.

  * •

We obtain that the second term on the right hand of (17) converges to 0 by
using [29, Lemma 2.2.3].

  * •

The linear convergence factor $1-\xi$ is chosen such that there exists
$K_{0}\in\mathbb{N}$ with

$(1+\eta_{k})\Big{(}\frac{1-\frac{\gamma_{B}}{C}}{1+\frac{\gamma_{A}}{C}}\Big{)}\leq
1-\xi\,,$

for all $k\geq K_{0}$ holds.

  * •

The convergence rate for the strongly monotone setting can be influenced a lot
by the decay rate of $\alpha_{k}$.

    1. (i)

If $\alpha_{k}=O(q^{k})$ for $q=1-\xi$, we have convergence rate of
$O(kq^{k})$ for $k>K_{0}$ where $K_{0}$ is sufficiently large.

    2. (ii)

If $\alpha_{k}=O(\tfrac{1}{k^{2}})$, we have convergence rate of
$O(\tfrac{1}{k})$ for $k>K_{0}$ where $K_{0}$ is sufficiently large.

####  3.2.2 Algorithm 2: Quasi-Newton Forward-Backward Splitting with
Relaxation

At the cost of a relaxation instead of an inertial step, we can significantly
generalize the flexibility of choosing the metric. Note that without loss of
generality, we assume $\tilde{z}_{k}\neq z_{k}$ for all $k\in\mathbb{N}$,
since otherwise $\tilde{z}_{k}$ is already solves the inclusion problem after
a finite number of iterations. This method is inspired by [17].

######  Assumption 2.

Let $\sigma\in(0,+\infty)$. $(M_{k})_{k\in\mathbb{N}}$ is a sequence in
$\mathcal{S}_{\sigma}(\mathcal{H})$ and $(\epsilon_{k})_{k\in\mathbb{N}}$ is a
sequence in $\mathcal{H}$ such that:

  1. (i)

For all $k\in\mathbb{N}$, we have
$(M_{k}-\tfrac{1}{\beta}I)\in\mathcal{S}_{c}(\mathcal{H})$, for some $c>0$,

  2. (ii)

$\mathrm{sup}_{k\in\mathbb{N}}||M_{k}||<\infty$,

  3. (iii)

$\sum_{k\in\mathbb{N}}\|\epsilon_{k}\|<+\infty$.

######  Theorem 3.2.

Consider Problem (10) and let the sequence $(z_{k})_{k\in\mathbb{N}}$ be
generated by Algorithm 2 where Assumption 2 holds. Then
$(\|z_{k}-z^{*}\|)_{k\in\mathbb{N}}$ is bounded for any $z^{*}\in zer(A+B)$
and $(z_{k})_{k\in\mathbb{N}}$ weakly converges to some $z^{*}\in zer(A+B)$,
i.e. $z_{k}\rightharpoonup z^{*}$ as $k\to\infty$.

Moreover, if $\epsilon_{k}\equiv 0$, $\|z_{k}-z^{*}\|$ decreases for any
$z^{*}\in zer(A+B)$ as $k\to\infty$. Furthermore, if $\gamma_{A}>0$ or
$\gamma_{B}>0$ and
$M_{k}-\tfrac{1}{\beta}I\in\mathcal{S}_{\epsilon}(\mathcal{H})$ for all
$k\in\mathbb{N}$, then $z_{k}$ converges linearly: there exist some
$\xi\in(0,1)$ such that

$\|z_{k}-z^{*}\|^{2}\leq(1-\xi)^{k}\|z_{0}-z^{*}\|^{2}\,.$ (18)

###### Proof.

See Appendix A.4. ∎

######  Remark 2.

The linear convergence factor is given by

$\xi=\tfrac{1}{2}\min\\{2(\gamma_{A}+\gamma_{B})\delta,c\delta\\}\,,$

where $\delta=\tfrac{c}{2(C+\tfrac{1}{\beta})^{2}}$.

######  Remark 3.

It is worth mentioning that Algorithm 2 can be interpreted as a closed loop
system which uses the previous iterates to update the relaxation parameter
$\alpha_{k}$.

##  4 Resolvent Calculus for Low-Rank Perturbed Metric

In this section, we extend the proximal calculus of [7] to the setting of
resolvent operators $J^{V}_{T}$ with a metric $V=M\pm
Q\in\mathcal{S}_{++}(\mathcal{H})$, where $M\in\mathcal{S}_{++}(\mathcal{H})$
and $Q\in\mathcal{S}_{0}(\mathcal{H})$. This is a key result for making our
quasi-Newton methods efficiently applicable for solving monotone inclusion
problemss. Computing the resolvent operator $J^{V}_{T}(z)$ can be reduced to
two simple problems consisting of evaluating $J^{M}_{T}$ at a shifted point
$z-M^{-1}v^{*}\in\mathcal{H}$ for some $v^{*}\in\mathcal{H}$ that is derived
from an $r$-dimensional root finding problem which can be solved by a semi-
smooth Newton method (Algorithm 3) or bisection method (Algorithm 4). In
conclusion, if $J^{M}_{T}$ can be computed efficiently, the same is true for
$J_{T}^{V}$. The result crucially relies on the abstract duality principle of
Attouch-Théra [4] (see Lemma 2.6).

###  4.1 General case

We first present a general result.

######  Theorem 4.1.

Let
$T\colon\mathcal{H}\mathrel{\vbox{\halign{#\cr\nointerlineskip\cr$\rightarrow$\crcr\nointerlineskip\cr$\rightarrow$\crcr}}}{\mathcal{H}}$
be a maximally monotone operator, $V\coloneqq M\pm Q$ be a symmetric positive
definite metric with bounded linear operators
$M,Q\colon\mathcal{H}\to\mathcal{H}$ such that $M$ is positive definite and Q
is positive semi-definite. Then, the resolvent operator $J^{V}_{T}$ can be
computed as follows:

$x^{*}=J^{V}_{T}(z)\quad\iff\quad\begin{cases}x^{*}=J^{M}_{T}(z-M^{-1}v^{*})\,\mathrm{and}\\\
v^{*}\,\mathrm{solves}\,0\in\pm Q^{+}v^{*}+z-J^{M}_{T}(z-M^{-1}v^{*})\\\
\mathrm{s.t.}\,v^{*}\in\mathrm{im}(\pm Q)\,.\end{cases}$ (19)

The solution $v^{*}=U\alpha^{*}$ where $\alpha^{*}$ is the unique root of

$l(\alpha)\coloneqq U^{*}Q^{+}U\alpha+U^{*}(z-J^{M}_{T}(z\mp
M^{-1}U\alpha))\,,$ (20)

where $Q^{+}$ is the pseudo-inverse of $Q$ and
$U\colon\mathbb{R}^{r}\to\mathrm{im}(Q),\,\alpha\mapsto
U\alpha\coloneqq\sum^{r}_{i=1}\alpha_{i}u_{i}$ is an isomorphism defined by
$r$ linearly independent $u_{1},...,u_{r}\in\mathrm{im}(Q)$. The function $l$
is Lipschitz continuous with constant $||U^{*}Q^{+}U||+||M^{-1/2}U||^{2}$ and
strictly monotone.

Now, we give the proof of Theorem 4.1. For convenience, we define translation
operator $\tau_{p}\colon\mathcal{H}\to\mathcal{H}$ by $\tau_{p}(x)=x-p$ with
inverse $\tau_{p}^{-1}=\tau_{-p}$.

###### Proof.

See Appendix A.5. ∎

######  Corollary 4.1.1.

Consider the assumptions of Theorem 4.1 with $Q$ given by $Q=UU^{*}$ for some
$U\in\mathcal{B}(\mathbb{R}^{r},\mathcal{H})$. In this case,

$x^{*}=J^{V}_{T}(z)\quad\iff\quad\begin{cases}x^{*}=J^{M}_{T}(z-M^{-1}U\alpha^{*})\,\mathrm{and}\\\
\alpha^{*}\,\mathrm{solves}\,l(\alpha)\\\ \end{cases}$ (21)

where $\alpha^{*}$ is the unique root of

$l(\alpha)\coloneqq\alpha+\left\langle U,z-J^{M}_{T}(z\mp
M^{-1}U\alpha)\right\rangle\,.$ (22)

The function $l$ is Lipschitz continuous with constant $1+||M^{-1/2}U||^{2}$
and strictly monotone.

###### Proof.

See Appendix A.6. ∎

###  4.2 Solving the Root-Finding Problem

The efficiency of the reduction in Theorem 4.1 relies also on the solution of
a root finding problem which we discuss thoroughly in this subsection. We
present a semi-smooth Newton approach, for which we prove local convergence.
In order to narrow down the neighborhood of the sought root, we complement the
semi-smooth Newton strategy by a bisection method in Section 4.2.2.

####  4.2.1 Semi-smooth Newton Methods

In order to solve $l(\alpha)=0$ in (20) efficiently, we turn to a semi-smooth
Newton method. A locally Lipschitz function is called semi-smooth if its
Clarke Jacobian defines a Newton approximation scheme [15, Definition 7.4.2].
If $l(\alpha)$ is semi-smooth and any element of the Clarke Jacobian
$\partial^{c}l(\alpha^{*})$ is non-singular, then we can apply the inexact
semi-smooth Newton method detailed in Algorithm 3, analog to [7].

Algorithm 3 Semi-smooth Newton method to solve $l(\alpha)=0$

A point $\alpha_{0}\in\mathbb{R}^{n}$. $N$ is the maximal number of
iterations.

Update for $k=0,\cdots,N$:

if $l(\alpha_{k})=0$ then

stop

else

Select $G_{k}\in\partial^{c}l(\alpha_{k})$, compute $\alpha_{k+1}$ such that
$l(\alpha_{k})+G_{k}(\alpha_{k+1}-\alpha_{k})=e_{k}\,,$ and
$e_{k}\in\mathbb{R}^{r}$ is an error term satisfying
$\|e_{k}\|\leq\eta_{k}\|G_{k}\|$ and $\eta_{k}\geq 0$.

end if

End

Semi-smoothness may seem very restrictive. However we will observe there
exists a broad class of functions $l(\alpha)$ satisfying that $l(\alpha)$ is
semi-smooth. As shown in [8], a tame locally Lipschitz function is semi-
smooth. Therefore, it is sufficient to ensure $l(\alpha)$ is tame. We refer to
[8] for the definition of tameness. If the monotone operator $T$ in
$l(\alpha)$ is a tame map, then $l(\alpha)$ is also tame. In this case, the
convergence result for Algorithm 3 can be adapted from [7].

######  Proposition 4.2.

Let $l(\alpha)$ be defined as in Theorem 4.1, where $T$ is a set-valued tame
mapping. Then $l(\alpha)$ is semi-smooth and all elements of
$\partial^{c}l(\alpha^{*})$ are non-singular where $\alpha^{*}$ is the unique
root of $l(\alpha)$ from (20). In turn there exists $\bar{\eta}$ such that if
$\eta_{k}\leq\bar{\eta}$ for every $k$, there exists a neighborhood of
$\alpha^{*}$ such that for all $\alpha_{0}$ in that neighborhood, the sequence
generated by Algorithm 3 is well defined and converges to $\alpha^{*}$
linearly. If $\eta_{k}\to 0$, the convergence is superlinear.

###### Proof.

See Appendix A.7. ∎

######  Example 1.

If $f$ is a tame function and locally Lipschitz, then by [18, Proposition
3.1], $\partial f$ is a tame map.

######  Example 2.

The assumption that $T$ is a tame mapping is not restrictive. For example, in
PDHG setting we have a set-valued operator $T=\begin{pmatrix}\partial
g&K^{*}\\\ -K&\partial f\end{pmatrix}$ as defined by (33). If $g$ and $f$ are
both tame functions, then $\partial g$ and $\partial f$ are tame as well [18].
As a result, $T$ is a tame mapping.

####  4.2.2 Bisection

We are going to solve a root finding problem for $l(\alpha)=0$ via the
bisection method in the Algorithm 4 in the special case when $\alpha$ is one-
dimensional. We can still have similar bound on the range of $\alpha^{*}$ as
in [7]. The key tool is a bound on the values of $\alpha$ given by the
following proposition.

######  Proposition 4.3.

For $r=1$, the root $\alpha^{*}$ of $l(\alpha)=0$ lies in the set
$[-\zeta,\zeta]$ where

$\zeta=\|u\|(2\|z\|+\|J^{V}_{T}(0)\|)\,.$ (23)

###### Proof.

See Appendix A.8. ∎

Algorithm 4 Bisection method to solve $l(\alpha)=0$ when $r=1$

Tolerance $\epsilon\geq 0$, number of iterates $N$

Compute the bound $\beta$ from (23), and set $k=0$.

Set $\alpha_{-}=-\beta$ and $\alpha_{+}=\beta$.

Update for $k=0,\cdots,N$:

Set $\alpha_{k}=\tfrac{1}{2}(\alpha_{-}+\alpha_{+})$.

if $l(\alpha_{k})>0$ then

$\alpha_{+}\leftarrow\alpha_{k}$,

else

$\alpha_{-}\leftarrow\alpha_{k}$.

end if

if $k>1$ and $|\alpha_{k}-\alpha_{k-1}|<\epsilon$ then

return $\alpha_{k}$

end if

End

Furthermore, we can combine the bisection method and the semi-smooth Newton
method. Since the non-smooth Newton method is locally convergent, it requires
a starting point in a sufficiently near neighborhood of the solution in order
to guarantee convergence. By using bisection, we can obtain a sequence of
points approaching the solution. When those points reach the neighborhood
required for convergence of the semi-smooth Newton Method, we turn to use the
semi-smooth Newton method to obtain faster convergence. In Proposition 4.2,
$\alpha_{0}$ is required to belong to a neighborhood of $\alpha^{*}$, which
can be achieved by Algorithm 4, i.e., we can assert to find a point $\alpha$
such that $|\alpha-\alpha^{*}|<\delta$ in $\log_{2}((2\zeta/\delta))$ steps,
where $\zeta$ is as in (23).

###  4.3 Implementation of the quasi-Newton Forward-Backward Step

In this section, we consider a fixed iteration $k$ and for convenience, we
omit this index in this subsection. We consider low-rank perturbed metric
$V=M\pm Q$ where $Q=UU^{*}$. Algorithm 1 and Algorithm 2 both have similar
forward-backward steps. And those forward-backward steps require a kind of
calculation in the form of $J^{V}_{A}(z-V^{-1}Bz)$ for some $z$ which can be
significantly simplified by applying Theorem 4.1. Instead of evaluating
$J^{V}_{A}(z-V^{-1}Bz)$ directly, we split the whole task into two subproblems
consisting of a root finding problem and an evaluation of resolvent with
respect to $M$ which is assumed to be simple to calculate. For simplicity, we
denote $\breve{z}=z-V^{-1}Bz$.

######  Proposition 4.4.

The forward-backward step from some $\breve{z}$ to some $\hat{z}$ in Algorithm
1 and Algorithm 2 can be solved by

$\hat{z}=J^{V}_{A}(\breve{z})=J_{A}^{M}(\breve{z}\mp M^{-1}U\alpha^{*})\,.$
(24)

Here, $\alpha^{*}\in\mathbb{R}^{r}$ is the unique root of
$\mathcal{L}\colon\mathbb{R}^{r}\to\mathbb{R}^{r}$,

$\mathcal{L}(\alpha^{*})\coloneqq U^{*}(\breve{z}-J_{A}^{M}(\breve{z}\mp
M^{-1}U\alpha^{*}))+\alpha^{*}=0\,.$ (25)

We notice that the calculation of $V^{-1}$ is expensive when the size of $V$
is large. We can use Sherman-Morrison formula. However the inverse of $M$ is
inevitable and in algorithms that we propose in Section 6 it can be
computationally hard. So, we introduce Proposition 4.5 to avoid computing
$M_{k}^{-1}$.

######  Proposition 4.5.

The forward-backward step from some $\breve{z}$ to some $\hat{z}$ in Algorithm
1 and Algorithm 2 can be solved by

$\hat{z}=J^{V}_{A}(\breve{z})=J_{A}^{M}(z-M^{-1}Bz\mp M^{-1}U\xi^{*}).$ (26)

Here, $\xi^{*}\in\mathbb{R}^{r}$ is the unique zero of
$\mathcal{J}\colon\mathbb{R}^{r}\to\mathbb{R}^{r}$,

$\mathcal{J}(\xi^{*})\coloneqq U^{*}_{k}(z_{k}-J_{A}^{M}(z-M^{-1}Bz\mp
M^{-1}U\xi^{*}))+\xi^{*}=0\,.$ (27)

The function $\mathcal{J}$ is Lipschitz continuous with constant
$1+||M^{-1/2}U||^{2}$ and strictly monotone.

###### Proof.

See Appendix A.9. ∎

##  5 A General 0SR1 Inertial Quasi-Newton Method for Monotone Inclusion

###  5.1 General algorithm with generic metric $M\pm Q_{k}$

We note that if we set in the inclusion problem (10) $A\equiv 0$, and
$B=\nabla f$ of some convex smooth function $f$, the update step consisting of
(12) and (11) reduces to Gradient Descent when $M_{k}\equiv I$ and to the
classic Newton method when $M_{k}=\nabla^{2}f(z_{k})$. Motivated by classic
Newton and quasi-Newton methods, we construct $M_{k}$ as an approximation of
the differential of $Bz_{k}$ at $z_{k}$. We generalize the quasi-Newton method
0SR1 from differentiable $\nabla f$ (SR1 method with 0-memory) to a cocoercive
operator $B$. The approximation $M_{k}$ shall satisfy the modified secant
condition:

$M_{k}s_{k}=y_{k},\quad\mathrm{where}\quad y_{k}=Bz_{k}-Bz_{k-1},\quad
s_{k}=z_{k}-z_{k-1}\,.$ (28)

Choose $M_{0}\in\mathcal{S}_{++}(\mathcal{H})$ which is positive definite. The
update is:

$M_{k}=M_{0}\pm\gamma_{k}u_{k}u_{k}^{*},\quad
u_{k}=(y_{k}-M_{0}s_{k})/\sqrt{|\left\langle
y_{k}-M_{0}s_{k},s_{k}\right\rangle|},$ (29)

where $\gamma_{k}\in[0,+\infty)$ needs to be selected such that $M_{k}$ is
positive definite and $\pm$ depends on the sign of $\left\langle
y_{k}-M_{0}s_{k},s_{k}\right\rangle$. If $\left\langle
y_{k}-M_{0}s_{k},s_{k}\right\rangle>0$, we use

$M_{k}=M_{0}+\gamma_{k}u_{k}u_{k}^{*},$ (30)

and, if $\left\langle y_{k}-M_{0}s_{k},s_{k}\right\rangle<0$, we use

$M_{k}=M_{0}-\gamma_{k}u_{k}u_{k}^{*}.$ (31)

In the rest part of this section, we investigate the conditions for $M_{k}$
defined above to meet Assumption 1 and Assumption 2.

######  Lemma 5.1.

Let $M_{0}$ be symmetric positive definite, $\tfrac{1}{\beta}$ be Lipschitz
constant of operator $B$.

  1. (i)

If $M_{k}=M_{0}+\gamma_{k}u_{k}u_{k}^{*}$ and $\gamma_{k}>0,\,\forall
k\in\mathbb{N}$, $M_{k}$ is always positive definite and symmetric and
Assumption 2 is guaranteed.

  2. (ii)

If $M_{k}=M_{0}-\gamma_{k}u_{k}u_{k}^{*}$ and $\gamma_{k}\leq
c\frac{\rho_{min}(M_{0}-\tfrac{1}{\beta}I)}{\|u_{k}\|^{2}}$ where
$\rho_{min}(M_{0}-\tfrac{1}{\beta}I)>0$ is the smallest eigenvalue of matrix
$M_{0}-\tfrac{1}{\beta}$ and $c\in(0,1)$, then $M_{k}$ is positive definite
and Assumption 2 is satisfied.

Now, in order to guarantee Assumption 1 for $M_{k}\coloneqq M_{0}\pm Q_{k}$
with $Q_{k}=\gamma_{k}u_{k}u_{k}^{*}$, we introduce the following lemma.

######  Lemma 5.2.

Let $(\eta_{k})_{k}\in\ell_{+}^{1}(\mathbb{N})$, $M_{0}-\tfrac{1}{\beta}I$ be
positive definite, $M_{k}\coloneqq M_{0}\pm\gamma_{k}u_{k}u_{k}^{*}$ and
$\gamma_{k}=\frac{\eta_{k}}{\|u_{k}\|_{2}^{2}}$ with
$\eta_{k}\leq\rho_{min}(M_{0}-\tfrac{1}{\beta}I)$ for any $k\in\mathbb{N}$.
Then, Assumption 1 is satisfied.

##  6 Inertial Quasi-Newton PDHG for Saddle-point Problems

In this section, we consider a min-max problem as follows:

$\min_{x\in\mathcal{H}_{1}}\max_{y\in\mathcal{H}_{2}}g(x)+G(x)+\left\langle
Kx,y\right\rangle-f(y)-F(y)$ (32)

with a linear mapping $K\colon\mathcal{H}_{1}\to\mathcal{H}_{2}$, proper lower
semi-continuous convex functions
$g\colon\mathcal{H}_{1}\to\overline{\mathbb{R}}\coloneqq\mathbb{R}\cup\\{\pm\infty\\}$
and $f\colon\mathcal{H}_{2}\to\overline{\mathbb{R}}$, and continuously
differentiable convex functions $G\colon\mathcal{H}_{1}\to\mathbb{R}$ and
$F\colon\mathcal{H}_{2}\to\mathbb{R}$ with Lipschitz continuous gradient. This
problem can be expressed as a special monotone inclusion problem from which we
derive in Algorithm 5 an inertial quasi-Newton Primal-Dual Hybrid Gradient
Method (PDHG) and in Algorithm 6 a quasi-Newton PDHG with relaxation step as
special applications of Algorithm 1 and Algorithm 2 respectively (see
Proposition 6.1). By Fermat’s rule, the optimality condition for (32) is the
following inclusion problem:

$0\in T(z)+B(z)\quad\mathrm{with}\quad z=\begin{pmatrix}x\\\
y\end{pmatrix},\quad\mathrm{where}\,T(z)=\begin{pmatrix}\partial
g(x)+K^{*}y\\\ -Kx+\partial
f(y)\end{pmatrix}\mathrm{and}\,B(x)=\begin{pmatrix}\nabla G(x)\\\ \nabla
F(y)\end{pmatrix}\,.$ (33)

Algorithm 5 Inertial quasi-Newton PDHG

$N\geq 0$, $(\|\epsilon_{k}\|)_{k\in\mathbb{N}}\in\ell_{+}^{1}(\mathbb{N})$

Update for $k=0,\cdots,N$:

  1. (i)

Compute $M_{k}$ according to a quasi-Newton framework.

  2. (ii)

Compute the inertial step with parameter $\alpha_{k}$:

$\begin{split}\bar{x}_{k}&=x_{k}+\alpha_{k}(x_{k}-x_{k-1})\,,\\\
\bar{y}_{k}&=y_{k}+\alpha_{k}(y_{k}-y_{k-1})\,.\\\ \end{split}$ (34)

  3. (iii)

Compute the main quasi-Newton PDHG step:

$\begin{split}x_{k+1}&=\mathrm{prox}^{\mathcal{T}}_{g}(\bar{x}_{k}-\mathcal{T}\nabla
G(\bar{x}_{k})-\mathcal{T}K^{*}\bar{y}_{k}\mp\mathcal{T}U_{k,x}\xi_{k})+\epsilon_{k,x}\,,\\\
y_{k+1}&=\mathrm{prox}^{\Sigma}_{f}(\bar{y}_{k}-\Sigma\nabla
F(\bar{y}_{k})+\Sigma K(2x_{k+1}-\bar{x}_{k})\mp\Sigma
U_{k,y}\xi_{k})+\epsilon_{k,y},\\\ \end{split}$ (35)

where $\xi_{k}$ solves $\mathcal{J}(\xi_{k})=0$ and
$\epsilon_{k}=\begin{pmatrix}\epsilon_{k,x}\\\ \epsilon_{k,y}\end{pmatrix}$ is
the error caused by computation at the $k$-th iterate.

End

Algorithm 6 Quasi-Newton PDHG with relaxation

$N\geq 0$, $(\|\epsilon_{k}\|)_{k\in\mathbb{N}}\in\ell_{+}^{1}(\mathbb{N})$

Update for $k=0,\cdots,N$:

  1. (i)

Compute $M_{k}$ according to quasi-Newton framework.

  2. (ii)

Compute the main quasi-Newton PDHG step:

$\begin{split}\tilde{x}_{k}&=\mathrm{prox}^{\mathcal{T}}_{g}(x_{k}-\mathcal{T}\nabla
G(x_{k})-\mathcal{T}K^{*}y_{k}\mp\mathcal{T}U_{k,x}\xi_{k})+\epsilon_{k,x}\,,\\\
\tilde{y}_{k}&=\mathrm{prox}^{\Sigma}_{f}(y_{k}-\Sigma\nabla F(y_{k})+\Sigma
K(2\tilde{x}_{k}-x_{k})\mp\Sigma U_{k,y}\xi_{k})+\epsilon_{k,y}\,,\\\
\end{split}$ (36)

where $\xi_{k}$ solves $\mathcal{J}(\xi_{k})=0$ and
$\epsilon_{k}=\begin{pmatrix}\epsilon_{k,x}\\\ \epsilon_{k,y}\end{pmatrix}$ is
the error caused by computation at the $k$-th iterate.

  3. (iii)

Relaxation step:

$v_{k}\coloneqq M_{k}\bigg{(}\begin{pmatrix}x_{k}\\\
y_{k}\end{pmatrix}-\begin{pmatrix}\tilde{x}_{k}\\\
\tilde{y}_{k}\end{pmatrix}\bigg{)}+\begin{pmatrix}\nabla G(\tilde{x}_{k})\\\
\nabla F(\tilde{y}_{k})\end{pmatrix}-\begin{pmatrix}\nabla G(x_{k})\\\ \nabla
F(y_{k})\end{pmatrix}$ (37)

to compute the relaxation parameter $t_{k}$

$t_{k}=\frac{\left\langle\begin{pmatrix}x_{k}\\\
y_{k}\end{pmatrix}-\begin{pmatrix}\tilde{x}_{k}\\\
\tilde{y}_{k}\end{pmatrix},v_{k}\right\rangle}{2\|v_{k}\|^{2}}$ (38)

and update $x_{k}$ and $y_{k}$ as follows

$\begin{pmatrix}x_{k+1}\\\
y_{k+1}\end{pmatrix}\leftarrow\begin{pmatrix}x_{k}\\\
y_{k}\end{pmatrix}-t_{k}v_{k}\,.$ (39)

End

######  Proposition 6.1 (PDHG method as a special proximal point algorithm
[17]).

The update step of PDHG can be regarded as a proximal point algorithm with
special metric $M$:

$M(z_{k+1}-z_{k})+T(z_{k+1})\ni-B(z_{k})\quad\text{where}\quad
M=\begin{pmatrix}\mathcal{T}^{-1}&-K^{*}\\\ -K&\Sigma^{-1}\end{pmatrix}\,.$
(40)

Therefore, we are able to apply Algorithm 1 and Algorithm 2 directly to a
saddle point problem by regarding $T$ as $A$ in the inclusion problem (10),
which will give us Algorithm 5 and 6, respectively. We are going to use a
variable metric $M_{k}=M\pm
Q_{k}\in\mathcal{S}(\mathcal{H}_{1}\times\mathcal{H}_{2})$ instead of
$M\in\mathcal{S}(\mathcal{H}_{1}\times\mathcal{H}_{2})$ in (40) for the $k$-th
iterate. Here, we set $Q_{k}=U_{k}U_{k}^{*}$ where
$U_{k}\colon\mathbb{R}^{r}\to\mathrm{im}(Q),\,\alpha\mapsto
U_{k}(\alpha)\coloneqq\sum^{r}_{i=1}\alpha_{i}u_{k,i}$ is an isomorphism
defined by $r$ linearly independent $u_{k,1},...,u_{k,r}\in\mathrm{im}(Q)$ for
each $k$-th iterate. We obtain the update step: Find
$z_{k+1}\in\mathcal{H}_{1}\times\mathcal{H}_{2}$ such that

$M_{k}(z_{k+1}-z_{k})+T(z_{k+1})\ni-B(z_{k})\,.$ (41)

In practice, Proposition 4.5 which is derived from our main result Theorem 4.1
turns out to be more tractable. To show how to calculate the update step,
using Proposition 4.5 in quasi-Newton PDHG in Algorithm 5 and in Algorithm 6,
we introduce Proposition 6.2. Let
$U_{k}\in\mathcal{B}(\mathbb{R}^{r},\mathcal{H}_{1}\times\mathcal{H}_{2})$. We
set $U_{k}=\begin{pmatrix}U_{k,x}\\\ U_{k,y}\end{pmatrix}$ with
$U_{k,x}\in\mathcal{B}(\mathbb{R}^{r},\mathcal{H}_{1})$ and
$U_{k,y}\in\mathcal{B}(\mathbb{R}^{r},\mathcal{H}_{2})$. The validity of the
update step is verified in Proposition 6.2.

######  Proposition 6.2.

The update step from $\bar{z}_{k}$ ($z_{k}$) to $z_{k+1}$ ($\tilde{z}_{k}$) in
Algorithm 5 (Algorithm 6) reduces for the quasi-Newton PDHG update step to
compute

$\begin{cases}x_{k+1}&=\mathrm{prox}^{\mathcal{T}}_{g}(\bar{x}_{k}-\mathcal{T}\nabla
G(\bar{x}_{k})-\mathcal{T}K^{*}\bar{y}_{k}\mp\mathcal{T}U_{k,x}\xi_{k})\\\
y_{k+1}&=\mathrm{prox}^{\Sigma}_{f}(\bar{y}_{k}-\Sigma\nabla
F(\bar{y}_{k})+\Sigma K(2x_{k+1}-\bar{x}_{k})\mp\Sigma
U_{k,y}\xi_{k})\,.\end{cases}$ (42)

where, $\xi_{k}\in\mathbb{R}^{r}$ is the unique zero of
$\mathcal{J}\colon\mathbb{R}^{r}\to\mathbb{R}^{r}$:

$\begin{split}\mathcal{J}(\xi)&=(U_{k,x})^{*}[\underbrace{\mathrm{prox}_{g}^{\mathcal{T}}(\bar{x}_{k}-\mathcal{T}\nabla
G(\bar{x}_{k})-\mathcal{T}K^{*}\bar{y}_{k}\mp\mathcal{T}U_{k,x}\xi)}_{x_{k+1}(\xi)}-\bar{x}_{k}]\\\
&+(U_{k,y})^{*}[\mathrm{prox}_{f}^{\Sigma}(\bar{y}_{k}-\Sigma\nabla
F(\bar{y}_{k})+\Sigma K(2x_{k+1}(\xi)-\bar{x}_{k})\mp\Sigma
U_{k,y}\xi)-\bar{y}_{k}]-\xi\,.\\\ \end{split}$ (43)

######  Remark 4.

By switching $\bar{z}_{k}$ with $z_{k}$ and $z_{k+1}$ with $\tilde{z}_{k}$, we
obtain the update step for Algorithm 6.

We would like to emphasize that using Proposition 6.2, we can avoid the
computation of $M^{-1}$ in the primal and dual setting, which is a
computationally significant advantage. The convergence of Algorithms 5 and 6
is a direct consequence of Theorems 3.1 and 3.2.

######  Proposition 6.3 (Convergence of quasi-Newton PDHG method).

Let $M_{0}=M$ as defined in (40) and $M_{k}=M_{0}\pm U_{k}U_{k}^{*}$.

  1. (i)

If $(M_{k})_{k\in\mathbb{N}}$ satisfies Assumption 1, then $(x_{k},y_{k})$
generated by Algorithm 5 converges weakly to some solution of (2).
Furthermore, if $g$ and $f$ are both strongly convex (or $G$ and $F$ are both
strongly convex), then we obtain the same convergence rate as in Theorem 3.1.

  2. (ii)

If $(M_{k})_{k\in\mathbb{N}}$ satisfies Assumption 2. Thus, $(x_{k},y_{k})$
generated by Algorithm 6 converges weakly to some solution of (2).
Furthermore, if $g$ and $f$ are both strongly convex ( or $G$ and $F$ are both
strongly convex) , then we obtain linear convergence.

##  7 Numerical experiments

The algorithms that we analyze in the experiment are summarized in Table 1.

Algorithm | Algorithm Name | metric
---|---|---
Forward-Backward Splitting (FBS) | Foward-Backward Primal-Dual Hybrid Gradient Method |  $M$ fixed as in (40)
Inertial Forward-Backward Splitting (IFBS) | Inertial Primal-Dual Hybrid Gradient Method |  $M$ fixed as in (40)
quasi-Newton Forward-Backward Splitting (QN-FBS) |  | quasi-Newton
---
Primal-Dual Hybrid Gradient Method Gradient
Variable metric as in (29) with $M_{0}=M$ from (40)
Relaxation quasi-Newton Forward-Backward Splitting (RQN-FBS) |  | Primal-Dual Hybrid Gradient Method
---
with relaxation step
Variable metric as in (29) with $M_{0}=M$ from (40)
Inertial quasi-Newton Forward-Backward Splitting (IQN-FBS) |  | Inertial quasi-Newton
---
Primal-Dual Hybrid Gradient Method
Variable metric as in (29) with $M_{0}=M$ from (40)

Table 1: Summary of algorithms used in the numerical experiments. Details are
provided within each section.

Note that PDHG is used interchangeably as FBS in the later experiments since
PDHG is a specialization of FBS.

###  7.1 TV-$l_{2}$ deconvolution

In this experiment, we solve a problem that is used for an image deconvolution
[9]. Given a blurry and noisy image $b\in\mathbb{R}^{MN}$ (interpreted as a
vector by stacking the M column vectors of length N), we seek to find a clean
image $x\in\mathbb{R}^{MN}$ by solving the following optimization problem:

$\min_{0\leq x\leq 255}\frac{1}{2}\|Ax-b\|_{2}^{2}+\mu\|Dx\|_{2,1}\,,$ (44)

where $A\in\mathbb{R}^{MN\times MN}$ is a linear operator that acts as a
blurring operator and $\|Dx\|_{2,1}$ implements a discrete version of the
isotropic total variation norm of $x$ using simple forward differences in
horizontal and vertical direction. The parameter $\mu>0$ stresses as the
influence of the regularization term $\|Dx\|_{2,1}$ vs the data fidelity term
$\tfrac{1}{2}\|Ax-b\|^{2}$. In order to deal with the non-smoothness, we
rewrite the problem as a saddle point problem.

$\min_{x}\max_{y}\left\langle
Dx,y\right\rangle+\delta_{\Delta}(x)+\frac{1}{2}\|Ax-b\|_{2}^{2}-\delta_{\\{\|\cdot\|_{2,\infty}\leq\mu\\}}(y)\,,$
(45)

where $\Delta\coloneqq\\{x\in\mathbb{R}^{MN}|0\leq x_{i}\leq 255,\forall
i\\}$. We can cast this problem into the general class of problems (32) by
setting $K=D$, $f=\delta_{\\{\|\cdot\|_{2,\infty}\leq\lambda\\}}$,
$G(x)=\frac{1}{2}\|Ax-b\|_{2}^{2}$, $F(p)=0$ and $g=\delta_{V}$. Here, $G$ is
$L$-smooth with $L=1$ since we assume $\rho_{max}(A^{*}A)\leq 1$, which is the
largest eigenvalue of $A^{*}A$. Let $z_{k}=\begin{pmatrix}x_{k}\\\
y_{k}\end{pmatrix}$. We compute the low-rank part
$Q_{k}=\gamma_{k}u_{k}u_{k}^{*}$ by (29) with
$Bz_{k}=\begin{pmatrix}A^{*}Ax_{k}-A^{*}b\\\ 0\end{pmatrix}$, which leads to a
metric that just affects the primal update. In each iteration, we use the
bisection Algorithm 4 and the semi-smooth Newton method 3 to locate the root.

Figure 1: We compare convergence of the inertial quasi-Newton PDHG (IQN-FBS)
to other algorithms in the Table 1 with $\mu=0.001$, $\tau=0.09$, $\sigma=0.9$
and inertial parameter $\alpha_{0}=10$,
$\alpha_{k}=\frac{10}{k^{1.1}(\max\\{\|z_{k}-z_{k-1}\|,\|z_{k}-z_{k-1}\|^{2}\\})}$.
We also observe that the curve of RQN-FBS converges faster at the beginning.
Our three quasi-Newton type algorithm RQN-, QN-, and IQN-FBS clearly
outperform the original FBS and IFBS algorithm.

Figure 1 shows the primal gap where the optimal primal value was computed by
running the original PDHG method for 10000 iterates. For the variable metric
at iterate $k$, we fixed $\gamma_{k}=5/\|u_{k}\|_{2}^{2}$. Thus, the
Assumption 2 is satisfied and the convergence of RQN-FBS (Algorithm 2) is
guaranteed. However, with the same setting, Assumption 1 is not satisfied
since there is no sequence $(\eta_{k})_{k}\in\ell_{+}^{1}$ such that
$M_{k+1}\preceq(1+\eta_{k})M_{k}$ for any $k\in\mathbb{N}$. In this practical
problem, we still observe the convergence of IQN-FBS (Algorithm 1). We notice
that our quasi-Newton type algorithms IQN-FBS, RQN-FBS and QN-FBS are much
faster than original FBS algorithm and inertial FBS (IFBS). This can be
explained by the fact that the Hessian of $G(x)$ is not identity and then by
applying our quasi-Newton SR1 methods, we can adapt metric to the local
geometry of the objective. All three quasi-Newton type algorithms share
similar convergence rates in the end. The extrapolation can not speed up those
algorithms at the beginning as shown by IFBS and IQN-FBS since for first
several iterates $\|z_{k}-z_{k-1}\|$ are large and extrapolation parameter is
small. However it does help eventually.

###  7.2 TV-$l_{2}$ deconvolution with infimal convolution type regularization

A source of optimization problems that fits (32) is derived from the
following:

$\min_{x\in\mathbb{R}^{n}}g(x)+G(x)+(f\square h)(Dx),$ (46)

where $f\square h(\cdot)\coloneqq\inf_{v\in\mathbb{R}^{m}}f(v)+h(\cdot-v)$
denotes the infimal convolution of $f$ and $h$. As a prototypical image
processing problem, we define a regularization term as infimal convolution
between the total variation norm and a weighted squared norm, i.e. $g=0$,
$G(x)=\frac{1}{2}\|Ax-b\|^{2}$, $h(\cdot)=\frac{1}{2}\|W\cdot\|^{2}$ and
$f(\cdot)=\mu\|\cdot\|_{2,1}$. This yields the problem:

$\min_{x}\frac{1}{2}\|Ax-b\|^{2}+\mu R(Dx)$ (47)

where
$R(\cdot)\coloneqq\inf_{v\in\mathbb{R}^{m}}\|v\|_{2,1}+\frac{1}{2\mu}\|W(\cdot-v)\|^{2}$,
$W$ is a diagonal matrix of weights which is given to favor discontinuities
along image edges and $A$, $b$, $D$ are defined as in the first experiment. In
practice, $W$ can be computed by additional edge finding steps or by extra
information. Here, we select $W$ such that $\tfrac{1}{2}\leq\|W\|\leq 1$. The
optimization problem (47) given in primal type can be converted into the
saddle point problem:

$\min_{x}\max_{y}\left\langle
Dx,y\right\rangle+\frac{1}{2}\|Ax-b\|^{2}-\delta_{\\{\|\cdot\|_{2,+\infty}\leq\mu\\}}(y)-\frac{1}{2}\|W^{-1}y\|^{2}\,.$
(48)

We compute the low-rank part $Q_{k}$ by (29) with
$Bz_{k}=\begin{pmatrix}A^{*}Ax_{k}-A^{*}b\\\
(W^{-1})^{*}W^{-1}y_{k}\end{pmatrix}$ which leads to metric that affects both
primal and dual update. Here, $\beta=1/4$. In each iteration, we combine the
bisection (Algorithm 4) and the semi-smooth Newton method (Algorithm 3) to
locate the root.

Figure 2 also shows the primal gap where the optimal primal value was still
computed by running original PDHG for 10000 iterates. For the variable metric
at iterate $k$, we fixed $\gamma_{k}=2/\|u_{k}\|^{2}$. Assumption 2 is always
satisfied. But Assumption 1 is not satisfied. We can observe from Figure 2:
though RQN-FBS turned out to be the slowest one for this experiment, IQN-FBS
is still the fast one. We observe that two quasi-Newton type methods (IQN-FBS
and QN-FBS) converge faster than IFBS and FBS. Inertial methods (IQN-FBS,
IFBS) are slightly faster respectively than QN-FBS and FBS.

Figure 2: We compare convergence of the inertial quasi-Newton PDHG (IQN-FBS)
to other algorithms in the Table 1 with $\mu=0.01$, $\tau=0.1$, $\sigma=0.1$
and the extrapolation parameter $\alpha_{0}=1$,
$\alpha_{k}=\max\\{\frac{10}{k^{1.1}(\max\\{\|z_{k}-z_{k-1}\|,\|z_{k}-z_{k-1}\|^{2}\\})},1\\}$.
Our quasi-Newton type algorithm IQN-FBS can converge faster than the original
FBS and IFBS algorithm.

####  7.2.1 Image denoising

Let $A=I$ in (48). We then obtain an image denoising problem with special norm
defined by infimal convolution of total variation and weighted norm, which has
strong convexity for both primal part and dual part. Besides, due to the
simple formula, we obtain dual problem explicitly which means we can calculate
primal and dual gap. The dual problem reads:

$\max_{\|y\|_{2,\infty}\leq
1}-\frac{1}{2}\|K^{*}y-b\|^{2}-\frac{1}{2}\|W^{-1}y\|^{2}\,.$ (49)

Figure 3: We compare convergence of the inertial quasi-Newton PDHG (IQN-FBS)
to other algorithms in the Table 1 with $\mu=0.1$, $\tau=0.1$, $\sigma=0.1$
and extrapolation parameter $\alpha_{0}=10$,
$\alpha_{k}=\frac{10}{\max\\{k^{2},k^{2}\|z_{k}-z_{k-1}\|^{2}\\}}$. All
algorithms converge linearly and faster than $O(\frac{1}{1.1^{k}})$. Our
quasi-Newton type algorithm IQN-FBS can converge faster than the original FBS
and IFBS algorithm.

Figure 3 shows the primal dual gap and it will decrease to zero by using any
algorithm from Table 1. To construct $Q_{k}$, we use
$\gamma_{k}=\frac{2}{\|u_{k}\|^{2}}$. As we can observe from Figure 3, we have
linear convergence for quasi-Newton type methods as what we expected in
Theorem 3.1 and 3.2.

###  7.3 Conclusion

In this paper, we extended the framework of [7] for variable metrics to the
setting of resolvent operators, solving efficiently the monotone inclusion
problem (10) consisting of a set-valued operator $A$ and a single-valued
operator $B$. We proposed two variants of quasi-Newton Forward-Backward
Splitting. We develop a general efficient Resolvent calculus that perfectly
applies to this quasi-Newton setting. The convergence of the variant with
relaxation requires mild assumptions on the metric which are easier to
satisfy, whereas the other variant implements an inertial feature and is
therefore often fast. As a special case of this framework, we developed an
inertial quasi-Newton primal-dual algorithm that can be flexibly applied to a
large class of saddle point problems.

##  Appendix A Appendix

###  A.1 Proof of Lemma 2.3

###### Proof.

We assume $y=M^{-1/2}\circ J_{M^{-1/2}TM^{-1/2}}\circ M^{1/2}(z)$. Then, we
obtain

$\begin{split}y&=M^{-1/2}\circ J_{M^{-1/2}TM^{-1/2}}\circ M^{1/2}(z)\\\
M^{1/2}y&=J_{M^{-1/2}TM^{-1/2}}\circ M^{1/2}(z)\\\
M^{1/2}y&=(I+M^{-1/2}TM^{-1/2})^{-1}(M^{1/2}z)\\\
(M^{1/2}z)&\in(I+M^{-1/2}TM^{-1/2})M^{1/2}y\\\
(M^{1/2}z)&\in(M^{1/2}+M^{-1/2}T)y\\\ Mz&\in(M+T)y\\\ y&=(M+T)^{-1}(Mz)\\\
y&=J^{M}_{T}(z)\,.\end{split}$ (50)

∎

###  A.2 Proof of Lemma 2.4

###### Proof.

This proof is adapted from [5]. Let $(u,v)=(J^{M}_{A}(x),J^{M}_{A}(y))$ for
some $x,y\in\mathcal{H}$. By the definition of resolvent operator $J^{M}_{A}$,
we obtain

$u=J^{M}_{A}(x)\iff M(x-u)\in Au\,.$

Similarly, we obtain $M(y-v)\in Av$. Then $\gamma_{A}$-strong monotonicity of
$A$ yields

$\begin{split}\left\langle
M(x-u)-M(y-v),u-v\right\rangle&\geq\gamma_{A}\|u-v\|^{2}\\\ \left\langle
M(x-y),u-v\right\rangle-\left\langle
M(u-v),u-v\right\rangle&\geq\gamma_{A}\|u-v\|^{2}\\\ \left\langle
M(x-y),u-v\right\rangle&\geq\gamma_{A}\|u-v\|^{2}+\|u-v\|_{M}^{2}\,.\\\
\end{split}$

Since $\|M\|$ is bounded by $C$, we obtain

$\left\langle
M(x-y),u-v\right\rangle\geq\gamma_{A}\|u-v\|^{2}+\|u-v\|_{M}^{2}\geq(1+\tfrac{\gamma_{A}}{C})\|u-v\|_{M}^{2}\,.$
(51)

Consequently, $J^{M}_{A}$ is $(1+\tfrac{\gamma_{A}}{C})$-cocoercive and
Lipschitz continuous with constant $1/(1+\tfrac{\gamma_{A}}{C})$ with respect
to the norm $\|\cdot\|_{M}$. ∎

###  A.3 Proof of Theorem 3.1

###### Proof.

For convenience, we set $B_{k}\coloneqq M_{k}^{-1}B$. Fix $z\in zer(A+B)$
which is equivalent to $z=J_{A}^{M_{k}}(z-B_{k}z)$. We set
$\hat{z}_{k+1}\coloneqq J_{A}^{M_{k}}(\bar{z}_{k}-B_{k}\bar{z}_{k})$, i.e.,
$z_{k+1}=\hat{z}_{k+1}+\epsilon_{k}$.
Boundedness:
First, we are going to show that $\|z_{k}-z\|_{M_{k}}$ is bounded. The
following are several useful estimations we will use later. Assumption 1
yields

$\|z_{k+1}-z\|_{M_{k+1}}^{2}\leq(1+\eta_{k})\|z_{k+1}-z\|_{M_{k}}^{2}\,.$ (52)

Since $z_{k+1}=\hat{z}_{k+1}+\epsilon_{k}$, it follows that

$\|z_{k+1}-z\|_{M_{k}}^{2}=\|\hat{z}_{k+1}-z+\epsilon_{k}\|_{M_{k}}^{2}\leq\|\hat{z}_{k+1}-z\|_{M_{k}}^{2}+2\|\epsilon_{k}\|_{M_{k}}\|\hat{z}_{k+1}-z\|_{M_{k}}+\|\epsilon_{k}\|_{M_{k}}^{2}\,.$
(53)

The assumption that $B$ is $\beta$-cocoercive yields that

$\left\langle\bar{z}_{k}-z,B_{k}\bar{z}_{k}-B_{k}z\right\rangle_{M_{k}}=\left\langle\bar{z}_{k}-z,B\bar{z}_{k}-Bz\right\rangle\geq\beta\|B\bar{z}_{k}-Bz\|^{2}\,.$
(54)

The assumption that
$M_{k}-\tfrac{1}{2\beta}I\in\mathcal{S}_{\epsilon}(\mathcal{H})$ yields that
$2\beta I-M_{k}^{-1}\in\mathcal{S}_{++}(\mathcal{H})$.
The cocoercivity of $B$ and monotonicity of $B$ yield the following
estimations $\mathrm{(I)}$ and $\mathrm{(II)}$ respectively:

$\displaystyle\|(\bar{z}_{k}-B_{k}\bar{z}_{k})-(z-B_{k}z)\|_{M_{k}}^{2}$
$\displaystyle=\|\bar{z}_{k}-z\|_{M_{k}}^{2}-2\left\langle\bar{z}_{k}-z,B_{k}\bar{z}_{k}-B_{k}z\right\rangle_{M_{k}}+\|B_{k}\bar{z}_{k}-B_{k}z\|_{M_{k}}^{2}$
$\displaystyle\overset{\mathrm{(i)}}{=}\|\bar{z}_{k}-z\|_{M_{k}}^{2}-2\left\langle\bar{z}_{k}-z,B\bar{z}_{k}-Bz\right\rangle+\|B\bar{z}_{k}-Bz\|_{M_{k}^{-1}}^{2}$
$\displaystyle\overset{\mathrm{(ii)}}{\leq}\|\bar{z}_{k}-z\|_{M_{k}}^{2}-\|B\bar{z}_{k}-Bz\|_{2\beta-
M_{k}^{-1}}^{2}$ (I)
$\displaystyle\overset{\mathrm{(iii)}}{\leq}\|\bar{z}_{k}-z\|_{M_{k}}^{2}-\|B\bar{z}_{k}-Bz\|_{\beta-
M_{k}^{-1}}^{2}-\gamma_{B}\|\bar{z}_{k}-z\|^{2}\,,$ (II)

where $\mathrm{(i)}$ uses
$\|B_{k}\bar{z}_{k}-B_{k}z\|_{M_{k}}^{2}=\|B\bar{z}_{k}-Bz\|_{M_{k}^{-1}}^{2}$,
$\mathrm{(ii)}$ uses (54) and $\mathrm{(iii)}$ uses monotonicity of $B$. Note
that we use the shorthand $\beta-M_{k}^{-1}$ for $\beta I-M_{k}^{-1}$. The
fact that $J^{M_{k}}_{A}$ is firmly non-expansive since $A$ is maximally
monotone with respect to $M_{k}$ implies that

$\begin{split}\|\hat{z}_{k+1}-z\|_{M_{k}}^{2}&=\|J^{M_{k}}_{A}(\bar{z}_{k}-B_{k}\bar{z}_{k})-J^{M_{k}}_{A}(z-B_{k}z)\|_{M_{k}}^{2}\\\
&\leq\|(\bar{z}_{k}-B_{k}\bar{z}_{k})-(z-B_{k}z)\|_{M_{k}}^{2}\\\
&\quad-\|(I-J^{M_{k}}_{A})(\bar{z}_{k}-B_{k}\bar{z}_{k})-(I-J^{M_{k}}_{A})(z-B_{k}z)\|_{M_{k}}^{2}\\\
&\leq\|\bar{z}_{k}-z\|_{M_{k}}^{2}-\|B\bar{z}_{k}-Bz\|_{2\beta-
M_{k}^{-1}}^{2}-\|(\bar{z}_{k}-\hat{z}_{k+1})-(B_{k}\bar{z}_{k}-B_{k}z)\|_{M_{k}}^{2}\,,\\\
\end{split}$ (55)

where the last inequality uses (I). It follows from Assumption 1 that the term
$\|B\bar{z}_{k}-Bz\|_{2\beta-M_{k}^{-1}}^{2}\geq 0$. We continue to bound the
first term on the right hand side of (55). Using [5, Lemma 2.14] and the
definition of $\bar{z}_{k}$, we obtain the following:

$\|\bar{z}_{k}-z\|_{M_{k}}^{2}=(1+\alpha_{k})\|z_{k}-z\|_{M_{k}}^{2}-\alpha_{k}\|z_{k-1}-z\|_{M_{k}}^{2}+(1+\alpha_{k})\alpha_{k}\|z_{k}-z_{k-1}\|_{M_{k}}^{2}\,,$
(56)

and by using the triangle inequality, we also obtain another estimation:

$\|\bar{z}_{k}-z\|_{M_{k}}\leq(1+\alpha_{k})\|z_{k}-z\|_{M_{k}}+\alpha_{k}\|z_{k-1}-z\|_{M_{k}}\,.$
(57)

In order to address complete update step, we make the following estimation:

$\begin{split}\|z_{k+1}-z\|_{M_{k+1}}^{2}&\leq(1+\eta_{k})\|z_{k+1}-z\|_{M_{k}}^{2}\\\
&\overset{\mathrm{(i)}}{\leq}(1+\eta_{k})\big{(}\|\hat{z}_{k+1}-z\|_{M_{k}}^{2}+2\|\epsilon_{k}\|_{M_{k}}\|\hat{z}_{k+1}-z\|_{M_{k}}+\|\epsilon_{k}\|_{M_{k}}^{2}\big{)}\\\
&\overset{\mathrm{(ii)}}{\leq}(1+\eta_{k})\big{(}\|\bar{z}_{k}-z\|_{M_{k}}^{2}+2\|\epsilon_{k}\|_{M_{k}}\|\bar{z}_{k}-z\|_{M_{k}}+\|\epsilon_{k}\|_{M_{k}}^{2}\\\
&\quad-\|B\bar{z}_{k}-Bz\|_{2\beta-
M_{k}^{-1}}^{2}-\|(\bar{z}_{k}-\hat{z}_{k+1})-(B_{k}\bar{z}_{k}-B_{k}z)\|_{M_{k}}^{2}\big{)}\\\
&\overset{\mathrm{(iii)}}{\leq}(1+\eta_{k})\big{(}(1+\alpha_{k})\|z_{k}-z\|_{M_{k}}^{2}-\alpha_{k}\|z_{k-1}-z\|_{M_{k}}^{2}\\\
&\quad+(1+\alpha_{k})\alpha_{k}\|z_{k}-z_{k-1}\|_{M_{k}}^{2}\\\
&\quad+2\|\epsilon_{k}\|_{M_{k}}((1+\alpha_{k})\|z_{k}-z\|_{M_{k}}+\alpha_{k}\|z_{k-1}-z\|_{M_{k}})\\\
&\quad+\|\epsilon_{k}\|_{M_{k}}^{2}-\|B\bar{z}_{k}-Bz\|_{2\beta-
M_{k}^{-1}}^{2}-\|(\bar{z}_{k}-\hat{z}_{k+1})-(B_{k}\bar{z}_{k}-B_{k}z)\|_{M_{k}}^{2}\big{)}\\\
&\overset{\mathrm{(iv)}}{\leq}(1+\eta_{k})\big{(}\|z_{k}-z\|_{M_{k}}^{2}\\\
&\quad+\alpha_{k}(\|z_{k}-z\|_{M_{k}}-\|z_{k-1}-z\|_{M_{k}})(\|z_{k}-z\|_{M_{k}}+\|z_{k-1}-z\|_{M_{k}})\\\
&\quad+(1+\alpha_{k})\alpha_{k}\|z_{k}-z_{k-1}\|_{M_{k}}^{2}\\\
&\quad+2\|\epsilon_{k}\|_{M_{k}}((1+\alpha_{k})\|z_{k}-z\|_{M_{k}}+\alpha_{k}\|z_{k-1}-z\|_{M_{k}})\\\
&\quad+\|\epsilon_{k}\|_{M_{k}}^{2}-\|B\bar{z}_{k}-Bz\|_{2\beta-
M_{k}^{-1}}^{2}-\|(\bar{z}_{k}-\hat{z}_{k+1})-(B_{k}\bar{z}_{k}-B_{k}z)\|_{M_{k}}^{2}\big{)}\\\
&\overset{\mathrm{(v)}}{\leq}(1+\eta_{k})\big{(}\|z_{k}-z\|_{M_{k}}^{2}+\alpha_{k}\|z_{k}-z_{k-1}\|_{M_{k}}(\|z_{k}-z\|_{M_{k}}+\|z_{k-1}-z\|_{M_{k}})\\\
&\quad+(1+\alpha_{k})\alpha_{k}\|z_{k}-z_{k-1}\|_{M_{k}}^{2}+2\|\epsilon_{k}\|_{M_{k}}((1+\alpha_{k})\|z_{k}-z\|_{M_{k}}+\alpha_{k}\|z_{k-1}-z\|_{M_{k}})\\\
&\quad+\|\epsilon_{k}\|_{M_{k}}^{2}-\|B\bar{z}_{k}-Bz\|_{2\beta-
M_{k}^{-1}}^{2}-\|(\bar{z}_{k}-\hat{z}_{k+1})-(B_{k}\bar{z}_{k}-B_{k}z)\|_{M_{k}}^{2}\big{)}\,.\\\
\end{split}$ (58)

where $\mathrm{(i)}$ uses (53), $\mathrm{(ii)}$ uses (55), $\mathrm{(iii)}$
uses (56) and (57), $\mathrm{(iv)}$ uses factorization of the quadratic, and
$\mathrm{(v)}$ uses the triangle inequality to obtain the bound
$\|z_{k}-z\|_{M_{k}}-\|z_{k-1}-z\|_{M_{k}}\leq\|z_{k}-z_{k-1}\|_{M_{k}}$.

Now, our goal is to conclude boundedness using Lemma 2.7. For simplicity, we
set:

$\begin{cases}e_{k}&\coloneqq\alpha_{k}\|z_{k}-z_{k-1}\|_{M_{k}}^{2}\\\
r_{k}&\coloneqq\alpha_{k}\|z_{k}-z_{k-1}\|_{M_{k}}\\\
\theta_{k}&\coloneqq\|z_{k}-z\|_{M_{k}}\\\
m_{k}&\coloneqq\|z_{k-1}-z\|_{M_{k}}\\\
p_{k}&\coloneqq\|B\bar{z}_{k}-Bz\|_{2\beta-M_{k}^{-1}}^{2}\\\
q_{k}&\coloneqq\|(\bar{z}_{k}-\hat{z}_{k+1})-(B_{k}\bar{z}_{k}-B_{k}z)\|_{M_{k}}^{2}\,.\\\
\end{cases}$

By Assumption 1, we have $m_{k}\leq(1+\eta_{k-1})\theta_{k-1}$. Without loss
of generality, we can assume $0<\eta_{k}<1$ for any $k\in\mathbb{N}$.
Replacing each term in (58) with new corresponding notations, we obtain:

$\begin{split}\theta_{k+1}^{2}&\leq(1+\eta_{k})\big{(}\theta^{2}_{k}+r_{k}(\theta_{k}+m_{k})\\\
&\quad+(1+\alpha_{k})e_{k}+2\|\epsilon_{k}\|_{M_{k}}((1+\alpha_{k})\theta_{k}+\alpha_{k}m_{k})+\|\epsilon_{k}\|_{M_{k}}^{2}-p_{k}-q_{k}\big{)}\\\
&\leq(1+\eta_{k})\big{(}\theta^{2}_{k}+r_{k}(\theta_{k}+(1+\eta_{k-1})\theta_{k-1})\\\
&\quad+(1+\alpha_{k})e_{k}+2\|\epsilon_{k}\|_{M_{k}}((1+\alpha_{k})\theta_{k}+\alpha_{k}(1+\eta_{k-1})\theta_{k-1})+\|\epsilon_{k}\|_{M_{k}}^{2}\big{)}\\\
&\leq(1+\eta_{k})\big{(}\theta^{2}_{k}+r_{k}(\theta_{k}+2\theta_{k-1})+(1+\Lambda)e_{k}+2\|\epsilon_{k}\|_{M_{k}}((1+\Lambda)\theta_{k}+2\Lambda\theta_{k-1})+\|\epsilon_{k}\|_{M_{k}}^{2}\big{)}\,,\\\
\end{split}$ (59)

where the last inequality uses $0<\alpha_{k}\leq\Lambda$ and $0<\eta_{k}<1$.
Now, we claim that $\theta_{k}$ is bounded in two steps. We introduce an
auxiliary bounded sequence $(C_{k})_{k\in\mathbb{N}}$ (step 1) such that
$\theta_{k}\leq C_{k}$ for any $k\in\mathbb{N}$ (step 2). The boundedness of
$\theta_{k}$ follows from that of $C_{k}$.
Step1: We construct a sequence $C_{k}$ as the following:

$\begin{cases}C_{0}&=\max\\{\theta_{0},1\\}\,,\\\
C_{k}&=(1+\eta_{k})C_{k}+\nu_{k}\,,\\\ \end{cases}$ (60)

where
$\nu_{k}=(1+\eta_{k})((1+\Lambda)e_{k}+2r_{k}+(1+3\Lambda)\|\epsilon_{k}\|_{M_{k}})$.
From our assumptions, it holds that $(M_{k})_{k\in\mathbb{N}}$ is bounded from
above,
$(r_{k})_{k\in\mathbb{N}}\in\ell_{+}^{1}(\mathbb{N}),(e_{k})_{k\in\mathbb{N}}\in\ell_{+}^{1}(\mathbb{N})$
and $(\|\epsilon_{k}\|)_{k\in\mathbb{N}}\in\ell_{+}^{1}(\mathbb{N})$, which
implies $(\nu_{k})_{k\in\mathbb{N}}\in\ell_{+}^{1}(\mathbb{N})$. Using Lemma
2.7, we obtain the convergence of $C_{k}$ to some $C_{\infty}<+\infty$.
Step2: From the update step of (60), we observe that $(C_{k})_{k}$ is a non-
decreasing sequence and $C_{k}\geq 1$ for any $k\in\mathbb{N}$. We claim that
for each $k$, $\theta_{k}\leq C_{k}$. We argument by induction. Clearly, we
have $\theta_{0}\leq C_{0}$. Assume $\theta_{i}\leq C_{i}$ holds true for
$i\leq k$. Then, (59) yields that

$\begin{split}\theta_{k+1}^{2}&\leq(1+\eta_{k})\big{(}C_{k}^{2}+r_{k}(C_{k}+2C_{k-1})+(1+\Lambda)e_{k}+2\|\epsilon_{k}\|_{M_{k}}((1+\Lambda)C_{k}+2\Lambda
C_{k-1})+\|\epsilon_{k}\|_{M_{k}}^{2}\big{)}\\\
&\overset{\mathrm{(*)}}{\leq}(1+\eta_{k})\big{(}C_{k}^{2}+4r_{k}C_{k}+2(1+\Lambda)e_{k}C_{k}+2(1+3\Lambda)\|\epsilon_{k}\|_{M_{k}}C_{k}+\|\epsilon_{k}\|_{M_{k}}^{2}\big{)}\\\
&\leq(1+\eta_{k})\big{(}C_{k}+2r_{k}+(1+\Lambda)e_{k}+(1+3\Lambda)\|\epsilon_{k}\|_{M_{k}}\big{)}^{2}\,,\\\
\end{split}$ (61)

where $\mathrm{(*)}$ uses $C_{k}\geq C_{k-1}\geq 1$ and $r_{k}>0$. By the
definition of $C_{k+1}$, we obtain

$\begin{split}\theta_{k+1}&\overset{\mathrm{(i)}}{\leq}\sqrt{1+\eta_{k}}\big{(}C_{k}+2r_{k}+(1+\Lambda)e_{k}+(1+3\Lambda)\|\epsilon_{k}\|_{M_{k}}\big{)}\\\
&\overset{\mathrm{(ii)}}{\leq}(1+\eta_{k})C_{k}+(1+\eta_{k})(2r_{k}+(1+\Lambda)e_{k}+(1+3\Lambda)\|\epsilon_{k}\|_{M_{k}})\\\
&\overset{\mathrm{(iii)}}{\leq}(1+\eta_{k})C_{k}+\nu_{k}\\\ &=C_{k+1}\,,\\\
\end{split}$ (62)

where $\mathrm{(i)}$ uses (61), $\mathrm{(ii)}$ holds true since
$(1+\eta_{k})>1$ and $\mathrm{(iii)}$ uses definition of $\nu_{k}$. This
concludes the induction, and we deduce that $\theta_{k}$ is bounded and
therefore, $z_{k}$ and $\bar{z}_{k}$ are both bounded.
Weak convergence:
This part of the proof is adapted from the one for [11, Theorem 4.1]. Since
$\theta_{k}$ is bounded, we set
$\zeta\coloneqq\sup_{k\in\mathbb{N}}\theta_{k}$. The last inequality in (58)
implies that

$\begin{split}\theta_{k+1}^{2}&\leq(1+\eta_{k})\big{(}\theta_{k}^{2}+r_{k}(\zeta+2\zeta)+(1+\Lambda)e_{k}+2(1+3\Lambda)\|\epsilon_{k}\|_{M_{k}}\zeta+\|\epsilon_{k}\|_{M_{k}}^{2}-p_{k}-q_{k}\big{)}\\\
&\leq(1+\eta_{k})\big{(}\theta_{k}^{2}+3r_{k}\zeta+(1+\Lambda)e_{k}+2(1+3\Lambda)\|\epsilon_{k}\|_{M_{k}}\zeta+\|\epsilon_{k}\|_{M_{k}}^{2}-p_{k}-q_{k}\big{)}\\\
&\leq\theta_{k}^{2}+\eta_{k}\theta_{k}^{2}+(1+\eta_{k})\big{(}3r_{k}\zeta+(1+\Lambda)e_{k}+(2+6\Lambda)\|\epsilon_{k}\|_{M_{k}}\zeta+\|\epsilon_{k}\|_{M_{k}}^{2}\big{)}-p_{k}-q_{k}\\\
&\leq\theta_{k}^{2}+\underbrace{\eta_{k}\zeta^{2}+2(3r_{k}\zeta+(1+\Lambda)e_{k}+(2+6\Lambda)\|\epsilon_{k}\|_{M_{k}}\zeta+\|\epsilon_{k}\|_{M_{k}}^{2})}_{{\delta_{k}}}-p_{k}-q_{k}\,.\\\
\end{split}$ (63)

We set
$\delta_{k}\coloneqq\eta_{k}\zeta^{2}+2(3r_{k}\zeta+(1+\Lambda)e_{k}+(2+6\Lambda)\|\epsilon_{k}\|_{M_{k}}\zeta+\|\epsilon_{k}\|_{M_{k}}^{2})$
and observe that $(\delta_{k})_{k\in\mathbb{N}}\in\ell_{+}^{1}(\mathbb{N})$.
Now, (63) yields that

$\theta_{k+1}^{2}\leq\theta_{k}^{2}+\delta_{k}\,.$ (64)

Using (64) and Lemma 2.7, we obtain the convergence of
$\theta_{k}^{2}=\|z_{k}-z\|_{M_{k}}^{2}$ for any $z\in zer(A+B)$. Rearranging
(63) to $p_{k}\leq\theta_{k}^{2}-\theta_{k+1}^{2}+\delta_{k}$, using
Assumption 1 and summing it for $k=0,\cdots,N$ , we obtain

$\begin{split}\epsilon\sum_{k=0}^{N}\|B\bar{z}_{k}-Bz\|^{2}&\leq\sum_{k=0}^{N}\|B\bar{z}_{k}-Bz\|_{2\beta-
M_{k}^{-1}}^{2}=\sum_{k=0}^{N}p_{k}\leq\theta_{0}^{2}-\theta_{N}^{2}+\sum_{k=0}^{N}\delta_{k}\leq\zeta^{2}+\sum_{k=0}^{N}\delta_{k}\,.\end{split}$
(65)

Since $(\delta_{k})_{k\in\mathbb{N}}\in\ell^{1}_{+}(\mathbb{N})$, by taking
limit as $N\to+\infty$, we obtain

$\sum_{k\in\mathbb{N}}\|B\bar{z}_{k}-Bz\|^{2}\leq\frac{1}{\epsilon^{2}}(\zeta^{2}+\sum_{k\in\mathbb{N}}\delta_{k})<+\infty\,.$
(66)

Similarly, we obtain from (63) using
$q_{k}\leq\theta_{k}^{2}-\theta_{k+1}^{2}+\delta_{k}$ that

$\sum_{k\in\mathbb{N}}\|(\bar{z}_{k}-\hat{z}_{k+1})-(B_{k}\bar{z}_{k}-B_{k}z)\|_{M_{k}}^{2}<+\infty\,.$
(67)

Set $z^{*}$ as an arbitrary weak sequential cluster point of
$(z_{k})_{k\in\mathbb{N}}$, namely, a subsequence $z_{k_{n}}\rightharpoonup
z^{*}$ as $n\to\infty$.

In order to obtain weak convergence of $z_{k}$, by Proposition 2.5 with
$\varphi(t)=t^{2}$ and (64) and Assumption 1, it suffices to show that
$z^{*}\in zer(A+B)$. It follows from the selection of $\alpha_{k}$ that:

$\begin{split}\|\bar{z}_{k}-z_{k}\|&\leq\alpha_{k}\|z_{k}-z_{k-1}\|\to 0\,.\\\
\end{split}$ (68)

Thus, (68) yields $\bar{z}_{k_{n}}\rightharpoonup z^{*}$. From (66), we obtain
that $B\bar{z}_{k_{n}}\to Bz$ as $n\to\infty$. Since $B$ is cocoercive, it is
maximally monotone and we can use the weak strong graph closedness of $B$ in
Proposition 2.2 to infer that $(z^{*},Bz)\in\mathrm{Graph}\,B$, i.e. $Bz\in
Bz^{*}$. However, since $B$ is single valued, we obtain $Bz^{*}=Bz$ and hence
$B\bar{z}_{k_{n}}\to Bz^{*}$. Setting $u_{k}\coloneqq
M_{k}(\bar{z}_{k}-\hat{z}_{k+1})-B\bar{z}_{k}$, by definition of the resolvent
$J^{M_{k}}_{A}$, we have $u_{k}\in A(\hat{z}_{k+1})$ for all $k\in\mathbb{N}$.
From (67), we obtain as $k\to+\infty$,

$\begin{split}\|u_{k}+Bz^{*}\|&=\|M_{k}(\bar{z}_{k}-\hat{z}_{k+1}-B_{k}\bar{z}_{k}+B_{k}z^{*})\|\\\
&\leq C\|\bar{z}_{k}-\hat{z}_{k+1}-B_{k}\bar{z}_{k}+B_{k}z^{*}\|\\\
&\leq\frac{C}{\sigma}\|\bar{z}_{k}-\hat{z}_{k+1}-B_{k}\bar{z}_{k}+B_{k}z^{*}\|_{M_{k}}\to
0\,.\\\ \end{split}$ (69)

Furthermore, from (66) and (67), we have

$\begin{split}\|\bar{z}_{k}-\hat{z}_{k+1}\|&\leq\|\bar{z}_{k}-\hat{z}_{k+1}-B_{k}\bar{z}_{k}+B_{k}z^{*}\|+\|B_{k}\bar{z}_{k}-B_{k}z^{*}\|\\\
&\leq\|\bar{z}_{k}-\hat{z}_{k+1}-B_{k}\bar{z}_{k}+B_{k}z^{*}\|+\frac{1}{\sigma}\|B\bar{z}_{k}-Bz^{*}\|\to
0\,.\\\ \end{split}$ (70)

Therefore, together with (70), $\bar{z}_{k_{n}}\rightharpoonup z^{*}$ implies
$\hat{z}_{k_{n}+1}\rightharpoonup z^{*}$ as $n\to\infty$. Now we already have
$u_{k_{n}}\to-Bz^{*}$ as $n\to\infty$ and

$(\forall
k\in\mathbb{N})\colon\quad(\hat{z}_{k_{n}+1},u_{k_{n}})\in\mathrm{Graph}\,A\,.$
(71)

Since $A$ is maximally monotone and using Proposition 2.2, we infer that
$-Bz^{*}\in Az^{*}$, hence $z^{*}\in zer(A+B)$. As mentioned above, the result
follows from Proposition 2.5 with $\varphi(t)=t^{2}$.
Convergence rate:
In the following part, we are going to show the convergence rate of Algorithm
1: Assume $\epsilon_{k}\equiv 0$ for $k\in\mathbb{N}$ and either
$\gamma_{A}>0$ or $\gamma_{B}>0$. Because of (52), Assumption 1 and Lipschitz
continuity of $J_{A}^{M_{k}}$, we obtain for any $z\in zer(A+B)$ that

$\begin{split}\|z_{k+1}-z\|_{M_{k+1}}^{2}&\leq(1+\eta_{k})\|z_{k+1}-z\|_{M_{k}}^{2}\\\
&\leq(1+\eta_{k})\Big{(}\frac{1}{1+\frac{\gamma_{A}}{C}}\Big{)}^{2}\|(\bar{z}_{k}-B_{k}\bar{z}_{k})-(z-B_{k}z)\|_{M_{k}}^{2}\\\
&\overset{\mathrm{(i)}}{\leq}(1+\eta_{k})\Big{(}\frac{1}{1+\frac{\gamma_{A}}{C}}\Big{)}^{2}\big{(}\|\bar{z}_{k}-z\|_{M_{k}}^{2}-\|B\bar{z}_{k}-Bz\|_{\beta-
M_{k}^{-1}}^{2}-\gamma_{B}\|\bar{z}_{k}-z\|^{2}\big{)}\\\
&\overset{\mathrm{(ii)}}{\leq}(1+\eta_{k})\Big{(}\frac{1}{1+\frac{\gamma_{A}}{C}}\Big{)}^{2}(1-\frac{\gamma_{B}}{C})\big{(}\|\bar{z}_{k}-z\|_{M_{k}}^{2}\big{)}\\\
&\overset{\mathrm{(iii)}}{=}(1+\eta_{k})\frac{(1-\frac{\gamma_{B}}{C})}{(1+\frac{\gamma_{A}}{C})^{2}}\big{(}(1+\alpha_{k})\|z_{k}-z\|_{M_{k}}^{2}-\alpha_{k}\|z_{k-1}-z\|_{M_{k}}^{2}\\\
&\quad+(1+\alpha_{k})\alpha_{k}\|z_{k}-z_{k-1}\|_{M_{k}}^{2}\big{)}\\\
&=(1+\eta_{k})\frac{(1-\frac{\gamma_{B}}{C})}{(1+\frac{\gamma_{A}}{C})^{2}}\big{(}\|z_{k}-z\|_{M_{k}}^{2}+\alpha_{k}(\|z_{k}-z\|_{M_{k}}^{2}-\|z_{k-1}-z\|_{M_{k}}^{2})\\\
&\quad+(1+\alpha_{k})\alpha_{k}\|z_{k}-z_{k-1}\|_{M_{k}}^{2}\big{)}\\\
&\overset{\mathrm{(iv)}}{=}(1+\eta_{k})\frac{(1-\frac{\gamma_{B}}{C})}{(1+\frac{\gamma_{A}}{C})^{2}}\big{(}\|z_{k}-z\|_{M_{k}}^{2}+\alpha_{k}\|z_{k}-z_{k-1}\|_{M_{k}}(\|z_{k}-z\|_{M_{k}}\\\
&\quad+\|z_{k-1}-z\|_{M_{k}})+(1+\alpha_{k})\alpha_{k}\|z_{k}-z_{k-1}\|_{M_{k}}^{2}\big{)}\\\
&\overset{\mathrm{(v)}}{=}(1+\eta_{k})\frac{(1-\frac{\gamma_{B}}{C})}{(1+\frac{\gamma_{A}}{C})^{2}}\big{(}\|z_{k}-z\|_{M_{k}}^{2}+\alpha_{k}\|z_{k}-z_{k-1}\|_{M_{k}}(\|z_{k}-z\|_{M_{k}}\\\
&\quad+(1+\eta_{k-1})\|z_{k-1}-z\|_{M_{k-1}})+(1+\alpha_{k})\alpha_{k}\|z_{k}-z_{k-1}\|_{M_{k}}^{2}\big{)}\\\
&\overset{\mathrm{(vi)}}{=}(1+\eta_{k})\frac{(1-\frac{\gamma_{B}}{C})}{(1+\frac{\gamma_{A}}{C})^{2}}\big{(}\|z_{k}-z\|_{M_{k}}^{2}+3\alpha_{k}\zeta\|z_{k}-z_{k-1}\|_{M_{k}}\\\
&\quad+(1+\Lambda)\alpha_{k}\|z_{k}-z_{k-1}\|_{M_{k}}^{2}\big{)}\\\
\end{split}$ (72)

where $\mathrm{(i)}$ uses (II), $\mathrm{(ii)}$ uses the fact $M_{k}$ is
bounded uniformly and the assumption that
$M_{k}-\tfrac{1}{\beta}I\in\mathcal{S}_{\epsilon}(\mathcal{H})$,
$\mathrm{(iii)}$ uses (56), $\mathrm{(iv)}$ uses factorization of the
quadratic and uses the triangle inequality to obtain the bound
$\|z_{k}-z\|_{M_{k}}-\|z_{k-1}-z\|_{M_{k}}\leq\|z_{k}-z_{k-1}\|_{M_{k}}$,
$\mathrm{(v)}$ uses Assumption 1 and $\mathrm{(vi)}$ uses boundedness of
$\alpha_{k}$ and $\|z_{k}-z\|_{M_{k}}$. Since either $\gamma_{A}>0$ or
$\gamma_{B}>0$,
$\frac{1-\frac{\gamma_{B}}{C}}{(1+\frac{\gamma_{A}}{C})^{2}}<1$. Then there
exists sufficient large $K_{0}>0$ such that for any $k>K_{0}$,
$(1+\eta_{k})\Big{(}\frac{1-\frac{\gamma_{B}}{C}}{(1+\frac{\gamma_{A}}{C})^{2}}\Big{)}<1-\xi<1$
for some $\xi\in(0,1)$. Thus, we infer that for any $k>K_{0}$:

$\|z_{k}-z\|_{M_{k}}^{2}\leq(1-\xi)^{k-K_{0}}\|z_{K_{0}}-z\|_{M_{K_{0}}}^{2}+\sum_{i=K_{0}}^{k-1}(1-\xi)^{k-i}\alpha_{i}(3\zeta\|z_{i}-z_{i-1}\|_{M_{i}}+(1+\Lambda)\|z_{i}-z_{i-1}\|_{M_{i}}^{2})\,.$
(73)

Let $\Theta=3\zeta+(1+\Lambda)$. Therefore (73) can be simplified as the
following:

$\|z_{k}-z\|_{M_{k}}^{2}\leq(1-\xi)^{k-K_{0}}\|z_{K_{0}}-z\|_{M_{K_{0}}}^{2}+\sum_{i=K_{0}}^{k-1}\Theta(1-\xi)^{k-i}\alpha_{i}\max\\{\|z_{i}-z_{i-1}\|_{M_{i}},\|z_{i}-z_{i-1}\|_{M_{i}}^{2}\\}\,.$
(74)

Since $z_{k}$ is bounded and $M_{k}\in S_{\sigma}(\mathcal{H})$, it follows
that

$\|z_{k}-z\|^{2}\leq\tfrac{1}{\sigma}(1-\xi)^{k-K_{0}}\|z_{K_{0}}-z\|_{M_{K_{0}}}^{2}+O(\sum_{i=K_{0}}^{k-1}(1-\xi)^{k-i}\alpha_{i})\,.$
(75)

Furthermore, if $\alpha_{i}\equiv 0$, for $k>K_{0}$ and for any $z\in
zer(A+B)$, we obtain linear convergence:

$\|z_{k+1}-z\|^{2}\leq\frac{1}{\sigma}(1-\xi)^{k-K_{0}}\|z_{K_{0}}-z\|_{M_{K_{0}}}^{2}\,.$
(76)

If $\alpha_{k}\neq 0$, $\alpha_{k}=O(\tfrac{1}{k^{2}})$ and $K_{0}$ large
enough, then $\|z_{k+1}-z\|^{2}$ converges in the rate of $O(\tfrac{1}{k})$
for $k>K_{0}$ according to [29, Lemma 2.2.4 (Chung)]; if $\alpha_{k}\neq 0$
and $\alpha_{k}=O(q^{k})$ for $q=1-\xi$ and $k>K_{0}$, then
$\|z_{k+1}-z\|^{2}$ converges in the rate of $O(kq^{k})$ for $k>K_{0}$ since
(75).

∎

###  A.4 Proof of Theorem 3.2

###### Proof.

For simplicity, we set

$\begin{cases}\hat{z}_{k}\coloneqq
J_{A}^{M_{k}}(z_{k}-M_{k}^{-1}Bz_{k})\mathrm{,i.e.,}\tilde{z}_{k}=\hat{z}+\epsilon_{k}\,,\\\
\delta_{k}\coloneqq\max\\{1,\rho\\}\|\epsilon_{k}\|\,,\end{cases}$ (77)

where $\rho=\sqrt{\tfrac{C}{\sigma}}\tfrac{(C+\tfrac{1}{\beta})}{2c}$ and
$(\delta_{k})_{k\in\mathbb{N}}\in\ell_{+}^{1}(\mathbb{N})$. We set $z^{*}$
such that $-Bz^{*}\in Az^{*}$.
Boundedness:
We claim that by choosing proper $t_{k}$ for each $k\in\mathbb{N}$, we have

$\|z_{k+1}-z^{*}\|\leq\|z_{k}-z^{*}\|+\delta_{k}.$ (78)

Note that if (78) is satisfied, it follows from Lemma 2.7 that
$\|z_{k}-z^{*}\|$ is bounded and converges.

To stress the relation between $z_{k+1}$ and $t_{k}$, we define
$z(t_{k})\coloneqq z_{k}-t_{k}[(M_{k}-B)(z_{k}-\tilde{z}_{k})]$ and we will
use $z_{k+1}$ and $z(t_{k})$ interchangeably. We also set
$\gamma_{k}(t_{k})\coloneqq(\delta_{k}+||z_{k}-z^{*}||)^{2}-||z(t_{k})-z^{*}||^{2}$.

In order to prove (78), it is sufficient to show that for proper $t_{k}$ at
each iterate, $\gamma_{k}(t_{k})>0$. It results from the definition of
$\gamma_{k}(t_{k})$ that

$\begin{split}\gamma_{k}(t_{k})&=(||z_{k}-z^{*}||+\delta_{k})^{2}-||z(t_{k})-z^{*}||^{2}\\\
&=\left\langle
z_{k}-z^{*}+z(t_{k})-z^{*},z_{k}-z(t_{k})\right\rangle+2\delta_{k}||z_{k}-z^{*}||+\delta_{k}^{2}\\\
&=\underbrace{2t_{k}\left\langle
z_{k}-z^{*},(M_{k}-B)(z_{k}-\tilde{z}_{k})\right\rangle}_{\mathrm{(I)}}-\underbrace{t_{k}^{2}||(M_{k}-B)(z_{k}-\tilde{z}_{k})||^{2}}_{\mathrm{(II)}}\\\
&\quad+2\delta_{k}||z_{k}-z^{*}||+\delta_{k}^{2}\,.\end{split}$ (79)

We need several useful properties to estimate $\mathrm{(I)}$ and
$\mathrm{(II)}$.

It follows from the definition of $\hat{z}_{k}$ that
$-M_{k}(\hat{z}_{k}-z_{k})-Bz_{k}\in A\hat{z}_{k}$, and $\gamma_{A}$-strong
monotonicity of $A$ implies that

$\begin{split}\gamma_{A}\|\hat{z}_{k}-z^{*}\|^{2}&\leq\left\langle\hat{z}_{k}-z^{*},-M_{k}(\hat{z}_{k}-z_{k})-Bz_{k}+Bz^{*}\right\rangle\\\
&=\left\langle\hat{z}_{k}-z^{*},M_{k}(z_{k}-\hat{z}_{k})-Bz_{k}+B\hat{z}_{k}-B\hat{z}_{k}+Bz^{*}\right\rangle\\\
&=\left\langle\hat{z}_{k}-z^{*},M_{k}(z_{k}-\hat{z}_{k})-Bz_{k}+B\hat{z}_{k}\right\rangle-\left\langle\hat{z}_{k}-z^{*},B\hat{z}_{k}-Bz^{*}\right\rangle\,.\\\
\end{split}$ (80)

We deduce by (80) and $\gamma_{B}$-strong monotonicity of $B$ that

$\begin{split}\left\langle\hat{z}_{k}-z^{*},M_{k}(z_{k}-\hat{z}_{k})-Bz_{k}+B\hat{z}_{k}\right\rangle&\geq\left\langle\hat{z}_{k}-z^{*},B\hat{z}_{k}-Bz^{*}\right\rangle+\gamma_{A}\|\hat{z}_{k}-z^{*}\|^{2}\\\
&\geq\gamma_{B}\|\hat{z}_{k}-z^{*}\|^{2}+\gamma_{A}\|\hat{z}_{k}-z^{*}\|^{2}\,.\\\
\end{split}$ (81)

$J^{M_{k}}_{A}$ is Lipschitz since $A$ is monotone. Therefore,

$\begin{split}\|\hat{z}_{k}-z^{*}\|_{M_{k}}^{2}&\leq\|(z_{k}-B_{k}z_{k})-(z^{*}-B_{k}z^{*})\|_{M_{k}}^{2}\\\
&=\|z_{k}-z^{*}\|_{M_{k}}^{2}-2\left\langle
z_{k}-z^{*},B_{k}z_{k}-B_{k}z^{*}\right\rangle_{M_{k}}+\|B_{k}z_{k}-B_{k}z^{*}\|_{M_{k}}^{2}\\\
&\leq\|z_{k}-z^{*}\|_{M_{k}}^{2}-\|Bz_{k}-Bz\|_{2\beta-M^{-1}_{k}}\\\
&\overset{(*)}{\leq}\|z_{k}-z^{*}\|_{M_{k}}^{2}\,,\end{split}$ (82)

where $(*)$ uses Assumption 2 that
$M_{k}-\frac{1}{\beta}I\in\mathcal{S}_{c}(\mathcal{H})$. Using the assumption
that $M_{k}-\frac{1}{\beta}I\in\mathcal{S}_{c}(\mathcal{H})$ again, we obtain

$\left\langle
z_{k}-\tilde{z}_{k},(M_{k}-B)(z_{k}-\tilde{z}_{k})\right\rangle\geq\|z_{k}-\tilde{z}_{k}\|_{M_{k}-\tfrac{1}{\beta}}^{2}>0\,.$
(83)

Combining (81), (82), the first term $\mathrm{(I)}$ in (79) can be estimated
by the following:

$\begin{split}\mathrm{(I)}&=2t_{k}\left\langle
z_{k}-z^{*},(M_{k}-B)(z_{k}-\tilde{z}_{k})\right\rangle\\\
&=2t_{k}\left\langle
z_{k}-\tilde{z}_{k},(M_{k}-B)(z_{k}-\tilde{z}_{k})\right\rangle+2t_{k}\left\langle\tilde{z}_{k}-\hat{z}_{k},(M_{k}-B)(z_{k}-\tilde{z}_{k})\right\rangle\\\
&\quad+2t_{k}\left\langle\hat{z}_{k}-z^{*},(M_{k}-B)(z_{k}-\hat{z}_{k})\right\rangle+2t_{k}\left\langle\hat{z}_{k}-z^{*},(M_{k}-B)(\hat{z}_{k}-\tilde{z}_{k})\right\rangle\\\
&\overset{\mathrm{(i)}}{\geq}2t_{k}\left\langle
z_{k}-\tilde{z}_{k},(M_{k}-B)(z_{k}-\tilde{z}_{k})\right\rangle+2t_{k}\left\langle\tilde{z}_{k}-\hat{z}_{k},(M_{k}-B)(z_{k}-\tilde{z}_{k})\right\rangle\\\
&\quad+2t_{k}(\gamma_{A}+\gamma_{B})\|\hat{z}_{k}-z^{*}\|^{2}+2t_{k}\left\langle\hat{z}_{k}-z^{*},(M_{k}-B)(\hat{z}_{k}-\tilde{z}_{k})\right\rangle\\\
&\overset{\mathrm{(ii)}}{\geq}2t_{k}\left\langle
z_{k}-\tilde{z}_{k},(M_{k}-B)(z_{k}-\tilde{z}_{k})\right\rangle-2t_{k}\|\epsilon_{k}\|\|(M_{k}-B)(z_{k}-\tilde{z}_{k})\|\\\
&\quad+2t_{k}(\gamma_{A}+\gamma_{B})\|\hat{z}_{k}-z^{*}\|^{2}-2t_{k}\|(M_{k}-B)\epsilon_{k}\|_{M_{k}^{-1}}\|\hat{z}_{k}-z^{*}\|_{M_{k}}\\\
&\overset{\mathrm{(iii)}}{\geq}2t_{k}\left\langle
z_{k}-\tilde{z}_{k},(M_{k}-B)(z_{k}-\tilde{z}_{k})\right\rangle-\|\epsilon_{k}\|^{2}-t_{k}^{2}\|(M_{k}-B)(z_{k}-\tilde{z}_{k})\|^{2}\\\
&\quad+2t_{k}(\gamma_{A}+\gamma_{B})\|\hat{z}_{k}-z^{*}\|^{2}-2\tfrac{1}{\sqrt{\sigma}}(C+\tfrac{1}{\beta})t_{k}\|\epsilon_{k}\|\|z_{k}-z^{*}\|_{M_{k}}\\\
&\overset{\mathrm{(iv)}}{\geq}2t_{k}\left\langle
z_{k}-\tilde{z}_{k},(M_{k}-B)(z_{k}-\tilde{z}_{k})\right\rangle-\|\epsilon_{k}\|^{2}-t_{k}^{2}\|(M_{k}-B)(z_{k}-\tilde{z}_{k})\|^{2}\\\
&\quad+2t_{k}(\gamma_{A}+\gamma_{B})\|\hat{z}_{k}-z^{*}\|^{2}-2\sqrt{\tfrac{C}{\sigma}}(C+\tfrac{1}{\beta})t_{k}\|\epsilon_{k}\|\|z_{k}-z^{*}\|\,,\\\
\end{split}$ (84)

where $\mathrm{(i)}$ uses (81), $\mathrm{(ii)}$ uses Cauchy inequality and
$\mathrm{(iii)}$ uses (82), $2ab\leq a^{2}+b^{2}$ and Assumption 2. We set
$b_{k}\coloneqq\left\langle
z_{k}-\tilde{z}_{k},(M_{k}-B)(z_{k}-\tilde{z}_{k})\right\rangle$ and
$a_{k}\coloneqq\|(M_{k}-B)(z_{k}-\tilde{z}_{k})\|^{2}$. The definition of
$\delta_{k}$ yields that $\delta_{k}\geq\|\epsilon_{k}\|^{2}$ and it follows
from (84) that:

$\begin{split}\gamma_{k}(t_{k})&\geq
2t_{k}b_{k}-2t_{k}^{2}a_{k}+2\delta_{k}\|z_{k}-z^{*}\|+\delta_{k}^{2}-\|\epsilon_{k}\|^{2}\\\
&\quad-2\sqrt{\tfrac{C}{\sigma}}(C+\tfrac{1}{\beta})t_{k}\|\epsilon_{k}\|\|z_{k}-z^{*}\|+2t_{k}(\gamma_{A}+\gamma_{B})\|\hat{z}_{k}-z^{*}\|^{2}\\\
&\geq\underbrace{2t_{k}b_{k}-2t_{k}^{2}a_{k}}_{\mathrm{(III)}}+2t_{k}(\gamma_{A}+\gamma_{B})\|\hat{z}_{k}-z^{*}\|^{2}+\underbrace{2(\delta_{k}-\sqrt{\tfrac{C}{\sigma}}(C+\tfrac{1}{\beta})t_{k})\|\epsilon_{k}\|\|z_{k}-z^{*}\|}_{\mathrm{(IV)}}\,.\end{split}$
(85)

We continue to find a proper $t_{k}$ such that $\gamma_{k}(t_{k})>0$. This
goal boils down to ensuring both $\mathrm{(III)}$ and $\mathrm{(IV)}$ are
positive.

Let $t_{k}=\tfrac{b_{k}}{2a_{k}}$. We first show $t_{k}=\tfrac{b_{k}}{2a_{k}}$
is the proper value to make sure $\mathrm{(III)}$ is positive. From (83), we
observe that $b_{k}>0$. Since $a_{k}>0$ and $b_{k}>0$, the quadratic term
$\mathrm{(III)}$ in (85) will be zero for $t_{k}=0$ or
$t_{k}=\tfrac{b_{k}}{a_{k}}$ and will be strictly positive for any
$t_{k}\in(0,b_{k}/a_{k})$ with the maximum value obtained at
$t_{k}=\tfrac{b_{k}}{2a_{k}}$. As a result, $\mathrm{(III)}$ is strictly
positive.

Second, we will show $\mathrm{(IV)}$ is positive when
$t_{k}=\tfrac{b_{k}}{2a_{k}}$. We observe that
$0<t_{k}=\tfrac{b_{k}}{2a_{k}}<\tfrac{1}{2c}$ for Assumption 2. Thus, the
definition of $\delta_{k}$ and $t_{k}=\tfrac{b_{k}}{2a_{k}}$ imply that

$\begin{split}\mathrm{(IV)}&=2(\delta_{k}-\sqrt{\tfrac{C}{\sigma}}(C+\tfrac{1}{\beta})t_{k})\|\epsilon_{k}\|\|z_{k}-z^{*}\|\\\
&\geq(2\delta_{k}-\sqrt{\tfrac{C}{\sigma}}\tfrac{(C+\tfrac{1}{\beta})}{c})\|\epsilon_{k}\|\|z_{k}-z^{*}\|\\\
&\geq 0\,.\end{split}$ (86)

Since $\mathrm{(III)}$ and $\mathrm{(IV)}$ both are positive when
$t_{k}=\tfrac{b_{k}}{2a_{k}}$, (85) yields:

$\begin{split}\gamma_{k}(t_{k})&\geq
2t_{k}b_{k}-2t_{k}^{2}a+2t_{k}(\gamma_{A}+\gamma_{B})\|\hat{z}_{k}-z^{*}\|^{2}+(2\delta_{k}-\sqrt{\tfrac{C}{\sigma}}\tfrac{(C+\tfrac{1}{\beta})}{c})\|\epsilon_{k}\|\|z_{k}-z^{*}\|\\\
&\geq\frac{b_{k}^{2}}{2a_{k}}+2t_{k}(\gamma_{A}+\gamma_{B})\|\hat{z}_{k}-z^{*}\|^{2}\\\
&>0\,,\\\ \end{split}$ (87)

It results from (87) and the definition of $\gamma_{k}(t_{k})$ that for each
$k\in\mathbb{N}$,
$(||z_{k}-z^{*}||+\delta_{k})\geq||z(t_{k})-z^{*}||=||z_{k+1}-z^{*}||$. We
conclude that if $t_{k}=\tfrac{b_{k}}{2a_{k}}$, then $\gamma_{k}(t_{k})>0\
\text{for all}\ k\in\mathbb{N}$ and the sequence $\|z_{k}-z^{*}\|$ is bounded
and converges as $k\to+\infty$ by using Lemma 2.7.
Weak convergence:
The sequence $(z_{k})_{k\in\mathbb{N}}$ generated by Algorithm 2 is bounded
and $\|z_{k}-z\|$ converges as $k\to\infty$ and $\gamma_{k}(t_{k})$ converges
to zero as $k\to\infty$ for all $z\in zer(A+B)$. Set $z^{*}$ as an arbitrary
weak sequential cluster point of $(z_{k})_{k\in\mathbb{N}}$ and there exists a
subsequence $(z_{k_{n}})_{n\in\mathbb{N}}$ such that $z_{k_{n}}\rightharpoonup
z^{*}$.

In order to obtain weak convergence of $z_{k}$, by Proposition 2.5 with
$\varphi(t)=t$ and fixed metric $M_{k}=I$ and (78), it suffices to show that
$z^{*}\in zer(A+B)$.

Using Assumption 2, $(M_{k})_{k\in\mathbb{N}}$ is bounded uniformly by $C$.
Together with boundedness of operator $B$, we obtain $a_{k}$ is bounded by
$(C+\tfrac{1}{\beta})^{2}\|z_{k}-\tilde{z}_{k}\|^{2}$ for each
$k\in\mathbb{N}$. Using Assumption 2, we obtain that
$b_{k}\geq\|z_{k}-\tilde{z}_{k}\|_{M_{k}-\tfrac{1}{\beta}I}^{2}\geq
c\|z_{k}-\tilde{z}_{k}\|^{2}$. By definition of $t_{k}$ and (87), we have
$\gamma_{k}(t_{k})\geq\tfrac{c}{2(C+\tfrac{1}{\beta})^{2}}\|z_{k}-\tilde{z}_{k}\|^{2}$.
Since $\gamma_{k}(t_{k})\to 0$, $\|z_{k}-\tilde{z}_{k}\|\to 0$ as
$k\to\infty$. Moreover, since $\epsilon_{k}\to 0$, $\|z_{k}-\hat{z}_{k}\|\to
0$. We set $u_{k}\coloneqq M_{k}(z_{k}-\hat{z}_{k})+B\hat{z}_{k}-Bz_{k}$.
Therefore, we obtain that $u_{k}\to 0$ as $k\to+\infty$ and
$\hat{z}_{k_{n}}\rightharpoonup z^{*}$ as $n\to+\infty$. We observe that
$u_{k}\in A\hat{z}_{k}+B\hat{z}_{k}$. Then, by using Proposition 2.2 and the
fact that $A+B$ is maximally monotone, we conclude that $0\in Az^{*}+Bz^{*}$
and $\|z_{k}-z^{*}\|$ decreases since $z^{*}\in zer(A+B)$. As mentioned above,
the result follows from Proposition 2.5 with $\varphi(t)=t$.
Linear convergence rate:
If we assume $\epsilon_{k}\equiv 0$, then $\hat{z}_{k}=\tilde{z}_{k}$ and
$\delta_{k}\equiv 0$. Therefore, from (87) we can obtain an estimation for
$\gamma_{k}(t_{k})$ when $t_{k}=\tfrac{b_{k}}{2a_{k}}$:

$\begin{split}\gamma_{k}(t_{k})=\|z_{k}-z^{*}\|^{2}-\|z_{k+1}-z^{*}\|^{2}\geq\frac{b_{k}^{2}}{2a_{k}}+2t_{k}(\gamma_{A}+\gamma_{B})\|\hat{z}_{k}-z^{*}\|^{2}\,.\\\
\end{split}$ (88)

The following part is to derive linear convergence for the case that either
$\gamma_{A}>0$ or $\gamma_{B}>0$. The definition of $b_{k}$ and that of
$a_{k}$ yield the following estimation for $t_{k}$:

$t_{k}=\frac{b_{k}}{2a_{k}}=\frac{\left\langle
z_{k}-\hat{z}_{k},(M_{k}-B)(z_{k}-\hat{z}_{k})\right\rangle}{2\|(M_{k}-B)(z_{k}-\hat{z}_{k})\|^{2}}\overset{(\mathrm{i})}{\geq}\frac{\|z_{k}-\hat{z}_{k}\|_{M_{k}-\tfrac{1}{\beta}}^{2}}{2(C+\tfrac{1}{\beta})^{2}\|z_{k}-\hat{z}_{k}\|^{2}}\overset{(\mathrm{ii})}{>}\frac{c}{2(C+\tfrac{1}{\beta})^{2}},$
(89)

where both $(\mathrm{i})$ and $(\mathrm{ii})$ use Assumption 2. For
convenience, we denote $\tfrac{c}{2(C+\tfrac{1}{\beta})^{2}}$ by $\delta$.
Using Assumption 2 again, we have the estimation for the first term at the
right hand side of (88):

$\tfrac{b_{k}^{2}}{2a_{k}}\geq\delta\|z_{k}-\hat{z}_{k}\|_{M_{k}-\tfrac{1}{\beta}I}^{2}>c\delta\|z_{k}-\hat{z}_{k}\|^{2}\,.$
(90)

Furthermore, combining (88) with (89), (90) and the definition
$\gamma_{k}(t_{k}):=\|z_{k}-z^{*}\|^{2}-\|z_{k+1}-z^{*}\|^{2}$, we obtain

$\begin{split}\|z_{k+1}-z^{*}\|^{2}&\leq\|z_{k}-z^{*}\|^{2}-2(\gamma_{A}+\gamma_{B})\delta\|\hat{z}_{k}-z^{*}\|^{2}-c\delta\|z_{k}-\hat{z}_{k}\|^{2}\\\
&\leq\|z_{k}-z^{*}\|^{2}-\tfrac{1}{2}\min\\{2(\gamma_{A}+\gamma_{B})\delta,c\delta\\}(2\|\hat{z}_{k}-z^{*}\|^{2}+2\|z_{k}-\hat{z}_{k}\|^{2})\\\
&\overset{(\mathrm{i})}{\leq}\|z_{k}-z^{*}\|^{2}-\tfrac{1}{2}\min\\{2(\gamma_{A}+\gamma_{B})\delta,c\delta\\}(\|\hat{z}_{k}-z^{*}\|+\|z_{k}-\hat{z}_{k}\|)^{2}\\\
&\overset{(\mathrm{ii})}{\leq}\|z_{k}-z^{*}\|^{2}-\tfrac{1}{2}\min\\{2(\gamma_{A}+\gamma_{B})\delta,c\delta\\}\|z_{k}-z^{*}\|^{2}\\\
&\leq(1-\tfrac{1}{2}\min\\{2(\gamma_{A}+\gamma_{B})\delta,c\delta\\})\|z_{k}-z^{*}\|^{2}\,,\\\
\end{split}$ (91)

where $(\mathrm{i})$ uses inequality $2x^{2}+2y^{2}\geq(x+y)^{2}$ and
$(\mathrm{ii})$ uses triangle inequality. Consequently, we obtain linear
convergence if $(\gamma_{A}+\gamma_{B})>0$:

$\|z_{k}-z^{*}\|^{2}\leq(1-\xi)^{k}\|z_{0}-z^{*}\|^{2}\,,$ (92)

where $\xi=\tfrac{1}{2}\min\\{2(\gamma_{A}+\gamma_{B})\delta,c\delta\\}>0$. ∎

###  A.5 Proof of Theorem 4.1

###### Proof.

Computing the resolvent operator shows the following equivalences

$\begin{split}&x^{*}=J^{V}_{T}(z)=(I+V^{-1}T)^{-1}(z)\\\
\iff&Vz\in(V+T)(x^{*})\\\ \iff&Mz\in(M+T)(x^{*})\pm Q(x^{*}-z)\\\
[y^{*}=M^{1/2}x^{*}]\iff&Mz\in(M+T)(M^{-1/2}y^{*})\pm Q(M^{-1/2}y^{*}-z)\\\
\iff&Mz\in(M+T)(M^{-1/2}y^{*})\pm QM^{-1/2}(y^{*}-M^{1/2}z)\\\
\iff&Mz\in(M^{1/2}+TM^{-1/2})(y^{*})\pm QM^{-1/2}(y^{*}-M^{1/2}z)\\\
[W=M^{-1/2}QM^{-1/2}]\iff&M^{1/2}z\in(I+M^{-1/2}TM^{-1/2})(y^{*})\pm
W(y^{*}-M^{1/2}z)\,.\\\ \end{split}$ (93)

Since existence of $x^{*}$ is guaranteed by the properties of $J_{T}^{V}$,
Lemma 2.6 yields the existence of a unique primal-dual pair $(x^{*},u^{*})$
that satisfies the equivalent relations in Lemma 2.6 with $B\coloneqq\pm
W\circ\tau_{M^{1/2}z}$ and
$A\coloneqq\tau_{M^{1/2}z}\circ(I+M^{-1/2}TM^{-1/2})$. The mapping $B$ is
single-valued and , as $T$ is a maximally monotone operator and $M$ is
positive definite, $A^{-1}=J_{M^{-1/2}TM^{-1/2}}\circ\tau_{-M^{1/2}z}$ is
single valued. Therefore, the solution of $J_{T}^{V}$ can be computed by
finding $u^{*}$ such that

$0\in B^{-1}u^{*}-A^{-1}(-u^{*})=[(\pm
W)^{-1}+M^{1/2}z]-J_{M^{-1/2}TM^{-1/2}}(M^{1/2}z-u^{*})$ (94)

and the using

$x^{*}=M^{-1/2}y^{*}\quad\mathrm{and}\quad
y^{*}=A^{-1}(-u^{*})=J_{M^{-1/2}TM^{=1/2}}(M^{1/2}z-u^{*})\,.$ (95)

Substituting $u^{*}=M^{-1/2}v^{*}$ in both problems, multiplying the former
one from left with $M^{-1/2}$, and using $M^{-1/2}W^{-1}M^{-1/2}=Q^{+}$ leads
to

$\begin{cases}0\in\pm Q^{+}v^{*}+z-M^{-1/2}\circ J_{M^{-1/2}TM^{-1/2}}\circ
M^{1/2}(z-M^{-1}v^{*})\\\ x^{*}=M^{-1/2}\circ J_{M^{-1/2}TM^{-1/2}}\circ
M^{1/2}(z-M^{-1}v^{*})\,.\end{cases}$ (96)

Since a solution to $J^{V}_{T}$ exists, there exists $v^{*}\in\mathrm{im}(Q)$
that satisfies the inclusion. Given the isomorphism
$U\colon\mathbb{R}^{r}\to\mathrm{im}(Q)$, which can be realized using $r$
linearly independent $u_{1},\cdots,u_{r}\in\mathcal{H}$ by
$\alpha\to\sum_{i=1}^{r}\alpha_{i}u_{i}$, the inclusion problem is equivalent
to finding the unique root $\alpha^{*}\in\mathbb{R}^{r}$ of $l(\alpha)$, where
$U^{*}$ denotes the adjoint of $U$ and for the case “$-$” of “$\pm$” we
substituted $\alpha$ by $-\alpha$.

The following shows that $l(\alpha)$ is Lipschitz continuous with constant
$\|U^{*}Q^{+}U\|+\|M^{-1/2}U\|^{2}$:(We abbreviate $J_{M^{-1/2}TM^{-1/2}}$ by
$J$ in the following)

$\begin{split}&\left\langle l(\alpha)-l(\beta),\alpha-\beta\right\rangle\\\
=&\|\alpha-\beta\|_{U^{*}Q^{+}U}^{2}-\left\langle J(M^{1/2}z\mp
M^{-1/2}U\alpha)-J(M^{1/2}z\mp
M^{-1/2}U\beta),M^{-1/2}U(\alpha-\beta)\right\rangle\\\
\leq&\|U^{*}Q^{+}U\|\|\alpha-\beta\|^{2}+\|M^{-1/2}U\|^{2}\|\alpha-\beta\|^{2}\,,\\\
\end{split}$ (97)

where, in the last line, we use the 1-Lipschitz continuity (non-expansive) of
$J$.

The following shows strict monotonicity of $l$. We rewrite $l(\alpha)$ as
follows:

$\begin{split}l(\alpha)&=U^{*}Q^{+}U\alpha+U^{*}M^{-1/2}(M^{1/2}z-J_{M^{-1/2}TM^{-1/2}}(M^{1/2}z\mp
M^{-1/2}U\alpha))\\\ &=U^{*}Q^{+}U\alpha\pm
U^{*}M^{-1}U\alpha+U^{*}M^{-1/2}(I-J_{M^{-1/2}TM^{-1/2}})(M^{1/2}z\mp
M^{-1/2}U\alpha)\\\ &=U^{*}(Q^{+}\pm
M^{-1})U\alpha+U^{*}M^{-1/2}J_{M^{-1/2}T^{-1}M^{-1/2}}(M^{1/2}z\mp
M^{-1/2}U\alpha)\,.\\\ \end{split}$ (98)

Using the 1-cocoercivity of $J_{M^{1/2}T^{-1}M^{1/2}}$, the function
$l(\alpha)$ can be seen to be strictly monotone if $\alpha\mapsto
U^{*}(Q^{+}\pm M^{-1})U\alpha$ is strictly monotone. This fact is clear for
the case “$+$” of “$\pm$”. Therefore, in the remainder, we show strictly
monotonicity of $\alpha\mapsto
U^{*}(Q^{+}-M^{-1})U\alpha=U^{*}M^{-1/2}(M^{1/2}Q^{+}M^{1/2}-I)M^{-1/2}U\alpha$.
We observe $M-Q\in\mathcal{S}_{0}(\mathcal{H})$ implies that
$\|M^{-1/2}QM^{-1/2}\|<1$ and by $1\leq\|AA^{-1}\|\leq\|A\|\|A^{-1}\|$, we
conclude that $\|M^{1/2}Q^{+}M^{1/2}\|_{\mathrm{im}(M^{-1/2}Q)}$ for the
restriction of the operator norm to $\mathrm{im}(M^{-1/2}Q)$, hence,
$Q^{+}-M^{-1}\in\mathcal{S}_{++}(\mathcal{H})$.

According to Lemma 2.3, we can replace $M^{-1/2}\circ
J_{M^{-1/2}TM^{-1/2}}\circ M^{1/2}$ with $J^{M}_{T}$. Then we obtain the
formula in the statement of Theorem 4.1. ∎

###  A.6 Proof of Theorem 4.1.1

###### Proof.

Since $Q=UU^{*}$ for some $U\in\mathcal{B}(\mathbb{R}^{r},\mathcal{H})$, we
have

$\mathrm{im}(Q)=\\{UU^{*}v|v\in\mathcal{H}\\}=\\{U\alpha|\alpha\in\mathbb{R}^{r}\\}\,.$
(99)

Since $Q^{+}$ is pseudo-inverse of $Q$ on $im(Q)$, the following holds for
arbitrary $v\in\mathcal{H}$:

$QQ^{+}Qv=Qv\iff UU^{*}Q^{+}UU^{*}v=UU^{*}v\iff UU^{*}Q^{+}U\alpha=U\alpha\,.$
(100)

Since the column vectors $\\{u_{i}\\}_{i=1,\cdots,r}$ of $U$ are independent
with each other, $UU^{*}Q^{+}U\alpha=U\alpha$ yields that
$U^{*}Q^{+}U\alpha=\alpha$. Therefore, the root finding problem in Theorem 4.1
simplifies to (22). ∎

###  A.7 Proof of Proposition 4.2

###### Proof.

Our proof relies on the convergence result [15, Theorem 7.5.5]. By the same
argument as the one in Appendix B.5 paper [7], we obtain
$\partial^{c}l(\alpha^{*})$ is non-singular. If $l(\alpha)$ is tame, then by
[8, Theorem 1], $l(\alpha)$ is semi-smooth. In order to apply this result it
remains to show that $l(\alpha)$ is tame. The property of definable functions
is preserved by operations including the sum, composition by a linear
operator, derivation and canonical projection ([34], [14]). Since $T$ is a
tame mapping, $I+M^{-1/2}TM^{-1/2}$ is tame as well as its graph. Here $I$ is
identity. Then the resolvent
$J_{M^{-1/2}TM^{-1/2}}=(I+M^{-1/2}TM^{-1/2})^{-1}$ which is defined by the
inverse of the same graph is tame ([18]) and single-valued. By the stability
of the sum and composition by linear operator, we obtain that $l(\alpha)$ is
tame. ∎

###  A.8 Proof of Proposition 4.3

###### Proof.

Let $p=J^{V}_{T}(z)$. Since resolvent operator is non-expansive,
$\|p\|=\|J^{V}_{T}(z)\|\leq\|z\|+\|J_{T}^{V}(0)\|$. By duality, optimal
$\alpha^{*}$ will satisfy $\alpha^{*}=u^{*}(p-z)$. Then,

$\begin{split}|\alpha^{*}|&=|u^{*}(p-z)|\\\
&\leq\|u\|(2\|z\|+\|J^{V}_{T}(0)\|)\,.\quad\qed\qed\\\ \end{split}$ (101)

###  A.9 Proof of Proposition 4.5

###### Proof.

$\mathcal{L}(\alpha)$ is as defined in Proposition 4.4. Substitute
$\alpha=\xi+U^{*}V^{-1}Bz$ in $\mathcal{L}(\alpha)$. We obtain
$\mathcal{J}(\xi)=\mathcal{L}(\alpha)$. It yields
$\alpha^{*}=\xi^{*}+U^{*}V^{-1}Bz$. In (24), we do the same substitution.

$\begin{split}\hat{z}&=J_{T}^{M}(\breve{z}\mp M^{-1}U\alpha^{*})\\\
&=J_{A}^{M}(\breve{z}\mp M^{-1}U\xi^{*}\mp M^{-1}UU^{*}V^{-1}Bz)\\\
&=J_{A}^{M}(z-M^{-1}MV^{-1}Bz\mp M^{-1}UU^{*}V^{-1}Bz\mp M^{-1}U\xi^{*})\\\
&=J_{A}^{M}(z-M^{-1}\underbrace{(M\pm UU^{*})}_{=V}V^{-1}Bz\mp
M^{-1}U\xi^{*})\\\ &=J_{A}^{M}(z-M^{-1}Bz\mp M^{-1}U\xi^{*})\,.\\\
\end{split}$ (102)

Since $\mathcal{L}(\alpha)$ is Lipschitz with constant $1+||M^{-1/2}U||^{2}$
and strongly monotone, the composition of translation preserves the same
properties. ∎

## REFERENCES

  * [1] B. Abbas, H. Attouch, and B. F. Svaiter, Newton-like dynamics and forward-backward methods for structured monotone inclusions in Hilbert spaces, Journal of Optimization Theory and Applications, 161 (2014), pp. 331–360.
  * [2] F. Alvarez, H. Attouch, J. Bolte, and P. Redont, A second-order gradient-like dissipative dynamical system with Hessian-driven damping.: Application to optimization and mechanics, Journal de mathématiques pures et appliquées, 81 (2002), pp. 747–779.
  * [3] H. Attouch and B. F. Svaiter, A continuous dynamical Newton-like approach to solving monotone inclusions, SIAM Journal on Control and Optimization, 49 (2011), pp. 574–598.
  * [4] H. Attouch and M. Théra, A general duality principle for the sum of two operators, Journal of Convex Analysis, 3 (1996), pp. 1–24.
  * [5] H. Bauschke and P. L. Combettes, Convex analysis and monotone operator theory in Hilbert Spaces, Springer, 2011.
  * [6] S. Becker and J. Fadili, A quasi-Newton proximal splitting method, Advances in Neural Information Processing Systems, (2012).
  * [7] S. Becker, J. Fadili, and P. Ochs, On quasi-Newton forward-backward splitting: Proximal calculus and convergence, SIAM Journal on Optimization, 29 (2019), pp. 2445–2481.
  * [8] J. Bolte, A. Daniilidis, and A. Lewis, Tame functions are semismooth, Mathematical Programming, 117 (2009), pp. 5–19.
  * [9] K. Bredies, D. Lorenz, et al., Mathematical image processing, Springer, 2018.
  * [10] A. Chambolle and T. Pock, A first-order primal-dual algorithm for convex problems with applications to imaging, Journal of mathematical imaging and vision, 40 (2011), pp. 120–145.
  * [11] P. Combettes, L. Condat, J. Pesquet, and B. Vu, A forward-backward view of some primal-dual optimization methods in image recovery, IEEE International Conference on Image Processing, (2014).
  * [12] P. L. Combettes and B. C. Vũ, Variable metric quasi-Fejér monotonicity, Nonlinear Analysis: Theory, Methods and Applications, 78 (2013), pp. 17–31.
  * [13] , Variable metric forward–backward splitting with applications to monotone inclusions in duality, Optimization, 63 (2014), pp. 1289–1318.
  * [14] M. Coste, An Introduction to o-minimal Geometry, Istituti editoriali e poligrafici internazionali Pisa, 2000.
  * [15] F. Facchinei and J. Pang, Finite-dimensional variational inequalities and complementarity problems, Springer, 2003.
  * [16] T. Goldstein, M. Li, X. Yuan, E. Esser, and R. Baraniuk, Adaptive primal-dual hybrid gradient methods for saddle-point problems, arXiv preprint arXiv:1305.0546, (2013).
  * [17] B. He and X. Yuan, Convergence analysis of primal-dual algorithms for a saddle-point problem: From contraction perspective, SIAM Journal on Imaging Sciences, 5 (2012), pp. 119–149.
  * [18] A. D. Ioffe, An invitation to tame optimization, SIAM Journal on Optimization, 19 (2009), pp. 1894–1917.
  * [19] C. Kanzow and T. Lechner, Globalized inexact proximal Newton-type methods for nonconvex composite functions, Computational Optimization and Applications, 78 (2021), pp. 377–410.
  * [20] C. Kanzow and T. Lechner, Efficient regularized proximal quasi-Newton methods for large-scale nonconvex composite optimization problems, tech. rep., University of Würzburg, Institute of Mathematics, January 2022.
  * [21] S. Karimi and S. Vavasis, IMRO: A proximal quasi-Newton method for solving $\ell\\_1$-regularized least squares problems, SIAM Journal on Optimization, 27 (2017), pp. 583–615.
  * [22] P. D. Khanh, B. Mordukhovich, V. T. Phat, and D. B. Tran, Generalized damped Newton algorithms in nonsmooth optimization with applications to lasso problems, arXiv preprint arXiv:2101.10555, (2021).
  * [23] J. D. Lee, Y. Sun, and M. A. Saunders, Proximal Newton-type methods for minimizing composite functions, SIAM Journal on Optimization, 24 (2014), pp. 1420–1443.
  * [24] C. Liu and L. Luo, Quasi-Newton methods for saddle point problems, arXiv preprint arXiv:2111.02708, (2021).
  * [25] Y. Liu, Y. Xu, and W. Yin, Acceleration of primal–dual methods by preconditioning and simple subproblem procedures, Journal of Scientific Computing, 86 (2021), pp. 1–34.
  * [26] D. A. Lorenz and T. Pock, An inertial forward-backward algorithm for monotone inclusions, Journal of Mathematical Imaging and Vision, 51 (2015), pp. 311–325.
  * [27] P. Patrinos, L. Stella, and A. Bemporad, Forward-backward truncated Newton methods for convex composite optimization, arXiv preprint arXiv:1402.6655, (2014).
  * [28] T. Pock and A. Chambolle, Diagonal preconditioning for first order primal-dual algorithms in convex optimization, IEEE International Conference on Computer Vision, (2011).
  * [29] B. Polyak, Introduction to optimization, Optimization Software, 1987\.
  * [30] A. Rodomanov and Y. Nesterov, Greedy quasi-Newton methods with explicit superlinear convergence, SIAM Journal on Optimization, 31 (2021), pp. 785–811.
  * [31] M. Schmidt, E. Berg, M. Friedlander, and K. Murphy, Optimizing costly functions with simple constraints: A limited-memory projected quasi-Newton algorithm, in AISTATS, 2009.
  * [32] M. Schmidt, D. Kim, and S. Sra, Projected Newton-type methods in machine learning, Optimization for Machine Learning, (2012).
  * [33] L. Stella, A. Themelis, and P. Patrinos, Forward–backward quasi-Newton methods for nonsmooth optimization problems, Computational Optimization and Applications, 67 (2017), pp. 443–487.
  * [34] L. Van den Dries and C. Miller, Geometric categories and o-minimal structures, Duke Mathematical Journal, 84 (1996), pp. 497–540.
  * [35] S. Wright and J. Nocedal, Numerical optimization, Springer Science, (1999).