Stochastic Primal-Dual Hybrid Gradient Algorithm with Arbitrary Sampling and Imaging Applications

TL;DR

This paper introduces a stochastic primal-dual hybrid gradient algorithm with arbitrary sampling, extending previous deterministic methods and demonstrating superior performance in imaging applications.

Contribution

It presents a novel stochastic extension of the primal-dual hybrid gradient algorithm that handles arbitrary sampling and improves imaging task results.

Findings

Outperforms deterministic variants in imaging tasks

Handles general convex-concave saddle point problems

Supports arbitrary sampling of dual variables

Abstract

We propose a stochastic extension of the primal-dual hybrid gradient algorithm studied by Chambolle and Pock in 2011 to solve saddle point problems that are separable in the dual variable. The analysis is carried out for general convex-concave saddle point problems and problems that are either partially smooth / strongly convex or fully smooth / strongly convex. We perform the analysis for arbitrary samplings of dual variables, and obtain known deterministic results as a special case. Several variants of our stochastic method significantly outperform the deterministic variant on a variety of imaging tasks.

Equations

$$\underset{x\in \mathbb{X}}{\min }\left\{\Phi (x):=\sum _{i=1}^n{f}_i({\mathbf{A}}_{i}x)+g(x)\right\}.$$

$$\underset{x\in \mathbb{X}}{\mathrm{minp}}\underset{y\in \mathbb{Y}}{\sup }\left\{\Psi (x,y):=\sum _{i=1}^n⟨{\mathbf{A}}_{i}x,y_i⟩-f_i^{\ast }(y_i)+g(x)\right\}.$$

$${\mathbf{A}}_i{x}^{♯}\in \partial f_i^{\ast }(y_i^{♯})i=1,\dots ,n,-{\mathbf{A}}^{\ast }{y}^{♯}\in \partial g(x^{♯}).$$

$$x^{(k+1)}$$

$$y^{(k+1)}$$

$${\overline{y}}^{(k+1)}$$

$${\mathrm{prox}}_f^{\tau }\left(y\right):=\arg \underset{x\in \mathbb{X}}{\min }\left\{\frac{1}{2}\Vert x-y{\Vert }_{{\tau }^{-1}}^2+f(x)\right\}$$

$$x^{(k+1)}$$

$$y_i^{(k+1)}$$

$${\overline{y}}^{(k+1)}$$

$${\mathbf{A}}^{\ast }{\overline{y}}^{(k+1)}$$

$${\mathbb{E}}_{\mathbb{S}}{\Vert \sum _{i\in \mathbb{S}}{\mathbf{C}}_i^{\ast }{z}_i\Vert }^2\le \sum _{i=1}^n{p}_i{v}_i\Vert z_i{\Vert }^2.$$

$$2{\mathbb{E}}_{\mathbb{S}}⟨\mathbf{Q}\mathbf{A}x,y^{+}-y⟩\ge -{\mathbb{E}}_{\mathbb{S}}\left\{\frac{1}{c}\Vert x{\Vert }_{{\mathbf{T}}^{-1}}^2+c\underset{i}{\max }\frac{v_i}{p_i}\Vert y^{+}-y{\Vert }_{\mathbf{Q}{\mathbf{S}}^{-1}}^2\right\}.$$

$$D_f^q(x,y):=f(x)-f(y)-⟨q,x-y⟩.$$

$$G_{{\mathbb{B}}_1\times {\mathbb{B}}_2}(x,y):=\underset{\tilde{y}\in {\mathbb{B}}_2}{\sup }\Psi (x,\tilde{y})-\underset{\tilde{x}\in {\mathbb{B}}_1}{\inf }\Psi (\tilde{x},y).$$

$$G_{\mathbb{B}}(w)\ge \Psi (x,y^{♯})-\Psi (x^{♯},y)=D_g^{-{\mathbf{A}}^{\ast }{y}^{♯}}(x,x^{♯})+D_{f^{\ast }}^{\mathbf{A}{x}^{♯}}(y,y^{♯})=D_h^q(w,w^{♯})\ge 0$$

$${\mathcal{F}}_i(y_i|\tilde{x},{\tilde{y}}_i):=f_i^{\ast }(y_i)-f_i^{\ast }({\tilde{y}}_i)-⟨{\mathbf{A}}_i\tilde{x},y_i-{\tilde{y}}_i⟩$$

$${\mathcal{F}}^p(y|\tilde{w}):=\sum _{i=1}^n\left(\frac{1}{p_i}-1\right){\mathcal{F}}_i(y_i|\tilde{x},{\tilde{y}}_i)$$

$$G_{\mathbb{B}}(w)=\underset{\tilde{w}\in \mathbb{B}}{\sup }\mathcal{G}(x|\tilde{w})+\mathcal{F}(y|\tilde{w})=\underset{\tilde{w}\in \mathbb{B}}{\sup }\mathcal{H}(w|\tilde{w}).$$

$$v_i<p_{i}i=1,\dots ,n.$$

$$D_h^q(w^{(k)},w^{♯})\to 0\text{a.s.}$$

$$\mathbb{E}{G}_{\mathbb{B}}(w_{(K)})\le \frac{C_{\mathbb{B}}}{K}$$

$$C_{\mathbb{B}}:=\underset{x\in {\mathbb{B}}_1}{\sup }\frac{1}{2}\Vert x^{(0)}-x{\Vert }_{{\mathbf{T}}^{-1}}^2+\underset{y\in {\mathbb{B}}_2}{\sup }\frac{1}{2}\Vert y^{(0)}-y{\Vert }_{\mathbf{Q}{\mathbf{S}}^{-1}}^2+\underset{w\in \mathbb{B}}{\sup }{\mathcal{F}}^p(y^{(0)}|w).$$

$$x^{(k+1)}$$

$${\hat{y}}_i^{(k+1)}$$

$$\Vert x^{(k)}-x{\Vert }_{{\mathbf{T}}_{(k)}^{-1}}^2+\Vert y^{(k)}-y{\Vert }_{{\mathbf{S}}_{(k)}^{-1}}^2$$

$$\ge \Vert x^{(k+1)}-x{\Vert }_{{\mathbf{T}}_{(k)}^{-1}+{\mu }_g\mathbf{I}}^2+\Vert {\hat{y}}^{(k+1)}-y{\Vert }_{{\mathbf{S}}_{(k)}^{-1}+\mathbf{M}}^2$$

$$+2\left(\mathcal{G}(x^{(k+1)}|w)+\mathcal{F}({\hat{y}}^{(k+1)}|w)\right)-2⟨\mathbf{A}(x^{(k+1)}-x),{\hat{y}}^{(k+1)}-{\overline{y}}^{(k)}⟩$$

$$+\Vert x^{(k+1)}-x^{(k)}{\Vert }_{{\mathbf{T}}_{(k)}^{-1}}^2+\Vert {\hat{y}}^{(k+1)}-y^{(k)}{\Vert }_{{\mathbf{S}}_{(k)}^{-1}}^2.$$

$$\phi ({\hat{y}}_i^{(k+1)})$$

$$=\left(\frac{1}{p_i}-1\right){\mathbb{E}}^{(k+1)}\phi (y_i^{(k+1)})-\left(\frac{1}{p_i}-1\right)\phi (y_i^{(k)})+{\mathbb{E}}^{(k+1)}\phi (y_i^{(k+1)}).$$

$$\mathcal{F}({\hat{y}}^{(k+1)}|w)$$

$$\Vert {\hat{y}}^{(k+1)}-\cdot {\Vert }_{\mathbf{B}}^2$$

