Randomized Bregman Coordinate Descent Methods for Non-Lipschitz Optimization

I Introduction

In this paper, we consider a composite optimization problem in the following form

where $r$ has $n$ separated blocks. More specifically, we have

where ${\mathbf{x}}_{i}$ denotes a subvector of ${\mathbf{x}}$ with dimension $N_{i}$ such that $\sum_{i=1}^{n}N_{i}=N$, and each $r_{i}$ is a (possibly nonsmooth) convex function.

Due to the block separable structure, Problem (1) can be solved by (block) coordinate descent (CD) methods and/or their variants, especially in the large scale optimization problems. Roughly speaking, these methods are based on the strategy of selecting one coordinate/block of variables at each iteration using some index selection procedure (e.g., cyclic, greedy, randomized). This often dramatically reduces the computational complexity of the algorithms per iteration as well as memory storage, making these methods simple and salable. See for instance and references therein and a short summary in Table I, as well as the recent comprehensive review paper for the up-to-date materials.

A widely used assumption in showing the convergence of CD methods in the literature is that the (partial) gradient of $f$ is globally Lipschitz-continuous. However, this could be a restrictive assumption violated in diverse applications in practice, such as matrix factorization , tensor decomposition , matrix/tensor completion , Poisson likelihood models , etc. Although this assumption may be relaxed by adopting conventional line search methods, the efficiency and computational complexity of the first-order method are unavoidably distorted, especially when the size of the problem is large. In fact, this longstanding issue also appears in the classical proximal gradient descent (PGD) method. Fortunately, this issue is solved in . They develop a new framework called Bregman proximal gradient (BPG) method that adapts the geometry of $f$ by the Bregman distance. In such a way, the decrease of the objective value can be still quantified. As a result, they are able to characterize the convergence behavior of BPG for minimizing convex composite problems without assuming globally Lipschitz-continuous gradient of the objective function. Further, this framework has been extended to the case of nonconvex optimization in .

Despite the crucial issue is solved in PGD-type methods, there are only few results on CD-type methods. A cyclic Bregman coordinate descent (CBCD) method has been proposed in , but no rates are given. In , the authors provide the convergence rate result using randomized (block) coordinate selection strategy in a special case where $F$ is smooth convex and $r\equiv 0$. To the best of our knowledge, how to deal with this crucial issue is still an open problem, when using CD methods to solve a nonsmooth and convex/nonconvex Problem (1). Furthermore, the accelerated version of the RBCD method has not been proposed yet, and its iteration complexity analysis is still open as well. In this paper, we bridge these gaps by proposing a randomized Bregman (block) coordinate descent (RBCD) method and its accelerated variant. The comprehensive convergence analyses are established. The main contributions are highlighted as follows.

We propose a randomized Bregman (block) coordinate descent (RBCD) method to solve the composite problem where the smooth part does not have the global Lipschitz-continuous (partial) gradient property.

By adapting the relative smoothness framework, we establish a rigorous convergence rate analysis of the RBCD method, showing that the convergence rate to an stationary point is $\mathcal{O}(n\varepsilon^{-2})$ if $F$ is nonconvex, where $k$ is the number of iterations.

If $F$ is convex, RBCD achieves the global sublinear convergence rate of $\mathcal{O}(n\varepsilon^{-1})$. The global linear convergence rate is obtained if $f$ is (relative) strongly convex.

The RBCD method can also be accelerated in the relative smoothness setting. The iteration complexity of $\mathcal{O}(n\varepsilon^{-1/\gamma})$ can be obtained through the notion of generalized translation variant (explained in the latter section) of the Bregman distance.

II Preliminaries

Notation. Throughout this paper, we use bold upper case letters denote matrices (e.g.. ${\mathbf{X}}$), bold lower case letters denote vectors (e.g., ${\mathbf{x}}$), and Calligraphic letters (e.g., $\mathcal{X}$) are used to denote sets. We use $\|\cdot\|$ to denote the Euclidean norm. $\delta_{\mathcal{X}}({\mathbf{x}})$ represents the indicator function: $\delta_{\mathcal{X}}({\mathbf{x}})=0$ if ${\mathbf{x}}\in\mathcal{X}$; otherwise, $\delta_{\mathcal{X}}({\mathbf{x}})=\infty$. If \mathcal{X}={\mbox{\mathbf{R}}}^{N}_{+}, the indicator function becomes $\delta_{+}({\mathbf{x}})$. For a function $f$, $\nabla f({\mathbf{x}})$ denotes its the gradient, while $\nabla_{i}f({\mathbf{x}})$ is the partial gradient with respect to the $i$-th block. Let $f_{i}({\mathbf{x}}_{i})$ be the function with respect to the $i$-th block, while the rest of blocks are fixed. Clearly, we have $\nabla_{i}f({\mathbf{x}})=\nabla f_{i}({\mathbf{x}}_{i})$. If $f$ is not differentiable, $\partial f$ denotes the subdifferential of $f$.

Given a convex function $\phi$, the Bregman proximal mapping of $\phi$ at a point ${\mathbf{x}}$ is defined as

