Acceleration of the Pdhgm on Strongly Convex Subspaces

We propose several variants of the primal-dual method due to Chambolle and Pock. Without requiring full strong convexity of the objective functions, our methods are accelerated on subspaces with strong convexity. This yields mixed rates, O (1/N 2) with respect to initialisation and O (1/N) with respect to the dual sequence, and the residual part of the primal sequence. We demonstrate the eecacy of the proposed methods on image processing problems lacking strong convexity, such as total generalised variation denoising and total variation deblurring.


Introduction
Let G : X → R and F : Y → R be convex, proper, and lower semicontinuous functionals on Hilbert spaces X and Y , possibly infinite dimensional. Also let K ∈ L(X ; Y ) be a bounded linear operator. We then wish to solve the problem This can under mild conditions on F (see, for example, [1,2]) also be written with the help of the convex conjugate F * in the minimax form min x∈X max y∈Y G(x) + K x, y − F * (y).
One possibility for the numerical solution of the latter form is the primal-dual algorithm of Chambolle and Pock [3], a type of proximal point or extragradient method, also classified as the 'modified primal-dual hybrid gradient method' or PDHGM by Esser et al. [4]. If either G or F * is strongly convex, the method can be accelerated to O(1/N 2 ) convergence rates of the iterates and an ergodic duality gap [3]. But what if we have only partial strong convexity? For example, what if for a projection operator P to a subspace X 0 ⊂ X , and strongly convex G 0 : X 0 → R? This kind of structure is common in many applications in image processing and data science, as we will more closely review in Sect. 5. Under such partial strong convexity, can we obtain a method that would give an accelerated rate of convergence at least for Px?
We provide a partially positive answer: we can obtain mixed rates, O(1/N 2 ) with respect to initialisation, and O(1/N ) with respect to bounds on the 'residual variables' y and (I − P)x. In this respect, our results are similar to the 'optimal' algorithm of Chen et al. [5]. Instead of strong convexity, they assume smoothness of G to derive a primal-dual algorithm based on backward-forward steps, instead of the backward-backward steps of [3].
The derivation of our algorithms is based, firstly, on replacing simple step length parameters by a variety of abstract step length operators and, secondly, a type of abstract partial strong monotonicity property the full details of which we provide in Sect. 2. Here T is an auxiliary step length operator. Our factor of strong convexity is a positive semidefinite operator ≥ 0; however, to make our algorithms work, we need to introduce additional artificial strong convexity through another operator , which may not satisfy 0 ≤ ≤ . This introduces the penalty term in (1). The exact procedure can be seen as a type of smoothing, famously studied by Nesterov [6], and more recently, for instance, by Beck and Teboulle [7]. In these approaches, one computes a priori a level of smoothing-comparable to -needed to achieve prescribed solution quality. One then solves a smoothed problem, which can be done at O(1/N 2 ) rate. However, to obtain a solution with higher quality than the a priori prescribed one, one needs to solve a new problem from scratch, as the smoothing alters the problem being solved. One can also employ restarting strategies, to take some advantage of the previous solution, see, for example, [8]. Our approach does not depend on restarting and a priori chosen solution qualities: the method will converge to an optimal solution to the original non-smooth problem. Indeed, the introduced additional strong convexity is controlled automatically.
The 'fast dual proximal gradient method', or FDPG [9], also possesses different type of mixed rates, O(1/N ) for the primal, and O(1/N 2 ) for the dual. This is, however, under standard strong convexity assumptions. Other than that, our work is related to various further developments from the PDHGM, such as variants for nonlinear K [10,11] and non-convex G [12]. The PDHGM has been the basis for inertial methods for monotone inclusions [13] and primal-dual stochastic coordinate descent methods without separability requirements [14]. Finally, the FISTA [15,16] can be seen as a primal-only relative of the PDHGM. Not attempting to do full justice here to the large family of closely related methods, we point to [4,17,18] for further references.
The contributions of our paper are twofold: firstly, to paint a bigger picture of what is possible, we derive a very general version of the PDHGM. This algorithm, useful as a basis for deriving other new algorithms besides ours, is the content of Sect. 2. In this section, we provide an abstract bound on the iterates of the algorithm, later used to derive convergence rates. In Sect. 3, we extend the bound to include an ergodic duality gap under stricter conditions on the acceleration scheme and the step length operators. A by-product of this work is the shortest convergence rate proof for the accelerated PDHGM known to us. Afterwards, in Sect. 4, we derive from the general algorithm two efficient mixed-rate algorithms for problems exhibiting strong convexity only on subspaces. The first one employs the penalty or smoothing ψ on both the primal and the dual. The second one only employs the penalty on the dual. We finish the study with numerical experiments in Sect. 5. The main results of interest for readers wishing to apply our work are Algorithms 3 and 4 along with the respective convergence results, Theorems 4.1 and 4.2.