$${\hat{y}}^{(k+1)}$$

$$\Vert x^{(k)}-x{\Vert }_{{\mathbf{T}}^{-1}}^2+\Vert y^{(k)}-y{\Vert }_{\mathbf{Q}{\mathbf{S}}^{-1}}^2+2{\mathcal{F}}^p(y^{(k)}|w)$$

Full text

{AMS}

65D18, 65K10, 74S60, 90C25, 90C15, 92C55, 94A08

1 Introduction

Many modern problems in a variety of disciplines (imaging, machine learning, statistics, etc.) can be formulated as convex optimization problems. Instead of solving the optimization problems directly, it is often advantageous to reformulate the problem as a saddle point problem. A very popular algorithm to solve such saddle point problems is the primal-dual hybrid gradient (PDHG)111We follow the terminology of [13] and call the algorithm simply PDHG. It corresponds to PDHGMu and PDHGMp in [20]. algorithm [36, 20, 12, 35, 13, 14]. It has been used to solve a vast amount of state-of-the-art problems—to name a few examples in imaging: image denoising with the structure tensor [21], total generalized variation denoising [10], dynamic regularization [6], multi-modal medical imaging [26], multi-spectral medical imaging [42], computation of non-linear eigenfunctions [25], regularization with directional total generalized variation [28]. Its popularity stems from two facts: First, it is very simple and therefore easy to implement. Second, it involves only simple operations like matrix-vector multiplications and evaluations of proximal operators which are for many problems of interest simple and in closed-form or easy to compute iteratively, cf. e.g. [32]. However, for large problems that are encountered in many real world applications, even these simple operations might be still too costly to perform too often.

We propose a stochastic extension of the PDHG for saddle point problems that are separable in the dual variable (cf. e.g. [17, 51, 53, 33]) where not all but only a few of these operations are performed in every iteration. Moreover, as in incremental optimization algorithms [46, 30, 9, 8, 7, 44, 18] over the course of the iterations we continuously build up information from previous iterations which reduces variance and thereby negative effects of stochasticity. Non-uniform samplings [39, 37, 51, 38, 2] have been proven very efficient for stochastic optimization. In this work we use the expected separable overapproximation framework of [37, 38, 40] to prove all statements for all non-trivial and iteration-independent samplings.

Related Work

The proposed algorithm can be seen as a generalization of the algorithm of [17, 53, 51] to arbitrary blocks and a much wider class of samplings. Moreover, in contrast to their results, our results generalize the deterministic case considered in [36, 12, 35, 14]. Fercoq and Bianchi [22] proposed a stochastic primal-dual algorithm with explicit gradient steps that allows for larger step sizes by averaging over previous iterates, however, this comes at the cost of prohibitively large memory requirements. Similar memory issues are encountered by a primal-dual algorithm of [3]. It is related to forward-backward splitting [29] and averaged gradient descent [9, 19] and therefore suffers the same memory issues as the averaged gradient descent. Moreover, Valkonen proposed a stochastic primal-dual algorithm that can exploit partial strong convexity of the saddle point functional [47]. Randomized versions of the alternating direction method of multipliers are discussed for instance in [52, 24]. In contrast to other works on stochastic primal-dual algorithms [34, 50], our analysis is not based on Fejér monotonicity [15]. We therefore do not prove almost sure convergence of the sequence but prove a variety of convergence rates depending on strong convexity assumptions instead.

As a word of warning, our contribution should not be mistaken by other \enquotestochastic primal-dual algorithms, where errors in the computation of matrix-vector products and evaluation of proximal operators are modeled by random variables, cf. e.g. [34, 15, 43]. In our work we deliberately choose to compute only a subset of a whole iteration to save computational cost. These two notations are related but are certainly not the same.

1.1 Contributions

We briefly mention the main contributions of our work.

Generalization of Deterministic Case

The proposed stochastic algorithm is a direct generalization of the deterministic setting [36, 12, 35, 13, 14]. In the degenerate case where in every iteration all computations are performed, our algorithm coincides with the original deterministic algorithm. Moreover, the same holds true for our analysis of the stochastic algorithm where we recover almost all deterministic statements [12, 35] in this degenerate case. Therefore, the theorems for both the deterministic and the stochastic case can be combined by a single proof.

Better Rates

Our analysis extends the simple setting of [51] such that the strong convexity assumptions and the sampling do not have to be uniform. Even in the special case of uniform strong convexity and uniform sampling, the proven convergence rates are slightly better than the ones proven in [51].

Arbitrary Sampling

The proposed algorithm is guaranteed to converge under a very general class of samplings [37, 38, 40] and thereby generalizes also the algorithm of [51] which has only been analyzed for two specific samplings. As long as the sampling is independent and identically distributed over the iterations and all computations have non-zero probability to be carried out, the theory holds and the algorithm will converge with the proven convergence rates.

Acceleration

We propose an acceleration of the stochastic primal-dual algorithm which accelerates the convergence from $\mathcal{O}(1/K)$ to $\mathcal{O}(1/K^{2})$ if parts of the saddle point functional are strongly convex and thereby results in a significantly faster algorithm.

Scaling Invariance

In the strongly convex case, we propose parameters for several serial samplings (uniform, importance, optimal), all based on the condition numbers of the problem and thereby independent of scaling.

2 General Problem

Let ${\mathbb{X}},{\mathbb{Y}}_{i},i=1,\ldots,n$ be real Hilbert spaces of any dimension and define the product space ${\mathbb{Y}}\coloneqq\prod_{i=1}^{n}{\mathbb{Y}}_{i}$. For $y\in{\mathbb{Y}}$, we shall write $y=(y_{1},y_{2},\dots,y_{n})$, where $y_{i}\in{\mathbb{Y}}_{i}$. Further, we consider the natural inner product on the product space ${\mathbb{Y}}$ given by $\langle y,z\rangle=\sum_{i=1}^{n}\langle y_{i},z_{i}\rangle$, where $y_{i},z_{i}\in{\mathbb{Y}}_{i}$. This inner product induces the norm $\|y\|^{2}=\sum_{i=1}^{n}\|y_{i}\|^{2}$. Moreover, for simplicity we will consider the space ${\mathbb{W}}\coloneqq{\mathbb{X}}\times{\mathbb{Y}}$ that combines both primal and dual variables.

Let ${\bf A}:{\mathbb{X}}\to{\mathbb{Y}}$ be a bounded linear operator. Due to the product space nature of ${\mathbb{Y}}$, we have $({\bf A}x)_{i}={\bf A}_{i}x$, where ${\bf A}_{i}:{\mathbb{X}}\to{\mathbb{Y}}_{i}$ are linear operators. The adjoint of ${\bf A}$ is given by ${\bf A}^{\ast}y=\sum_{i=1}^{n}{\bf A}^{\ast}_{i}y_{i}$. Moreover, let $f:{\mathbb{Y}}\to{\mathbb{R}}_{\infty}\coloneqq{\mathbb{R}}\cup\{+\infty\}$ and $g:{\mathbb{X}}\to{\mathbb{R}}_{\infty}$ be convex functions. In particular, we assume that $f$ is separable, i.e. $f(y)=\sum_{i=1}^{n}f_{i}(y_{i})$.

Given the setup described above, we consider the optimization problem

[TABLE]

