Bregman Finito/MISO for nonconvex regularized finite sum minimization without Lipschitz gradient continuity

Puya Latafat, Andreas Themelis, Masoud Ahookhosh, Panagiotis Patrinos

Introduction

We consider the following regularized finite sum minimization

where $\overline{C}$ denotes the closure of $C\coloneqq\bigcap_{i=1}^{N}\operatorname{int}\operatorname{dom}h_{i}$ , for some convex functions $h_{i}$ , $i\in[N]\coloneqq{\mathopen{}\left\{1,\dots,N{}\mathrel{\mid}{}{}\right\}\mathclose{}}$ . Our goal in this paper is to study such problems without imposing convexity assumptions on $f_{i}$ and $g$ , and in a setting where $f_{i}$ are differentiable but their gradients need not be Lipschitz continuous. Our full setting is formalized in I.

To relax the Lipschitz differentiability assumption, we adopt the notion of smoothness relative to a distance-generating function , and following we will use the terminology of relative smoothness. Despite the lack of Lipschitz differentiability, in many applications the involved functions satisfy a descent property where the usual quadratic upper bound is replaced by a Bregman distance (cf. Items 1and 2.1). Owing to this property, Bregman extensions for many classical schemes have been proposed .

In the setting of finite sum minimization, the incremental aggregated algorithm PLIAG was proposed recently as a Bregman variant of the incremental aggregated gradient method . The analysis of PLIAG is limited to the convex case and requires restrictive assumptions for the Bregman kernel [69, Thm. 1, Assump. 8]. Stochastic mirror descent (SMD) is another relevant algorithm which can tackle more general stochastic optimization problems. SMD may be viewed as a Bregman extension of the stochastic (sub)gradient method and has long been studied . More recently, studied SMD for convex and relatively smooth formulations, and (sub)gradient versions have been analyzed under relative continuity in a convex setting , as well as relative weak convexity .

Motivated by these recent works, we propose a Bregman extension of the popular Finito/MISO algorithm in a fully nonconvex setting and with very general sampling strategies that will be made percise shortly after. In a nutshell, our analysis revolves around the fact that, regardless of the index selection strategy, the function $\mathcal{L}:\varmathbb R^{n}\times\varmathbb R^{nN}\rightarrow\overline{\varmathbb}R$ defined as

where $\hat{h}_{i}^{\ast}$ denotes the convex conjugate of $\hat{h}_{i}\coloneqq\nicefrac{{h_{i}}}{{\gamma_{i}}}-\nicefrac{{f_{i}}}{{N}}$ , monotonically decreases along the iterates $(z^{k},\bm{s}^{k})_{k\in\varmathbb N}$ generated by Algorithm 1 (see I for the requirements on $h_{i},f_{i}$ ). Our methodology leverages an interpretation of Finito/MISO as a block-coordinate algorithm that was observed in in the Euclidean setting. In fact, the analysis is here further simplified after noticing that the smooth function can be “hidden” in the distance-generating function, resulting in a Lyapunov function $\mathcal{L}$ that can be expressed as a Bregman Moreau envelope (cf. 3.2).

We cover a wide range of sampling strategies for the index set $\mathcal{I}^{k+1}$ at step 5, which we can summarize into the following two umbrella categories:

The randomized setting ( $\mathcal{S}_{1}$ ), in which $\mathcal{P}_{k}$ denotes the probability conditional to the knowledge at iteration $k$ , covers, for instance, a mini-batch strategy of size $b$ . Another notable case is when each index $i$ is selected at random with probability $p_{i}$ independently of other indices.

The essentially cyclic rule ( $\mathcal{S}_{2}$ ) is also very general and has been considered by many authors . Two notable special cases of single index selection rules complying with ( $\mathcal{S}_{2}$ ) are the cyclic and shuffled cyclic sampling strategies:

We remark that, in the cyclic case, our algorithm generalizes DIAG for smooth strongly convex problems, which itself may be seen as a cyclic variant of Finito/MISO.

One iteration of Algorithm 1 involves the computation of $z^{k}$ at step 4 and that of the gradients ${\nabla}(\nicefrac{{h_{i}}}{{\gamma_{i}}}-\nicefrac{{f_{i}}}{{N}})$ , $i\in\mathcal{I}^{k+1}$ , at step 6. Consequently, the overall complexity of each iteration is independent of the number $N$ of functions appearing in problem (P), and is instead proportional to the number of sampled indices, which the user is allowed to upper bound to any integer between 1 and $N$ . As is the case for all incremental gradient methods, the low iteration cost comes at the price of having to store in memory a table $\bm{s}^{k}$ of $N$ many $\varmathbb R^{n}$ vectors, which can become problematic when $N$ grows large. Other incremental algorithms for convex optimization such as IAG , IUG , SAG , and SAGA , can considerably reduce memory allocation from $O(nN)$ to $O(n)$ in applications such as logistic regression and lasso where the gradients ${\nabla}f_{i}$ can be expressed as scaled versions of the data vectors. Despite the favorable performance of the Finito/MISO algorithm on such problems as observed in , this memory reduction trick can not be employed due to the fact that the vectors $s_{i}$ stored in the table depend not only on the gradients, but also on the vectors ${\nabla}h_{i}(z^{k})$ . Nevertheless, inspired by the popular stochastic methods SVRG and SARAH , by suitably interleaving incremental and full gradient evaluations it is possible to completely waive the need of a memory table and match the $O(n)$ storage requirement.

Since full gradient updates correspond to selecting all indices, Algorithm 2 may be viewed as Algorithm 1 with an essentially cyclic sampling rule of period $N$ , a claim that will be formally shown in 4.12. In fact, not only does it naturally inherit all the convergence results, but its particular sampling strategy also allows us to waive convexity requirements on $g$ that are necessary for more general essentially cyclic rules.

2 Contributions

As a means to informally summarize the content of the paper, in Table 1 we synopsize the convergence results of the two algorithms.

To the best of our knowledge, this is the first analysis of an incremental aggregated method in a fully nonconvex setting and without Lipschitz differentiability assumptions. Our analysis, surprisingly simple and yet covering randomized and essentially cyclic samplings altogether, relies on a sure descent property on the Bregman Moreau envelope (cf. 4.2).

We propose a novel low-memory variant of the (Bregman) Finito/MISO algorithm, that in the spirit of SVRG and SARAH alternates between incremental steps and a full proximal gradient step. It is highly interesting even in the Euclidean case, as it can accommodate fully nonconvex formulations while maintaining a $O(n)$ memory requirement.

Linear convergence of Algorithm 1 in the randomized case is established when the cost function $\varphi$ is strongly convex, yet with no convexity requirement on $f_{i}$ or $g$ . To the best of our knownledge, this is a novelty even in the Euclidean case, for all available results are bound to strong convexity of each term $f_{i}$ in the sum; see e.g., . This type of assumption has also been considered in for the case of SVRG. Although in practice it is hardly ever the case that $f$ is strongly convex without each $f_{i}$ also being (strongly) convex, our analysis being agnostic to the decomposition leads to tighter and more general results.

We leverage the Kurdyka-Łojasiewicz (KL) property to establish global (as apposed to subsequential) convergence as well as linear convergence, for Algorithm 1 with (essentially) cyclic sampling and for the low-memory Algorithm 2.

3 Organization

We conclude this section by introducing some notational conventions. The problem setting is formally described in Section 2 together with a list of related definitions and known facts involving Bregman distances, relative smoothness, and proximal mapping. Section 3 offers an alternative interpretation of Algorithm 1 as the block-coordinate Bregman proximal point Algorithm 3, which majorly simplifies the analysis, addressed in Section 4. Some auxiliary results are deferred to Appendix A. Section 5 applies the proposed algorithms to sparse phase retrieval problems, and Section 6 concludes the paper.

4 Notation

The set of real and extended-real numbers are $\varmathbb R\coloneqq(-\infty,\infty)$ and $\overline{\varmathbb}R\coloneqq\varmathbb R\cup{\mathopen{}\left\{\infty{}\mathrel{\mid}{}{}\right\}\mathclose{}}$ , while the positive and strictly positive reals are $\varmathbb R_{+}\coloneqq[0,\infty)$ and $\varmathbb R_{++}\coloneqq(0,\infty)$ . With ${\rm id}$ we indicate the identity function $x\mapsto x$ defined on a suitable space. We denote by ${\mathopen{}\left\langle{}{}\cdot{}{},{}{}\cdot{}{}\right\rangle\mathclose{}}$ and $\|{}\cdot{}\|$ the standard Euclidean inner product and the induced norm. For a vector $\bm{w}=(w_{1},\ldots,w_{r})\in\varmathbb R^{\sum_{i}n_{i}}$ , $w_{i}\in\varmathbb R^{n_{i}}$ is used to denote its $i$ -th block coordinate. $\operatorname{int}E$ and $\operatorname{bdry}E$ respectively denote the interior and boundary of a set $E$ , and for a sequence $(x^{k})_{k\in\varmathbb N}$ we write $(x^{k})_{k\in\varmathbb N}\subseteq E$ to indicate that $x^{k}\in E$ for all $k\in\varmathbb N$ . We say that $(x^{k})_{k\in\varmathbb N}$ converges at $Q$ -linear rate (resp. $R$ -linear rate) to a point $x$ if there exists $c\in(0,1)$ such that $\|x^{k+1}-x\|\leq c\|x^{k}-x\|$ (resp. $\|x^{k}-x\|\leq\rho c^{k}$ for some $\rho>0$ ) holds for all $k\in\varmathbb N$ .

We use the notation $H:\varmathbb R^{n}\rightrightarrows\varmathbb R^{m}$ to indicate a mapping from each point $x\in\varmathbb R^{n}$ to a subset $H(x)$ of $\varmathbb R^{m}$ . The graph of $H$ is the set $\operatorname{gph}H{}\coloneqq{}{\mathopen{}\left\{(x,y)\in\varmathbb R^{n}\times\varmathbb R^{m}{}\mathrel{\mid}{}y\in H(x)\right\}\mathclose{}}$ . We say that $H$ is outer semicontinuous (osc) if $\operatorname{gph}H$ is a closed subset of $\varmathbb R^{n}\times\varmathbb R^{m}$ , and locally bounded if for every bounded $U\subset\varmathbb R^{n}$ the set $\bigcup_{x\in U}H(x)$ is bounded.