Notation
To make the notation definite, we denote by L(X ; Y ) the space of bounded linear operators between Hilbert spaces X and Y . For T, S ∈ L(X ; X ), the notation T ≥ S means that T − S is positive semidefinite; in particular, T ≥ 0 means that T is positive semidefinite. In this case, we also denote The identity operator is denoted by I , as is standard.
For 0 ≤ M ∈ L(X ; X ), which can possibly not be selfadjoint, we employ the notation a, b M := Ma, b , and a M := a, a M . ( We also use the notation T −1, * := (T −1 ) * .

Background
As in the introduction, let us be given convex, proper, lower semicontinuous functionals G : X → R and F * : Y → R on Hilbert spaces X and Y , as well as a bounded linear operator K ∈ L(X ; Y ). We then wish to solve the minimax problem assuming the existence of a solution u = ( x, y) satisfying the optimality conditions − K * y ∈ ∂G( x), and K x ∈ ∂ F * ( y  [19] for several sufficient conditions. The primal-dual method of Chambolle and Pock [3] for the solving (P) consists of iterating In the basic version of the algorithm, ω i = 1, τ i ≡ τ 0 , and σ i ≡ σ 0 , assuming that the step length parameters satisfy τ 0 σ 0 K 2 < 1. The method has O(1/N ) rate for the ergodic duality gap [3]. If G is strongly convex with factor γ , we may use the acceleration scheme [3] to achieve O(1/N 2 ) convergence rates of the iterates and an ergodic duality gap, defined in [3]. To motivate our choices later on, observe that σ 0 is never used expect to calculate σ 1 . We may therefore equivalently parametrise the algorithm by δ = 1 − K 2 τ 0 σ 0 > 0. We note that the order of the steps in (3) is different from the original ordering in [3]. This is because with the present order, the method (3) may also be written in the proximal point form. This formulation, first observed in [20] and later utilised in [10,11,21], is also what we will use to streamline our analysis. Introducing the general variable splitting notation, for the monotone operator and the preconditioning or step length operator M basic,i := We note that the optimality conditions (OC) can also be encoded as 0 ∈ H ( u).

Abstract Partial Monotonicity
Our plan now is to formulate a general version of (3), replacing τ i and σ i by operators T i ∈ L(X ; X ) and i ∈ L(Y ; Y ). In fact, we will need two additional operators T i ∈ L(X ; X ) andT i ∈ L(Y ; Y ) to help communicate change in T i to i . They replace ω i in (3b) and (7), operating asT i+1 K T −1 i ≈ ω i K from both sides of K . The role of T i is to split the original primal step length τ i in the space X into the two parts T i and T i with potentially different rates. The role ofT i is to transfer T i into the space Y , to eventually control the dual step length i . In the basic algorithm (3), we would simply To start the algorithm derivation, we now formulate abstract forms of partially strong monotonicity. As a first step, we take subsets of invertible operators as well as subsets of positive semidefinite operators We assume T andT closed with respect to composition: We use the sets K andK as follows. We suppose that ∂G is partially strongly (ψ, T , K)-monotone, meaning that for all x, x ∈ X, T ∈ T , ∈ [0, ] + K holds for some family of functionals {ψ T : X → R}, and a linear operator 0 ≤ ∈ L(X ; X ) which models partial strong monotonicity. The inequality in (G-PM), and all such set inequalities in the remainder of this paper, is understood to hold for all elements of the sets ∂G(x ) and ∂G(x). The operator T ∈ T acts as a testing operator, and the operator ∈ K as introduced strong monotonicity. The functional ψ T −1, * ( − ) is a penalty corresponding to the test and the introduced strong monotonicity. The role of testing will become more apparent in Sect. 2.4. Similarly to (G-PM), we assume that ∂ F * is (φ,T ,K)monotone with respect toT in the sense that for all y, y ∈ Y , T ∈T , R ∈K holds for some family of functionals {φ T : Y → R}. Again, the inequality in (F * -PM) is understood to hold for all elements of the sets ∂ F * (y ) and ∂ F * (y).
In our general analysis, we do not set any conditions on ψ and φ, as their role is simply symbolic transfer of dissatisfaction of strong monotonicity into a penalty in our abstract convergence results.
Let us next look at a few examples on how (G-PM) or (F * -PM) might be satisfied. First we have the very wellbehaved case of quadratic functions.
Example 2.2 An indicator function ι A of a convex bounded set A satisfies the conditions of Lemma 2.1. This is generally what we will use and need.

A General Algorithm and the Idea of Testing
The only change we make to the proximal point formulation (5) of the method (3) is to replace the basic step length or preconditioning operator M basic,i by the operator As we have remarked, the operatorsT i+1 and T i play the role of ω i , acting from both sides of K . Our proposed algorithm can thus be characterised as solving on each iteration i ∈ N for the next iterate u i+1 the preconditioned proximal point problem To study the convergence properties of (PP), we define the testing operator It will turn out that multiplying or 'testing' (PP) by this operator will allow us to derive convergence rates. The testing of (PP) by S i is why we introduced testing into the monotonicity conditions (G-PM) and (F * -PM). If we only tested (PP) with S i = I , we could at most obtain ergodic convergence of the duality gap for the unaccelerated method. But by testing with something appropriate and faster increasing, such as (11), we are able to extract better convergence rates from (PP). We also set for some i ∈ [0, ] + K and R i+1 ∈K. We will see in Sect. 2.6 that¯ i is a factor of partial strong monotonicity for H with respect to testing by S i . With this, taking a fixed δ > 0, the properties will turn out to be the crucial defining properties for the convergence rates of the iteration (PP). The method resulting from the combination of (PP), (C1), and (C2) can also be expressed as Algorithm 1. The main steps in developing practical algorithms based on Algorithm 1 will be in the choice of the various step length operators. This will be the content of Sects. 3 and 4. Before this, we expand the conditions (C1) and (C2) to see how they might be satisfied and study abstract convergence results.

A Simplified Condition
We expand as well as and .
We observe that if S, T ∈ L(X ; Y ), then for arbitrary invertible Z ∈ L(Y ; Y ) a type of Cauchy (or Young) inequality holds, namely The inequality here can be verified using the basic Cauchy inequality 2 x, y ≤ x 2 + y 2 . Applying (14) in (12), we see that (C2) is satisfied when for some invertible Z i ∈ L(X ; X ). The second condition in (15) is satisfied as an equality if By the spectral theorem for self-adjoint operators on Hilbert spaces (e.g. [22,Chapter 12]), we can make the choice (16) if Equivalently, by the same spectral theorem, T −1 i T i ∈ Q. Therefore, we see from (15) that (C2) holds when Also, (C1) can be rewritten as

Basic Convergence Result
Our main result on Algorithm 1 is the following theorem, providing some general convergence estimates. It is, however, important to note that the theorem does not yet directly prove convergence, as its estimates depend on the rate of decrease in T N T * N , as well as the rate of increase in the penalty sum coming from the dissatisfaction of strong convexity. Deriving these rates in special cases will be the topic of Sect. 4. Theorem 2.1 Let us be given K ∈ L(X ; Y ), and convex, proper, lower semicontinuous functionals G : X → R and F * : Y → R on Hilbert spaces X and Y , satisfying (G-PM) and (F * -PM). Pick δ ∈ (0, 1), and suppose (C1) and (C2) are satisfied for each i ∈ N for some invertible T i ∈ L(X ; X ), Then, the iterates of Algorithm 1 satisfy for Remark 2.1 The term D i+1 , coming from the dissatisfaction of strong convexity, penalises the basic convergence, which is on the right-hand side of (17) presented by the constant C 0 . If T N T N is of the order O(1/N 2 ), at least on a subspace, and we can bound the penalty D i+1 ≤ C for some constant C, then we clearly obtain mixed O(1/N 2 ) + O(1/N ) convergence rates on the subspace. If we can assume that D i+1 actually converges to zero at some rate, then it will even be possible to obtain improved convergence rates. Since typically T i ,T i+1 0 reduce to scalar factors within D i+1 , this would require prior knowledge of the rates of convergence x i → x and y i → y. Boundedness of the iterates , we can, however, usually ensure.
Recalling the definition of S i from (11), and of H from (6), it follows An application of (G-PM) and (F * -PM) consequently gives Using the expression (13) for S i¯ i , and (18) for D i+1 , we thus deduce For We observe that S i M i in (12) is self-adjoint as we have assumed that are self-adjoint. In consequence, using (20) we obtain Combining (19) and (21) through (PP), we thus obtain Summing (22) over i = 1, . . . , N − 1, and applying (C2) to estimate S N M N from below, we obtain (17).

Scalar Off-diagonal Updates and the Ergodic Duality Gap
One relatively easy way to satisfy (G-PM), (F * -PM), (C1) and (C2) is to take the 'off-diagonal' step length operatorsT i and T i as equal scalars. Another good starting point would be to choose T i = T i . We, however, do not explore this route in the present work. Instead, we now specialise Theorem 2.1 to the scalar case. We then explore ways to add estimates of the ergodic duality gap into (17). While this would be possible in the general framework through convexity notions analogous to (G-PM) and (F * -PM), the resulting gap would not be particularly meaningful. We therefore concentrate on the scalar off-diagonal updates to derive estimates on the ergodic duality gap.
Observe also that M i is under this setup self-adjoint if T i and i+1 are. For simplicity, we now assume φ and ψ to satisfy the identities The monotonicity conditions (G-PM) and (F * -PM) then simplify into holding for all x, x ∈ X , and ∈ [0, ] + K; and holding for all y, y ∈ Y , and R ∈K.

The Ergodic Duality Gap and Convergence
To study the convergence of an ergodic duality gap, we now introduce convexity notions analogous to (G-pm) and (F*-pm). Namely, we assume to hold for all x, x ∈ X and ∈ [0, ] + K and to hold for all y, y ∈ Y and R ∈K. Clearly these imply (G-pm) and (F*-pm).
To define an ergodic duality gap, we set and define the weighted averages With these, the ergodic duality gap at iteration N is defined as the duality gap for (x N , Y N ), namely and we have the following convergence result.
Thus, q N is of the order (N 2 ), while τ N T N = τ 2 N I is of the order O(1/N 2 ). Therefore, (26) shows O(1/N 2 ) convergence of the squared distance to solution. For O(1/N 2 ) convergence of the ergodic duality gap, we need to slow down (4) to

Remark 3.2
The result (28) can be improved to estimate Lemma 2] we know the following for the rule ω i = 1/ √ 1 + 2γ τ i : given λ > 0 and N ≥ 0 with γ τ N ≤ λ, for any ≥ 0 holds If we pick N = 0 and λ = γ τ 0 , this says Proof (Theorem 3.1) The non-gap estimate in the last paragraph of the theorem statement, where λ = 1, we modify G N := 0, is a direct consequence of Theorem 2.1. We therefore concentrate on the estimate that includes the gap, and fix λ = 1/2. We begin by expanding Since then i ∈ ([0, ] + K)/2, and R i+1 ∈K/2, we may take = 2 i and R = 2R i+1 in (G-pc) and (F*-pc). It follows Using (2) and (24), we can make all of the factors '2' and '1/2' in this expression annihilate each other. With D i+1 as in (27) and λ = 1/2, we therefore have A little bit of reorganisation and referral to (13) for the expansion of S i¯ i thus gives Let us write Observing here the switches between the indices i + 1 and i of the step length parameters in comparison with the last step of (29), we thus obtain We note that S i M i in (12) is self-adjoint as we have assumed T i and i+1 to be, and taken T i andT i+1 to be scalars times the identity. We therefore deduce from the proof of Theorem 2.1 that (21) holds. Using (PP) to combine (21) and (30), we thus deduce Summing this for i = 0, . . . , N − 1 gives with C 0 from (27) the estimate We want to estimate the sum of the gaps G i + in (31). Using the convexity of G and F * , we observe Also, by (25) and simple reorganisation All of (32)-(34) together give Another use of (25) gives Thus, where the remainder At a solution u = ( x, y) to (OC), we have Kx ∈ ∂ F * (ŷ), so r N ≥ 0 provided q N ≤q N . But q N −q N = τ −1 0 − τ −1 N , so this is guaranteed by our assumption (C3 ). Using (35) in (31) therefore gives A referral to (C2) to estimate S N M N from below shows (26), concluding the proof.

Convergence Rates in Special Cases
To derive a practical algorithm, we need to satisfy the update rules (C1) and (C2), as well as the partial monotonicity conditions (G-PM) and (F * -PM). As we have already discussed in Sect. 3, this can be done when for some τ i > 0 we set The result of these choices is Algorithm 2, whose convergence we studied in Theorem 3.1. Our task now is to verify its conditions, in particular (G-pc) and (F*-pc) [alternatively (F*-pm) and (G-pm)], as well as (C1 ), (C2 ), and (C3 ) for of the projection form γ P.

An Approach to Updating
We have not yet defined an explicit update rule for i+1 , merely requiring that it has to satisfy (C2 ) and (C1 ). The former in particular requires Hiring the help of some linear operator F ∈ L(L(Y ; Y ); L(Y ; Y )) satisfying our approach is to define Then, (C2 ) is satisfied provided , the condition (C1 ) reduces into the satisfaction for each i ∈ N of To apply Theorem 3.1, all that remains is to verify in special cases the conditions (40) together with (C3 ) and the partial strong convexity conditions (G-pc) and (F*-pc).

When is a Multiple of a Projection
We now take =γ P for someγ > 0, and a projection operator P ∈ L(X ; X ): idempotent, P 2 = P, and self-adjoint, P * = P. We let P ⊥ := I − P. Then, P ⊥ P = P P ⊥ = 0. With this, we assume that K is such that for someγ ⊥ > 0 holds To unify our analysis for gap and non-gap estimates of Theorem 3.1, we now pick λ = 1/2 in the former case, and λ = 1 in the latter. We then pick 0 ≤ γ ≤ λγ , and 0 ≤ γ ⊥ i ≤ λγ ⊥ , and set With this, τ i , τ ⊥ i > 0 guarantee T i ∈ Q. Moreover, T i is selfadjoint. Moreover, i ∈ λ([0, ] + K), exactly as required in both the gap and the non-gap cases of Theorem 3.1. Since we are encouraged to take Observe that (43) satisfies (38). Inserting (43) into (39), we obtain Since i+1 is now equivalent to a scalar, (40b), we also take R i+1 = ρ i+1 I , assuming for someρ > 0 that Setting we thus expand (40) as We are almost ready to state a general convergence result for projective . However, we want to make one more thing more explicit. Since the choices (42) satisfy we suppose for simplicity that for some ψ ⊥ : P ⊥ X → R and φ : Y → R. The conditions (G-pc) and (F*-pc) reduce in this case to the satisfaction for someγ ,γ ⊥ ,ρ > 0 of for all x, x ∈ X and 0 ≤ γ ⊥ ≤γ ⊥ , as well as of for all y, y ∈ Y and 0 ≤ ρ ≤ρ. Analogues of (G-pm) and (F*-pm) can be formed.
To summarise the findings of this section, we state the following proposition.
Observe that presently Proof As we have assumed through (47), or otherwise already verified its conditions, we may apply Theorem 3.1. Multiplying (26) by τ N τ N , we obtain Now, observe that solving (45a) exactly gives Therefore, we have the estimate With this, (50) yields (48).

Primal and Dual Penalties with Projective
We now study conditions that guarantee the convergence of the sum τ N τ N N −1 i=0 D i+1 in (48). Indeed, the right-hand sides of (45b) and (45c) relate to D i+1 . In most practical cases, which we study below, φ and ψ transfer these righthand side penalties into simple linear factors within D i+1 . Optimal rates are therefore obtained by solving (45b) and (45c) as equalities, with the right-hand sides proportional to each other. Since η i ≥ 0, and it will be the case that η i = 0 for large i, we, however, replace (45c) by the simpler condition Then, we try to make the left-hand sides of (45b) and (53) proportional with only τ ⊥ i+1 as a free variable. That is, for some proportionality constant ζ > 0, we solve Multiplying both sides of (54) by ζ −1 τ i+1 τ ⊥ i+1 , gives on τ ⊥ i+1 the quadratic condition Thus, Solving (45b) and (53) as equalities, (54) and (55) give Note that this quantity is non-negative exactly when ω ⊥ i ≥ ω i . We have This quickly yields ω ⊥ i ≥ ω i if ω i ≤ 1. In particular, (56) is non-negative when ω i ≤ 1.
The next lemma summarises these results for the standard choice of ω i .

4: Update
(59c) 6: until a stopping criterion is fulfilled. exploit with specific φ and ψ the telescoping property stemming from the non-negativity of the last term of (56). There is still, however, one matter to take care of. We need ρ i ≤ λρ and γ ⊥ i ≤ λγ ⊥ , although in many cases of practical interest, the upper bounds are infinite and hence inconsequential. We calculate from (55) and (57) that Therefore, we need to choose ζ and τ 0 to satisfy 2ζ γ τ 0 ≤ (λγ ⊥ ) 2 . Likewise, we calculate from (56), (57), and (60) that This tells us to choose τ 0 and ζ to satisfy 2 Overall, we obtain for τ 0 and ζ the condition This can always be satisfied through suitable choices of τ 0 and ζ .
If now φ ≡ C φ and ψ ≡ C ⊥ ψ , using the non-negativity of (56), we calculate Similarly Using these expression to expand (49), we obtain the following convergence result. (61). Then, Algorithm 3 satisfies for some C 0 , C τ > 0 the estimate
Proof During the course of the derivation of Algorithm 3, we have verified (45), solving (45a) as an equality. Moreover, Lemma 4.1 and (61) guarantee (47). We may therefore apply Proposition 4.1. Inserting (62) and (63) into (48) and (49) gives The condition ζ ≤ (τ ⊥ 0 ) −2 now guarantees τ ⊥ N ≤ ζ −1/2 through (58). Now we note that τ i is not used in Algorithm 3, so it only affects the convergence rate estimates. We therefore simply take τ 0 = τ 0 , so that τ N = τ N for all N ∈ N. With this and the bound τ N ≤ C τ /N from Remark 3.2, (64) follows by simple estimation of (65).   N ) rate, similarly to that derived in [5] for a type of forward-backward splitting algorithm for smooth G. Ours is of course backward-backward type algorithm. It is interesting to note that using the differentiability properties of infimal convolutions [23,Proposition 18.7], and the presentation of a smooth G as an infimal convolution, it is formally possible to derive a forward-backward algorithm from Algorithm 3. The difficulties lie in combining this conversion trick with conditions on the step lengths.