Instead of solving Eq. 1 directly, it is often desirable to reformulate the problem as a saddle point problem with the help of the Fenchel conjugate. If $f$ is proper, convex, and lower semi-continuous, then $f(y)=f^{\ast\ast}(y)=\sup_{z\in{\mathbb{Y}}}\langle z,y\rangle-f^{\ast}(z)$ where $f^{\ast}:{\mathbb{Y}}\to{\mathbb{R}}\cup\{-\infty,+\infty\}$, $f^{\ast}(y)=\sum_{i=1}^{n}f_{i}^{\ast}(y_{i})$ is the Fenchel conjugate of $f$ (and $f^{\ast\ast}$ its biconjugate defined as the conjugate of the conjugate). Then solving Eq. 1 is equivalent to finding the primal part $x$ of a solution to the saddle point problem (called a saddle point)

[TABLE]

We will assume that the saddle point problem Eq. 2 has a solution. For conditions for existence and uniqueness, we refer the reader to [4]. A saddle point $w^{\sharp}=(x^{\sharp},y^{\sharp})=(x^{\sharp},y^{\sharp}_{1},\ldots,y^{\sharp}_{n})\in{\mathbb{W}}$ satisfies the optimality conditions

[TABLE]

An important notion in this work is strong convexity. A functional $g$ is called $\mu_{g}$-convex if $g-\frac{\mu_{g}}{2}\|\cdot\|^{2}$ is convex. In general, we assume that $g$ is $\mu_{g}$-convex, $f_{i}^{\ast}$ is $\mu_{i}$-convex with nonnegative strong convexity parameters $\mu_{g},\mu_{i}\geq 0$. The convergence results in this contribution cover three different cases of regularity: i) no strong convexity $\mu_{g},\mu_{i}=0$, ii) semi strong convexity $\mu_{g}>0$ or $\mu_{i}>0$ and iii) full strong convexity $\mu_{g},\mu_{i}>0$. For notational convenience we make use of the operator ${\bf M}\coloneqq\operatorname{diag}(\mu_{1}{\bf I},\ldots,\mu_{n}{\bf I})$.

A very popular algorithm to solve the saddle point problem Eq. 2 is the Primal-Dual Hybrid Gradient (PDHG) algorithm [36, 20, 12, 35, 13, 14]. It reads (with extrapolation on $y$)

[TABLE]

where the proximal operator (or proximity / resolvent operator) is defined as

[TABLE]

and the weighted norm by $\|x\|_{\tau^{-1}}^{2}=\langle\tau^{-1}x,x\rangle$. Its convergence is guaranteed if the step size parameters $\sigma,\tau$ are positive and satisfy $\sigma\tau\|{\bf A}\|^{2}<1,\theta=1$ [12]. Note that the definition of the proximal operator is well-defined for an operator-valued step size $\tau$. In the case of a separable function $f$ and with operator-valued step sizes the PDHG algorithm takes the form

[TABLE]

Here the step size parameters ${\bf S}=\operatorname{diag}({\bf S}_{1},\ldots,{\bf S}_{n})$ (a block diagonal operator), ${\bf S}_{1},\ldots,{\bf S}_{n}$ and ${\bf T}$ are symmetric and positive definite. The algorithm is guaranteed to converge if $\|{\bf S}^{1/2}{\bf A}{\bf T}^{1/2}\|<1$ and $\theta=1$ [35].

3 Algorithm

In this work we extend the PDHG algorithm to a stochastic setting where in each iteration we update a random subset ${\mathbb{S}}$ of the dual variables Eq. 3b. This subset is sampled in an i.i.d. fashion from a fixed but otherwise arbitrary distribution, whence the name \enquotearbitrary sampling. In order to guarantee convergence, it is necessary to assume that the sampling is \enquoteproper [41, 38]. A sampling is proper if for each dual variable $i$ we have $i\in{\mathbb{S}}$ with a positive probability $p_{i}>0$. Examples of proper samplings include the full sampling where ${\mathbb{S}}=\{1,\ldots,n\}$ with probability 1 and serial sampling where ${\mathbb{S}}=\{i\}$ is chosen with probability $p_{i}$. It is important to note that also other samplings are admissible. For instance for $n=3$, consider the sampling that selects ${\mathbb{S}}=\{1,2\}$ with probability $1/3$ and $\mathbb{S}=\{2,3\}$ with probability $2/3$. Then the probabilities for the three blocks are $p_{1}=1/3$, $p_{2}=1$ and $p_{3}=2/3$ which makes it a proper sampling. However, if only ${\mathbb{S}}=\{1,2\}$ is chosen with probability 1, then this sampling is not proper as the probability for the third block is zero: $p_{3}=0$.

The algorithm we propose is formalized as Algorithm 1. As in the original PDHG, the step size parameters ${\bf T},{\bf S}_{i}$ have to be self-adjoint and positive definite operators for the updates to be well-defined. The extrapolation is performed with a scalar $\theta>0$ and an operator ${\bf Q}\coloneqq\operatorname{diag}(p_{1}^{-1}{\bf I},\ldots,p_{n}^{-1}{\bf I})$ of probabilities $p_{i}$ that an index is selected in each iteration.

Remark 1**.**

*Both, the primal and dual iterates $x^{(k)}$ and $y^{(k)}$ are random variables but only the dual iterate $y^{(k)}$ depends on the sampling ${\mathbb{S}}^{(k)}$. However, $x^{(k)}$ depends of course on all previous samplings ${\mathbb{S}}^{(i)},i<k$. *

Remark 2**.**

Due to the sampling each iteration requires both ${\bf A}_{i}$ and ${\bf A}^{\ast}_{i}$ to be evaluated only for each selected index $i\in{\mathbb{S}}^{(k+1)}$. To see this, note that

[TABLE]

*where ${\bf A}^{\ast}y^{(k)}$ can be stored from the previous iteration (needs the same memory as the primal variable $x$) and the operators ${\bf A}^{\ast}_{i}$ are evaluated only for $i\in{\mathbb{S}}^{(k+1)}$. *

4 General Convex Case

We first analyze the convergence of Algorithm 1 in the general convex case without making use of any strong convexity or smoothness assumptions. In order to analyze the convergence for the large class of samplings described in the previous section we make use of the expected separable overapproximation (ESO) inequality [38].

Definition 4.1** (Expected Separable Overapproximation (ESO) [38]).**

Let ${\mathbb{S}}\subset\{1,\ldots,n\}$ be a random set and $p_{i}\coloneqq\mathbb{P}(i\in{\mathbb{S}})$ the probability that $i$ is in ${\mathbb{S}}$. We say that $\{v_{i}\}\subset{\mathbb{R}}^{n}$ fulfill the ESO inequality (depending on ${\bf C}=[{\bf C}_{1};\ldots;{\bf C}_{n}],{\bf C}^{\ast}z=\sum_{i=1}^{n}C_{i}^{\ast}z_{i}$ and ${\mathbb{S}}$) if for all $z\in{\mathbb{Y}}$ it holds that

[TABLE]

*Such parameters $\{v_{i}\}$ are called ESO parameters. *

Remark 4.2**.**

*Note that for any bounded linear operator ${\bf C}$ such parameters always exist but are obviously not unique. For the efficiency of the algorithm it is desirable to find ESO parameters such that Eq. 4 is as tight as possible; i.e., we want the parameters $\{v_{i}\}$ be small. As we shall see, the ESO parameters influence the choice of the extrapolation parameter $\theta$ in the strongly convex case. *

The ESO inequality was first proposed by Richtárik and Takáč [41] to study parallel coordinate descent methods in the context of uniform samplings, which are samplings for which $p_{i}=p_{j}$ for all $i,j$. Improved bounds for ESO parameters were obtained in [23] and used in the context of accelerated coordinate descent. Qu et al. [38] perform an in-depth study of ESO parameters. The ESO inequality is also critical in the study mini-batch stochastic gradient descent with [27] or without [45] variance reduction.