where $D_{h}({\mathbf{u}},{\mathbf{x}})=h({\mathbf{u}})-h({\mathbf{x}})-\langle\nabla h({\mathbf{x}}),{\mathbf{u}}-{\mathbf{x}}\rangle$ is the Bregman distance with the reference convex function $h$. This mapping is well-defined since the functions $\phi$ and $h$ are convex. The convexity of $h$ also implies $D_{h}({\mathbf{x}},{\mathbf{y}})\geq 0,\forall{\mathbf{x}},{\mathbf{y}}$. If, in addition, $h$ is strictly convex, $D_{h}({\mathbf{x}},{\mathbf{y}})=0$ if and only if ${\mathbf{x}}={\mathbf{y}}$. In the rest of this paper, we assume $h$ is strictly convex. Note that $D_{h}({\mathbf{x}},{\mathbf{y}})$ is not symmetric in general. Therefore, we use symmetric coefficient $\theta(h)$, defined by

to measure the symmetry. When $\phi=\delta_{\mathcal{{\mathbf{x}}}}$, the Bregman proximal mapping reduces to the Bregman projection

Problem Formulation. Our goal is to solve the following composite optimization problem

where the following assumptions are made throughout this paper.

$r$ is convex, block separable, proper and loser semi-continuous.

$F^{*}=\inf_{\mathbf{x}}F({\mathbf{x}})>-\infty$.

An estimate ${\mathbf{x}}$ is said to be a stationary point of $F$ if it satisfies

Note that the objective function $F$ could be convex or nonconvex since we don’t make the convexity assumption of $f$, which is the case in . In addition, the function $r$ could be an indicator function of a closed convex set, so that the problem formulation in (6) includes the case where minimizing a nonsmooth objective function over a closed convex set.

III Randomized Bregman Coordinate Descent

In this section, we introduce the randomized Bregman (block) coordinate descent (RBCD) method for solving problem (6). Given the current estimate ${\mathbf{x}}$, the $i$-th block of coordinates is selected uniformly at random, then the new estimate ${\mathbf{x}}^{+}$ is updated as follows

where, for some stepsize $\alpha$, the vector $T_{i}({\mathbf{x}})$ is defined as

Note that we drop the index $i$ in $D_{h_{i}}$ to simplify the notation. The algorithm is summarized in Algorithm 1.

Here the stepsize $\alpha$ can be determined by a conventional line search method and the global convergence results can be established. However, line search methods are usually expensive since this subroutine requires to evaluate the objective function multiple times to ensure the sufficient descent in the objective value. To establish convergence results for a CD-type method with a constant stepsize, the common assumption is that $\nabla f({\mathbf{x}})$ (or $\nabla_{i}f({\mathbf{x}})$) is globally Lipschitz-continuous . However, this assumption may be restrict to some modern optimization problems. See for instances and reference therein. In the following section, we review the notion of relative smoothness introduced in . This notion allows us to establish the convergence results for RBCD method without the assumption of global Lipschitz-continuous gradient.

IV Convergence Analyses of RBCD

We start with the definition of relative smoothness , by which a new descent lemma is obtained without the assumption of the global Lipschitz-continuity of (partial) gradient.

[21, Definition 1.1] A pair of functions $(g,h)$ are said to be relatively smooth if $h$ is convex and there exists a scalar $L>0$ such that $Lh-g$ is convex.

Moreover, the relative smoothness nicely translates the Bregman distance to produce a non-Lipschitz descent lemma .

[20, Lemma 1] The pair of functions $(g,h)$ is relatively smooth if and only if for all ${\mathbf{x}}$ and ${\mathbf{y}}$, it holds that

When $h=\frac{1}{2}\|\cdot\|^{2}$, the classical descent lemma is recovered, i.e., $g({\mathbf{y}})-g({\mathbf{x}})-\langle\nabla g({\mathbf{x}}),{\mathbf{y}}-{\mathbf{x}}\rangle\leq\frac{L}{2}\|{\mathbf{y}}-{\mathbf{x}}\|^{2}.$

To use Lemma 1, we additionally make the following assumptions for the rest of this paper.

The functions $(f_{i},h_{i})$ are relatively smooth with constants $L_{i}>0,\forall i$.

With the relative smoothness between $(f_{i},h_{i})$, the following result shows the basic descent property of the proposed method.

For any ${\mathbf{x}}$, and any $i\in\{1,2,\cdots,n\}$, let ${\mathbf{x}}^{+}$ to be defined as in E.q. (8). Then we have

where $\theta_{i}=\theta(h_{i})$. In particular, with $0<\alpha<\frac{1+\theta_{i}}{L_{i}}$, a sufficient descent in the objective value of $F$ is guaranteed.

Maximizing the function $g(\alpha)=(1+\theta_{i}-L_{i}\alpha)\alpha$ with respect to $\alpha$ yields the stepsize $\alpha^{*}=\frac{1+\theta_{i}}{2L_{i}}$. Substituting the obtained stepsize into (11) yields the following result.

For any ${\mathbf{x}}$, let ${\mathbf{x}}^{+}$ to be defined as in E.q. (8). With stepsize $\alpha=\frac{1+\theta_{i}}{2L_{i}}$, we have

