Variable smoothing for convex optimization problems using stochastic gradients

TL;DR

This paper introduces novel stochastic gradient algorithms for structured convex optimization problems with nonsmooth functions, leveraging variable smoothing via the Moreau envelope, and demonstrates their effectiveness in image denoising and deblurring.

Contribution

It proposes new stochastic algorithms using variable smoothing for convex problems with Lipschitz continuous nonsmooth functions, extending primal-dual methods with regularization techniques.

Findings

Effective in large-scale problems

Applicable to total variational denoising

Improves convergence with stochastic smoothing

Abstract

We aim to solve a structured convex optimization problem, where a nonsmooth function is composed with a linear operator. When opting for full splitting schemes, usually, primal-dual type methods are employed as they are effective and also well studied. However, under the additional assumption of Lipschitz continuity of the nonsmooth function which is composed with the linear operator we can derive novel algorithms through regularization via the Moreau envelope. Furthermore, we tackle large scale problems by means of stochastic oracle calls, very similar to stochastic gradient techniques. Applications to total variational denoising and deblurring are provided.

Equations

$$\underset{x\in \mathcal{H}}{\min }f(x)+g(Kx),$$

$$\underset{x\in \mathcal{H}}{\min }f(x)+\sum _{i=1}^m{g}_i(K_{i}x),$$

$$g^{\ast }(x):=\underset{p\in \mathcal{H}}{\sup }\left\{⟨x,p⟩-g(p)\right\}\forall x\in \mathcal{H}.$$

$$\text{leftidx}{}^{\mu }g(\cdot ):={\left(g^{\ast }+\frac{\mu }{2}\Vert \cdot {\Vert }^2\right)}^{\ast }(\cdot )=\underset{p\in \mathcal{H}}{\sup }\left\{⟨\cdot ,p⟩-g^{\ast }(p)-\frac{\mu }{2}\Vert p{\Vert }^2\right\}.$$

$$\nabla (\text{leftidx}{}^{\mu }g)(x)=\frac{1}{\mu }\left(x-{\mathrm{prox}}_{\mu g}\left(x\right)\right)={\mathrm{prox}}_{\frac{1}{\mu }{g}^{\ast }}\left(\frac{x}{\mu }\right)\forall x\in \mathcal{H}$$

$$x-\mu \nabla (\text{leftidx}{}^{\mu }g)(x)={\mathrm{prox}}_{\mu g}\left(x\right)\forall x\in \mathcal{H}.$$

$$\nabla \left(\text{leftidx}{}^{\mu }g\circ K\right)(x)=K^{\ast }\nabla \left(\text{leftidx}{}^{\mu }g\right)(Kx)=\frac{1}{\mu }{K}^{\ast }\left(Kx-{\mathrm{prox}}_{\mu g}\left(Kx\right)\right)=K^{\ast }{\mathrm{prox}}_{\frac{1}{\mu }{g}^{\ast }}\left(\frac{Kx}{\mu }\right),$$

$$\text{dom}{g}^{\ast }\subseteq B(0,L_g),$$

$$g^{{\mu }_2}(x)\le g^{{\mu }_1}(x)\le g^{{\mu }_2}(x)+({\mu }_2-{\mu }_1)\frac{L_g^2}{2}.$$

$$\text{leftidx}{}^{{\mu }_1}g(x)=$$

$$\le$$

$$\text{leftidx}{}^{\mu }g(x)\le g(x)\le \text{leftidx}{}^{\mu }g(x)+\mu \frac{L_g^2}{2}.$$

$$\underset{p\in \mathcal{H}}{\mathrm{argmax}}\left\{⟨\cdot ,p⟩-g^{\ast }(p)-\frac{\mu }{2}\Vert p{\Vert }^2\right\}=\nabla \text{leftidx}{}^{\mu }g(\cdot ).$$

$$\underset{p\in \mathcal{H}}{\mathrm{argmax}}\left\{⟨x,p⟩-g^{\ast }(p)-\frac{\mu }{2}\Vert p{\Vert }^2\right\}=$$

$$=$$

$$\mu \mapsto \text{leftidx}{}^{\mu }g(x).$$

$$\frac{\partial }{\partial \mu }\text{leftidx}{}^{\mu }g(x)=-\frac{1}{2}\Vert \nabla \text{leftidx}{}^{\mu }g(x){\Vert }^2\forall x\in \mathcal{H}\forall \mu \in (0,+\infty ).$$

$$\frac{\partial }{\partial \mu }\text{leftidx}{}^{\mu }g(x)=$$

$$=$$

$$\text{leftidx}{}^{{\mu }_1}g(x)\le \text{leftidx}{}^{{\mu }_2}g(x)+({\mu }_2-{\mu }_1)\frac{1}{2}\Vert \nabla \text{leftidx}{}^{{\mu }_1}g(x){\Vert }^2.$$

$$\text{leftidx}{}^{{\mu }_2}g(x)\ge$$

$$=$$

$$\text{leftidx}{}^{\mu }g(x)+⟨\nabla \text{leftidx}{}^{\mu }g(x),y-x⟩\le g(y)-\frac{\mu }{2}\Vert \nabla \text{leftidx}{}^{\mu }g(x){\Vert }^2.$$