We will frequently need to estimate the expected value of inner products which we will do by means of ESO parameters. Recall that we defined weighted norms as $\|x\|_{{\bf T}^{-1}}^{2}\coloneqq\langle{\bf T}^{-1}x,x\rangle$. The proof of this lemma can be found in the appendix.

Lemma 4.3**.**

Let ${\mathbb{S}}\subset\{1,\ldots,n\}$ be a random set and $y^{+}_{i}=\hat{y}_{i}$ if $i\in{\mathbb{S}}$ and $y_{i}$ otherwise. Moreover, let $\{v_{i}\}$ be some ESO parameters of ${\bf S}^{1/2}{\bf A}{\bf T}^{1/2}$ and $p_{i}={\mathbb{P}}(i\in{\mathbb{S}})$. Then for any $x\in{\mathbb{X}}$ and $c>0$

[TABLE]

Example 4.4** (Full Sampling).**

*Let ${\mathbb{S}}=\{1,\ldots,n\}$ with probability 1 such that $p_{i}=\mathbb{P}(i\in{\mathbb{S}})=1$. Then $v_{i}=\|{\bf S}^{1/2}{\bf A}{\bf T}^{1/2}\|^{2}$ are ESO parameters of ${\bf S}^{1/2}{\bf A}{\bf T}^{1/2}$. Thus, the deterministic condition on convergence, $\|{\bf S}^{1/2}{\bf A}{\bf T}^{1/2}\|^{2}<1$, implies a bound on some ESO parameters $v_{i}<p_{i}$. *

Example 4.5** (Serial Sampling).**

*Let ${\mathbb{S}}=\{i\}$ be chosen with probability $p_{i}>0$. Then $v_{i}=\|{\bf S}_{i}^{1/2}{\bf A}_{i}{\bf T}^{1/2}\|^{2}$ are ESO parameters of ${\bf S}^{1/2}{\bf A}{\bf T}^{1/2}$. *

The analysis for the general convex case will use the notation of Bregman distance which is defined for any function $f:{\mathbb{X}}\to{\mathbb{R}}_{\infty}$, $x,y\in{\mathbb{X}}$ and $q\in\partial f(y)$ in the subdifferential of $f$ at $y$ as

[TABLE]

Next to Bregman distances, one can measure optimality by the partial primal-dual gap. Let ${\mathbb{B}}_{1}\times{\mathbb{B}}_{2}\subset{\mathbb{W}}={\mathbb{X}}\times{\mathbb{Y}}$, then we define the partial primal-dual gap as

[TABLE]

It is convenient to define ${\mathbb{B}}\coloneqq{\mathbb{B}}_{1}\times{\mathbb{B}}_{2}\subset{\mathbb{W}}$ and to denote the gap as $G_{{\mathbb{B}}}(w)\coloneqq G_{{\mathbb{B}}_{1}\times{\mathbb{B}}_{2}}(x,y).$ Note that if ${\mathbb{B}}$ contains a saddle point $w^{\sharp}=(x^{\sharp},y^{\sharp})$, then we have that

[TABLE]

where the first equality is obtained by adding a zero and we used $h(w)\coloneqq g(x)+f^{\ast}(y)$ and $q\coloneqq(-{\bf A}^{\ast}y^{\sharp},{\bf A}x^{\sharp})\in\partial h(w^{\sharp})$ for the last equality. The non-negativity stems from the fact that Bregman distances of convex functionals are non-negative and $h$ is convex indeed.

We will make frequent use of the following \enquotedistance functions

[TABLE]

and ${\mathcal{F}}(y|\tilde{w})\coloneqq\sum_{i=1}^{n}\mathcal{F}_{i}(y_{i}|\tilde{x},\tilde{y}_{i})\,.$ Note that these are strongly related to Bregman distances; if $w^{\sharp}$ is a saddle point, then $\mathcal{F}(y|w^{\sharp})=D_{f^{\ast}}^{{\bf A}x^{\sharp}}(y,y^{\sharp})$ is the Bregman distance of $f^{\ast}$ between $y$ and $y^{\sharp}$. Similarly, we make use of the weighted distance

[TABLE]

and the distance for the primal functional ${\mathcal{G}}(x|\tilde{w})\coloneqq g(x)-g(\tilde{x})-\langle-{\bf A}^{\ast}\tilde{y},x-\tilde{x}\rangle.$ We note that these distances are also related to the partial primal-dual gap as

[TABLE]

Theorem 4.6**.**

Let $g,f_{i},i=1,\ldots,n$ be convex, $\theta=1$ and ${\bf T},{\bf S}$ be chosen so that there exist ESO parameters $\{v_{i}\}$ of ${\bf S}^{1/2}{\bf A}{\bf T}^{1/2}$ with

[TABLE]

Then, the Bregman distance between iterates of Algorithm 1 $w^{(k)}=(x^{(k)},y^{(k)})\in{\mathbb{W}}$ and any saddle point $w^{\sharp}\in{\mathbb{W}}$ converges to zero almost surely,

[TABLE]

Moreover, the ergodic sequence $w_{(K)}\coloneqq\frac{1}{K}\sum_{k=1}^{K}w^{(k)}$ converges with rate $1/K$ in an expected partial primal-dual gap sense, i.e. for any set ${\mathbb{B}}\coloneqq{\mathbb{B}}_{1}\times{\mathbb{B}}_{2}\subset{\mathbb{W}}$ it holds

[TABLE]

where the constant is given by

[TABLE]

*The same rate holds for the expected Bregman distance, $\mathbb{E}D_{h}^{q}(w_{(K)},w^{\sharp})\leq C_{\{w^{\sharp}\}}/K$. *

Remark 4.7**.**

*The convergence Eq. 6 in Bregman distance implies convergence in norm as soon as $h$ is strictly convex. If $h$ is not strictly convex, then the convergence has to be seen in a more generalized sense. For example, if $h$ is a $\ell^{1}$-norm (and thus not strictly convex), then the Bregman distance between $w^{(k)}$ and $w^{\sharp}$ is zero if and only if they have the same support and sign. Thus, the convergence statement is related to the support and sign of $w^{\sharp}$. In the extreme case $h\equiv 0$, then $D_{h}^{q}(\cdot,w^{\sharp})\equiv 0$ and the convergence statement has no meaning. *

The proof of this theorem utilizes a standard inequality for which we provide the proof in the appendix for completeness.

Lemma 4.8**.**

Consider the deterministic updates

[TABLE]

with iteration varying step sizes ${\bf T}_{(k)}$ and ${\bf S}_{(k)}=\operatorname{diag}([{\bf S}_{(k)}]_{1},\ldots,[{\bf S}_{(k)}]_{n})$. Then for any $(x,y)\in{\mathbb{W}}$ it holds that

[TABLE]

Proof 4.9** (Proof of Theorem 4.6).**

The result of Lemma 4.8 (with constant step sizes) has to be adapted to the stochastic setting as the dual iterate is updated only with a certain probability. First, a trivial observation is that for any mapping $\varphi$ it holds that

[TABLE]

Thus, for the generalized distance of $f^{\ast}$ we arrive at

[TABLE]

and for any block diagonal matrix ${\bf B}=\operatorname{diag}({\bf B}_{1},\ldots,{\bf B}_{n})$

[TABLE]

Using Eqs. 10, 11 and 12, we can rewrite the estimate of Lemma 4.8 as