With the stepsize $\alpha=\frac{1+\theta_{i}}{2L_{i}}$, Corollary 1 quantifies the descent in the objective value. Therefore, the stepsize $\alpha^{k}=\frac{1+\theta_{i_{k}}}{2L_{i_{k}}}$ is an appropriate choice for Algorithm 1.

Since only one block is selected and updated per iteration, the quantity $D_{h}({\mathbf{x}}^{+},{\mathbf{x}})$ introduced in cannot be used to measure the optimality of the RBCD method. Given an estimate ${\mathbf{x}}$, we introduce the reference function $H$ and the corresponding Bregman mapping as follows:

Based on this mapping, the following result shows that the quantity $D_{H}(T({\mathbf{x}}),{\mathbf{x}})$ can be used to measure the optimality of $F$.

A vector ${\mathbf{x}}$ is a stationary point of $F$ if and only if $D_{H}(T({\mathbf{x}}),{\mathbf{x}})=0$.

Clearly, when $F$ is convex, then the current estimate ${\mathbf{x}}$ is a global minimum if $D_{H}(T({\mathbf{x}}),{\mathbf{x}})=0$.

Instead of using the classical convexity definition, we here use the relative strongly convexity introduced in , which is similar to the relative smoothness.

[21, Definition 1.2.] A function $g$ is $\mu$-strongly convex relative to $h$ if for any ${\mathbf{x}}$ and ${\mathbf{y}}$, there exists a scalar $\mu\geq 0$ such that

Note that if $\mu=0$, the classical convexity for a smooth function $g$ is recovered. Moreover, when $h=\frac{1}{n}\|\cdot\|$, the classical strongly convexity is recovered. In the rest of this subsection, we assume $f$ is strongly convex relative to $H$.

$f$ is $\mu$-strongly convex relative to $H$, i.e., there exists a scalar $\mu\geq 0$ such that for every ${\mathbf{y}}$ and ${\mathbf{x}}$

Since $r$ is assumed to be convex, the function $F$ is also $\mu$-strongly convex relative to $H$, i.e.,

for some ${\mathbf{v}}\in\partial F({\mathbf{x}})$. Moreover, by Assumption 2, we have

Substituting ${\mathbf{y}}=T_{i}({\mathbf{x}})$ in E.q. (17) and combing it with the inequality (19), we immediately obtain that $\mu\leq 1$.

The following lemma provides the key inequalities used to prove the convergence results of the RBCD method.

For any vector ${\mathbf{x}}$, let ${\mathbf{x}}^{+}$ to be defined as in E.q. (8) by picking up $i\in\{1,2,\cdots,n\}$ uniformly at random. Set stepsize $\alpha=\frac{1+\theta_{i}}{2L_{i}}$. For any vector ${\mathbf{u}}$, the expectation of $F({\mathbf{x}}^{+})$ satisfies

and the expectation of $D_{H}({\mathbf{x}}^{+},{\mathbf{x}})$ satisfies

By applying Lemma 4, the main convergence results are established in Theorem 1. Note that this result generalizes [2, Theorem 1] through replacing the proximal mapping by the Bregman proximal mapping so that the assumption of global Lipschitz-continues (partial) gradient is not necessary.

Let $\{{\mathbf{x}}^{k}\}$ be the sequence generated by Algorithm 1. Then for any $k\geq 0$, the iterates ${\mathbf{x}}^{k}$ satisfies

Further, if $f$ is $\mu$-strongly convex relative to $H$, then

where $\theta=\underset{i}{\min}\{\theta_{i}\}$.

Therefore, if $F$ is convex, the sequence $\{{\mathbf{x}}^{k}\}$ needs at most $\mathcal{O}(n\varepsilon^{-1})$ to converge to an $\varepsilon$-solution. Further, the classical linear convergence rate is obtained if $f$ is strongly convex (relative to $H$).

IV-B Nonconvex case

In this subsection, we establish the convergence results for the case where $F$ is nonconvex. Since $r$ is convex, $f$ is nonconvex. Due to the nonconvexity, it is of interest to find a stationary point. Lemma 3 implies that $D_{H}(T({\mathbf{x}}),{\mathbf{x}})$ can be used to measure the optimality. The following result shows the descent property of the proposed method in terms of the optimality gap $D_{H}(T({\mathbf{x}}),{\mathbf{x}})$.

For any ${\mathbf{x}}$, let ${\mathbf{x}}^{+}$ to be defined as in E.q.(8) by picking up the index $i$ uniformly at random. Let $\alpha=\frac{1+\theta_{i}}{2L_{i}}$. Then the following inequality holds:

Using Lemma 5, we can establish the convergence results of the RBCD method for nonconvex $F$.

Let $\{{\mathbf{x}}^{k}\}$ to be the sequence generated by Algorithm 1. Let stepsize $\alpha^{k}=\frac{1+\theta_{i_{k}}}{2L_{i_{k}}}$, then

The sequence $\{F({\mathbf{x}}^{k})\}$ is non-increasing.

where $F^{*}=\inf F({\mathbf{x}})>-\infty$.