Dual Penalty Only with Projective
Continuing with the projective setup of Sect. 4.2, we now study the case K = {0}, that is, when only the dual penalty φ is available with ψ ≡ 0. To use Proposition 4.1, we need to satisfy (47) and (45), with (45a) holding as an equality.
With respect to τ ⊥ i+1 , the left-hand side of (45c) is maximised (and the penalty on the right-hand side minimised) when (66) is minimised. Thus, we solve (66) exactly, which gives In consequence ω ⊥ i = ω −1 i , and (45c) becomes In order to simultaneously satisfy (45a), this suggests for some, yet undetermined, a i > 0, to choose Since η i ≥ 0, (67) is satisfied with the choice (68) if we take To use Proposition 4.1, we need to satisfy ρ i+1 ≤ λρ. Since (68) implies that { τ i } ∞ i=0 is non-increasing, we can satisfy this for large enough i if a i 0. To ensure satisfaction for all i ∈ N, it suffices to take {a i } ∞ i=0 non-increasing, and satisfy the initial condition The rule τ i+1 = ω i τ i and (68) give τ −2 i+1 = τ −2 i + a i . We therefore see that Assuming φ to have the structure (46), moreover, Thus, the rate (48) in Proposition 4.1 states The convergence rate is thus completely determined by μ N 0 and μ N 1 .
For a more generally applicable algorithm, suppose φ(y i+1 − y) ≡ C φ as in Theorem 4.1. We need to choose a i . One possibility is to pick some q ∈ (0, 1] and The concavity of i → q i for q ∈ (0, 1] easily shows that {a i } ∞ i=0 is non-increasing. With the choice (71), we then compute
Proof We apply Proposition 4.1 whose assumptions we have verified during the course of the present section. In particular, τ i ≤ τ 0 through the choice (68) that forces ω i ≤ 1. Also, have already derived the rate (70) from (48). Inserting (72) into (70), noting that the former is only valid for N ≥ 2, immediately gives (74).