$$\text{leftidx}{}^{\mu }g(x)+⟨\nabla \text{leftidx}{}^{\mu }g(x),y-x⟩=$$

$$=$$

$$\le$$

$$=$$

$$\text{leftidx}{}^{\mu }g(Kx)+⟨(\nabla \text{leftidx}{}^{\mu }g\circ K)(x),y-x⟩\le g(Ky)-\frac{\mu }{2}\Vert \nabla \text{leftidx}{}^{\mu }g(Kx){\Vert }^2\forall x,y\in \mathcal{H}.$$

$$h(x)+⟨\nabla h(x),y-x⟩\le h(y)-\frac{1}{2L_h}\Vert \nabla h(x)-\nabla h(y){\Vert }^2\forall x,y\in \mathcal{H}.$$

$$\text{leftidx}{}^{\mu }g(Kx)+⟨\nabla (\text{leftidx}{}^{\mu }g\circ K)(x),y-x⟩\le \text{leftidx}{}^{\mu }g(Ky)-\frac{1}{2\mu }\Vert \nabla \text{leftidx}{}^{\mu }g(Kx)-\nabla \text{leftidx}{}^{\mu }g(Ky){\Vert }^2\forall x,y\in \mathcal{H}.$$

$$(1-\alpha )\Vert x-y{\Vert }^2+\alpha \Vert y{\Vert }^2\ge \alpha (1-\alpha )\Vert x{\Vert }^2.$$

$$\underset{x\in \mathcal{H}}{\min }F(x):=f(x)+g(Kx),$$

Full text

Keywords. structured convex optimization problem, variable smoothing algorithm, convergence rate, stochastic gradients

AMS Subject Classification. 90C25, 90C15, 65Y20

1 Introduction

The problem at hand is the following structured convex optimization problem

[TABLE]

for real Hilbert spaces $\mathcal{H}$ and $\mathcal{G}$, $f:\mathcal{H}\to\mkern 1.5mu\overline{\mkern-1.5mu\mathbb{R}\mkern-1.5mu}\mkern 1.5mu$ a proper, convex and lower semicontinuous function, $g:\mathcal{G}\to\mathbb{R}$ a, possibly nonsmooth, convex and Lipschitz continuous function, and $K:\mathcal{H}\to\mathcal{G}$ a linear continuous operator.

Our aim will be to devise an algorithm for solving (1) following the full splitting paradigm (see[15, 9, 5, 25, 16, 8, 6]). In other words, we allow only proximal evaluations for simple nonsmooth functions, but no proximal evaluations for compositions with linear continuous operators, like, for instance, for $g\circ K$.

We will accomplish this feat by the means of a smoothing strategy, which, for the purpose of this paper, means, making use of the Moreau-Yosida approximation (see [11, 10, 7, 19, 17]). The approach can be described as follows: we “smooth” $g$, i.e. we replace it by its Moreau envelope, and solve the resulting optimization problem by an accelerated proximal-gradient algorithm (see [13, 3, 18]).

The only other family of methods able to solve problems of type (1) are the so called primal-dual algorithms, first and foremost the primal-dual hybrid gradient (PDHG) introduced in [15]. In comparison, this method does not need the Lipschitz continuity of $g$ in order to proof convergence. However, in this very general case, convergence rates can only be shown for the so-called restricted primal-dual gap function. In order to derive from here convergence rates for the primal objective function, either Lipschitz continuity of $g$ or finite dimensionality of the problem plus the condition that $g$ must have full domain are necessary (see, for instance, [5, Theorem 9]). This means, that for infinite dimensional problems the assumptions required by both, PDHG and our method, for deriving convergence rates for the primal objective function are in fact equal, but for finite dimensional problems the assumption of PDHG are weaker. In either case, however, we are able to proof these rates for the sequence of iterates ${(x_{k})}_{k\geq 1}$ itself whereas PDHG only has them for the sequence of so-called ergodic iterates, i.e. ${(\frac{1}{k}\sum_{i=1}^{k}x_{i})}_{k\geq 1}$, which is naturally undesirable as the averaging slows the convergence down.

Furthermore, we will also consider the case where only a stochastic oracle of the proximal operator of $g$ is available to us. This setup corresponds e.g. to the case where the objective function is given as

[TABLE]

where, for $i=1,\dots,m$, $\mathcal{G}_{i}$ are real Hilbert spaces, $g_{i}:\mathcal{G}_{i}\rightarrow\mathbb{R}$ are convex and Lipschitz continuous functions and $K_{i}:\mathcal{H}\rightarrow\mathcal{G}_{i}$ are linear continuous operators, but the number of summands being large we wish to not compute all proximal operators of all $g_{i},i=1,\dots,m$, for purpose of making iterations cheaper to compute.

For the finite sum case (2), there exist algorithms of similar spirit such those in [14, 21]. Some algorithms do in fact deal with a similar setup of stochastic gradient like evaluations, see [23], but only for smooth terms in the objective function.