Every limit point of $\{{\mathbf{x}}^{k}\}$ is a stationary point.

Suppose $H$ is $\sigma$-strongly convex with respect to the Euclidean norm $\|\cdot\|$. Then we have $D_{H}({\mathbf{y}},{\mathbf{x}})\geq\frac{\sigma}{2}\|{\mathbf{y}}-{\mathbf{x}}\|^{2}$. Combining the strongly convexity of $H$ with Theorem 2, we immediately obtain the following convergence rate result

Therefore, the sequence $\{{\mathbf{x}}^{k}\}$ converges to a stationary point at the rate of $\mathcal{O}(\frac{\sqrt{n}}{\sqrt{k}})$. In another word, to obtain an $\varepsilon$-stationary point, i.e., $\|T({\mathbf{x}})-{\mathbf{x}}\|\leq\varepsilon$, the RBCD method needs to run $\mathcal{O}(n\varepsilon^{-2})$ iterations.

V Accelerated Randomized Bregman Coordinate Descent

In this section, we restrict ourselves to the unconstrained smooth minimization problem as follows

where $f$ is convex and satisfies Assumption 1. The closed convex set ${\mathcal{X}}$ satisfies ${\mathcal{X}}={\mathcal{X}}_{1}\times\cdots\times{\mathcal{X}}_{n}$ such that ${\mathbf{x}}_{i}\in{\mathcal{X}}_{i}$ $\forall i$. It is equivalent to consider $r_{i}$ as an indicator function of the closed convex set ${\mathcal{X}}_{i}$.

The accelerated randomized Bregman coordinate descent (ARBCD) method is given as Algorithm 2. At the $k$-th iteration, the ARBCD method selects a coordinate $i_{k}$ uniformly at random, and generates the three vectors ${\mathbf{y}}^{k}$, ${\mathbf{z}}^{k+1}$, and ${\mathbf{x}}^{k+1}$, where the vectors ${\mathbf{y}}^{k}$ and ${\mathbf{x}}^{k+1}$ are the affine combinations of ${\mathbf{x}}^{k}$ and ${\mathbf{z}}^{k}$, and ${\mathbf{y}}^{k}$, ${\mathbf{z}}^{k}$, and ${\mathbf{z}}^{k+1}$, respectively, and the vector ${\mathbf{z}}^{k+1}$ is obtained as follows

Note that Step 1 and 3 of Algorithm 2 need $\mathcal{O}(N)$ operations, while $\mathcal{O}(1)$ operations are usually expected in a general coordinate descent method. In the latter section, we will show an efficient implementation of the ARBCD method so that the ARBCD method only needs $\mathcal{O}(1)$ operations at each iteration.

VI Convergence Analysis of ARBCD

which is the full-dimensional update version of ${\mathbf{z}}^{k+1}_{i_{k}}$ in E.q. (28). Therefore, the vector ${\mathbf{z}}^{k+1}$ can be computed by

It follows from the definition of ${\mathbf{x}}^{k+1}$ in Step 3 of Algorithm 2 that we have

Clearly, the vector ${\mathbf{x}}^{k+1}$ and ${\mathbf{y}}^{k}$ are only one coordinate part from each other, which satisfies the relative smoothness property in Assumption 2.

One of the challenges to establish the convergence results is from the nature of Bregman distances. Since a Bregman distance is in general not a norm, it does not hold the homogeneous translation invariant, i.e.,

To handle this issue, introduces the notion of triangle scaling property (TSP).

[22, Definition 2] The Bregman distance defined with a convex reference function $h$ has the triangle scaling property if there exists some scalar $\gamma>0$ such that for all ${\mathbf{u}},{\mathbf{v}},{\mathbf{w}}$,

In contrast, we introduce the more general notion of the generalized translation invariant (GTI) in the following definition, and show it is equivalent to triangle scaling property, when restricting $\theta\in$.

[Generalized Translation Invariant] The Bregman distance defined with a convex reference function $h$ has the generalized translation invariant property if there exists some scalar $\gamma\geq 0$ such that for all ${\mathbf{u}},{\mathbf{v}},{\mathbf{w}}$

The Bregman distance has the generalized translation invariant with $\theta\in$ if and only if it holds the triangle scaling property.

Here we gives three examples to show the existences of GNI in some Bregman divergences, while the proof is included in Appendix.

The norms. Let $\|\cdot\|_{\mathbf{A}}$ be a norm, ${\mathbf{A}}$ be a positive define matrix, $h({\mathbf{x}})=(1/2)\|{\mathbf{x}}\|^{2}_{\mathbf{A}}$, and $D_{h}({\mathbf{x}},{\mathbf{y}})=(1/2)\|{\mathbf{x}}-{\mathbf{y}}\|_{\mathbf{A}}^{2}=(1/2){\mathbf{x}}^{T}{\mathbf{A}}{\mathbf{y}}$. It is easy to see that $\gamma=2$.

The Kullback-Leibler (KL) divergence. Let $h$ be the negative Boltzmann-Shannon entropy: $h({\mathbf{x}})=\sum_{i=1}^{N}{\mathbf{x}}_{i}\log{\mathbf{x}}_{i}$ defined over {\mbox{\mathbf{R}}}_{+}^{N}. The Bregman distance is given by