Examples from Image Processing and the Data Sciences
We now consider several applications of our algorithms. We generally have to consider discretisations, since many interesting infinite-dimensional problems necessitate Banach spaces. Using Bregman distances, it would be possible to generalise our work form Hilbert spaces to Banach spaces, as was done in [24] for the original method of [3]. This is, however, outside the scope of the present work.

Regularised Least Squares
A large range of interesting application problems can be written in the Tikhonov regularisation or empirical loss minimisation form Here α > 0 is a regularisation parameter, G 0 : Z → R typically convex and smooth fidelity term with data f ∈ Z . The forward operator A ∈ L(X ; Z )-which can often also be data-maps our unknown to the space of data. The operator K ∈ L(X ; Y ) and the typically non-smooth and convex F : Y → R act as a regulariser. We are particularly interested in strongly convex G 0 and A with a non-trivial null-space. Examples include, for example, Lasso-a type of regularised regression-with G 0 = x 2 2 /2, K = I , and F(x) = x 1 , on finite-dimensional spaces. If the data of the Lasso is 'sparse', in the sense that A has a non-trivial null-space, then, based on accelerating the strongly convex part of the variable, our algorithm can provide improved convergence rates compared to standard non-accelerated methods.
In image processing examples abound, we refer to [25] for an overview. In total variation (TV) regularisation, we still take F(x) = x 1 , but is K = ∇ the gradient operator. Strictly speaking, this has to be formulated in the Banach space BV( ), but we will consider the discretised setting to avoid this problem. For denoising of Gaussian noise with TV regularisation, we take A = I , and again G 0 = x 2 2 /2. This problem is not so interesting to us, as it is fully strongly convex. In a simple form of TV inpainting-filling in missing regions of an image-we take A as a subsampling operator S mapping an image x ∈ L 2 ( ) to one in L 2 ( \ d ), for d ⊂ the defect region that we want to recreate. Observe that in this case, = S * S is directly a projection operator. This is therefore a problem for our algorithms! Related problems include reconstruction from subsampled magnetic resonance imaging (MRI) data (see, for example, [11,26]), where we take A = SF for F the Fourier transform. Still, A * A is a projection operator, so the problem perfectly suits our algorithms.
Another related problem is total variation deblurring, where A is a convolution kernel. This problem is slightly more complicated to handle, as A * A is not a projection operator. Assuming periodic boundary conditions on a box we can write A = F * â F, multiplying the Fourier transform by someâ ∈ L 2 ( ). If |â| ≥ γ on a subdomain, we obtain a projection form (it would also be possible to extend our theory to non-constant γ , but we have decided not to extend the length of the paper by doing so. Dualisation likewise provides a further alternative).