[TABLE]

where we have used the identity

[TABLE]

to simplify the expression. With the extrapolation $\overline{y}^{(k)}=y^{(k)}+{\bf Q}(y^{(k)}-y^{(k-1)})$, the inner product term can be reformulated as

[TABLE]

and with Lemma 4.3 and $\gamma^{2}\coloneqq\max_{i}v_{i}/p_{i}$ it holds that

[TABLE]

Taking expectations with respect to ${\mathbb{S}}^{1},\ldots,{\mathbb{S}}^{K}$ (denoting this by $\mathbb{E}$) on Eq. 13, using the estimates Eqs. 15 and 16 and denoting

[TABLE]

leads with $\gamma<1$ (follows directly from Eq. 5) to

[TABLE]

Summing Eq. 17 over $k=0,\ldots,K-1$ (note that $y^{(-1)}=y^{(0)}$) and using the estimate (follows directly from Lemma 4.3)

[TABLE]

yields

[TABLE]

All assertions of the theorem follow from inequality Eq. 18. Inserting a saddle point $w=w^{\sharp}$ and taking the limit $K\to\infty$, it follows from Eq. 18 that $\mathbb{E}\sum_{k=1}^{\infty}D_{h}^{q}(w^{(k)},w^{\sharp})<\infty$ which implies almost surely $\sum_{k=1}^{\infty}D_{h}^{q}(w^{(k)},w^{\sharp})<\infty$ and thus Eq. 6.

To see Eq. 7, note first that

[TABLE]

and $\frac{\Delta^{(0)}}{2}-{\mathcal{F}}^{p}(y^{(K)}|w)\leq C_{{\mathbb{B}}}$ if $w\in{\mathbb{B}}$ with $C_{{\mathbb{B}}}$ as defined in Eq. 8. Moreover, the generalized distance ${\mathcal{H}}(\cdot|w)$ is convex, thus, dividing Eq. 18 by $K$ yields

[TABLE]

*for any $w\in{\mathbb{B}}$. Taking the supremum over $w\in{\mathbb{B}}$ yields Eq. 7. Noting that $D_{h}^{q}(w,w^{\sharp})=G_{\{w^{\sharp}\}}(w)$ completes the proof. *

5 Semi-Strongly Convex Case

In this section we propose two algorithms that converge as $\mathcal{O}(1/K^{2})$ if either $f_{i}^{\ast}$ or $g$ is strongly convex. For simplicity we restrict ourselves from now on to scalar-valued step sizes, i.e. ${\bf T}=\tau{\bf I}$ and ${\bf S}_{i}=\sigma_{i}{\bf I}$. However, large parts of what follows holds true for operator-valued step sizes, too.

Theorem 5.1** (Dual Strong Convexity).**

Let $f_{i}^{\ast}$ be strongly convex with constants $\mu_{i}>0,i=1,\ldots,n$. Consider Algorithm 3 and let the initial step sizes $\tilde{\sigma}_{(0)},\tau_{(0)}$ be chosen such that

[TABLE]

and for the ESO parameters $\{v_{i}\}$ of ${\bf S}_{(0)}^{1/2}{\bf A}\tau_{(0)}^{1/2}$ it holds that

[TABLE]

with $[{\bf S}_{(k)}]_{i}=\sigma^{(k)}_{i}{\bf I}$ and

[TABLE]

Then there exists $\tilde{K}\in{\mathbb{N}}$ such that for all $K\geq\tilde{K}$ it holds

[TABLE]

*where the metric on ${\mathbb{Y}}$ is defined by ${\bf Y}_{(k)}\coloneqq{\bf Q}{\bf S}_{(k)}^{-1}+2{\bf M}({\bf Q}-{\bf I})$. *

Remark 5.2**.**

*As already noted in [12], $\tilde{K}$ is usually fairly small so that the estimate in Theorem 5.1 has practical relevance. *

Remark 5.3**.**

For serial sampling the condition on the ESO parameters Eq. 20 is equivalent to

[TABLE]

*In particular, it implies condition Eq. 19 on $\tilde{\sigma}_{(0)}$. *

Remark 5.4**.**

*The convergence of Algorithm 2 with acceleration on the primal variable is similar to the deterministic case, cf. Appendix C.2 of [13], and omitted here for brevity. It converges with rate $\mathcal{O}(1/K^{2})$ if the ESO parameters satisfy $v_{i}<p_{i}$. *

This theorem requires an estimate on the expected contraction similar to the proof of Theorem 4.6 and shown in the appendix.

Lemma 5.5**.**

Let $x^{(k+1)},\hat{y}^{(k+1)}$ be defined as in Lemma 4.8 and $y^{(k+1)}_{i}=\hat{y}^{(k+1)}_{i}$ with probability $p_{i}>0$ and unchanged else. Moreover, let

[TABLE]

and $\{v_{i}\}$ be some ESO parameters of ${\bf S}_{(k)}^{1/2}{\bf A}\tau_{(k)}^{1/2}$. Then with $\gamma^{2}=\max_{i}\frac{v_{i}}{p_{i}}$ it holds

[TABLE]

Proof 5.6** (Proof of Theorem 5.1).**

The update on the step sizes in Algorithm 3 imply that

[TABLE]

for all $i=1,\ldots,n$ and therefore

[TABLE]

To see Eq. 23, the auxiliary sequence $\tilde{\sigma}_{(k)}$ satisfies

[TABLE]

such that Eq. 23 is satisfied as soon as

[TABLE]

Note that the transformation from $\tilde{\sigma}_{(k)}$ to $\sigma^{(k)}_{i}$ is well-defined if $\tilde{\sigma}_{(k)}<\min_{i}\frac{p_{i}}{2(1-p_{i})}$ which is the case as $\tilde{\sigma}_{(k)}$ is monotonically non-increasing and $\tilde{\sigma}_{(0)}$ satisfies the condition. By construction of the sequence $\tilde{\sigma}_{(k+1)}=\theta_{(k)}\tilde{\sigma}_{(k)}$, Eq. 26 is solved with equality by $\theta_{(k)}=(1+2\tilde{\sigma}_{(k)})^{-1/2}$. Moreover, the sequence $\sigma^{(k)}_{i}$ is also non-increasing as

[TABLE]

thus, with Eq. 20 we see that the ESO parameters of ${\bf S}_{(k)}^{1/2}{\bf A}\tau_{(k)}^{1/2}$ are also bounded by $p_{i}$.

For the actual proof of the theorem, note that the inequalities Eqs. 24 and 25 imply

[TABLE]

with

[TABLE]

Thus, combining Lemma 5.5 ($\mu_{g}=0$) and Eq. 27 yields

[TABLE]

With $\gamma\theta_{(k-1)}\leq 1,{\bf S}_{(k+1)}\leq\theta_{(k)}{\bf S}_{(k)}$ and \bar{\Delta}^{(k)}\coloneqq\mathbb{E}\Bigl{\{}\Delta^{(k)}+\|y^{(k)}-y^{(k-1)}\|^{2}_{{\bf Q}{\bf S}_{(k)}^{-1}}\Bigr{\}} we derive the recursion

[TABLE]

Using this inequality recursively, $y^{(-1)}=y^{(0)}$, we arrive at

[TABLE]

where the second estimated follows directly from Lemma 4.3 and the third inequality from $\gamma\leq 1$ which holds by assumption Eq. 20.