The Itakura-Saito (IS) distance. Let $h$ be the Burg’s entropy: $h({\mathbf{x}})=-\sum_{i=1}^{N}\log{\mathbf{x}}_{i}$ on {\mbox{\mathbf{R}}}_{++}^{N}. The Bregman distance associated with $h$ is given by

To satisfy the definition of GNI, we must have $\gamma=0$. Similar to TSP, however, $\gamma=0$ is the uniform value for $D_{\text{IS}}$, and the intrinsic $\gamma$ value can be $2$ if the three points are close to each other [22, Theorem 1].

Note that the GTI is more general since TSP needs $\theta\in$, but GTI holds for all \theta\in{\mbox{\mathbf{R}}}.

To use the notion of GTI, we make the following assumption.

The Bregman distances $D_{h}(\cdot,\cdot)$ have the generalized translation invariant with the constant $\gamma>0$, $\forall i$.

Using the notion of GTI, we will show that the ARBCD method converges with a sublinear rate of $\mathcal{O}(n\varepsilon^{-1/\gamma})$. We start with recalling the critical lemma [33, Lemma 3.2] for a Bregman proximal mapping.

[33, Lemma 3.2] For a convex function $\phi$ and a vector $x$, if the Bregman proximal mapping is defined as

The key relationship between two consecutive iterates in Algorithm 2 is established in the following lemma.

Suppose Assumptions 1, 2, and 4 holds. For any vector ${\mathbf{u}}$, the sequences generated by Algorithm 2 satisfy, for all $k\geq 0$,

The following lemma introduces a sequence $\{\beta_{k}\}$ that satisfies the condition in Step 4 of Algorithm 2.

[22, Lemma 3] The sequence $\beta_{k}=\frac{\gamma}{k+\gamma}$ satisfies

Combing Lemma 8 with Lemma 9, the main convergence results for the ARBCD are established in the following theorem.

Suppose Assumptions 1, 2, and 4 hold. If $\beta_{k}=\frac{\gamma}{k+\gamma}$ for all $k\geq 0$, then the following inequality holds, for any vector ${\mathbf{u}}$,

Note that due to the affine combinations in Step 1 and 3 of Algorithm 2, the current implementation requires $\mathcal{O}(N)$ operations. In the next section, we introduce an efficient implementation so that only $\mathcal{O}(1)$ operations are needed at each iteration.

VII Efficient implementation

In order to avoid full-dimensional vector operations, the previous works propose a strategy that changes the variables for the accelerated coordinated descent methods in the global Lipschitz-continuous (partial) gradient setting. Here we show this scheme can be adapted so that the full-dimensional operations can be avoided in the relative smoothness setting, which is given as Algorithm 3. Instead of computing the vector ${\mathbf{z}}^{k+1}$, a search direction ${\mathbf{d}}^{k}_{i_{k}}$ is computed in Algorithm 3 as follows

The sequences $\{{\mathbf{x}}^{k},{\mathbf{y}}^{k},{\mathbf{z}}^{k}\}$ and $\{{\mathbf{u}}^{k},{\mathbf{v}}^{k}\}$ generated from Algorithm 2 and 3, respectively, satisfy

for all $k\geq 1$. That is, these two algorithms are equivalent.

Note that in Algorithm 3, only a single block coordinates of the vectors ${\mathbf{u}}^{k}$ and ${\mathbf{v}}^{k}$ are updated at each iteration, which cost $\mathcal{O}(N_{i})$ operations. Although computing the partial gradient in E.q. (42) may still cost full-dimensional operations in general, the previous works introduce a number of optimization problems where the partial gradient can be computed cheaply without actually forming ${\mathbf{y}}^{k}$.

VIII Numerical Experiments

To showcase the strength of the proposed methods, we consider two applications of relatively smooth convex optimization: Poisson inverse problem, and relative-entropy nonnegative regression.

A large number of problems in nuclear medicine, night vision, astronomy and hyperspectral imaging can be described as inverse problems where data measurements are collected according to a Poisson process whose underling intensity function is indirectly related to an object of interest through a linear system. This class of problems have been studied intensively in the literature. See for instance and references therein, as well as a more recent comprehensive review for the up-to-date references.

Formally, in a Poisson inversion problem we are given a nonnegative observation matrix {\mathbf{A}}\in{\mbox{\mathbf{R}}}_{+}^{M\times N}, a noisy measurement vector {\mathbf{b}}\in{\mbox{\mathbf{R}}}_{+}^{M}, and the goal is to recover the signal or image of interest {\mathbf{x}}\in{\mbox{\mathbf{R}}}_{+}^{N}. Under the Poisson assumption, we can rewrite the observation model as follows

Therefore, a natural and widely used measure of proximity of two nonnegative vectors is based on the KL divergence. Particularly, minimizing the KL-divergence $D_{\text{KL}}({\mathbf{b}},{\mathbf{A}}{\mathbf{x}})$ is equivalent to maximize the Poisson log-likelihood function. The optimization problem can be formulated as follows