In Section 2 we will cover the preliminaries about the Moreau-Yosida envelope as well as useful identities and estimates connected to it. In Section 3 we will deal with the deterministic case and prove a convergence rate of $\mathcal{O}(\frac{1}{k})$ for the function values at the iterates. Next up, in Section 4, we will consider the stochastic case as described above and prove a convergence rate of $\mathcal{O}\left(\frac{\log(k)}{\sqrt{k}}\right)$. Last but not least, we will look at some numerical examples in image processing in Section 5.

It is important to note that the proof for the deterministic setting differs surprisingly from the one for the stochastic setting. The technique for the stochastic setting is less refined in the sense that there is no coupling between the smoothing parameter and the extrapolation parameter. Where as this technique works also works for the deterministic setting it gives a worse convergence rate of $\mathcal{O}\left(\frac{\log{k}}{k}\right)$. The tight coupling of the two sequences of parameters, however does not work in the proof of the stochastic algorithm as it does not allow for the particular choice of the smoothing parameters needed there.

2 Preliminaries

Definition 2.1**.**

For a proper, convex and lower semicontinuous function $g:\mathcal{H}\rightarrow\mkern 1.5mu\overline{\mkern-1.5mu\mathbb{R}\mkern-1.5mu}\mkern 1.5mu$, its convex conjugate is denoted by $g^{*}$ defined as a function from $\mathcal{H}$ to $\mkern 1.5mu\overline{\mkern-1.5mu\mathbb{R}\mkern-1.5mu}\mkern 1.5mu$, given by

[TABLE]

As mentioned in the introduction, we want to smooth a nonsmooth function by considering its Moreau envelope. The next definition will clarify exactly what object we are talking about.

Definition 2.2**.**

For a proper, convex and lower semicontinuous function $g:\mathcal{H}\rightarrow\mkern 1.5mu\overline{\mkern-1.5mu\mathbb{R}\mkern-1.5mu}\mkern 1.5mu$, its Moreau envelope with the parameter $\mu\geq 0$ is defined as a function from $\mathcal{H}$ to $\mathbb{R}$, given by

[TABLE]

From this definition, however, it is not completely evident that the Moreau envelope indeed fulfills its purpose in being a smooth representation of the original function. The next lemma will remedy this fact.

Lemma 2.1** (see [2, Proposition 12.29]).**

Let $g:\mathcal{H}\to\mkern 1.5mu\overline{\mkern-1.5mu\mathbb{R}\mkern-1.5mu}\mkern 1.5mu$ be a proper, convex and lower semicontinuous function and $\mu>0$. Then its Moreau envelope is Fréchet differentiable on $\mathcal{H}$. In particular, the gradient itself is given by

[TABLE]

and is $\mu^{-1}$-Lipschitz continuous.

In particular, for all $\mu>0$, a gradient step with respect to the Moreau envelope corresponds to a proximal step

[TABLE]

The previous lemma establishes two things. Not only does it clarify the smoothness of the Moreau envelope, but it also gives a way of computing its gradient. Obviously, a smooth representation whose gradient we would not be able to compute would not be any good.

As mentioned in the introduction, we want to smooth the nonsmooth summand of the objective function which is composed with the linear operator as this can be considered the crux of problem (1). The function $g\circ K$ will be smoothed via considering instead $\leftidx{{}^{\mu}}g\circ K:\mathcal{H}\to\mathbb{R}$. Clearly, by the chain rule, this function is continuously differentiable with gradient given for every $x\in\mathcal{H}$ by

[TABLE]

and is thus Lipschitz continuous with Lipschitz constant $\frac{\lVert K\rVert^{2}}{\mu}$.

Lipschitz continuity will play an integral role in our investigations, as can be seen by the following lemmas.

Lemma 2.2** (see [4, Proposition 4.4.6]).**

Let $g:\mathcal{H}\to\mathbb{R}$ be a convex and $L_{g}$-Lipschitz continuous function. Then, the domain of its Fenchel conjugate is bounded, i.e.

[TABLE]

where $B(0,L_{g})$ denotes the open ball with radius $L_{g}$ around the origin.

Lemma 2.3**.**

Let $g:\mathcal{H}\to\mathbb{R}$ be a convex and $L_{g}$-Lipschitz continuous function. For $\mu_{2}\geq\mu_{1}\geq 0$ and every $x\in\mathcal{H}$ it holds

[TABLE]

Proof.

For $\mu_{2}\geq\mu_{1}\geq 0$ and every $x\in\mathcal{H}$ we have, by definition,

[TABLE]

where we used Lemma 2.2 in the last inequality. ∎

Lemma 2.4**.**

For $\mu\geq 0$ and every $x\in\mathcal{H}$ it holds that

[TABLE]

Proof.

The statement follows from Lemma 2.3 using the fact that $\leftidx{{}^{0}}g(x)=g(x)$. ∎

Lemma 2.5**.**

The maximizing argument in the definition of the Moreau-Yosida envelope is given by its gradient, i.e. for $\mu>0$ it holds that

[TABLE]

Proof.

Let $x\in\mathcal{H}$ be fixed. It holds

[TABLE]

and the conclusion follows by using Lemma 2.1. ∎

Lemma 2.6**.**