The domain and epigraph of an extended-real-valued function $h:\varmathbb R^{n}\rightarrow\overline{\varmathbb}R$ are the sets $\operatorname{dom}h{}\coloneqq{}{\mathopen{}\left\{x\in\varmathbb R^{n}{}\mathrel{\mid}{}h(x)<\infty\right\}\mathclose{}}$ and $\operatorname{epi}h{}\coloneqq{}{\mathopen{}\left\{(x,\alpha)\in\varmathbb R^{n}\times\varmathbb R{}\mathrel{\mid}{}h(x)\leq\alpha\right\}\mathclose{}}$ . Function $h$ is said to be proper if $\operatorname{dom}h\neq\emptyset$ , and lower semicontinuous (lsc) if $\operatorname{epi}h$ is a closed subset of $\varmathbb R^{n+1}$ . We say that $h$ is level bounded if its $\alpha$ -sublevel set $\operatorname{lev}_{\leq\alpha}h{}\coloneqq{}{\mathopen{}\left\{x\in\varmathbb R^{n}{}\mathrel{\mid}{}h(x)\leq\alpha\right\}\mathclose{}}$ is bounded for all $\alpha\in\varmathbb R$ . The conjugate of $h$ , is defined by $h^{\ast}(y){}\coloneqq{}\sup_{x\in\varmathbb R^{n}}{\mathopen{}\left\{{\mathopen{}\left\langle{}y{},{}x{}\right\rangle\mathclose{}}-h(x){}\mathrel{\mid}{}{}\right\}\mathclose{}}$ . The indicator function of a set $E\subseteq\varmathbb R^{n}$ is denoted by $\operatorname{\delta}_{E}$ , namely $\operatorname{\delta}_{E}(x)=0$ if $x\in E$ and $\infty$ otherwise.

We denote by $\hat{\partial}h:\varmathbb R^{n}\rightrightarrows\varmathbb R^{n}$ the regular subdifferential of $h$ , where

A necessary condition for local minimality of $x$ for $h$ is $0\in\hat{\partial}h(x)$ , see [53, Th. 10.1]. The (limiting) subdifferential of $h$ is $\partial h:\varmathbb R^{n}\rightrightarrows\varmathbb R^{n}$ , where $v\in\partial h(x)$ iff $x\in\operatorname{dom}h$ and there exists a sequence $(x^{k},v^{k})_{k\in\varmathbb N}\subseteq\operatorname{gph}\hat{\partial}h$ such that $(x^{k},h(x^{k}),v^{k}){}\to{}(x,h(x),v)$ as $k\to\infty$ . Finally, the set of $r$ times continuously differentiable functions from $X$ to $\varmathbb R$ is denoted by $\mathcal{C}^{r}(X)$ .

Problem setting and preliminaries

Throughout this paper, problem (P) is studied under the following assumptions.

$f_{i}:\varmathbb R^{n}\rightarrow\overline{\varmathbb}R$ are $L_{f_{i}}$ -smooth relative to Legendre kernels $h_{i}$ (s 2.2and 2.4);

$g:\varmathbb R^{n}\rightarrow\overline{\varmathbb}R$ is proper and lower semicontinuous (lsc);

a solution exists: $\operatorname*{arg\,min}{\mathopen{}\left\{\varphi(x){}\mathrel{\mid}{}x\in\overline{C}\right\}\mathclose{}}\neq\emptyset$ ;

for given $\gamma_{i}\in(0,\nicefrac{{N}}{{L_{f_{i}}}})$ , $i\in[N]$ , it holds that for any $s\in\varmathbb R^{n}$

As it will become clear in Section 3, the subproblem (2.1) is in fact a reformulation of a (Bregman) proximal mapping. Assumption 4 excludes boundary points from $\operatorname{range}T$ . This is a standard assumption that usually holds in practice , e.g., when $g$ is convex or when the intersection of $\operatorname{dom}h_{i}$ , $i\in[N]$ , is an open set.

For a convex function $h:\varmathbb R^{n}\rightarrow\overline{\varmathbb}R$ that is continuously differentiable on $\operatorname{int}\operatorname{dom}h\neq\emptyset$ , the Bregman distance $\operatorname{D}_{h}:\varmathbb R^{n}\times\varmathbb R^{n}\rightarrow\overline{\varmathbb}R$ is defined as

Function $h$ will be referred to as a distance-generating function.

A proper, lsc, and strictly convex function $h:\varmathbb R^{n}\rightarrow\overline{\varmathbb}R$ with $\operatorname{int}\operatorname{dom}h\neq\emptyset$ and such that $h\in\mathcal{C}^{1}(\operatorname{int}\operatorname{dom}h)$ is said to be a Legendre kernel if it is (i) $1$ -coercive, i.e., such that $\lim_{\|x\|\to\infty}\nicefrac{{h(x)}}{{\|x\|}}=\infty$ , and (ii) essentially smooth, i.e., if $\|{\nabla}h(x_{k})\|\to\infty$ for every sequence $(x_{k})_{k\in\varmathbb N}\subseteq\operatorname{int}\operatorname{dom}h$ converging to a boundary point of $\operatorname{dom}h$ .

Let $h:\varmathbb R^{n}\rightarrow\overline{\varmathbb}R$ be a Legendre kernel, $x\in\varmathbb R^{n}$ , and $y,z\in\operatorname{int}\operatorname{dom}h$ . Then:

$h^{\ast}\in\mathcal{C}^{1}(\varmathbb R^{n})$ is strictly convex and ${\nabla}h^{-1}={\nabla}h^{\ast}$ [52, Thm. 26.5 and Cor. 13.3.1].

$\operatorname{D}_{h}(x,z)=\operatorname{D}_{h}(x,y)+\operatorname{D}_{h}(y,z)+\langle{}x-y{},{}{\nabla}h(y)-{\nabla}h(z){}\rangle$ [22, Lem. 3.1].

$\operatorname{D}_{h}(y,z)=\operatorname{D}_{h^{\ast}}({\nabla}h(z),{\nabla}h(y))$ [8, Thm. 3.7(v)].

$\operatorname{D}_{h}({}\cdot{},z)$ and $\operatorname{D}_{h}(z,{}\cdot{})$ are level bounded [9, Lem. 7.3(v)-(viii)].

If $\operatorname{dom}h$ is closed and $\operatorname{D}_{h}(x^{k},y^{k})\to 0$ for some $x^{k}\in\operatorname{dom}h$ and $y^{k}\in\operatorname{int}\operatorname{dom}h$ , then $(x^{k})_{k\in\varmathbb N}$ converges to a point $x$ iff so does $(y^{k})_{k\in\varmathbb N}$ [56, Thm. 2.4].

Moreover, for any convex set $\mathcal{U}\subseteq\operatorname{int}\operatorname{dom}h$ and $u,v\in\mathcal{U}$ the following hold:

If $h$ is $\mu_{h,\mathcal{U}}$ -strongly convex on $\mathcal{U}$ , then $\frac{\mu_{h,\mathcal{U}}}{2}\|v-u\|^{2}{}\leq{}\operatorname{D}_{h}(v,u){}\leq{}\frac{1}{2\mu_{h,\mathcal{U}}}\|{\nabla}h(v)-{\nabla}h(u)\|^{2}$ .

We say that a proper, lsc function $f:\varmathbb R^{n}\rightarrow\overline{\varmathbb}R$ is smooth relative to a Legendre kernel $h:\varmathbb R^{n}\rightarrow\overline{\varmathbb}R$ if $\operatorname{dom}f\supseteq\operatorname{dom}h$ , and there exists $L_{f}\geq 0$ such that $L_{f}h\pm f$ are convex functions on $\operatorname{int}\operatorname{dom}h$ . We will simply say that $f$ is $L_{f}$ -smooth relative to $h$ to make the modulus $L_{f}$ explicit.

Let $f:\varmathbb R^{n}\rightarrow\overline{\varmathbb}R$ be $L_{f}$ -smooth relative to a Legendre kernel $h:\varmathbb R^{n}\rightarrow\overline{\varmathbb}R$ . Then, $f\in\mathcal{C}^{1}(\operatorname{int}\operatorname{dom}h)$ and the following hold:

$\bigl{|}f(y)-f(x)-\langle{}{\nabla}f(x){},{}y-x{}\rangle\bigr{|}{}\leq{}L_{f}\operatorname{D}_{h}(y,x)$ for all $x,y\in\operatorname{int}\operatorname{dom}h$ .

$-L_{f}\nabla^{2}h{}\preceq{}\nabla^{2}f{}\preceq{}L_{f}\nabla^{2}h$ on $\operatorname{int}\operatorname{dom}h$ , provided that $f,h\in\mathcal{C}^{2}(\operatorname{int}\operatorname{dom}h)$ .

Relative to a Legendre kernel $h:\varmathbb R^{n}\rightarrow\overline{\varmathbb}R$ , the Bregman proximal mapping of $\psi$ is the set-valued map $\operatorname{prox}_{\psi}^{h}:\operatorname{int}\operatorname{dom}h\rightrightarrows\varmathbb R^{n}$ given by

The following hold for a Legendre kernel $h:\varmathbb R^{n}\rightarrow\overline{\varmathbb}R$ and a proper, lsc, lower bounded function $\psi:\varmathbb R^{n}\rightarrow\overline{\varmathbb}R$ :

$\operatorname{prox}_{\psi}^{h}$ is locally bounded, compact-valued, and outer semicontinuous on $\operatorname{int}\operatorname{dom}h$ .

$\psi^{h}$ is real-valued and continuous on $\operatorname{int}\operatorname{dom}h$ ; in fact, it is locally Lipschitz if so is ${\nabla}h$ .

Let $h$ be a Legendre kernel and $\psi:\varmathbb R^{n}\rightarrow\overline{\varmathbb}R$ be proper, lsc, and lower bounded on $\overline{\operatorname{dom}h}$ . Then, for every $x\in\operatorname{int}\operatorname{dom}h$ , $y\in\operatorname{dom}h$ , and $\bar{x}\in\operatorname{prox}_{\psi}^{h}(x)$

$\psi^{h}(x){}\mathrel{{}\mathop{=}\limits^{\text{\clap{\tiny(def)}}}}{}\psi(\bar{x}){}+{}\operatorname{D}_{h}(\bar{x},x){}\leq{}\psi(y){}+{}\operatorname{D}_{h}(y,x)$ , and in particular $\psi^{h}(x)\leq\psi(x)$ ;

if $\psi$ is convex, then $\psi^{h}(x){}\leq{}\psi(y){}+{}\operatorname{D}_{h}(y,x){}-{}\operatorname{D}_{h}(y,\bar{x})$ [59, Lem. 3.1].

Moreover, if $\operatorname{range}\operatorname{prox}_{\psi}^{h}\subseteq\operatorname{int}\operatorname{dom}h$ , then the following also hold [1, Prop.3.3]:

$\operatorname*{inf}_{\overline{\operatorname{dom}h}}\psi{}\leq{}\operatorname*{inf}_{\operatorname{int}\operatorname{dom}h}\psi{}={}\operatorname*{inf}\psi^{h}$ and $\operatorname*{arg\,min}\psi^{h}{}={}\operatorname*{arg\,min}_{\operatorname{int}\operatorname{dom}h}\psi$ .

$\psi+\operatorname{\delta}_{\overline{\operatorname{dom}h}}$ is level bounded iff so is $\psi^{h}$ .

A block-coordinate interpretation

By introducing $N$ copies of $x$ , problem (P) can equivalently be written as