To apply the RBCD and ARBCD methods, we need to identify a series of adequate reference functions $h_{i}$. Here we use Burg’s entropy and the corresponding Bregman distance, i.e., the IS distance.

Let $f_{i}({\mathbf{x}}_{i})=D_{\text{KL}}({\mathbf{b}},{\mathbf{A}}{\mathbf{x}})$ and $h_{i}({\mathbf{x}}_{i})$ to be defined as

Then the functions $(f_{i},h_{i})$ are relatively smooth with any scalar $L_{i}$ satisfying

Equipped with Lemma 10, Theorem 1 is applicable and warrants the convergence. Since $\theta(h_{i})=0$, we can take the stepsize $\alpha_{k}=\frac{1}{2\|{\mathbf{b}}\|_{1}}$, $\forall k\geq 0$. To solve Poisson inverse problems, the E.q. (9) can be written as

It follows from [22, Theorem 1] that the intrinsic TSE of a Bregman distance is $2$, even the uniform TSE is not. In addition, numerically shows the convergence and efficiency of the Accelerated Bregman Proximal method (ABPG) with $\gamma=2$. Thus, we here also use $\gamma=2$ for the ARBCD method. As a result, E.q. (28) becomes

We compare the proposed algorithms RBCD and ARBCD with two state-of-the-art algorithms: Bregman Proximal Gradient (BPG) method and accelerated Bregman Proximal Gradient (ABPG) method. All algorithms are implemented in Matlab code.

Figure 1 shows the computational results for a randomly generated dataset with $M=500$ and $N=500$. The entries in ${\mathbf{A}}$ and ${\mathbf{b}}$ are generated randomly from a uniform distribution over the interval $$. Each algorithm starts with the same initial values. Note that the CD-type methods has a inner loop of$N$iterations as their computational complexity is$N$ times cheaper than the gradient-based methods. As a result, the computational complexity in each iteration is identical.

In Figure 1, we can see the RBCD method is only slightly better than the BPG method, because the RBCD method uses the most updated coordinate to update, and BPG and RBCD methods use the same stepsize $\alpha_{k}=\frac{1}{2\|b\|_{1}}$. Figure 1 also shows that the accelerated methods ABPG and ARBCD are both faster than their non-accelerated variants. We can also conclude that the ARBCD method is faster than the other methods. It is well-known that the accelerated (proximal) gradient method does not guarantee the descent in the objective values at each iteration. Instead, the number of ripples are on the traces of the objective values. This criteria can be found on the ABPG method as well in Figure 1. On the other hand, we does not find such ripples or bumps from the ARBCD method. Particularly, Figure 1 shows that the ARBCD method provides consistent descent in the objective values.

It is easy to check numerically that $D_{\text{IS}}$ does not hold GTI or TSP property for any scalar $\gamma>0.5$. We conduct another experiment to explore the impact of the parameter $\gamma$. Figure 2 shows the convergence behaviors of the ABPG and ARBCD methods with $\gamma=0.1,1.0$ and $2.0$. The larger $\gamma$ is, the more acceleration the ABPG method obtains. However, it seems the ARBCD method holds the opposite relationship with the $\gamma$ values. The ARBCD method achieves the maximum acceleration when the $\gamma$ is minimum.

VIII-B Relative-entropy nonnegative regression

Anther formulation to solve the nonnegative linear inverse problem introduced in Section VIII-A is to minimize $D_{\text{KL}}({\mathbf{A}}{\mathbf{x}},{\mathbf{b}})$, i.e.,

In this case, the following result shows that the function $f$ is relative smooth to the Boltzman-Shannon entropy defined by

Let $f_{i}({\mathbf{x}}_{i})=D_{\text{KL}}({\mathbf{A}}{\mathbf{x}},{\mathbf{b}})$ and $h_{i}({\mathbf{x}}_{i})$ to be defined as

Then the functions $(f_{i},h_{i})$ are relatively smooth with any scalar $L_{i}$ satisfying

where ${\mathbf{a}}_{ij}$ is the $(i,j)$-th entry of ${\mathbf{A}}$.

Figure 3-4 shows the computational results for a randomly generated dataset with $M=500$ and $N=500$. Figure 3 shows the almost identical convergence behaviors as in Figure 1, where the RBCD and ARBCD methods are slightly faster than the BPG and ABPG methods, respectively, and the ARBCD method is faster than the rest methods. As the $\gamma$ values increases, Figure 4 shows improved convergence for the ABPG method. However, the smallest value of $\gamma$, i.e., $\gamma=0.1$, causes the divergence of the ARBCD method. Therefore, the choice of the hyperparameter $\gamma$ has significant influence on the performance of the ARBCD method.

IX Conclusion

In this paper, we propose a randomized Bregman (block) coordinate descent (RBCD) method and its accelerated variant ARBCD method for minimizing a composite problems, where the smooth part of the objective function does not satisfies the global Lipschitz-continuous (partial) gradient property. By using the relative smoothness, we establish the iteration complexity of $\mathcal{O}(n\varepsilon^{-2})$ to obtain an $\varepsilon$-stationary point in the case where $F$ is nonconvex. Besides, the iteration complexity is improved to $\mathcal{O}(n\varepsilon^{-1})$ if $f$ is convex, and the global linear convergence rate can be achieved by RBCD if $f$ is strongly convex. We introduce the notion of generalized translation invariant. Thanks to this notion, we are able to establish the convergence result for the ARBCD method which uses the acceleration technique. Thus, the iteration complexity is further improved to $\mathcal{O}(n\varepsilon^{-1/\gamma})$ by the ARBCD method.