For a proper, convex and lower semicontinuous function $g:\mathcal{H}\to\mkern 1.5mu\overline{\mkern-1.5mu\mathbb{R}\mkern-1.5mu}\mkern 1.5mu$ and every $x\in\mathcal{H}$ we can consider the mapping from $(0,+\infty)$ to $\mathbb{R}$ given by

[TABLE]

This mapping is convex and differentiable and its derivative is given by

[TABLE]

Proof.

Let $x\in\mathcal{H}$ be fixed. From the definition of the Moreau-Yosida envelope we can see that the mapping given in (14) is a pointwise supremum of functions which are linear in $\mu$. It is therefore convex. Furthermore, since the objective function is strongly concave, this supremum is uniquely attained at $\nabla\leftidx{{}^{\mu}}g(x)=\operatorname*{arg\,max}_{p\in\mathcal{H}}\left\{\left\langle x,p\right\rangle-g^{*}(p)-\frac{\mu}{2}\lVert p\rVert^{2}\right\}$. According to the Danskin Theorem, the function $\mu\mapsto\leftidx{{}^{\mu}}g(x)$ is differentiable and its gradient is given by

[TABLE]

∎

Lemma 2.7**.**

For $\mu_{1},\mu_{2}>0$ and every $x\in\mathcal{H}$ it holds

[TABLE]

Proof.

Let $x\in\mathcal{H}$ be fixed. Via Lemma 2.6 we know that the map $\mu\mapsto\leftidx{{}^{\mu}}g(x)$ is convex and differentiable. We can therefore use the gradient inequality to deduce that

[TABLE]

which is exactly the statement of the lemma. ∎

Lemma 2.8**.**

For $\mu>0$ and every $x,y\in\mathcal{H}$ we have that

[TABLE]

Proof.

Using Lemma 2.5 and the definition of the Moreau-Yosida envelope we get that

[TABLE]

∎

In the convergence proof of Section 3 we will need the inequality in the above lemma at the points $x:=Kx$ and $y:=Ky$, namely

[TABLE]

The following lemma is a standard result for convex and Fréchet differentiable functions.

Lemma 2.9** (see [20]).**

For a convex and Fréchet differentiable function $h:\mathcal{H}\to\mathbb{R}$ with $L_{h}$-Lipschitz continuous gradient we have that

[TABLE]

By applying Lemma 2.9 with $h:=\leftidx{{}^{\mu}}g$, and $x:=Kx$ and $y:=Ky$, we obtain

[TABLE]

The following technical result will be used in the proof of the convergence statement.

Lemma 2.10**.**

For $\alpha\in(0,1)$ and every $x,y\in\mathcal{H}$ we have that

[TABLE]

3 Deterministic Method

Problem 3.1**.**

The problem at hand reads

[TABLE]

for a proper, convex and lower semicontinuous function $f:\mathcal{H}\to\mkern 1.5mu\overline{\mkern-1.5mu\mathbb{R}\mkern-1.5mu}\mkern 1.5mu$, a convex and $L_{g}$-Lipschitz continuous $(L_{g}>0)$ function $g:\mathcal{G}\to\mathbb{R}$, and a nonzero linear continuous operator $K:\mathcal{H}\to\mathcal{G}$.

The idea of the algorithm which we propose to solve (1) is to smooth $g$ and then to solve the resulting problem by means of an accelerated proximal-gradient.

Algorithm 3.1** (Variable Accelerated SmooThing (VAST)).**

Let $y_{0}=x_{0}\in\mathcal{H},{(\mu_{k})}_{k\geq 0}\!\subseteq(0,+\infty)$, and ${(t_{k})}_{k\geq 1}$ a sequence of real numbers with $t_{1}=1$ and $t_{k}\geq 1$ for every $k\geq 2$. Consider the following iterative scheme

[TABLE]

Remark 3.1*.*

The assumption $t_{1}=1$ can be removed but guarantees easier computation and is also in line with classical choices of ${(t_{k})}_{k\geq 1}$ in [18, 13].

Remark 3.2*.*

The sequence ${(u_{k})}_{k\geq 1}$ given by

[TABLE]

despite not appearing in the algorithm, will feature a prominent role in the convergence proof. Due to the convention $t_{1}=1$ we have that

[TABLE]

We also denote

[TABLE]

The next theorem is the main result of this section and it will play a fundamental role when proving a convergence rate of $\mathcal{O}(\frac{1}{k})$ for the sequence ${(F(x_{k}))}_{k\geq 0}$. Similar convergence rate results for smoothing algorithms have been obtained in [11, 10, 7]. The techniques used in this section, however are in flavor of [24], with the most notable difference of a much easier choice of the involved parameters.

Theorem 3.1**.**

Consider the setup of Problem 3.1 and let ${(x_{k})}_{k\geq 0}$ and ${(y_{k})}_{k\geq 0}$ be the sequences generated by Algorithm 3.1. Assume that for every $k\geq 1$

[TABLE]

and

[TABLE]

Then, for every optimal solution $x^{*}$ of Problem 3.1, it holds

[TABLE]

The proof of this result relies on several partial results which we will prove as follows.

Lemma 3.1**.**

The following statement holds for every $z\in\mathcal{H}$ and every $k\geq 0$