Satisfaction of convexity conditions
In all of the above examples, when written in the saddle point form (P), F * is a simple pointwise ball constraint. Lemma 2.1 thus guarantees (F * -pcr). If F(x) = x 1 and K = I , then clearly P ⊥x can be bounded in Z = L 1 forx the optimal solution to (75). Thus, for some M > 0, we can add to (75) the artificial constraint In finite dimensions, this gives a bound in L 2 . Lemma 2.1 gives (G-pcr) withγ ⊥ = ∞. In case of our total variation examples, F(x) = x 1 and K = ∇. Provided mean-zero functions are not in the kernel of A, one can through Poincar's inequality [27] on BV( ) and a two-dimensional connected domain ⊂ R 2 show that even the original infinite-dimensional problems have bounded solutions in L 2 ( ). We may therefore again add the artificial constraint (76) with Z = L 2 to (75).

Dynamic bounds and pseudo-duality gaps
We seldom know the exact bound M, but can derive conservative estimates. Nevertheless, adding such a bound to Algorithm 4 is a simple, easily implemented projection of P ⊥ (x i − T i K * y i ) into the constraint set. In practise, we do not use or need the projection, and update the bound M dynamically so as to ensure that the constraint (76) is never active. Indeed, A having a non-trivial nullspace also causes duality gaps for (P) to be numerically infinite. In [28], a 'pseudo-duality gap' was therefore introduced, based on dynamically updating M. We will also use this type of dynamic duality gaps in our reporting.