Appendix A Appendix

From the optimality of $T_{i}({\mathbf{x}})$ in (9), we have

for some ${\mathbf{v}}_{i}^{+}\in\partial r_{i}(T_{i}({\mathbf{x}}))$. The convexity of $r_{i}$ implies

where the second inequality is due to $D_{h}({\mathbf{x}}_{i},T_{i}({\mathbf{x}}))\geq\theta D_{h}(T_{i}({\mathbf{x}}),{\mathbf{x}}_{i})$. Since ${\mathbf{x}}_{j}^{+}={\mathbf{x}}_{j}$ $\forall i\neq j$, we obtain

A-B Proof of Lemma 3

$(\Longrightarrow)$. Suppose ${\mathbf{x}}$ is a stationary point. Then we have

for some ${\mathbf{v}}\in\partial r({\mathbf{x}})$. From the convexity of $r$, it follows that for any vector ${\mathbf{u}}$

for some ${\mathbf{v}}^{+}\in\partial r(T({\mathbf{x}}))$. It follows that

Let ${\mathbf{u}}=T({\mathbf{x}})$ and combine the equations (58) and (60). Then we obtain

Since $D_{H}({\mathbf{x}},T({\mathbf{x}})),D_{H}(T({\mathbf{x}}),{\mathbf{x}})\geq 0$, we obtain $D_{H}(T({\mathbf{x}}),{\mathbf{x}})=0$.

$(\Longleftarrow)$. Suppose $D_{H}(T({\mathbf{x}}),{\mathbf{x}})=0$. The (strict) convexity of $H$ implies $T({\mathbf{x}})={\mathbf{x}}$. From (59), we obtain

which indicates ${\mathbf{x}}$ is a stationary point. ∎

A-C Proof of Lemma 4

Since each block $i$ is selected uniformly at random, we have

where $(i)$ follows from the relative smoothness of $(f_{i},h_{i})$; $(ii)$ uses the fact of $T_{i}({\mathbf{x}})=T({\mathbf{x}})_{i}$; $(iii)$ is based on the convexity of $f$ and $r$; $(iv)$ uses the the fact of $\langle\nabla h({\mathbf{z}})-\nabla h({\mathbf{x}}),{\mathbf{y}}-{\mathbf{z}}\rangle=D_{h}({\mathbf{y}},{\mathbf{x}})-D_{h}({\mathbf{y}},{\mathbf{z}})-D_{h}({\mathbf{z}},{\mathbf{x}})$.

Taking the expectation of Eq.(61) with respect to $i$ yields

A-D Proof of Theorem 1

Combining (21) with (20), let ${\mathbf{u}}={\mathbf{x}}^{*}$, and we have

Taking the expectation of (63) with respect to $\{i_{0},i_{1},\cdots\}$ yields

where the last inequality is because $\{F({\mathbf{x}}^{l})\}$ is a descent sequence. Subtracting $F({\mathbf{x}}^{*})$ on both sides and rearrange yields

Dividing both sides by $\frac{n+k}{n}$ yields the desired result.

If $f$ is $\mu$-strongly convex relative to $H$, we have

Subtracting $F({\mathbf{x}}^{*})$ on the both sides and rearrange yields

The relative strongly convexity of $F$ implies

Clearly, we have $\beta\leq 1$ since $\mu\leq 1$. Then

Combining the inequality above with (64) yields

Taking the expectation with respect to $\{i_{0},i_{1},\cdots\}$ on the both sides of the relation above, we have

Dropping $D_{H}(x^{*},{\mathbf{x}}^{k})$ on the left hand yields the desired result. ∎

A-E Proof of Lemma 5

Taking the expectation of (12) with respect to $i$ yields

where $(i)$ is because $T_{i}({\mathbf{x}})=T({\mathbf{x}})_{i}$. ∎

A-F Proof of Theorem 2

$(i)$. The result is directly obtained from Lemma 5.

$(ii)$. Taking the expectation of (24) with respect to all variables and rearranging yields

Taking the telescopic sum of the above inequality for $l=0,1,\cdots,k$ gives us

Since $F$ is lower bounded, taking the limit $k\rightarrow\infty$ yields the desired result.

$(iii)$. The inequality (66) further implies that

Dividing $k+1$ on both sides gives us the desired result. $(iv)$. Let ${\mathbf{x}}^{*}$ to be a limit point of $\{{\mathbf{x}}^{k}\}$ and there exists a subsequence $\{{\mathbf{x}}^{k_{p}}\}$ such that ${\mathbf{x}}^{k_{p}}\rightarrow{\mathbf{x}}^{*}$ as $p\rightarrow\infty$.

Since the functions $r_{i}$ are lower semi-continuous, we have for all $i$,