where $\Delta{}\coloneqq{}{\mathopen{}\left\{\bm{x}=(x_{1},\ldots,x_{N})\in\varmathbb R^{nN}{}\mathrel{\mid}{}x_{1}=x_{2}=\cdots=x_{N}\right\}\mathclose{}}$ is the consensus set. The equivalence between (3.1) and the original problem (P) is formally established in A.1. Note that Assumption 1 implies that $F$ as in (3.1) is smooth with respect to the Legendre kernel

making Bregman forward-backward iterations $\bm{x}^{+}{}\in{}\operatorname*{arg\,min}\{\langle{}{\nabla}F(\bm{x}){},{}{}\cdot{}{}\rangle{}+{}G({}\cdot{}){}+{}\tfrac{1}{\gamma}\operatorname{D}_{H}({}\cdot{},\bm{x}){}\mathrel{\mid}{}{}\}$ for some stepsize $\gamma>0$ a suitable option to address problem (3.1). In fact, it can be easily verified that $L_{F}=\tfrac{1}{N}\max_{i=1\ldots N}L_{f_{i}}$ is a smoothness modulus of $F$ relative to $H$ , indicating that fixed point iterations $\bm{x}\leftarrow\bm{x}^{+}$ under I converge (in some sense to be made precise) to a stationary point of the problem whenever $\gamma\in(0,\nicefrac{{1}}{{L_{F}}})$ . Notice that a higher degree of flexibility can be granted by considering an $N$ -uple of individual stepsizes $\Gamma=(\gamma_{1},\ldots,\gamma_{N})$ , giving rise to the forward-backward operator $\operatorname{T}_{\Gamma}^{F,G}:\varmathbb R^{nN}\rightrightarrows\varmathbb R^{nN}$ in the Bregman metric $(\bm{z},\bm{x})\mapsto\sum_{i=1}^{N}\tfrac{1}{\gamma_{i}}\operatorname{D}_{h_{i}}(z_{i},x_{i})$ , namely

This intuition is validated in the next result, which asserts that whenever the stepsizes $\gamma_{i}$ are selected as in Algorithm 1 the operator $\operatorname{T}_{\Gamma}^{F,G}$ coincides with a proximal mapping on a suitable Legendre kernel function $\hat{H}$ . This observation leads to a much simpler analysis of Algorithm 1, which will be shown to be a block-coordinate variant of a Bregman proximal point method.

Suppose that Assumption 1 holds and let $\gamma_{i}\in(0,\nicefrac{{N}}{{L_{f_{i}}}})$ be selected as in Algorithm 1. Then, $\hat{h}_{i}\coloneqq\frac{1}{\gamma_{i}}h_{i}-\tfrac{1}{N}f_{i}$ (with the convention $\infty-\infty=\infty$ ) is a Legendre kernel with $\operatorname{dom}\hat{h}_{i}=\operatorname{dom}h_{i}$ , $i\in[N]$ , and thus so is the function

Moreover, for any $(\bm{z},\bm{x})\in\varmathbb R^{nN}\times\varmathbb R^{nN}$ it holds that

and in particular the forward-backward operator (3.3) satisfies

When I is satisfied, then the following also hold:

$\operatorname{D}_{\hat{H}}(\bm{z},\bm{x}){}\geq{}\sum_{i=1}^{N}(\tfrac{1}{\gamma_{i}}-\tfrac{L_{f_{i}}}{N})\operatorname{D}_{h_{i}}(z_{i},x_{i})$ .

$\operatorname{prox}_{\Phi}^{\hat{H}}(\bm{x}){}={}{\mathopen{}\left\{(z,\cdots,z){}\mathrel{\mid}{}z\in T(\sum_{i=1}^{N}{\nabla}\hat{h}_{i}(x_{i}))\right\}\mathclose{}}$ , with $T$ as in (2.1), is a nonempty and compact subset of $C\times\cdots\times C$ for any $\bm{x}\in\operatorname{int}\operatorname{dom}h_{1}\times\cdots\times\operatorname{int}\operatorname{dom}h_{N}$ .

If $\bm{z}\in\operatorname{prox}_{\Phi}^{\hat{H}}(\bm{x})$ , then ${\nabla}\hat{H}(\bm{x}){}-{}{\nabla}\hat{H}(\bm{z}){}\in{}\hat{\partial}\Phi(\bm{z})$ ; the converse also holds when $\Phi$ is convex.

If $h_{i}$ is $\mu_{h_{i},\mathcal{U}_{i}}$ -strongly convex on a convex set $\mathcal{U}_{i}\subseteq\operatorname{dom}h_{i}$ , then $\hat{h}_{i}$ is $\mu_{\hat{h}_{i},\mathcal{U}_{i}}$ -strongly convex on $\mathcal{U}_{i}$ with $\mu_{\hat{h}_{i},\mathcal{U}_{i}}{}\geq{}\big{(}\tfrac{1}{\gamma_{i}}-\tfrac{L_{f_{i}}}{N}\big{)}\mu_{h_{i},\mathcal{U}_{i}}$ .

The claims on $\hat{h}_{i}$ are shown in [1, Thm. 4.1], and (3.5) and (3.6) then easily follow.

1 This is an immediate consenquence of Item 1.

2 Let $\bm{x}$ be as in the statement, and observe that $\bm{x}\in\operatorname{int}\operatorname{dom}\hat{H}$ ; nonemptyness and compactness of $\operatorname{prox}_{\Phi}^{\hat{H}}$ then follows from Item 1. Let now $\bm{u}\in\operatorname{prox}_{\Phi}^{\hat{H}}(\bm{x})$ be fixed, and note that the consensus constraint encoded in $\Phi$ ensures that $u_{i}=u_{j}$ for all $i,j\in[N]$ . Thus,

where the inclusion follows from Assumption 4.

3 Observe first that necessarily $\bm{x}\in\operatorname{int}\operatorname{dom}h_{i}\times\cdots\times\operatorname{int}\operatorname{dom}h_{N}$ , for otherwise no such $\bm{z}$ exists. Moreover, from assertion 2 it follows that also $\bm{z}$ belongs to such open set, onto which $\hat{H}$ is continuously differentiable. The claim then follows from the necessary condition for optimality of $\bm{z}$ in the minimization problem (2.4) — which is also sufficient when $\Phi$ is convex, for so is $\Phi+\operatorname{D}_{\hat{H}}({}\cdot{},\bm{x})$ in this case — having

The last equality follows from [53, Ex. 8.8(c)], owing to smoothness of $\hat{H}$ at $\bm{z}$ .

Algorithm 3 presents a block coordinate (BC) proximal point algorithm with the distance generating function $\hat{H}$ . Note that in a departure from most of the existing literature on BC proximal methods that consider separable nonsmooth terms (see e.g., ), here the nonsmooth function $G$ in (3.1) is nonseparable. It is shown in the next lemma that this conceptual algorithm is equivalent to the Bregman Finito/MISO Algorithm 1.

$s_{i}^{k}=\frac{1}{\gamma_{i}}{\nabla}h_{i}(x_{i}^{k})-\tfrac{1}{N}{\nabla}f_{i}(x_{i}^{k})$ (or, equivalently, $x_{i}^{k}={\nabla}\hat{h}_{i}^{\ast}(s_{i}^{k})$ )

$\varphi(z^{k})=\Phi(\bm{u}^{k})=\Phi^{\hat{H}}(\bm{x}^{k})-\operatorname{D}_{\hat{H}}(\bm{u}^{k},\bm{x}^{k})$

$\Phi^{\hat{H}}(\bm{x}^{k})=\mathcal{L}(z^{k},\bm{s}^{k})$ , where $\mathcal{L}$ is as in (1.1).

Let the index sets $\mathcal{I}^{k+1}$ and $\mathcal{J}^{k+1}$ be chosen identically, $k\in\varmathbb N$ . It follows from Item 2 that $u_{i}^{k}=u_{j}^{k}$ for all $k\in\varmathbb N$ and $i,j\in[N]$ , with

and the last term is $\mathcal{L}(z^{k},\bm{s}^{k})$ (cf. Items 3and 1), yielding assertion 6. ∎

Convergence analysis

The block coordinate interpretation of Algorithm 1 presented in Section 3 plays a crucial role in the proposed methodology, and leads to a remarkably simple convergence analysis. In fact, many key facts can be established without confining the discussion to a particular sampling strategy. These preliminary results are presented in the next subsection and will be extensively referred to in the subsequent subsections that are instead devoted to a specific sampling strategy.

Unlike classical analyses of BC proximal methods that employ the cost as a Lyapunov function (see e.g., [10, §11]), here, the nonseparability of $G$ precludes this possibility. To address this challenge, we instead employ the Bregman Moreau envelope equipped with the distance generating function $\hat{H}$ (see (3.4)). Before showing its Lyapunov-type behavior for Algorithm 3, we list some of its properties and its relation with the original problem. The proof is a simple consequence of Items 2and 2.7 and the fact that $\hat{H}$ is a Legendre kernel with $\operatorname{dom}\hat{H}=\operatorname{dom}h_{1}\times\dots\times\operatorname{dom}h_{N}$ (cf. 3.1).

$\Phi^{\hat{H}}$ is continuous on $\operatorname{dom}\Phi^{\hat{H}}=\operatorname{int}\operatorname{dom}h_{1}\times\dots\times\operatorname{int}\operatorname{dom}h_{N}$ , in fact, locally Lipschitz if so is ${\nabla}h_{i}$ on $\operatorname{int}\operatorname{dom}h_{i}$ , $i\in[N]$ ;

$\min_{\overline{C}}\varphi{}\leq{}\operatorname*{inf}_{C}\varphi{}={}\operatorname*{inf}\Phi^{\hat{H}}$ and $\operatorname*{arg\,min}\Phi^{\hat{H}}{}={}{\mathopen{}\left\{(x^{\star},\dots,x^{\star}){}\mathrel{\mid}{}x^{\star}\in\operatorname*{arg\,min}_{C}\varphi\right\}\mathclose{}}$ ;

$\Phi^{\hat{H}}$ is level bounded iff so is $\varphi+\operatorname{\delta}_{\overline{C}}$ .

Suppose that I holds, and consider the iterates generated by Algorithm 3. Then, $\bm{u}^{k}=(u^{k},\ldots,u^{k})$ for some $u^{k}\in C$ and $\bm{x}^{k}\in C\times\dots\times C\subseteq\operatorname{int}\operatorname{dom}\hat{H}$ for every $k\in\varmathbb N$ , and the algorithm is thus well defined. Moreover:

$\Phi^{\hat{H}}(\bm{x}^{k+1}){}\leq{}\Phi^{\hat{H}}(\bm{x}^{k}){}-{}\operatorname{D}_{\hat{H}}(\bm{x}^{k+1},\bm{x}^{k}){}={}\Phi^{\hat{H}}(\bm{x}^{k}){}{}-{}\sum_{i\in\mathcal{J}^{k+1}}\operatorname{D}_{\hat{h}_{i}}(u^{k},x_{i}^{k})$ for every $k\in\varmathbb N$ ; when $\Phi$ is convex (i.e., when so is $\varphi$ ), then the inequality can be strengthened to $\Phi^{\hat{H}}(\bm{x}^{k+1}){}\leq{}\Phi^{\hat{H}}(\bm{x}^{k}){}-{}\operatorname{D}_{\hat{H}}(\bm{x}^{k+1},\bm{x}^{k}){}-{}\operatorname{D}_{\hat{H}}(\bm{u}^{k},\bm{u}^{k+1})$ .

$(\Phi^{\hat{H}}(\bm{x}^{k}))_{k\in\varmathbb N}$ monotonically decreases to a finite value $\varphi_{\star}\geq\operatorname*{inf}_{C}\varphi\geq\min_{\overline{C}}\varphi$ .

The sequence $(\operatorname{D}_{\hat{H}}(\bm{x}^{k+1},\bm{x}^{k}))_{k\in\varmathbb N}$ has finite sum (and in particular vanishes); the same holds also for $(\operatorname{D}_{\hat{H}}(\bm{u}^{k},\bm{u}^{k+1}))_{k\in\varmathbb N}$ when $\Phi$ is convex (i.e., when so is $\varphi$ ).

If $\varphi+\operatorname{\delta}_{\overline{C}}$ is level bounded, then $(\bm{x}^{k})_{k\in\varmathbb N}$ and $(\bm{u}^{k})_{k\in\varmathbb N}$ are bounded.

If $\operatorname{dom}h_{i}$ is closed, a subsequence $(x_{i}^{k})_{k\in K}$ converges to a point $x^{\star}$ iff so does $(x_{i}^{k+1})_{k\in K}$ .

If $C=\varmathbb R^{n}$ , then $\Phi^{\hat{H}}$ is constantly equal to $\varphi_{\star}$ as above on the limit set of $(\bm{x}^{k})_{k\in\varmathbb N}$ .

It follows from Item 2 that $u^{k}\in C$ holds for every $k\in\varmathbb N$ . Notice that for every $i\in[N]$ and $k\in\varmathbb N$ , either $x_{i}^{k}=x^{\rm init}\in C$ (by initialization), or there exists $k_{i}\leq k$ such that $x_{i}^{k}=z^{k_{i}}\in C$ . It readily follows that $\bm{x}^{k}\in C\times\dots\times C\subseteq\operatorname{int}\operatorname{dom}H=\operatorname{int}\operatorname{dom}\hat{H}$ , hence that $\operatorname{prox}_{\Phi}^{\hat{H}}(\bm{x}^{k})\neq\emptyset$ for all $k\in\varmathbb N$ by Item 2, whence the well definedness of the algorithm. We now show the numbered claims.

1 It follows from Items 1and 2 that $\Phi^{\hat{H}}(\bm{x}^{k+1}){}\leq{}\Phi(\bm{u}^{k})+\operatorname{D}_{\hat{H}}(\bm{u}^{k},\bm{x}^{k+1}){}-{}c_{k}$ , where $c_{k}\geq 0$ can be taken as $c_{k}=\operatorname{D}_{\hat{H}}(\bm{u}^{k},\bm{u}^{k+1})$ when $\Phi$ is convex. Therefore,

The claim follows by noting that the inner product is zero:

2 Monotonic decrease of $(\Phi^{\hat{H}}(\bm{x}^{k}))_{k\in\varmathbb N}$ follows from assertion 1. This ensures that the sequence converges to some value $\varphi_{\star}$ , bounded below by $\min_{\overline{C}}\varphi$ in light of Item 2.

owing to Item 2 and Assumption 3. When $\varphi$ is convex, the tighter bound in assertion 1 yields the similar claim for $(\operatorname{D}_{\hat{H}}(\bm{u}^{k},\bm{u}^{k+1}))_{k\in\varmathbb N}$ .

4 It follows from assertion 2 that the entire sequence $(\bm{x}^{k})_{k\in\varmathbb N}$ is contained in the sublevel set $\{\bm{w}{}\mathrel{\mid}{}\Phi^{\hat{H}}(\bm{w})\leq\Phi^{\hat{H}}(\bm{x}^{0})\}$ , which is bounded provided that $\varphi+\operatorname{\delta}_{\overline{C}}$ is level bounded as shown in Item 3. In turn, boundedness of $(\bm{u}^{k})_{k\in\varmathbb N}$ then follows from local boundedness of $\operatorname{T}_{\Gamma}^{F,G}=\operatorname{prox}_{\Phi}^{\hat{H}}$ , cf. (3.6) and Item 1.

5 Follows from Item 5, since $x_{i}^{k}\in\operatorname{int}\operatorname{dom}h_{i}=\operatorname{int}\operatorname{dom}\hat{h}_{i}$ for every $k$ (with equality owing to 3.1), and $\operatorname{D}_{\hat{h}_{i}}(x_{i}^{k+1},x_{i}^{k})\to 0$ by assertion 3.

6 Follows from assertion 2 and the continuity of $\Phi^{\hat{H}}$ , see Item 2. ∎

In conclusion of this subsection we provide an overview of the ingredients that are needed to show that the limit points of the sequence $(z^{k})_{k\in\varmathbb N}$ generated by Algorithm 1 are stationary for problem (P). As will be shown in 4.4, these amount to the vanishing of the residual $\operatorname{D}_{\hat{H}}(\bm{u}^{k},\bm{x}^{k})$ together with some assumptions on the distance-generating functions $h_{i}$ . For the iterates of Algorithm 1, this translates to $\operatorname{D}_{\hat{h}_{i}^{\ast}}(s_{i}^{k},{\nabla}\hat{h}_{i}(z^{k})){}\to{}0$ for all indices $i\in[N]$ , indicating that all vectors $s_{i}^{k+1}$ in the table should be good estimates of ${\nabla}\hat{h}_{i}(z^{k+1}){}={}\frac{1}{\gamma_{i}}{\nabla}h_{i}(z^{k+1}){}-{}\frac{1}{N}{\nabla}f_{i}(z^{k+1})$ , as opposed to $\frac{1}{\gamma_{i}}{\nabla}h_{i}(z^{k}){}-{}\frac{1}{N}{\nabla}f_{i}(z^{k})$ and for the indices in $\mathcal{I}^{k+1}$ only (cf. step 6). As a result, we may view this property as jointly having $z^{k}-z^{k+1}$ vanish, desirable if any convergence of $(z^{k})_{k\in\varmathbb N}$ is expected, and the fact that a consensus is eventually reached among the sampled blocks.

In line with any result in the literature we are aware of, a complete convergence analysis for nonconvex problems will ultimately require $C=\varmathbb R^{n}$ . For convex problems, that is, when the cost function $\varphi$ is convex without any among $f_{i}$ and $g$ being necessarily so, the following requirement will instead suffice to our purposes in the randomized sampling setting of ( $\mathcal{S}_{1}$ ).

For $i\in[N]$ , $\operatorname{dom}h_{i}$ is closed, and whenever $\operatorname{int}\operatorname{dom}h_{i}\ni z^{k}\to z\in\operatorname{bdry}\operatorname{dom}h_{i}$ it holds that $\operatorname{D}_{h_{i}}(z,z^{k})\to 0$ .

II is vacuously satisfied when $\operatorname{dom}h_{i}=\varmathbb R^{n}$ , having $\operatorname{bdry}\varmathbb R^{n}=\emptyset$ .

While II always holds on $\varmathbb R$ , it may fail in higher dimensions [8, Ex. 7.32].

For any $i\in[N]$ , function $h_{i}$ complies with II iff so does $\hat{h}_{i}$ , owing to the inequalities $\tfrac{N-\gamma_{i}L_{f_{i}}}{N\gamma_{i}}\operatorname{D}_{h_{i}}{}\leq{}\operatorname{D}_{\hat{h}_{i}}{}\leq{}\tfrac{N+\gamma_{i}L_{f_{i}}}{N\gamma_{i}}\operatorname{D}_{h_{i}}$ (cf. (3.7)). ∎

Suppose that I holds, and consider the iterates generated by Algorithm 1. Let $x_{i}^{k}={\nabla}\hat{h}_{i}^{\ast}(s_{i}^{k})$ and $z^{k}=u^{k}$ be the corresponding iterates generated by Algorithm 3 as in 3.2, and suppose that

$\operatorname{D}_{\hat{H}}(\bm{u}^{k},\bm{x}^{k})\to 0$ (or equivalently, $\operatorname{D}_{\hat{h}_{i}^{\ast}}(s_{i}^{k},{\nabla}\hat{h}_{i}(z^{k})){}\to{}0$ , $i\in[N]$ ).

Then, letting $\varphi_{\star}$ be as in Item 2, the following hold:

$\varphi(z^{k})=\Phi(\bm{u}^{k})\to\varphi_{\star}$ as $k\to\infty$ .

If $\operatorname{dom}h_{i}$ is closed, $i\in[N]$ , then having (a) $(z^{k})_{k\in K}\to z$ , (b) $(x_{i}^{k})_{k\in K}\to z$ $\exists i\in[N]$ , and (c) $(z^{k+1},x_{i}^{k+1})_{k\in K}\to(z,z)$ $\forall i\in[N]$ , are all equivalent conditions. In particular, if $(z^{k})_{k\in\varmathbb N}$ is bounded (e.g., when $\varphi+\operatorname{\delta}_{\overline{C}}$ is level bounded), then $\|z^{k+1}-z^{k}\|\to 0$ holds, and the set of its limit points, be it $\omega$ , is thus nonempty, compact, and connected.

Under II, $\varphi\equiv\varphi_{\star}$ on $\omega$ (the set of limit points of $(z^{k})_{k\in\varmathbb N}$ ).

If $C=\varmathbb R^{n}$ , then every $z^{\star}\in\omega$ is stationary for (P).

Assumption 1 can be written as $\operatorname{D}_{\hat{h}_{i}}(z^{k},{\nabla}\hat{h}_{i}^{\ast}(s_{i}^{k})){}\to{}0$ , $i\in[N]$ . In turn, by the conjugate identity in Item 3, the equivalent expression in terms of $s_{i}^{k}$ and $z^{k}$ is obtained.

1 As shown in Item 2, $(\Phi^{\hat{H}}(\bm{x}^{k}))_{k\in\varmathbb N}$ monotonically decreases to $\varphi_{\star}$ . In turn, Item 5 and Assumption 1 then imply that $\varphi(z^{k})=\Phi(\bm{u}^{k})$ converges to $\varphi_{\star}$ .

2 The equivalences owe to Items 5and 5 (as $\operatorname{dom}\hat{h}_{i}=\operatorname{dom}h_{i}$ ), and imply $\|z^{k+1}-z^{k}\|\to 0$ if $(z^{k})_{k\in\varmathbb N}$ is bounded. The claim on $\omega$ then follows from [19, Rem. 5].