TGV 2 Regularised Problems
So far, we have considered very simple regularisation terms. Total generalised variation, TGV, was introduced in [29] as a higher-order generalisation of TV. It avoids the unfortunate stair-casing effect of TV-large flat areas with sharp transitions-while preserving the critical edge preservation property that smooth regularisers lack. We concentrate on the second-order TGV 2 . In all of our image processing examples, we can replace TV by TGV 2 .
As with total variation, we have to consider discretised models due the original problem being set in the Banach space BV( ). For two parameters α, β > 0, the regularisation functional is written in the differentiation cascade form of [30] as Here E = (∇ T + ∇)/2 is the symmetrised gradient. With x = (v, w) and y = (y 1 , y 2 ), we may write the problem in the saddle point form (P) with If A = I , as is the case for denoising, we have = γ P for P = I 0 0 0 , perfectly uncoupling in both Algorithm 3 and Algorithm 4 the prox updates for G into ones for G 1 and G 2 . The condition (F * -pcr) withρ = ∞ is then immediate from Lemma 2.1.
Moreover, the Sobolev-Korn inequality [31] allows us to bound on a connected domain ⊂ R 2 an optimalŵ to (77) as for some constant C > 0. We may assume thatw = 0, as the affine part of w is not used in (77). Therefore we may again replace G 2 = 0 by the artificial constraint G 2 (w) = ι · L 2 ≤M (w). By Lemma 2.1, G will then satisfy (G-pcr) withγ ⊥ = ∞.