As $\bar{\Delta}^{(0)}=\|x^{(0)}-x^{\sharp}\|_{\tau_{(0)}^{-1}}^{2}+\|y^{(0)}-y^{\sharp}\|_{{\bf Y}_{(0)}}^{2}$, $\theta_{(k)}=\frac{\tilde{\sigma}_{(k+1)}}{\tilde{\sigma}_{(k)}}$ and

[TABLE]

which holds by the definition of $\tilde{\sigma}_{(k)}$, it holds that

[TABLE]

*Finally, the assertion follows by Corollary 1 of [12]. *

6 Strongly Convex Case

If both $f_{i}^{\ast}$ and $g$ are strongly convex, we may find step size parameters such that the Algorithm 1 converges linearly.

Theorem 6.1**.**

Let $(x^{\sharp},y^{\sharp})\in{\mathbb{W}}$ be a saddle point and $g,f_{i}^{\ast}$ be strongly convex with constants $\mu_{g},\mu_{i}>0,i=1,\ldots,n$. Let the step sizes $\tau,\sigma_{1},\ldots,\sigma_{n},0<\theta<1$ be chosen such that the ESO parameters $\{v_{i}\}$ of ${\bf S}^{1/2}{\bf A}\tau^{1/2}$ can be estimated as

[TABLE]

and the extrapolation $\theta$ satisfies the lower bounds

[TABLE]

Then the iterates of Algorithm 1 converge linearly to the saddle point, in particular

[TABLE]

holds where the metrics are given by ${\bf X}\coloneqq(\tau^{-1}+2\mu_{g}){\bf I}$, ${\bf Y}\coloneqq({\bf S}^{-1}+2{\bf M}){\bf Q}$ and $\gamma^{2}=\max_{i}v_{i}/p_{i}$.

Proof 6.2**.**

The requirements Eq. 29 on the step sizes $\tau,\sigma_{1},\ldots,\sigma_{n}$ and $\theta$ imply $\theta\|\cdot\|_{{\bf X}}^{2}\geq\|\cdot\|^{2}_{\tau^{-1}}$ and $\theta\|\cdot\|_{{\bf Y}}^{2}\geq\|\cdot\|^{2}_{{\bf Q}{\bf S}^{-1}+2{\bf M}({\bf Q}-{\bf I})}$. Thus, we directly get

[TABLE]

where we denoted

[TABLE]

Combining Eqs. 30 and 5.5 with constant step sizes yields

[TABLE]

Multiplying both sides by $\theta^{-(k+1)}$ and summing over $k=0,\ldots,K-1$ yields

[TABLE]

*where we used again Lemma 4.3 and the non-negativity of norms for the second inequality. Thus, the assertion is proven. *

6.1 Optimal Parameters for Serial Sampling

This analysis is to optimize the convergence rate $\theta$ of Theorem 6.1 for three different serial sampling options where exactly one block is chosen in each iteration. Other sampling strategies, including multi-block, parallel, etc. [38] will be subject of future work.

We will derive the rates and parameters in terms of the condition numbers $\kappa_{i}\coloneqq\|{\bf A}_{i}\|^{2}/(\mu_{g}\mu_{i})$ as these are scaling invariant, thus we cannot improve the rates by simple rescaling of the problem. This can be seen as follows. If we rewrite problem Eq. 2 in terms of the scaled variables $\overline{x}\coloneqq\alpha x$ and $\overline{y}_{i}\coloneqq\beta_{i}y_{i}$, then the corresponding operators $\overline{\bf A}_{i}\coloneqq{\bf A}_{i}/(\alpha\beta_{i})$ have norm $\|\overline{\bf A}_{i}\|=\|{\bf A}_{i}\|/(\alpha\beta_{i})$, the function $\overline{g}(\overline{x})\coloneqq g(\overline{x}/\alpha)$ is $\overline{\mu}_{g}\coloneqq\mu_{g}/\alpha^{2}$ strongly convex and the functions $\overline{f}^{\ast}_{i}(\overline{y}_{i})\coloneqq f^{\ast}_{i}(\overline{y}_{i}/\beta_{i})$ are $\overline{\mu}_{i}\coloneqq\mu_{i}/\beta_{i}^{2}$ strongly convex. Thus the condition numbers are scaling invariant as

[TABLE]

With $\bar{\sigma}_{i}\coloneqq\sigma_{i}\mu_{i}$ and $\bar{\tau}\coloneqq\tau\mu_{g}$, the conditions on the step sizes Eq. 29 become

[TABLE]

for some $\rho<1$. The last condition arises from the ESO parameters of serial sampling which are $v_{i}=\sigma_{i}\tau\|{\bf A}_{i}\|^{2}$, see Example 4.5. Finding optimal parameters is equivalent to equating the above inequalities. Note that the first two conditions (with equality) are equivalent to $\theta\bar{\tau}=(1-\theta)/2$ and $\bar{\sigma}_{i}=\frac{1-\theta}{2(p_{i}-(1-\theta))}$. With these choices, the third condition in Eq. 31 reads

[TABLE]

It follows from Eq. 32 that with $\tilde{\kappa}=1+\kappa/\rho^{2}$ it holds

[TABLE]

Example 6.3** (Serial uniform sampling).**

We first consider uniform sampling, i.e. every block is sampled with the same probability $p_{i}=1/n$. Then it is easy to see that the smallest achievable rate is given by

[TABLE]

and the step sizes become

[TABLE]

Example 6.4** (Serial importance sampling).**

Instead of uniform sampling we may sample \enquoteimportant blocks more often, i.e. we sample every block with a probability proportional to the square root of its condition number $p_{i}=\sqrt{\kappa_{i}}/\sum_{j}\sqrt{\kappa_{j}}$. Then the smallest rate that achieves Eq. 33 is given by

[TABLE]

with $\nu\coloneqq\min_{j}\sqrt{\kappa_{j}}/(1+\sqrt{\tilde{\kappa}_{j}})$ and the step sizes are

[TABLE]

Example 6.5** (Serial optimal sampling).**

Instead of a predefined probability we will seek for an \enquoteoptimal sampling that minimizes the linear convergence rate $\theta$. The optimal sampling can be found by equating condition Eq. 33 for $i=1,\ldots,n$

[TABLE]

Summing Eq. 36 from 1 to $n$ and using that for serial sampling $\sum_{i=1}^{n}p_{i}=1$ leads to

[TABLE]

with step size parameters

[TABLE]

and probabilities

[TABLE]

Remark 6.6** (Minibatches).**

*All arguments above can readily be extended to samplings where at each iteration not only one but a fixed number of blocks are chosen. *

Remark 6.7** (Better Sampling).**

It is easy to see that optimal sampling is better than uniform sampling: if all condition numbers are the same, then the rates for uniform sampling Eq. 34 and optimal sampling Eq. 37 are equal but if they are not, then the rate of optimal sampling is strictly smaller and thus better.

Moreover, optimal sampling is better than importance sampling. To see this, note that due to the monotonicity of $\sqrt{x}/(1+\sqrt{1+x})$ we get

[TABLE]

Remark 6.8** (Comparison to Zhang and Xiao [51]).**

The algorithm of Zhang and Xiao [51] is (almost222In contrast to our work, they have an extrapolation on both primal and dual variables. However, both extrapolations are related as our extrapolation factor is the product of their extrapolation factors.) a special case of the proposed algorithm where each block is picked with probability $p_{i}=1/n$. Here $m$ denotes the size of each block to be processed at every iteration and $n$ the number of blocks. Moreover, they only consider the strongly convex case where $g$ is $\mu_{g}$-strongly convex and all $f_{i}^{\ast}$ are $\mu_{f}$-strongly convex. Then with $R$ being the largest norm of the rows in ${\bf A}$ they achieve