3 Let $z^{\star}\in\omega$ be fixed, and let $(z^{k})_{k\in K}$ be a subsequence converging to $z^{\star}$ . Assertion 2 ensures that $(\bm{x}^{k})_{k\in K}\to\bm{z}^{\star}\coloneqq(z^{\star},\dots,z^{\star})$ , hence

4 Suppose that $C=\varmathbb R^{n}$ and $(z^{k})_{k\in K}\to z^{\star}$ for some infinite $K\subseteq\varmathbb N$ and $z^{\star}\in\varmathbb R^{n}$ , so that, by virtue of assertion 2, $(\bm{x}^{k},\bm{x}^{k+1})_{k\in K}\to(\bm{z}^{\star},\bm{z}^{\star})$ . Since $(z^{k},\ldots,z^{k})=\bm{u}^{k}\in\operatorname{prox}_{\Phi}^{\hat{H}}(\bm{x}^{k})$ , the osc of $\operatorname{prox}_{\Phi}^{\hat{H}}$ (Item 1) ensures that $\bm{z}^{\star}{}\in{}\operatorname{prox}_{\Phi}^{\hat{H}}(\bm{z}^{\star})$ , hence $0\in\hat{\partial}\Phi(\bm{z}^{\star})$ owing to Item 3. By invoking Item 4 we conclude that $z^{\star}$ is stationary for (P). ∎

The analysis for the randomized case dealt in this section will make use of the following result, known as the Robbins-Siegmund supermartingale theorem, and stated here in simplified form following [13, Prop. 2].

For $k\in\varmathbb N$ , let $\xi_{k}$ and $\eta_{k}$ be random variables, and $\mathcal{F}_{k}\subseteq\mathcal{F}_{k+1}$ be sets of random variables such that

$\mathcal{F}_{k}\subseteq\mathcal{F}_{k+1}$ ;

$0\leq\xi_{k},\eta_{k}$ are functions of the random variables in $\mathcal{F}_{k}$ ;

Then, almost surely, $\sum_{k\in\varmathbb N}\eta_{k}<\infty$ and $\xi_{k}$ converges to a (nonnegative) random variable.

Suppose that I holds. Then, denoting $p_{\rm min}=\min_{i}p_{i}$ , the iterates generated by Algorithm 1 with indices selected according to the randomized rule ( $\mathcal{S}_{1}$ ) satisfy

where $\mathcal{L}$ is as in (1.1) (and satisfies $\mathcal{L}(z^{k},\bm{s}^{k})=\Phi^{\hat{H}}(\bm{x}^{k})$ , cf. Item 6). Moreover, letting $\omega$ denote the set of limit points of $(z^{k})_{k\in\varmathbb N}$ , the following assertions hold almost surely:

The sequence $(\operatorname{D}_{\hat{h}^{\ast}_{i}}(s_{i}^{k},{\nabla}\hat{h}_{i}(z^{k})))_{k\in\varmathbb N}$ has finite sum (and in particular vanishes), $i\in[N]$ .

The sequence $(\varphi(z^{k}))_{k\in\varmathbb N}$ converges to the finite value $\varphi_{\star}\leq\varphi(x^{\rm init})$ of Item 2.

If II is satisfied, then $\varphi\equiv\varphi_{\star}$ on $\omega$ .

If $C=\varmathbb R^{n}$ , then $0\in\hat{\partial}\varphi(z^{\star})$ for every $z^{\star}\in\omega$ .

When $\varphi$ is convex (without $g$ or any $f_{i}$ necessarily being so) and $\operatorname{dom}h_{i}$ is closed, $i\in[N]$ ,

the limit points of $(z^{k})_{k\in\varmathbb N}$ all belong to $\operatorname*{arg\,min}_{\overline{C}}\varphi$ ;

if either II holds and $\varphi+\operatorname{\delta}_{\overline{C}}$ is level bounded, or $C=\varmathbb R^{n}$ , then $(z^{k})_{k\in\varmathbb N}$ and $({\nabla}\hat{h}_{i}^{*}(s_{i}^{k}))_{k\in\varmathbb N}$ , $i\in[N]$ , converge to the same point in $\operatorname*{arg\,min}_{\overline{C}}\varphi$ .

By 3.2, we will consider the simpler setting of Algorithm 3. We have

which is (4.1). We thus infer from 4.5 that $(\operatorname{D}_{\hat{H}}(\bm{u}^{k},\bm{x}^{k}))_{k\in\varmathbb N}$ has almost surely finite sum, and the proof of assertions 1–4 then follows from 4.4.

In what follows, suppose that $\varphi$ is convex and $\operatorname{dom}h_{i}$ is closed, $i\in[N]$ , so that $\overline{C}{}={}\bigcap_{i=1}^{N}\overline{\operatorname{int}\operatorname{dom}h_{i}}{}={}\bigcap_{i=1}^{N}\operatorname{dom}h_{i}$ [14, Prop. 1.3.8], and in particular $\operatorname{D}_{\hat{h}_{i}}(y,x)<\infty$ holds for any $(y,x)\in\operatorname{dom}h_{i}\times\operatorname{int}\operatorname{dom}h_{i}\supseteq\overline{C}\times C$ .

5 For any $x^{\star}\in\operatorname*{arg\,min}_{\overline{C}}\varphi$ , so that $\bm{x}^{\star}\coloneqq(x^{\star},\ldots,x^{\star})\in\operatorname*{arg\,min}_{\overline{C}\times\cdots\times\overline{C}}\Phi$ (cf. Item 5), the three-point identity (Item 2), convexity of $\Phi$ (Item 7) and the inclusion ${\nabla}\hat{H}(\bm{x}^{k})-{\nabla}\hat{H}(\bm{u}^{k}){}\in{}\hat{\partial}\Phi(\bm{u}^{k})$ (Item 3) give the bound

where $\bm{u}^{k}=(u^{k},\ldots,u^{k})$ . From 4.5 we conclude that

6 Suppose that $(z^{k})_{k\in K}\to z^{\star}$ . Then, $(\bm{u}^{k})_{k\in K}\to\bm{u}^{\star}$ for $\bm{u}^{k}=(z^{k},\ldots,z^{k})$ and $\bm{u}^{\star}=(z^{\star},\ldots,z^{\star})$ . Notice that $z^{\star}\in\overline{C}$ , since $z^{k}\in C$ for all $k$ (cf. 4.2). We have

Therefore, $\bm{u}^{\star}$ is a minimizer of $\Phi$ on $\overline{C}\times\cdots\times\overline{C}$ , and thus $z^{\star}$ is a minimizer of $\varphi$ on $\overline{C}$ by virtue of Item 5.

7 If $\varphi+\operatorname{\delta}_{\overline{C}}$ is level bounded, then by Item 4 $(\bm{x}^{k})_{k\in\varmathbb N}$ and $(\bm{u}^{k})_{k\in\varmathbb N}$ are bounded. Alternatively, if $C=\varmathbb R^{n}$ , then boundedness of the former sequence follows from Item 4, (4.5) and Assumption 3, and in turn that of the latter from (4.4). In either cases II holds, as discussed in Item 1. Boundedness of the sequences ensures the existence of $K\subseteq\varmathbb N$ , $z^{\star}$ and $\bm{u}^{\star}$ as in the proof of assertion 6. The vanishing of $\operatorname{D}_{\hat{H}}(\bm{u}^{k},\bm{x}^{k})$ shown in (4.4) implies through Item 2 that $(\bm{x}^{k})_{k\in\varmathbb N}$ and $(\bm{u}^{k})_{k\in\varmathbb N}$ have same limit points, and that $(\bm{x}^{k})_{k\in K}\to\bm{u}^{\star}$ . In turn, $(\sum_{i=1}^{N}p_{i}^{-1}\operatorname{D}_{\hat{h}_{i}}(u^{\star},x_{i}^{k}))_{k\in K}\to 0$ holds by II. Hence, since the entire sequence is convergent (by (4.5)) we have $(\sum_{i=1}^{N}p_{i}^{-1}\operatorname{D}_{\hat{h}_{i}}(u^{\star},x_{i}^{k}))_{k\in\varmathbb N}\to 0$ , which by Item 5 implies $(x_{i}^{k})_{k\in\varmathbb N}\to u^{\star}$ , $i\in[N]$ . As discussed above, this implies that $(\bm{u}^{k})_{k\in\varmathbb N}\to\bm{u}^{\star}$ , and the identity $\bm{u}^{k}=(z^{k},\ldots,z^{k})$ of Item 2 yields the claimed convergence.∎

In Item 7 the assumption that $\varphi+\operatorname{\delta}_{\overline{C}}$ is level bounded can be relaxed by instead requiring that for every $v\in\operatorname{dom}h_{i}$ and $\alpha\in\varmathbb R$ , the level set ${\mathopen{}\left\{w\in\operatorname{int}\operatorname{dom}h_{i}{}\mathrel{\mid}{}\operatorname{D}_{h_{i}}(v,w)\leq\alpha\right\}\mathclose{}}$ is bounded, as this would suffice to ensure boundedness of the sequences. In fact, together with the closed-domain requirement this is a standing assumption in many works dealing with Bregman distances, specifically those involving Bregman functions “with zone $S$ ” ( $S$ being the interior of the domain), see e.g., .

We conclude this subsection with an analysis of the strongly convex case, in which linear convergence (in expectation) will be shown. Remarkably, strong convexity of the cost function $\varphi$ alone will suffice, without imposing any such requirement on the individual terms $f_{i}$ or $g$ which, in fact, are even allowed to be nonconvex.

Consider the iterates of Algorithm 1. Additionally to I, suppose that

$\varphi$ is $\mu_{\varphi}$ -strongly convex;

$h_{i}$ has locally Lipschitz gradient on the whole space $\varmathbb R^{n}$ , $i\in[N]$ , (hence $C=\varmathbb R^{n}$ ).

$g$ is convex, and $f\coloneqq\frac{1}{N}\sum_{i}f_{i}$ is $\mu_{f}$ -strongly convex (yet each $f_{i}$ can be nonconvex);

The convergence results in this subsection require convexity of the nonsmooth term $g$ and local strong convexity and smoothness of $h_{i}$ (as is the case when $h_{i}\in\mathcal{C}^{2}(\varmathbb R^{n})$ with $\nabla^{2}h_{i}\succ 0$ ). The proof of subsequential convergence is an adaptation of that of [37, Thm. 2.8].

Additionally to I, assume that $g$ is convex, $C=\varmathbb R^{n}$ , and that either one of the following assumptions holds:

(either) each $h_{i}$ is strongly convex and Lipschitz differentiable,

(or) $\varphi$ is level bounded and each $h_{i}$ is locally strongly convex and locally Lipschitz differentiable.

Then, all the claims in Items 1, 2, 3and 4 hold surely.