Numerical Results
We demonstrate our algorithms on TGV 2 denoising and TV deblurring. Our tests are done on the photographs in Fig. 1, both at the original resolution of 768 × 512, and scaled down by a factor of 0.25 to 192 × 128 pixels. It is image #23 from the free Kodak image suite. Other images from the collection that we have experimented on give analogous computational results. For both of our example problems, we calculate a target solution by taking one million iterations of the basic PDHGM (3). We also tried interior point methods for this, but they are only practical for the smaller denoising problem. We evaluate Algorithms 3 and 4 against the standard unaccelerated PDHGM of [3], as well as (a) the mixedrate method of [5], denoted here C-L-O, (b) the relaxed PDHGM of [20,32], denoted here 'Relax', and (c) the adaptive PDHGM of [33], denoted here 'Adapt'. All of these methods are very closely linked and have comparable low costs for each step. This makes them straightforward to compare.
As we have discussed, for comparison and stopping purposes, we need to calculate a pseudo-duality gap as in [28], because the real duality gap is in practise infinite when A has a non-trivial nullspace. We do this dynamically; upgrading, the M in (76) every time, we compute the duality gap. For both of our example problems, we use for simplicity Fig. 2 Step length parameter evolution, both axes logarithmic. 'Alg.3' and 'Alg.4 q=1' have the same parameters as our numerical experiments for the respective algorithms, in particular ζ = τ ⊥,−2 0 for Algorithm 3, which yields constant τ ⊥ . 'Alg.3 ζ /100' uses the value ζ = τ ⊥,−2 0 /100, which causes τ ⊥ to increase for some iterations. 'Alg.4 q=0.1' uses the value q = 0.1 for Algorithm 4, everything else being kept equal In the calculation of the final duality gaps comparing each algorithm, we then take as M the maximum over all evaluations of all the algorithms. This makes the results fully comparable. We always report the duality gap in decibels 10 log 10 (gap 2 /gap 2 0 ) relative to the initial iterate. Similarly, we report the distance to the target solutionû in decibels 10 log 10 ( u i −û 2 / û 2 ), and the primal objective value val(x) := G(x) + F(K x) relative to the target as 10 log 10 (val(x) 2 /val(x) 2 ). Our computations were performed in MATLAB+C-MEX on a MacBook Pro with 16GB RAM and a 2.8 GHz Intel Core i5 CPU. TGV 2 denoising The noise in our high-resolution test image, with values in the range [0, 255], has standard deviation 29.6 or 12 dB. In the downscaled image, these become, respectively, 6.15 or 25.7 dB. As parameters (β, α) of the TGV 2 regularisation functional, we choose (4.4, 4) for the downscale image, and translate this to the original image by multiplying by the scaling vector (0.25 −2 , 0.25 −1 ) corresponding to the 0.25 downscaling factor. See [34] for a discussion about rescaling and regularisation factors, as well as for a justification of the β/α ratio.
For the PDHGM and our algorithms, we take γ = 0.5, corresponding to the gap convergence results. We choose δ = 0.01, and parametrise the PDHGM with σ 0 = 1.9/ K and τ * 0 = τ 0 ≈ 0.52/ K solved from τ 0 σ 0 = (1 − δ) K 2 . These are values that typically work well. For forwarddifferences discretisation of TGV 2 with cell width h = 1, we have K 2 ≤ 11.4 [28]. We use the same value of δ for Algorithm 3 and Algorithm 4, but choose τ ⊥ 0 = 3τ * 0 , and τ 0 = τ 0 = 80τ * 0 . We also take ζ = τ ⊥,−2 0 for Algorithm 3. These values have been found to work well by trial and error, while keeping δ comparable to the PDHGM. A similar choice of τ 0 with a corresponding modification of σ 0 would significantly reduce the performance of the PDHGM. For Algorithm 4, we take exponent q = 0.1 for the sequence {a i }. This gives in principle a mixed O(1/N 1.5 )+ O(1/N 0.5 ) rate, possibly improved by the convergence of the dual sequence. We plot the evolution of the step length for these and some other choices in Fig. 2. For the C-L-O, we use the detailed parametrisation from [35,Corollary 2.4], taking as Y the true L 2 -norm Bregman divergence of B(0, α)× B(0, β), and X = 10 · f 2 /2 as a conservative estimate of a ball containing the true solution. For 'Adapt', we use the exact choices of α 0 , η, and c from [33]. For 'Relax', we use the value 1.5 for the inertial ρ parameter of [32]. For both of these algorithms, we use the same choices of σ 0 and τ 0 as for the PDHGM.
We take fixed 20,000 iterations and initialise each algorithm with y 0 = 0 and x 0 = 0. To reduce computational overheads, we compute the duality gap and distance to target only every 10 iterations instead of at each iteration. The results are in Fig. 3 and Table 1. As we can see, Algorithm 3 performs extremely well for the low-resolution image, especially in its initial iterations. After about 700 or 200 iterations, depending on the criterion, the standard and relaxed PDHGM start to overtake. This is a general effect that we have seen in our tests: the standard PDHGM performs in practise very well asymptotically, although in principle all that exists is a O(1/N ) rate on the ergodic The plot is logarithmic, with the decibels calculated as in Sect. 5.3. The poor high-resolution results for 'Adapt' [33] have been omitted to avoid poor scaling of the plots. a Gap, low resolution, b target, low resolution, c value, low resolution, d gap, high resolution, e target, high resolution, f value, high resolution Adapt ------ The CPU time and number of iterations (at a resolution of 10) needed to reach given solution quality in terms of the duality gap, distance to target, or primal objective value duality gap. Algorithm 4, by contrast, does not perform asymptotically so well. It can be extremely fast on its initial iterations, but then quickly flattens out. The C-L-O surprisingly performs better on the high-resolution image than on the low-resolution image, where it does somewhat poorly in comparison with the other algorithms. The adaptive PDHGM performs very poorly for TGV 2 denoising, and we have indeed excluded the high-resolution results from our reports to keep the scaling of the plots informative. Overall, Algorithm 3 gives good results fast, although the basic and relaxed PDHGM seems to perform, in practise, better asymptotically. TV deblurring Our test image has now been distorted by Gaussian blur of standard deviation 4, which we intent to remove. We denote byâ the Fourier presentation of the blur operator as discussed in Sect. 5.1. For numerical stability of the pseudo-duality gap, we zero out small entries, replacing thisâ byâχ |â( · )|≥ â ∞ /1000 (ξ ). Note that this is only needed for the stable computation of G * for the pseudo-duality gap, to compare the algorithms; the algorithms themselves are stable without this modification. To construct the projection operator P, we then setp(ξ ) = χ |â( · )|≥0.3 â ∞ (ξ ), and P = F * p F.
We use TV parameter 2.55 for the high-resolution image and the scaled parameter 2.55 * 0.15 for the low-resolution image. We parametrise all the algorithms almost exactly as  The CPU time and number of iterations (at a resolution of 10) needed to reach given solution quality in terms of the duality gap, distance to target, or primal objective value TGV 2 denoising above, of course with appropriate U and K 2 ≤ 8 corresponding to K = ∇ [36]. The only difference in parameterisation is that we take q = 1 instead of q = 0.1 for Algorithm 4.
The results are in Fig. 4 and Table 2. It does not appear numerically feasible to go significantly below −100 or −80 dB gap. Our guess is that this is due to the numerical inaccuracies of the fast Fourier transform implementation in MATLAB. The C-L-O performs very well judged by the duality gap, although the images themselves and the primal objective value appear to take a little bit longer to converge. The relaxed PDHGM is again slightly improved from the standard PDHGM. The adaptive PDHGM performs very well, slightly outperforming Algorithm 3, although not Algorithm 4. This time Algorithm 4 performs remarkably well.

Conclusion
To conclude, overall, our algorithms are very competitive within the class of proposed variants of the PDHGM. Within our analysis, we have, moreover, proposed very streamlined derivations of convergence rates for even the standard PDHGM, based on the proximal point formulation and the idea of testing. Interesting continuations of this study include whether the conditionT i K = K T i can reasonably be relaxed such thatT i and T i would not have to be scalars, as well as the relation to block coordinate descent methods, in particular [14,37].