[TABLE]

Proof.

Let $k\geq 0$ be fixed. Since, by the definition of the proximal map, $x_{k+1}$ is the minimizer of a $\frac{1}{\gamma_{k+1}}$-strongly convex function we know that for every $z\in\mathcal{H}$

[TABLE]

Next we use the $L_{k+1}$-smoothness of $\leftidx{{}^{\mu_{k+1}}}g\circ K$ and the fact that $\frac{1}{\gamma_{k+1}}\geq L_{k+1}$ to deduce

[TABLE]

∎

Lemma 3.2**.**

Let $x^{*}$ be an optimal solution of Problem 3.1. Then it holds

[TABLE]

Proof.

We use the gradient inequality to deduce that for every $z\in\mathcal{H}$ and every $k\geq 0$

[TABLE]

and plug this into the statement of Lemma 3.1 to conclude that

[TABLE]

For $k=0$ we get that

[TABLE]

Now we us the fact that $u_{1}=x_{1}$ and $y_{0}=x_{0}$ to obtain the conclusion. ∎

Lemma 3.3**.**

Let $x^{*}$ be an optimal solution of Problem 3.1. The following descent-type inequality holds for every $k\geq 0$

[TABLE]

Proof.

Let be $k\geq 0$ fixed. We apply Lemma 3.1 with $z:=\left(1-\frac{1}{t_{k+1}}\right)x_{k}+\frac{1}{t_{k+1}}x^{*}$ to deduce that

[TABLE]

Using the convexity of $f$ gives

[TABLE]

Now, we use (21) to deduce that

[TABLE]

and (23) to conclude that

[TABLE]

Combining (38), (39) and (40) gives

[TABLE]

The first term on the right hand side is $\leftidx{{}^{\mu_{k+1}}}g(Kx_{k})$ but we would like it to be $\leftidx{{}^{\mu_{k}}}g(Kx_{k})$. Therefore we use Lemma 2.7 to deduce that

[TABLE]

Next up we want to estimate all the norms of gradients by using Lemma 2.10 which says that

[TABLE]

Combining (42) and (43) gives

[TABLE]

Now we combine the two terms containing $\lVert\nabla\leftidx{{}^{\mu_{k+1}}}g(Kx_{k})\rVert^{2}$ and get that

[TABLE]

By subtracting $F(x^{*})=f(x^{*})+g(Kx^{*})$ on both sides we finally obtain

[TABLE]

∎

Now we are in the position to prove Theorem 3.1.

Proof of Theorem 3.1.

We start with the statement of Lemma 3.3 and use the assumption that

[TABLE]

to make the last term in the inequality disappear for every $k\geq 0$

[TABLE]

Now we use the assumption that

[TABLE]

to get that for every $k\geq 0$

[TABLE]

Let be $N\geq 2$. Summing (50) from $k=1$ to $N-1$ and getting rid of the nonnegative term $\lVert u_{N}-x^{*}\rVert^{2}$ gives

[TABLE]

Since $t_{1}=1$, the above inequality is fulfilled also for $N=1$. Using Lemma 3.2 shows that

[TABLE]

The above inequality, however, is still in terms of the smoothed objective function. In order to go to the actual objective function we apply Lemma 2.4 and deduce that

[TABLE]

∎

Corollary 3.1**.**

There exists a choice of parameters ${(\mu_{k})}_{k\geq 1},{(t_{k})}_{k\geq 1},{(\gamma_{k})}_{k\geq 1}$ such that

[TABLE]

and

[TABLE]

are for every $k\geq 1$ fulfilled and are e.g. given by

[TABLE]

and for every $k\geq 1$

[TABLE]

For this choice of the parameters we have that

[TABLE]

Proof.

Since $\gamma_{k}$ and $\mu_{k}$ are a scalar multiple of each other (55) is equivalent to

[TABLE]

and further to (by taking into account that $t_{k+1}>1$ for every $k\geq 1$)

[TABLE]

Our update choice in (57) for the sequence ${(\mu_{k})}_{k\geq 1}$ is exactly chosen in such a way that it satisfies this. Plugging (60) into (54) gives for every $k\geq 1$ the condition

[TABLE]

which is equivalent to

[TABLE]

and further to

[TABLE]

Plugging in $t_{k+1}=\sqrt{t_{k}^{2}+2t_{k}}$ we get that this equivalent to

[TABLE]

which is evidently fulfilled. Thus, the choices in (57) are indeed feasible for our algorithm.

Now we want to prove the claimed convergence rates. Via induction we show that

[TABLE]

Evidently, this holds for $t_{1}=1$. Assuming that (65) holds for $k\geq 1$, we easily see that

[TABLE]

and, on the other hand,

[TABLE]

In the following we prove a similar estimate for the sequence ${(\mu_{k})}_{k\geq 1}$. To this end we show, again by induction, the following recursion for every $k\geq 2$

[TABLE]

For $k=2$ this follows from the definition (60). Assume now that (68) holds for $k\geq 2$. From here we have that

[TABLE]

Using (68) together with (65) we can check that for every $k\geq 1$

[TABLE]