We remark that linear convergence in the strongly convex case may be obtained in a similar fashion to [37, Thm. 2.9] and [12, Thm. 3.9]. This however results in a rate more conservative than the one obtained for the randomized case, cf. [37, Eq. (2.21)], which is not consistent with what is observed in practice. Although a refined analysis not relying on conservative triangle inequalities may be possible, this direction is not investigated here.

Our next goal is to establish global (and linear) convergence results without convexity assumptions on $f_{i}$ or their sum. To this end, we leverage the Kurdyka-Łojasiewicz property , which has become the standard tool in the analysis of nonconvex proximal methods, and most notably holds for the class of semialgebraic functions .

A proper lsc function $q:\varmathbb R^{n}\rightarrow\overline{\varmathbb}R$ has the Kurdyka-Łojasiewicz (KL) property with exponent $\theta\in(0,1)$ if for every $\bar{w}\in\operatorname{dom}\partial q$ there exist $\varepsilon,\eta,\varrho>0$ such that $\psi^{\prime}(q(w)-q(\bar{w}))\operatorname{dist}(0,\partial q(w))\geq 1$ holds for every $w$ satisfying $\|w-\bar{w}\|<\varepsilon$ and $q(\bar{w})<q(w)<q(\bar{w})+\eta$ , where $\psi(s)\coloneqq\varrho s^{1-\theta}$ .

As will be detailed in 4.11, global convergence is established when the model $\mathcal{M}:\varmathbb R^{nN}\times\varmathbb R^{nN}\rightarrow\overline{\varmathbb}R$ defined as $\mathcal{M}(\bm{w},\bm{x})\coloneqq\Phi(\bm{w})+\operatorname{D}_{\hat{H}}(\bm{w},\bm{x})$ has the KL property. The next assumption provides easily verifiable requirements in terms of $f_{i}$ , $h_{i}$ , and $\varphi$ .

for $i\in[N]$ , $f_{i},h_{i}\in\mathcal{C}^{2}(\varmathbb R^{n})$ (hence $C=\varmathbb R^{n}$ ) with $\nabla^{2}h_{i}\succ 0$ ;

$\varphi$ has the KL property with exponent $\theta\in(0,1)$ (e.g., when $f_{i}$ and $g$ are semialgebraic) and is level bounded.

Suppose that s Iand III are satisfied and that $g$ is convex. Then, the following hold for the iterates generated by Algorithm 1 with an essentially cyclic rule ( $\mathcal{S}_{2}$ ):

$(z^{k})_{k\in\varmathbb N}$ converges to a stationary point $z^{\star}$ for $\varphi$ .

If $\theta>\nicefrac{{1}}{{2}}$ , then there exists $c>0$ such that $\varphi(z^{k})-\varphi(z^{\star}){}\leq{}ck^{-\frac{1}{2\theta-1}}$ hold for all $k\in\varmathbb N$ .

If $\theta\in(0,\nicefrac{{1}}{{2}}]$ , then $(z^{k})_{k\in\varmathbb N}$ and $(\varphi(z^{k}))_{k\in\varmathbb N}$ converge at $R$ -linear rate.

Notice that Assumption 1 along with level boundedness of $\varphi$ in Assumption 2 ensures that the requirement in 2 is satisfied. By the the first claim of 4.9, $(\operatorname{D}_{\hat{h}^{\ast}_{i}}(s_{i}^{k},{\nabla}\hat{h}_{i}(z^{k})))_{k\in\varmathbb N}$ converges to zero, and thus we may invoke 4.4 to conclude that the set $\omega$ of limit points of $(z^{k})_{k\in\varmathbb N}$ is nonempty, compact, connected, and made of stationary points for $\varphi$ , with $\varphi{}\equiv{}\varphi_{\star}{}\coloneqq{}\lim_{k\to\infty}\Phi^{\hat{H}}(\bm{x}^{k})$ on $\omega$ . If $\Phi^{\hat{H}}(\bm{x}^{k})=\varphi_{\star}$ holds for some $k\in\varmathbb N$ , then it follows from Item 1 that $(\bm{x}^{k})_{k\in\varmathbb N}$ is asymptotically constant, and the assertion holds trivially. In what follows we thus assume that $\Phi^{\hat{H}}(\bm{x}^{k})>\varphi_{\star}$ holds for all $k$ . The assumptions together with Item 3 ensure that $\Phi$ enjoys the KL property with exponent $\theta$ . Since $\hat{H}$ is locally strongly convex, we may invoke [68, Lem. 5.1] to infer that the function $\mathcal{M}:\varmathbb R^{nN}\times\varmathbb R^{nN}\rightarrow\overline{\varmathbb}R$ defined as $\mathcal{M}(\bm{w},\bm{x})=\Phi(\bm{w})+\operatorname{D}_{\hat{H}}(\bm{w},\bm{x})$ has the KL property with exponent $\vartheta\coloneqq\max{\mathopen{}\left\{\theta,\nicefrac{{1}}{{2}}{}\mathrel{\mid}{}{}\right\}\mathclose{}}$ at every point of the compact set $\bm{\Omega}\coloneqq{\mathopen{}\left\{(\bm{z}^{\star},\bm{z}^{\star}){}\mathrel{\mid}{}\bm{z}^{\star}\in\bm{\omega}\right\}\mathclose{}}$ , where $\bm{\omega}\coloneqq{\mathopen{}\left\{\bm{z}=(z,\ldots,z){}\mathrel{\mid}{}z\in\omega\right\}\mathclose{}}$ .Consistently with the locality of the KL property and the compactness of $\bm{\omega}$ , the global strong convexity requirement in [68, Lem. 5.1] can clearly be replaced by local strong convexity. Similarly, if $\Phi$ is a KL function with exponent $\theta$ , then it is trivially a KL function with exponent $\vartheta$ , thus complying with the requirement in the reference. Notice that $\Phi^{\hat{H}}(\bm{x}^{k}){}={}\mathcal{M}(\bm{u}^{k},\bm{x}^{k})$ and $\mathcal{M}(\bm{z}^{\star},\bm{z}^{\star}){}={}\Phi^{\hat{H}}(\bm{z}^{\star}){}={}\varphi_{\star}$ hold for every $k\in\varmathbb N$ and $\bm{z}^{\star}\in\bm{\omega}$ (cf. 4.9), and that $\partial\mathcal{M}(\bm{w},\bm{x}){}={}{\mathopen{}\left(\partial\Phi(\bm{w})+{\nabla}\hat{H}(\bm{w})-{\nabla}\hat{H}(\bm{x}),\nabla^{2}\hat{H}(\bm{x})(\bm{x}-\bm{w})\right)\mathclose{}}$ . By Item 3 we have ${\nabla}\hat{H}(\bm{x}^{k})-{\nabla}\hat{H}(\bm{u}^{k})\in\partial\Phi(\bm{u}^{k})$ , which in turn implies

where $C=\sup_{k}\|\nabla^{2}\hat{H}(\bm{x}^{k})\|$ is finite due to boundedness of $(\bm{x}^{k})_{k\in\varmathbb N}$ and continuity of $\nabla^{2}\hat{H}$ . Let $\psi(t)\coloneqq\rho t^{1-\vartheta}$ be a desingularizing function for $\mathcal{M}$ on $\bm{\Omega}$ [3, Lem. 1(ii)], namely such that

holds for some $\varepsilon>0$ and all $(\bm{w},\bm{x})$ $\varepsilon$ -close to $\bm{\Omega}$ such that $0<\mathcal{M}(\bm{w},\bm{x})-\varphi_{\star}<\varepsilon$ . Since $\mathcal{M}(\bm{u}^{k},\bm{x}^{k})=\Phi^{\hat{H}}(\bm{x}^{k})\searrow\varphi_{\star}$ (cf. Item 2) and $(\bm{u}^{k},\bm{x}^{k})_{k\in\varmathbb N}$ is bounded and accumulates on $\bm{\Omega}$ , by discarding early iterates we may assume that the inequality above holds for $(\bm{w},\bm{x})=(\bm{u}^{k},\bm{x}^{k})$ , $k\in\varmathbb N$ , which combined with (4.10) results in

Let $\Delta_{k}{}\coloneqq{}\psi(\Phi^{\hat{H}}(\bm{x}^{k})-\varphi_{\star})$ , so that $\Phi^{\hat{H}}(\bm{x}^{k})-\varphi_{\star}{}={}(\Delta_{k}/\rho)^{\nicefrac{{1}}{{1-\theta}}}$ . By concavity of $\psi$ we have

for some constant $c>0$ . As argued in the proof of 4.9, by suitably shifting, we conclude that for all $k\in\varmathbb N$

The rest of the proof is standard (see [3, Thm. 2]), and is provided for completeness. Since $\Delta_{k}\geq 0$ , by telescoping we conclude that $(\|\bm{u}^{k}-\bm{x}^{k}\|)_{k\in\varmathbb N}$ has finite sum, and since $\|\bm{x}^{k+1}-\bm{x}^{k}\|\leq\|\bm{u}^{k}-\bm{x}^{k}\|$ for all $k\in\varmathbb N$ , $(\bm{x}^{k})_{k\in\varmathbb N}$ has finite length. Therefore, $(\bm{x}^{k})_{k\in\varmathbb N}$ and $(\bm{u}^{k})_{k\in\varmathbb N}$ converge to the same stationary point, be it $\bm{z}^{\star}$ , owing to 4.9.

We now show the convergence rates. It follows from (4.11) that

Combined with (4.13), it results in $c_{1}\Delta_{k}^{\nicefrac{{\vartheta}}{{1-\vartheta}}}{}\leq{}\Delta_{k}-\Delta_{k+T}$ for some $c_{1}>0$ . We may now invoke A.3 to infer that, for every $t\in[T]$ , $(\Delta_{t+\nu T})_{\nu\in\varmathbb N}$ converges $Q$ -linearly (to 0) if $\theta\leq\nicefrac{{1}}{{2}}$ (which corresponds to $\vartheta=\nicefrac{{1}}{{2}}$ ), and $\Delta_{t+\nu T}\leq c_{3}\nu^{-\frac{1-\theta}{2\theta-1}}$ for some $c_{3}>0$ otherwise (that is, if $\vartheta=\theta>\nicefrac{{1}}{{2}}$ ). Note that the former case implies that $(\Delta_{k})_{k\in\varmathbb N}$ converges $R$ -linearly, whereas the latter implies that $\Delta_{k}\leq c_{4}k^{-\frac{1-\theta}{2\theta-1}}$ holds for every $k$ and some $c_{4}\geq c_{3}$ . Recalling that $\Phi^{\hat{H}}(\bm{x}^{k})-\varphi_{\star}{}={}(\Delta_{k}/\rho)^{\nicefrac{{1}}{{1-\theta}}}$ , the claimed rates of for the cost function follow from Item 1. Similarly, when $(\Delta_{k})_{k\in\varmathbb N}$ converges $Q$ -linearly then so does $(\|\bm{x}^{k}-\bm{u}^{k}\|)_{k\in\varmathbb N}$ , as it follows from (4.13), and in turn so does $(\|\bm{x}^{k}-\bm{x}^{k+1}\|)_{k\in\varmathbb N}$ owing to the inequality $\|\bm{x}^{k}-\bm{x}^{k+1}\|\leq\|\bm{x}^{k}-\bm{u}^{k}\|$ . These two facts imply that $(\bm{x}^{k})_{k\in\varmathbb N}$ and $(\bm{u}^{k}=(z^{k},\dots,z^{k}))_{k\in\varmathbb N}$ are $R$ -linearly convergent. ∎