[TABLE]

If the minibatch size is $m=1$, the blocks are chosen to be single rows and the probabilities are uniform, then their rate is slightly worse than ours

[TABLE]

*for any $\rho\geq\frac{1}{2}$. For $m>1$, the rates differ even more as the condition numbers are conservatively estimated. Similarly, the rates can be improved by non-uniform sampling if the row norms are not equal. *

7 Numerical Results

All numerical examples are implemented in python using numpy and the operator discretization library (ODL) [1]. The python code and all example data will be made available on github upon acceptance of this manuscript.

7.1 Non-Strongly Convex PET Reconstruction

In this example we consider positron emission tomography (PET) reconstruction with a total variation (TV) prior. The goal in PET imaging is to reconstruct the distribution of a radioactive tracer from its line integrals [31]. Let ${\mathbb{X}}={\mathbb{R}}^{d_{1}\times d_{2}},d_{1}=d_{2}=250$ be the space of tracer distributions (images) and ${\mathbb{Y}}_{i}={\mathbb{R}}^{|{\mathbb{B}}_{i}|}$ the data spaces where ${\mathbb{B}}_{i}\subset\{1,\ldots,N\},N=200\cdot 250$ (200 views around the object) are subsets of indices with ${\mathbb{B}}_{i}\cap{\mathbb{B}}_{j}=\emptyset$ if $i\neq j$ and $\cup_{i=1}^{n}{\mathbb{B}}_{i}=\{1,\ldots,N\}$. All samplings in this example divide the views equidistantly. It is standard that PET reconstruction can be posed as the optimization problem Eq. 1 where the data fidelity term is given by the Kullback–Leibler divergence

[TABLE]

where it is convention that $0\log 0\coloneqq 0$. The operator ${\bf A}$ is a scaled X-ray transform where in each of 200 directions 250 line integrals are computed with the astra toolbox [49, 48]. The prior is the TV of $x$ with non-negativity constraint, i.e. $g(x)=\alpha\|\nabla x\|_{1,2}+\chi_{\geq 0}(x)$, with regularization parameter $\alpha=0.2$ and the gradient operator $\nabla x=(\nabla_{1}x,\nabla_{2}x)\in{\mathbb{R}}^{d_{1}\cdot d_{2}\times 2}$ is discretized by forward differences in horizontal and vertical direction, cf. [11] for details. The $1,2$-norm of these gradients is defined as $\|x\|_{1,2}\coloneqq\sum_{j}\sqrt{(\nabla_{1}x_{j})^{2}+(\nabla_{2}x_{j})^{2}}$. The Fenchel conjugate of the Kullback–Leibler divergence Eq. 38 is

[TABLE]

its proximal operator given by

[TABLE]

The proximal operator for $g$ is approximated with 20 iterations of the fast gradient projection method (FGP) [5] with a warm start applied to the dual problem.

Parameters

In this experiment we choose $\gamma=0.99$, $\theta=1$ and all samplings are uniform, i.e. $p_{i}=1/n$. The number of subsets varies between $n=1$ (deterministic case), 50 and 250. The other step size parameters are chosen as

• •

PDHG, Pesquet&Repetti [34]: $\sigma_{i}=\tau=\gamma/\|{\bf A}\|\approx 6.9\cdot 10^{-4}$

• •

SPDHG: $\sigma_{i}=\gamma/\|{\bf A}_{i}\|\approx 2.2\cdot 10^{-3}$, $\tau=\gamma/(n\max_{i}\|{\bf A}_{i}\|)\approx 2.2\cdot 10^{-4}$

Results

Fig. 2 on the left shows that the ergodic Bregman distance converges with rate $1/k$ as proven in Theorem 4.6. On the right we compare the deterministic PDHG with the randomized SPDHG and the algorithm of Pesquet&Repetti. It can be clearly seen that the proposed SPDHG converges much faster than both the algorithm of Pesquet&Repetti and the deterministic PDHG. Some example images are found in Fig. 2 after 5 epochs which again highlight the speed-up gained by randomization.

7.2 TV denoising with Gaussian Noise (Primal Acceleration)

In the second example we consider denoising of an image that is degraded by Gaussian noise with the help of the anisotropic TV. This can be achieved by solving Eq. 1 with ${\mathbb{X}}={\mathbb{R}}^{d_{1}\times d_{2}}$, $d_{1}=442$, $d_{2}=331$, the data fit $g(x)=1/(2\alpha)\|x-b\|^{2}_{2}$ is the squared Euclidean norm and the prior the (anisotropic) TV $f_{i}(y_{i})=\|y_{i}\|_{1},{\bf A}_{i}=\nabla_{i}$ and $n=2$. Instead of the isotropic TV as in the previous example we consider here the anisotropic version as it is separable in the direction of the gradient. The regularization parameter is chosen to be $\alpha=0.12$. See e.g. [12] for details on convex conjugates and proximal operators of these functionals.

Parameters

In this experiment we choose $\gamma=0.99$ and the sampling to be uniform, i.e. $p_{i}=1/n$. The number of subsets are either $n=1$ in the deterministic case or $n=2$ in the stochastic case. The (initial) step size parameters are

• •

PDHG, PA-PDHG, Pesquet&Repetti: $\sigma_{i}^{(0)}=\tau^{(0)}=\gamma/\|{\bf A}\|\approx 0.35$

• •

SPDHG, PA-SPDHG: $\sigma_{i}^{(0)}=\gamma/\|{\bf A}_{i}\|\approx 0.50$, $\tau^{(0)}=\gamma/(n\max_{i}\|{\bf A}_{i}\|)\approx 0.25$

The step sizes for acceleration vary with the iteration with the primal step size $\tau_{(k)}$ getting smaller and the dual step size $\sigma^{(k)}$ getting larger. The extrapolation factor $\theta$ is chosen to be 1 for non-accelerated and converging to 1 for accelerated algorithms.

Results

The quantitative results in Fig. 4 show that the accelerated algorithms are much faster than the non-accelerated versions. Moreover, it can be seen that the stochastic variant of the accelerated PA-PDHG is even faster than its deterministic variant. In addition, the results show that the accelerated SPDHG indeed converges as $1/K^{2}$ in the norm of the primal part. Visual assessment of the denoised images in Fig. 4 confirms these conclusions.

7.3 Huber-TV Deblurring (Dual Acceleration)

In the third example we consider deblurring with known convolution kernel where the forward operator ${\bf A}_{1}$ resembles the convolution of images in ${\mathbb{X}}={\mathbb{R}}^{d_{1}\times d_{2}},d_{1}=408,d_{2}=544$ with a motion blur of size $15\times 15$. The noise is modeled to be Poisson with a constant background of 200 compared to the approximate data mean of $694.3$. We further assume to have the knowledge that the reconstructed image should be non-negative and upper-bounded by 100. By the nature of the forward operator ${\bf A}x\geq 0$ whenever $x\geq 0$. Therefore the solution to Eq. 1 with the Kullback–Leibler divergence Eq. 38 remains the same if we replace the Kullback–Leibler divergence by the differentiable

[TABLE]

which has a $\left(\max_{i}b_{i}/r_{i}^{2}\right)$ Lipschitz continuous gradient. The Lipschitz constant is well-defined and non-zero as both the data $b_{i}$ as well as the background $r_{i}$ are positive. In our numerically example it is approximately $0.31$.

Prior smoothness information is represented by the anisotropic TV with Huberized norm

[TABLE]