where we used in the last step the fact that $t_{k+1}\leq t_{k}+1$.

The last thing to check is the fact that $\mu_{k}$ goes to zero like $\frac{1}{k}$. First we check that for every $k\geq 1$

[TABLE]

This can be seen via

[TABLE]

By bringing $t_{k+1}$ to the other side we get that

[TABLE]

from which we can deduce (71) by dividing by $t_{k+1}^{2}-t_{k+1}$.

We plug in the estimate (71) in (68) and get for every $k\geq 2$

[TABLE]

With the above inequalities we can to deduce the claimed convergence rates. First note that from Theorem 3.1 we have

[TABLE]

Now, in order to obtain the desired conclusion, we used the above estimates and deduce for every $N\geq 1$

[TABLE]

where we used that

[TABLE]

∎

Remark 3.3*.*

Consider the choice (see [18])

[TABLE]

and

[TABLE]

Since

[TABLE]

we see that in this setting we have to choose

[TABLE]

Thus, the sequence of optimal function values ${(F(x_{N}))}_{N\geq 1}$ approaches a $b\lVert K\rVert^{2}\frac{L_{g}}{2}$-approximation of the optimal objective value $F(x^{*})$ with a convergence rate of $\mathcal{O}(\frac{1}{N^{2}})$, i.e.

[TABLE]

4 Stochastic Method

Problem 4.1**.**

The problem is the same as in the deterministic case

[TABLE]

other than the fact that at each iteration we are only given a stochastic estimator of the quantity

[TABLE]

For the stochastic quantities arising in this section we will use the following notation. For every $k\geq 0$, we denote by $\sigma(x_{0},\dots,x_{k})$ the smallest $\sigma$-algebra generated by the family of random variables $\{x_{0},\dots,x_{k}\}$ and by $\mathbb{E}_{k}(\cdot):=\mathbb{E}(\cdot|\sigma(x_{0},\dots,x_{k}))$ the conditional expectation with respect to this $\sigma$-algebra.

Algorithm 4.1** (stochastic Variable Accelerated SmooThing (sVAST)).**

Let $y_{0}=x_{0}\in\mathcal{H},{(\mu_{k})}_{k\geq 1}$ a sequence of positive and nonincreasing real numbers, and ${(t_{k})}_{k\geq 1}$ a sequence of real numbers with $t_{1}=1$ and $t_{k}\geq 1$ for every $k\geq 2$. Consider the following iterative scheme

[TABLE]

where we make the standard assumptions about our gradient estimator of being unbiased, i.e.

[TABLE]

and having bounded variance

[TABLE]

for every $k\geq 0$.

Note that we use the same notations as in the deterministic case

[TABLE]

Lemma 4.1**.**

Let $g:\mathcal{H}\to\mathbb{R}$ be a convex and $L_{g}$-Lipschitz continuous function. Then its Moreau envelope $\leftidx{{}^{\mu}}g$ satisfies for every $x,y\in\mathcal{H}$

[TABLE]

Proof.

From Lemma 2.4 we deduce for every $x,y\in\mathcal{H}$

[TABLE]

and, analogously,

[TABLE]

The statement of the lemma follows. ∎

Lemma 4.2**.**

The following statement holds for every $z\in\mathcal{H}$ and every $k\geq 0$

[TABLE]

Proof.

Here we have to proceed a little bit different from Lemma 3.1. Namely, we have to treat the gradient step and the proximal step differently. For this purpose we define the auxiliary variable

[TABLE]

Let be $k\geq 1$ fixed. From the gradient step we get

[TABLE]

Taking the conditional expectation gives

[TABLE]

Using the gradient inequality we deduce

[TABLE]

and therefore

[TABLE]

Also from the smoothness of $(\leftidx{{}^{\mu_{k}}}g\circ K)$ we deduce via the Descent Lemma that

[TABLE]

Plugging in the definition of $z_{k}$ and using the fact that $L_{k}=\frac{1}{\gamma_{k}}$ we get

[TABLE]

Now we take the conditional expectation to deduce that

[TABLE]

Multiplying (97) by $\gamma_{k}$ and adding it to (94) gives

[TABLE]

Now we use the assumption about the bounded variance to deduce that

[TABLE]

Next up for the proximal step we deduce

[TABLE]

Taking the conditional expectation and combining (99) and (100) we get

[TABLE]

From here, using now Lemma 4.1, we get that

[TABLE]

Now we use

[TABLE]

to obtain that

[TABLE]

∎

Lemma 4.3**.**

Let $x^{*}$ be an optimal solution of Problem 4.1. Then it holds

[TABLE]

Proof.

Applying the previous lemma with $k=0$ and $z=x^{*}$, we get that

[TABLE]

Therefore, using the fact that $y_{0}=x_{0}$ and $u_{1}=x_{1}$,

[TABLE]

which finishes the proof. ∎

Theorem 4.1**.**

Consider the setup of Problem 4.1 and let ${(x_{k})}_{k\geq 0}$ and ${(y_{k})}_{k\geq 0}$ denote the sequences generated by Algorithm 4.1. Assume that for all $k\geq 1$

[TABLE]

Then, for every optimal solution $x^{*}$ of Problem 4.1, it holds