At the $k$-th iteration, suppose the index $i$ is selected, then the convexity of $r_{i}$ implies that

Let $\{{\mathbf{x}}^{k_{q}}\}$ be the subsequence of $\{{\mathbf{x}}^{k_{p}}\}$ such that the index $i$ is selected. Choosing $k=k_{q}-1$ in the above inequality, and letting $q\rightarrow$ yields

where we use the facts ${\mathbf{x}}^{k_{q}}\rightarrow{\mathbf{x}}^{*}$ as $q\rightarrow\infty$. Thus, combining (68) with (67), we have

Since $i$ is selected arbitrarily, we have

Furthermore, by the continuity of $f$, we obtain

From $(ii)$ and Lemma 3, it follows that ${\mathbf{x}}^{*}$ is a stationary point of $F$. ∎

A-G Proof of Lemma 6

$\Longrightarrow$. Suppose the Bregman distance $D_{h}(\cdot,\cdot)$ holds the generalized translation variant, and let ${\mathbf{u}}=(1-\theta){\mathbf{x}}+\theta{\mathbf{w}}$ for any ${\mathbf{x}}$. Then we have

Since the above inequality holds for all $\theta$, it must hold for $\theta\in$.

$\Longleftarrow$. Suppose the triangle scaling property holds. Let ${\mathbf{y}}=(1-\theta){\mathbf{u}}+\theta{\mathbf{w}}$, then we have

Therefore, the generalized translation invariant holds for $\theta\in$. ∎

A-H Proof of Remark 2

Without loss the generality, we assume $N=1$. Using the log sum inequality, we obtain

Without loss generality, we assume $N=1$. As the GNI property in Definition 4 is defined for all ${\mathbf{u}},{\mathbf{v}},{\mathbf{w}}$, we consider a special case of ${\mathbf{u}}=\theta{\mathbf{w}}$. Then, we have

To obtain $D_{\text{IS}}(\theta{\mathbf{v}},\theta{\mathbf{w}})\leq\left|\theta\right|^{\gamma}D_{\text{IS}}({\mathbf{v}},{\mathbf{w}})$ for all \theta\in{\mbox{\mathbf{R}}}, we must have $\gamma=0$, otherwise $1>\theta^{\gamma}$ for all $\theta\in(0,1)$.

A-I Proof of Lemma 8

Based on the relation in E.q. (31), we know ${\mathbf{x}}^{k+1}$ and ${\mathbf{y}}^{k}$ satisfy the relative smoothness property since they are only one coordinate difference from each other. Therefore, we obtain

where $(i)$ is using the generalized translation invariant, $(ii)$ is due to E.q. (70), and $(iii)$ is due to E.q. (30). Taking the expectation with respect to $i_{k}$ on both sides yields for all ${\mathbf{u}}$

Dividing $\beta_{k}^{\gamma}$ on both sides, we have

Multiplying both sides by $n^{\gamma}$, we obtain

Finally applying the condition in Step 4 of Algorithm 2 yields the desired result. ∎

A-J Proof of Theorem 3

Taking the expectation with respect to $\{i_{0},i_{1},\cdots,\}$ yields

The direct consequence of E.q. (74) is, for any ${\mathbf{u}}$,

Using $D_{H}({\mathbf{u}},{\mathbf{z}}^{k+1})\geq 0$, and the initialization $\beta_{0}=1$ and ${\mathbf{z}}^{0}={\mathbf{x}}^{0}$, we obtain

A-K Proof of Proposition 1

It is straightforward to see that ${\mathbf{x}}^{0}={\mathbf{y}}^{0}={\mathbf{z}}^{0}={\mathbf{v}}^{0}$. Suppose the recursive hypotheses hold for the $k$-th iteration. From the optimality of E.q. (42), we have

where $(i)$ is due to the optimality, and $(ii)$ is due to the recursive hypotheses. Similarly, from the optimality of E.q. (28), we obtain

where $(i)$ is due to the optimality, and $(ii)$ is due to the recursive hypotheses. Combing (75) and (76) yields

where $(i)$ is due to E.q. (77) and $(ii)$ is due to the recursive hypotheses.

where $(i)$ and $(iii)$ is due to recursive hypotheses, and $(ii)$ is due to Step 4 of Algorithm 3. ∎

A-L Proof of Lemma 10

where ${\mathbf{a}}_{i}$ is the $i$-th row of ${\mathbf{A}}$. The first- and second-order derivatives of $f_{j}$ are given by

It follows from the nonnegativity of ${\mathbf{A}}$ and ${\mathbf{x}}$ that we have

A-M Proof of Lemma 11

Fixing the $j$-th coordinate of ${\mathbf{x}}$, define $f_{j}({\mathbf{x}}_{j})$ as follows

Then the first- and second-derivatives of $f_{j}$ are given by

Using the nonnegativity of ${\mathbf{A}}$ and ${\mathbf{x}}$, we obtain ${\mathbf{a}}_{ij}{\mathbf{x}}_{j}\leq\langle{\mathbf{a}}_{i},{\mathbf{x}}\rangle$, which further implies

Invoking the inequality above, we obtain the desired result

References