for $i=2,3$ where $y_{j}=\nabla_{i-1}x_{j}$ are finite differences, $\eta=1$ and regularization parameter $\alpha=0.1$. The constraints on the image are enforced by the characteristic function $g=\imath_{{\mathbb{B}}}$ with ${\mathbb{B}}=\{x\in{\mathbb{X}}\,|\,0\leq x_{j}\leq 100\}$.

The convex conjugate of the modified Kullback–Leibler divergence Eq. 40 is

[TABLE]

which is $\left(\min_{i}\frac{r_{i}^{2}}{b_{i}}\right)$-strongly convex with proximal operator

[TABLE]

Parameters

In this experiment we choose $\gamma=0.99$ and consider uniform sampling, i.e. $p_{i}=1/n$. The number of subsets are either $n=1$ in the deterministic case or $n=3$ in the stochastic case. The (initial) step size parameters are chosen to be

• •

PDHG: $\sigma_{i}=\tau=\gamma/\|{\bf A}\|\approx 0.095$

• •

DA-PDHG: $\tilde{\sigma}^{(0)}_{i}=\mu_{f}/\|{\bf A}\|\approx 0.096$ , $\tau^{(0)}=\gamma/\|{\bf A}\|\approx 0.095$

• •

DA-SPDHG: $\tilde{\sigma}^{(0)}=\min_{i}\frac{\mu_{i}p_{i}^{2}}{\tau^{0}\|{\bf A}_{i}\|^{2}+2\mu_{i}p_{i}(1-p_{i})}\approx 0.073$,

DA-SPDHG: $\tau^{(0)}=1/(n\max_{i}\|{\bf A}_{i}\|)\approx 0.032$

Results

The quantitative results in Fig. 6 show that the algorithm converges indeed with rate $O(1/K^{2})$ as proven in Theorem 5.1. Moreover, they also show that randomization and acceleration can be used in conjunction for further speed-ups. Some example images are shown in Fig. 6 which show that randomization may lead to sharper images with the same number of epochs.

7.4 PET Reconstruction (Linear Rate)

For the final example we turn back to PET reconstruction but this time with linear convergence rate. This means we want to solve the same minimization problem as in the first example, but now we replace the Kullback–Leibler functional by its modified version as in the previous example. We note again that this does not change the solution of the minimization problem. Moreover, to make TV strongly convex we add another regularization term $\mu/2\|x\|_{2}^{2}$ to $g$. Note that the proximal operator of TV (indeed any functional) with added squared $\ell^{2}$-norm, i.e. $g(x)=\alpha\operatorname{TV}(x)+\mu/2\|x\|_{2}^{2}$, can be solved by means of the original proximal operator $\operatorname{prox}^{\sigma}_{g}\left(z\right)=\operatorname{prox}^{\sigma\alpha/(1+\sigma\mu)}_{\operatorname{TV}}\left(z/(1+\sigma\mu)\right)$. The regularization parameters are chosen as $\alpha=0.05$ and $\mu=0.5$.

Parameters

In this experiment we choose $\rho=0.99$ and the sampling to be uniform as the operators ${\bf A}_{i}$ all have similar norms. The step size parameters are chosen as derived in Section 6.1, in particular, we choose

• •

PDHG: $\sigma\approx 3.8\cdot 10^{-4},\tau\approx 4.8\cdot 10^{-3},\theta\approx 0.995$

• •

Pesquet&Repetti: $\sigma_{i}=\tau=\gamma/\|{\bf A}\|\approx 1.4\cdot 10^{-3}$

• •

SPDHG ($n=10$ subsets): $\sigma_{i}\approx 1.2\cdot 10^{-3}$, $\tau\approx 1.5\cdot 10^{-3}$, $\theta^{n}\approx 0.985$

• •

SPDHG ($n=50$): $\sigma_{i}\approx 2.4\cdot 10^{-3}$, $\tau\approx 5.8\cdot 10^{-4}$, $\theta^{n}\approx 0.971$

Note that the contraction rates of one epoch $\theta^{n}$ already indicate that SPDHG ($n=50$) may be faster than PDHG and SPDHG ($n=10$).

Results

The quantitative results in Fig. 8 in terms of both distance to saddle point and objective value show that randomization speeds up the convergence so that both SPDHG and the algorithm of Pesquet&Repetti are faster than the deterministic PDHG. Interestingly, while more subsets make SPDHG faster, this does not hold for the algorithm of Pesquet&Repetti where the speed seems to be constant with respect to the number of subsets. Moreover, the plot on the left confirms the linear convergence as proven in Theorem 6.1. The visual results in Fig. 8 confirm these observations as SPDHG with 50 subsets and 10 epochs is (in contrast to PDHG) visually already very close to the saddle point.

8 Conclusions and Future Work

We proposed a natural stochastic generalization of the deterministic PDHG algorithm to convex-concave saddle point problems that are separable in the dual variable. The analysis was carried out in the context of arbitrary samplings which enabled us to obtain known deterministic convergence results as special cases. We proposed optimal choices of the step size parameters with which the proposed algorithm showed superior empirical performance on a variety of optimization problems in imaging.

In the future, we would like to extend the analysis to include iteration dependent (adaptive) probabilities [16] and strong convexity parameters to further exploit the structure of many relevant problems. Moreover, the present optimal sampling strategies are only for scalar-valued step sizes and serial sampling. In the future, we wish to extend this to other sampling strategies such as multi-block or parallel sampling.

Appendix A Postponed Proofs

Proof A.1** (Proof of Lemma 4.3).**

With the definition of $y^{+}$ and ${\bf Q}$ we have by completing the norm for any $x\in{\mathbb{X}}$ that

[TABLE]

where we used $z_{i}\coloneqq p_{i}^{-1}{\bf S}_{i}^{-1/2}(\hat{y}_{i}-y_{i})$. Moreover, the expectation of the second term of the right hand side of Eq. 41 can be estimated as

[TABLE]

where the first inequality is due to the ESO inequality Eq. 4. Inserting $z$ leads to

[TABLE]

*where the last equation holds true by the definition of the expectation. Combining the expected value of inequality Eq. 41 with Eqs. 42 and 43 yields the assertion. *

Proof A.2** (Proof of Lemma 4.8).**

By the definition of the proximal operator, for any $(x,y)\in{\mathbb{W}}$ it holds that

[TABLE]

for $i=1,\ldots,n$. Summing twice all inequalities and exploiting the identity

[TABLE]

yields

[TABLE]

*where we used the definition of the inner product and the norm on the product space ${\mathbb{Y}}$. It now suffices to complete the generalized distances ${\mathcal{G}}(x^{(k+1)}|w)$ and ${\mathcal{F}}(\hat{y}^{(k+1)}|w)$. *

Proof A.3** (Proof of Lemma 5.5).**

We follow a similar line of arguments as in the proof of Theorem 4.6. Note that for any saddle point $w^{\sharp}=(x^{\sharp},y^{\sharp})$ we have

[TABLE]

such that the estimate of Lemma 4.8 can be written with $w=w^{\sharp}$ as

[TABLE]

using the rule Eq. 14. With Eqs. 11 and 12 and again Eq. 14 we arrive at

[TABLE]

Inserting the extrapolation Eq. 22 into the inner product yields

[TABLE]

The assertion is shown by taking the expectations $\mathbb{E}^{(k,k-1)}\coloneqq\mathbb{E}^{(k)}\mathbb{E}^{(k-1)}$ on Eq. 44, using Eq. 45 and estimating the last inner product by Lemma 4.3 as

[TABLE]