[TABLE]

Proof of Theorem 4.1.

Let be $k\geq 0$ fixed. Lemma 4.2 for $z:=\left(1-\frac{1}{t_{k+1}}\right)x_{k}+\frac{1}{t_{k+1}}x^{*}$ gives

[TABLE]

From here and from the convexity of $F^{k+1}$ follows

[TABLE]

Now, by multiplying both sides with by $t_{k+1}^{2}$, we deduce

[TABLE]

Next, by adding $t_{k}^{2}(F^{k+1}(x_{k})-F^{k+1}(x^{*}))$ on both sides of (109), gives

[TABLE]

Utilizing (9) together with the assumption that ${(\mu_{k})}_{k\geq 1}$ is nonincreasing leads to

[TABLE]

Now, using that $t_{k}^{2}\geq t_{k+1}^{2}-t_{k+1}$, we get

[TABLE]

Multiplying both sides with $\gamma_{k+1}$ and putting all terms on the correct sides yields

[TABLE]

At this point we would like to discard the term $\gamma_{k+1}\rho_{k+1}(F^{k+1}(x_{k})-F^{k+1}(x^{*}))$ which we currently cannot as the positivity of $F^{k+1}(x_{k})-F^{k+1}(x^{*})$ is not ensured. So we add $\gamma_{k+1}\rho_{k+1}\mu_{k+1}\frac{L_{g}^{2}}{2}$ on both sides of (110) and get

[TABLE]

Using again (9) to deduce that

[TABLE]

we can now discard said term from (111), giving

[TABLE]

Last but not least we use the that $F^{k}(x_{k})-F^{k}(x^{*})+\mu_{k}\frac{L_{g}^{2}}{2}\geq F(x_{k})-F(x^{*})\geq 0$ and $\gamma_{k+1}\leq\gamma_{k}$ to follow that

[TABLE]

Combining (113) and (114) yields

[TABLE]

Let be $N\geq 2$. We take the expected value on both sides (115) and sum from $k=1$ to $N-1$. Getting rid of the non-negative terms $\lVert u_{N}-x^{*}\rVert^{2}$ gives

[TABLE]

Since $t_{1}=1$, the above inequality holds also for $N=1$. Now, using Lemma 4.3 we get that for every $N\geq 1$

[TABLE]

From Lemma 2.4 we follow that

[TABLE]

therefore, for every $N\geq 1$

[TABLE]

By using the fact that $\gamma_{k}=\frac{1}{L_{k}}=\frac{\mu_{k}}{\lVert K\rVert^{2}}$ for every $k\geq 1$ gives

[TABLE]

Thus,

[TABLE]

∎

Corollary 4.1**.**

Let

[TABLE]

and, for $b>0$,

[TABLE]

Then,

[TABLE]

Furthermore, we have that $F(x_{N})$ converges almost surely to $F(x^{*})$ as $N\rightarrow+\infty$.

Proof.

First we notice that the choice of $t_{k+1}=\frac{1+\sqrt{1+4t_{k}^{2}}}{2}$ fulfills that

[TABLE]

Now we derive the stated convergence result by first showing via induction that

[TABLE]

Assuming that this holds for $k\geq 1$, we have that

[TABLE]

and

[TABLE]

Furthermore, for every $N\geq 1$ we have that

[TABLE]

The statement of the convergence rate in expectation follows now by plugging in our parameter choices into the statement of Theorem 4.1, using the estimate (117) and checking that

[TABLE]

The almost sure convergence of ${(F(x_{N}))}_{N\geq 1}$ can be deduced by looking at (115) and dividing by $\gamma_{k+1}t_{k+1}^{2}$, which gives for every $k\geq 0$

[TABLE]

Using the fact that $\gamma_{k+1}t_{k+1}^{2}\geq\gamma_{k}t_{k}^{2}$ we deduce that for every $k\geq 0$

[TABLE]

Plugging in our choice of parameters gives for every $k\geq 0$

[TABLE]

where $C>0$.

Thus, by the famous Robbins-Siegmund Theorem (see [22, Theorem 1]) we get that ${(F^{k+1}(x_{k+1})-F^{k+1}(x^{*})+\mu_{k+1}\frac{L_{g}^{2}}{2})}_{k\geq 0}$ converges almost surely. In particular, from the convergence to [math] in expectation we know that the almost sure limit must also be the constant zero. ∎

The formulation of the previous section can be used to deal e.g. with problems of the form

[TABLE]

for $f:\mathcal{H}\to\mkern 1.5mu\overline{\mkern-1.5mu\mathbb{R}\mkern-1.5mu}\mkern 1.5mu$ a proper, convex and lower semicontinuous function, $g_{i}:\mathcal{G}_{i}\to\mathbb{R}$ convex and $L_{g_{i}}$-Lipschitz continuous functions and $K_{i}:\mathcal{H}\to\mathcal{G}_{i}$ linear continuous operators for $i=1,\dots,m$.

Clearly one could consider

[TABLE]

and

[TABLE]

in order to reformulate the problem as

[TABLE]

and use Algorithm 3.1 together with the parameter choices described in Corollary 3.1 on this. This results in the following algorithm.