4 Low-memory variant

We now analyze Algorithm 2, which, as shown next, is simply a particular implementation of Algorithm 1.

As a consequence of 4.12, Algorithm 2 inherits all the convergence results of Section 4.3. In addition, here, the convexity requirement of $g$ in s 4.9and 4.11 can be lifted thanks to the periodic full sampling of the indices.

if II is satisfied, then $\varphi\equiv\varphi_{\star}$ on $\omega$ ;

if $C=\varmathbb R^{n}$ , then $0\in\hat{\partial}\varphi(z^{\star})$ for every $z^{\star}\in\omega$ .

whence assertion 1 follows. Assertions 2 and 3 follow by patterning the arguments of Items 3and 4. ∎

In the next theorem global convergence results are provided under III in a fully nonconvex setting. Moreover, in the spirit of 4.11, linear and sublinear convergence rates are obtained according to the KL exponent.

Suppose that s Iand III hold. Then, the following hold for the iterates generated by Algorithm 2:

Notice that III entails local strong convexity and Lipschitz differentiability of each $h_{i}$ . Thus, as discussed in the proof of 4.9, $\hat{H}$ is $\mu_{\hat{H},\bm{\mathcal{U}}}$ -strongly convex on a convex compact set $\mathcal{\bm{\mathcal{U}}}$ containing the iterates. Let the indexing subsequence $(k_{r})_{r\in\varmathbb N}$ be as in the proof of 4.13, and observe that $k_{r}\leq(N+1)r$ holds for every $r\in\varmathbb N$ . The assertions are established with the same arguments as in 4.11 with the difference that here we examine the generated iterates at subindices $k_{r}$ . That is, (4.12) is replaced by

Application to phase retrieval and numerical simulations

In this section we study two examples related to the phase retrieval problem, which consists of recovering a signal based on intensity measurements, and arises in many important applications including X-ray crystallography, speech processing, electron microscopy, astronomy, and optical imaging; see, e.g., . Here, we consider phase retrieval problems with real-valued data, that is, given $a_{i}\in\varmathbb R^{n}\setminus{\mathopen{}\left\{0{}\mathrel{\mid}{}{}\right\}\mathclose{}}$ and scalars $b_{i}\in\varmathbb R_{+}$ , $i\in[N]$ , the goal is to find $x\in\varmathbb R^{n}$ such that

accounting for the fact that in real-world applications the recorded intensities are likely corrupted by noise, and may involve outliers due to measurement errors. To tackle such problems, we consider the following sparse phase retrieval formulation:

where $\mathscr{L}$ is a loss function, and $g$ is a sparsity inducing function (e.g., $l_{1}$ - or $l_{0}$ -norm). In particular, we study the case of squared loss $\mathscr{L}(y,z)=\tfrac{1}{4}(y-z)^{2}$ , and Poisson loss $\mathscr{L}(y,z)=z-y\log(z)$ , suitable when measurements follow the Poisson model ( $b_{i}\approx\operatorname{Poisson}({\mathopen{}\left\langle{}a_{i}{},{}x{}\right\rangle\mathclose{}}^{2})$ ). Other formulations with $l_{1}$ -loss have been studied in the literature .

Let $f_{i}$ and $h_{i}$ be as in (5.3). Then, $f_{i}$ is $L_{f_{i}}$ -smooth relative to $h_{i}$ with $L_{f_{i}}=(3\|a_{i}\|^{4}+\|a_{i}\|^{2}|b_{i}|)$ . Moreover, denoting $\bar{\gamma}{}={}(\sum_{i=1}^{N}\nicefrac{{1}}{{\gamma_{i}}})^{-1}$ , for any $y(s){}\in{}\operatorname{prox}_{\bar{\gamma}g}(\bar{\gamma}s)$ the operator $T$ as defined in (2.1) may be computed as follows:

If $g=\lambda\|{}\cdot{}\|_{1}$ , $\lambda\geq 0$ , then $T(s)=t^{\star}y(s)$ , where $t^{\star}$ is the real positive root of the equation $\|y(s)\|^{2}t^{3}+t-1=0$ .Nonnegative real roots of the cubic equation $t^{3}+pt+q=0$ for some $p>0$ and $q\leq 0$ are given by Cardano’s formula $t^{\star}{}={}(c-\nicefrac{{q}}{{2}})^{\nicefrac{{1}}{{3}}}{}-{}\big{(}c+\nicefrac{{q}}{{2}})^{\nicefrac{{1}}{{3}}}$ , where $c{}={}(\nicefrac{{q^{2}}}{{4}}+\nicefrac{{p^{3}}}{{27}})^{\nicefrac{{1}}{{2}}}$ ; see, e.g., .

2 Sparse phase retrieval with Poisson loss

We now assume that the recorded intensities follow the Poisson model ( $b_{i}\sim\operatorname{Poisson}({\mathopen{}\left\langle{}a_{i}{},{}x{}\right\rangle\mathclose{}}^{2})$ ). In this setting we adapt the Poisson loss $\mathscr{L}(y,z)=z-y\log(z)$ and consider the regularized problem (5.2) with $g=\lambda\|{}\cdot{}\|_{1}$ , $\lambda\geq 0$ . This problem may be written in the form of (P) by setting

As shown next, the nonconvex function $f_{i}$ is smooth relative to $h_{i}$ , and the operator $T$ as in (2.1) is easily computable.

Let $f_{i}$ and $h_{i}$ be as in (5.4), with $a_{i}\in\varmathbb R^{n}_{+}\setminus{\mathopen{}\left\{0{}\mathrel{\mid}{}{}\right\}\mathclose{}}$ and $b=(b_{1},\ldots,b_{N})\in\varmathbb R^{n}_{+}\setminus{\mathopen{}\left\{0{}\mathrel{\mid}{}{}\right\}\mathclose{}}$ . Then, $h_{i}$ is a $2\|a_{i}\|^{2}$ -strongly convex Legendre kernel (with $\operatorname{dom}h_{i}=\varmathbb R_{++}^{n}$ ), and $f_{i}$ is $L_{f_{i}}$ -smooth relative to $h_{i}$ with $L_{f_{i}}=1$ . Moreover, denoting $c_{a}{}={}\sum_{i=1}^{N}\tfrac{4}{\gamma_{i}}\|a_{i}\|^{2}$ and $c_{b}{}={}\sum_{i=1}^{N}\tfrac{4}{\gamma_{i}}b_{i}$ , the operator $T$ as defined in (2.1) with $g=\lambda\|{}\cdot{}\|_{1}$ , $\lambda\geq 0$ , is given by

where the first inequality follows by the Cauchy-Schwarz inequality, the second one from Jensen’s inequality $\langle a,y\rangle^{2}=\sum_{j=1}^{n}\big{(}a_{j}x_{j}(\tfrac{y_{j}}{x_{j}})\big{)}^{2}\leq\sum_{i=1}^{n}a_{i}x_{i}\sum_{j=1}^{n}a_{j}x_{j}(\tfrac{y_{j}}{x_{j}})^{2}$ , and the third one from the fact that $\sum_{i}\alpha_{i}\sum_{j}\beta_{j}\geq\sum_{i}\alpha_{i}\beta_{i}$ for every $\alpha_{i},\beta_{i}\geq 0$ , $i\in[n]$ . The closed-form solution for the proximal mapping $T(s)$ follows directly from its first-order optimality conditions. ∎

3 Experimental setup

For Algorithms 1and 2, we always use $\gamma_{i}=\nicefrac{{0.99N}}{{L_{f_{i}}}}$ . For SMD we used the popular square-summable stepsize $\gamma_{k}=\nicefrac{{\alpha}}{{(L_{f}k)}}$ , where $k$ is the iteration counter, $L_{f}$ is the smoothness modulus of $f=\nicefrac{{1}}{{N}}\sum_{i=1}^{N}f_{i}$ relative to a suitable Bregman kernel $h$ , and $\alpha>0$ is tuned for performance. In particular, for SMD in the problems described above, $L_{f}=\sum_{i=1}^{N}\frac{1}{N}{L_{f_{i}}}$ and $h(x)=\frac{1}{4}\|x\|^{4}+\frac{1}{2}\|x\|^{2}$ for simulations related to Section 5.1, and $h(x)=\frac{1}{N}\sum_{i=1}^{N}\|a_{i}\|^{2}\|x\|^{2}-\frac{2}{N}\sum_{i=1}^{N}b_{i}\sum_{j=1}^{n}\log(x_{j})$ for those related to Section 5.2 are used.

As a measure of suboptimality, we consider

In the first set of simulations we consider $16\times 16$ gray-scale images from a digits dataset https://web.stanford.edu/~hastie/ElemStatLearn/data.html. and a QR code dataset . The images are vectorized resulting in the signal $x\in\varmathbb R^{n}$ with $n=256$ . The data matrix $A\in\varmathbb R^{N\times n}$ ( $a_{i}$ being the $i$ -th row) with $N=nd$ , $d=5$ is generated following the procedure described in [29, §6.3]. Let $M\in\varmathbb R^{n\times n}$ be a normalized Hadamard matrix. We generate $d$ many i.i.d. diagonal sign matrices $S_{i}$ with diagonal elements in ${\mathopen{}\left\{-1,1{}\mathrel{\mid}{}{}\right\}\mathclose{}}$ selected uniformly at random, and set $A{}={}{[MS_{1},\dots,MS_{d}]}$ . Typically $d\geq 3$ is sufficient for near complete recovery on noiseless data. In our simulations, we corrupted a fraction of the measurements $b_{i}=\langle{}a_{i}{},{}x{}\rangle^{2}$ independently by setting $b_{i}=0$ with probability $p_{\rm c}=\nicefrac{{1}}{{50}}$ . All of the plotted algorithms are initialized using the initialization scheme described in [29, §3].