Algorithm 4.2**.**

Let $y_{0}=x_{0}\in\mathcal{H},\mu_{1}=b\lVert\bm{K}\rVert$, for $b>0$, and $t_{1}=1$. Consider the following iterative scheme

[TABLE]

However, problem (122) also lends itself to be tackled via the stochastic version of our method, Algorithm 4.1, by randomly choosing a subset of the summands. Together with the parameter choices described in Corollary 4.1 which results in the following scheme.

Algorithm 4.3**.**

Let $y_{0}=x_{0}\in\mathcal{H},b>0$, and $t_{1}=1$. Consider the following iterative scheme

[TABLE]

with $\epsilon_{k}:=(\epsilon_{1,k},\epsilon_{2,k},\dots,\epsilon_{m,k})$ a sequence of i.i.d., ${\{0,1\}}^{m}$ random variables and $p_{i}=\mathbb{P}[\epsilon_{i,1}=1]$.

Remark 4.1*.*

In theory Algorithm 4.1 could be used to treat more general stochastic problems than finite sums like (122), but in this case it is not clear anymore how a gradient estimator can be found, so we do not discuss it here.

5 Numerical Examples

We will focus our numerical experiments on image processing problems. The examples are implemented in python using the operator discretization library (ODL) [1]. We define the discrete gradient operators $D_{1}$ and $D_{2}$ representing the discretized derivative in the first and second coordinate respectively, which we will need for the numerical examples. Both map from $\mathbb{R}^{m\times n}$ to $\mathbb{R}^{m\times n}$ and are defined by

[TABLE]

and

[TABLE]

The operator norm of $D_{1}$ and $D_{2}$, respectively, is $2$ (where we equipped $\mathbb{R}^{m\times n}$ with the Frobenius norm). This yields an operator norm of $\sqrt{8}$ for the total gradient $D:=D_{1}\times D_{2}$ as a map from $\mathbb{R}^{m\times n}$ to $\mathbb{R}^{m\times n}\times\mathbb{R}^{m\times n}$, see also [12].

We will compare our methods, i.e. the Variable Accelerated SmooThing (VAST) and its stochastic counterpart (sVAST) to the Primal Dual Hybrid Gradient (PDHG) of [15] as well as its stochastic version (sPDHG) from [14]. Furthermore, we will illustrate another competitor, the method by Pesquet and Repetti, see [21], which is another a stochastic version of PDHG (see also [25]).

In all examples we choose the parameters in accordance with [14]:

• •

for PDHG and Pesquet&Repetti: $\tau=\sigma_{i}=\frac{\gamma}{\lVert K\rVert}$

• •

for sPDHG: $\sigma_{i}=\frac{\gamma}{\lVert K\rVert}$ and $\tau=\frac{\gamma}{n\max_{i}\lVert K_{i}\rVert}$,

where $\gamma=0.99$.

5.1 Total Variation Denoising

The task at hand is to reconstruct an image from its noisy observation. We do this by solving

[TABLE]

with $\alpha>0$ as regularization parameter, in the following setting: $f=\alpha\lVert\cdot-b\rVert_{2},g_{1}=g_{2}=\lVert\cdot\rVert_{1},K_{1}=D_{1},K_{2}=D_{2}$.

Figure 1 illustrates the images (of dimension $m=442$ and $n=331$) used in for this example. These include the groundtruth, i.e. the uncorrupted image, as well as the data for the optimization problem $b$, which visualizes the level of noise. In Figure 2 we can see that for the deterministic setting our method is as good as PDHG. For the objective function values, Subfigure 2(b), this is not too surprising as both algorithms share the same convergence rate. For the distance to a solution however we completely lack a convergence result. Nevertheless in Subfigure 2(a) we can see that our method performs also well with respect to this measure.

In the stochastic setting we can see in Figure 2 that, while sPDHG provides some benefit over its deterministic counterpart, the stochastic version of our method, although significantly increasing the variance, provides great benefit, at least for the objective function values.

Furthermore, Figure 3, shows the reconstructions of sPDHG and our method which are, despite the different objective function values, quite comparable.

5.2 Total Variation Deblurring

For this example we want to reconstruct an image from a blurred and noisy image. We assume to know the blurring operator $C:\mathbb{R}^{m\times n}\rightarrow\mathbb{R}^{m\times n}$. This is done by solving

[TABLE]

for $\alpha>0$ as regularization parameter, in the following setting: $f=0,g_{1}=\alpha\lVert\cdot-b\rVert_{2},g_{2}=g_{3}=\lVert\cdot\rVert_{1},K_{1}=C,K_{2}=D_{1},K_{2}=D_{2}$.

Figure 4 shows the images used to set up the optimization problem (131), in particular Subfigure 4(b) which corresponds to $b$ in said problem.

In Figure 5 we see that while PDGH performs better in the deterministic setting, in particular in the later iteration, the stochastic variable smoothing method provides a significant improvement where sPDHG method seems not to converge. It is interesting to note that in this setting even the deterministic version of our algorithm exhibits a slightly chaotic behaviour. Although neither of the two methods is monotone in the primal objective function PDHG seems here much more stable.