In the last set of simulations we consider synthetic data. We generate a standard random Gaussian matrix $A\in\varmathbb R^{N\times n}$ with $n=200$ , $N\in{\mathopen{}\left\{400,1000{}\mathrel{\mid}{}{}\right\}\mathclose{}}$ . The data vector $a_{i}$ , $i\in[N]$ , is set equal to the absolute value of the $i$ -th row of $A$ . We also drew a random vector from $\mathcal{N}(0,I_{n})$ and set the signal $x$ equal to its absolute value. We generated the measurements according to the Poisson model $b_{i}\sim\operatorname{Poisson}(\langle{}a_{i}{},{}x{}\rangle^{2})$ , $i\in[N]$ , and further corrupted the measurements $b_{i}$ by setting them equal to the nearest integer to the absolute value of $\|x\|^{2}\mathcal{N}(0,1)$ with probability $p_{\rm c}=\nicefrac{{1}}{{10}}$ . All methods were initialized at the same random point. We ran simulations with regularization parameter $\lambda\in{\mathopen{}\left\{\nicefrac{{0.01}}{{N}},\nicefrac{{0.1}}{{N}},\nicefrac{{1}}{{N}}{}\mathrel{\mid}{}{}\right\}\mathclose{}}$ . We only report the results for $\lambda=\nicefrac{{0.1}}{{N}}$ due to space limitations, nevertheless remarking that for other values similar plots were observed. The results are illustrated in Fig. 3. Similar to the previous experiments, SMD performed the worst, while the best results are observed for Algorithm 1 with cyclic rule ( $\mathcal{S}_{2}^{\text{\sc cycl}}$ ). The low-memory Algorithm 2 usually outperforms the randomized variant of Algorithm 1.

Conclusions

A Bregman incremental aggregated method was developed that extends Finito/MISO to non-Lipschitz and nonconvex settings. The basic algorithm was studied under randomized and essentially cyclic sampling strategies. Furthermore, a variant with $O(n)$ memory requirements is developed that is novel even in the Euclidean case. A sure descent property established on a Bregman Moreau envelope leads to a surprisingly simple convergence analysis. As one particularly interesting result, in the randomized setting linear convergence is established under strong convexity of the cost without requiring convexity of the individual functions $f_{i}$ or $g$ . Future research directions include extending the analysis to the framework of the Douglas-Rachford splitting, momentum-type schemes, as well as applications to nonconvex distributed asynchronous optimization.

Appendix A Auxiliary results

Let $\Phi$ and $\Delta$ be as in (3.1). Then:

cost: $\Phi(\bm{x})=\varphi(x)$ if $\bm{x}=(x,\ldots,x)$ , and $\Phi(\bm{x})=\infty$ otherwise.

subdifferential: $\hat{\partial}(\Phi+\operatorname{\delta}_{\overline{C}\times\cdots\times\overline{C}})(\bm{x}){}={}{\mathopen{}\left\{\bm{v}=(v_{1},\ldots,v_{N})\in\varmathbb R^{nN}{}\mathrel{\mid}{}\sum_{i=1}^{n}v_{i}{}\in{}\hat{\partial}(\varphi+\operatorname{\delta}_{\overline{C}})(x)\right\}\mathclose{}}$ if $\bm{x}=(x,\cdots,x)$ for some $x\in\varmathbb R^{n}$ , and is empty otherwise; the same relation still holds if the regular subdifferential $\hat{\partial}$ is replaced by the limiting subdifferential $\partial$ .

KL property: $\varphi$ has the KL property at $x$ iff so does $\Phi$ at $\bm{x}=(x,\dots,x)$ , in which case the desingularizing functions are the same up to a positive scaling.

stationary points: a point $\bm{x}^{\star}$ is stationary for problem (3.1) iff $\bm{x}^{\star}=(x^{\star},\ldots,x^{\star})$ for some $x^{\star}\in\varmathbb R^{n}$ which is stationary for problem (P).

minimizers: $\bm{x}^{\star}$ is a (local) minimizer of problem (3.1) iff $\bm{x}^{\star}=(x^{\star},\ldots,x^{\star})$ for some $x^{\star}\in\varmathbb R^{n}$ which is a (local) minimizer for problem (P); in fact, $\operatorname*{inf}_{\overline{C}\times\cdots\times\overline{C}}\Phi=\operatorname*{inf}_{\overline{C}}\varphi$ .

level boundedness: $\varphi+\operatorname{\delta}_{\overline{C}}$ is level bounded iff so is $\Phi+\operatorname{\delta}_{\overline{C}\times\cdots\times\overline{C}}$ .

convexity: $\varphi:\varmathbb R^{n}\rightarrow\overline{\varmathbb}R$ is convex iff so is $\Phi:\varmathbb R^{Nn}\rightarrow\overline{\varmathbb}R$ .

For notational convenience, up to possibly replacing $g$ with $g+\operatorname{\delta}_{\overline{C}}$ we may assume without loss of generality that $C=\varmathbb R^{n}$ .

1 Trivial consequence of the fact that $\operatorname{dom}\Phi\subseteq\Delta$ (the consensus set, cf. (3.1)).

2 In light of the previous point, having $\bm{x}=(x,\ldots,x)$ for some $x\in\varmathbb R^{n}$ is necessary for the nonemptiness of $\hat{\partial}\Phi(\bm{x})$ . Let $\bm{x}=(x,\ldots,x)$ and $\bm{v}\in\hat{\partial}\Phi(\bm{x})$ be fixed. Then,

where the equality comes from the fact that $\operatorname{dom}\Phi\subseteq\Delta$ together with assertion 1. This shows that $\sum_{i}v_{i}\in\hat{\partial}\varphi(x)$ . Conversely, let $u\in\hat{\partial}\varphi(x)$ and $\bm{v}\in\varmathbb R^{nN}$ be such that $\sum_{i}v_{i}=u$ . By reading (A.1) from right to left we obtain that $\bm{v}\in\hat{\partial}\Phi(\bm{x})$ . Having shown the identity of the regular subdifferential, the same claim with the limiting subdifferential follows by definition.

where the first and last equalities are due to assertion 2.

4–7 Directly follow from assertions 1 and 2. ∎

$\operatorname{prox}_{\Phi}^{\hat{H}}$ is Lipschitz continuous on $\bm{\mathcal{U}}$ .

If in addition $f_{i}$ and $h_{i}$ are twice continuously differentiable on $\mathcal{U}_{i}$ , $i\in[N]$ , then

$\Phi^{\hat{H}}$ is continuously differentiable on $\bm{\mathcal{U}}$ with ${\nabla}\Phi^{\hat{H}}=\nabla^{2}\hat{H}\circ({\rm id}-\operatorname{prox}_{\Phi}^{\hat{H}})$ ;

Let $(\alpha_{k})_{k\in\varmathbb N}\subset\varmathbb R_{+}$ be a sequence, and suppose that there exist $c>0$ and $\delta\in[1,\infty)$ such that $\alpha_{k+1}^{\delta}{}\leq{}c(\alpha_{k}-\alpha_{k+1})$ holds for every $k\in\varmathbb N$ .

If $\delta=1$ , then $(\alpha_{k})_{k\in\varmathbb N}$ is $Q$ -linearly convergent (to ).

If $\delta\in(1,\infty)$ , then there exists $c^{\prime}>0$ such that $\alpha_{k}{}\leq{}c^{\prime}k^{-\frac{1}{\delta-1}}$ holds for all $k\in\varmathbb N$ .

We will use the equivalent BC-reformulation of Algorithm 3, through the identities shown in 3.2. We start by observing that $x_{i}^{k}\in\{x^{\rm init},z^{k}{}\mathrel{\mid}{}k\in\varmathbb N\}\subseteq\mathcal{U}$ holds for any $k\in\varmathbb N$ and $i\in[N]$ , as it follows from the $\bm{x}$ -update at step 6 and the fact that $u_{i}^{k}=z^{k}$ , cf. Item 2. Let $\bm{x}^{\star}=(x^{\star},\ldots,x^{\star})$ be the unique minimizer of $\Phi$ (cf. Item 5). As shown in (4.2), denoting $\bm{v}^{k}{}\coloneqq{}{\nabla}\hat{H}(\bm{x}^{k})-{\nabla}\hat{H}(\bm{u}^{k}){}\in{}\hat{\partial}\Phi(\bm{u}^{k})$ we have

where the inequality follows from strong convexity of $\varphi$ . For any $\varepsilon_{i}>0$ , $i\in[N]$ , one has

Plugged in (B.1) with $\varepsilon_{i}>0$ such that $\sum_{i=1}^{N}\varepsilon_{i}=\mu_{\varphi}$ , so as to cancel the square norm therein,

Combining this with (4.1) (recall the equivalences in 3.2) yields

Since all indices are updated at least once every $T$ iterations, one has that

is well defined for each index $i\in[N]$ and $\nu\in\varmathbb N$ . In other words, since $i$ is sampled at iteration $\nu T+t_{\nu}(i)-1$ and not in any one between $\nu T$ and $\nu T+t_{\nu}(i)-2$ , it holds that

recalling $\bm{u}^{k}=(u^{k},\ldots,u^{k})$ . We now proceed to establish a descent inequality for $\Phi^{\hat{H}}$ holding every interval of $T$ iterations. First,

holds for all $t\in[T]$ . Next, for every $i\in[N]$ it holds that

where the first inequality uses the $\lambda$ -Lipschitz continuity of $\operatorname{prox}_{\Phi}^{\hat{H}}$ , the second one the triangular inequality, and the last one the fact that $t_{\nu}(i)\leq T$ . For all $i\in[N]$ , it follows from Item 7 and the triangular inequality that

By squaring and summing over $i\in[N]$ we obtain

Since by Eq. $\mathcal{S}_{2}$ in any interval of length $T$ every index is updated at least once, by suitably shifting, for every $t\in[T]$ the same holds for the sequences $(\bm{x}^{\nu T+t})_{\nu\in\varmathbb N}$ and $(\bm{u}^{\nu T+t})_{\nu\in\varmathbb N}$ . Thus,

By telescoping the inequality and using the fact that the envelope is lower bounded (Item 2 and Assumption 3), all the assertions follow from 4.4. ∎

where the first equality follows by step 6, the second one from the induction hypothesis, and the last one from the fact that $s_{i}^{k+1}=s_{i}^{k}$ for $i\notin\mathcal{I}^{k+1}$ as in step 6.

In what follows, let $\mathcal{N}_{\rm full}{}\coloneqq{}{\mathopen{}\left\{k{}\mathrel{\mid}{}\mathcal{K}^{k}=\emptyset\right\}\mathclose{}}$ , and observe that $k\in\mathcal{N}_{\rm full}$ iff the if statement at step 5 is true. In particular, it follows from step 6 that

The claim is true for $k=0$ ; suppose it holds up to iteration $k\geq 0$ . We consider two cases:

Case 1: $k\in\mathcal{N}_{\rm full}$ . We have

where the first equality follows by step 8, the second one by induction, the third one by the fact that $\mathcal{I}^{k+1}=\mathcal{I}_{\text{\sc lm}}^{k+1}=[N]$ (cf. (B.11) and (B.12)), and the last one from (B.10). It follows that the minimization problems defining $z^{k+1}$ and $z_{\text{\sc lm}}^{k+1}$ (at step 4 and step 4, respectively) coincide, thus ensuring that $z^{k+1}=z_{\text{\sc lm}}^{k+1}$ is a feasible update for Algorithm 1.