Accelerated Bregman Proximal Gradient Methods for Relatively Smooth Convex Optimization

convex optimization, relative smoothness, Bregman divergence, proximal gradient methods, accelerated gradient methods.

Introduction

where $L_{k}>0$ for all $k\geq 0$. Here, we use the gradient $\nabla f(x_{k})$ to construct a local quadratic approximation of $f$ around $x_{k}$ while leaving $\Psi$ untouched. Our assumption that $C$ and $\Psi$ are simple means that the minimization problem in (2) can be solved efficiently, especially if it admits a closed-form solution.

This smoothness assumption implies (see, e.g., [25, Lemma 1.2.3])

See, e.g., , and [5, Chapter 10]. Under the same assumption, accelerated proximal gradient methods () can achieve a faster $O(k^{-2})$ convergence rate:

which is optimal (up to a constant factor) for this class of convex optimization problems .

While the uniform smoothness condition (3) is central in the development and analysis of first-order methods, there are many applications where the objective function does not have this property, despite being convex and differentiable. For example, in D-optimal experiment design (e.g., ) and Poisson inverse problems (e.g., ), the objective functions involve the logarithm in the form of log-determinant or relative entropy, whose gradients may blow up towards the boundary of the feasible region. In order to develop efficient first-order algorithms for solving such problems, the notion of relative smoothness was introduced by several recent works .

The function $f$ is called $L$-smooth relative to $h$ on $C$ if there is an $L>0$ such that

As shown in and , this notion of relative smoothness is equivalent to the following statements:

The definition of relative smoothness in (7) gives an upper approximation of $f$ that is similar to (4). In fact, (4) is a special case of (7) with $h=(1/2)\|x\|^{2}$ and $D_{h}(x,y)=(1/2)\|x-y\|^{2}$. Therefore it is natural to consider a more general algorithm by replacing the squared Euclidean distance in (2) with a Bregman distance:

Here, our assumption that $C$ and $\Psi$ are simple means that the minimization problem in (8) can be solved efficiently. Similar to the proximal gradient method (2), this algorithm can also be interpreted through operator splitting mechanism: it is the composition of a Bregman proximal step and a Bregman gradient step (see details in [3, Section 3.1]). Therefore, it is called the Bregman proximal gradient (BPG) method .

This is a generalization of (5). The same convergence rate for the general case (with nontrivial $\Psi$) is obtained in , where the authors also discussed the effect of a symmetry measure for the Bregman distance. Similar results are also obtained in and . In addition, introduced the notion of relative strong convexity and obtained linear convergence of the BPG method when both relative smoothness and relative strong convexity hold. More recently, studied stochastic gradient descent and randomized coordinate descent methods in the relatively smooth setting, and extended this framework to minimize relatively continuous convex functions.

A natural question is whether the $O(k^{-1})$ rate can be improved with first-order methods under the relative smoothness assumption, especially whether the accelerated $O(k^{-2})$ rate can be achieved . Very recently, it is shown by Dragomir et al. that the $O(k^{-1})$ rate is optimal for the class of relatively smooth functions, thus cannot be improved in general. However, we note that the class of relatively smooth functions is very broad, containing differentiable functions whose gradients has arbitrarily large Lipschitz constants. Indeed, the worst-case function constructed in to prove the lower bound is obtained by smoothing a nonsmooth function, which demonstrate pathological nonsmooth behavior. This is in sharp contrast to the situation under the uniform Lipschitz assumption, which uses a fixed quadratic function as the relatively smooth measure.

Ideally, it would be most informative to derive both upper and lower bounds on the convergence rate of first-order methods for every fixed function $h$ in the relatively smooth setting, or at least for the popular ones that are frequently encountered in application (such as the KL divergence). It is plausible that the achievable convergence rates for particular functions $h$ (more likely particular combinations of $f$ and $h$) can be better than $O(k^{-1})$ in theory or at least in practice. A full spectrum investigation is beyond the scope of this paper. Instead, we study a structural property of general Bregman divergences called triangle scaling and develop adaptive first-order methods that, although without a priori guarantee, often demonstrate the $O(k^{-2})$ convergence rate empirically in many applications. Moreover, these methods produce simple numerical certificates of the fast rates whenever they happen.

2 Contributions and outline

In Section 3, we propose a basic accelerated Bregman proximal gradient (ABPG) method that attains an $O(k^{-\gamma})$ convergence rate, where $\gamma\leq 2$ is the TSE of the Bregman divergence. More specifically, under the assumption (7), the basic ABPG method produces a sequence $\{x_{k}\}$ satisfying

The exact value of $\gamma$ depends on a triangle scaling property of the Bregman distance. For $D_{h}(x,y)=(1/2)\|x-y\|^{2}$, we have $\gamma=2$ and $L=L_{f}$, hence the result in (9) recovers that in (6). We also give an adaptive variant that can automatically search for the largest possible $\gamma$ for which the convergence rate in (9) holds for finite $k$ even though $\gamma$ is larger than the TSE.

In Section 5, we present an accelerated Bregman dual-averaging algorithm that has similar convergence rates as the basic ABPG method, but omit discussions of its adaptive variants.

Finally, in Section 6, we present numerical experiments with three applications: the D-optimal experiment design problem, a Poisson linear inverse problem, and relative-entropy nonnegative regression. In all experiments, the ABPG methods, especially the adaptive variants, demonstrate superior performance compared with the BPG method. Moreover, we obtain numerical certificates for the empirical $O(k^{-2})$ rate in all our experiments.

The relative smoothness condition directly extends the upper approximation property (4) with more general Bregman distances. Nesterov took an alternative approach by extending the Lipschitz condition (3). Specifically, he considered functions with Hölder continuous gradients with a parameter $\nu\in$:

and obtained $O(k^{-(1+\nu)/2})$ rate with a universal gradient method and $O(k^{-(1+3\nu)/2})$ rate with accelerated schemes. These methods are called “universal” because they do not assume the knowledge of $\nu$ and automatically ensure the best possible rate of convergence. The accelerated $O(k^{-(1+3\nu)/2})$ rate interpolates between $O(k^{-1/2})$ and $O(k^{-2})$ with $\nu\in$. There seems to be no simple connection or correspondence between the Hölder smoothness property and the combination of relative smoothness and the triangle scaling property studied in this paper.

Gutman and Peña studied iteration complexity of first-order methods using a general framework of perturbed Fenchel duality. Their framework provides alternative derivations of the convergence rates of Bregman proximal gradient methods under the relative smooth setting and the ones under Hölder continuity assumption.

Development and analysis of optimization methods in the relatively smooth setting require some delicate assumptions in order to cover many interesting applications without loss of rigor. Here we adopt the same assumptions made in regarding problem (1).

$\inf_{x\in C}\{f(x)+\Psi(x)\}>-\infty$, i.e., problem (1) is bounded below.

Sufficient conditions for the well-posedness of (8) are given in [3, Lemma 2]. The same conditions also ensure that our proposed accelerated methods are well-posed.

Triangle scaling of Bregman distance

In this section, we define the triangle scaling property for Bregman distances and discuss two different notions of triangle scaling exponent (TSE).

We call $\gamma$ a uniform triangle scaling exponent (TSE) of $D_{h}$.

Figure 1 gives a geometric illustration of the points involved in the above definition.

If $D_{h}(x,y)$ is jointly convex in $(x,y)$, then the inequality (10) holds with $\gamma=1$ because

The squared Euclidean distance. Let $h(x)=(1/2)\|x\|_{2}^{2}$ and $D_{h}(x,y)=(1/2)\|x-y\|_{2}^{2}$. Obviously, here $D_{h}$ is jointly convex in its two arguments. But it is also easy to see that

Therefore the squared Euclidean distance has a uniform TSE $\gamma=2$, which is much larger than $1$ obtained by following the jointly convex argument.

Bregman distance induced by strongly convex and smooth functions. If $h$ is $\mu$-strongly convex and $L$-smooth over its domain, then the inequality (10) would hold with $\gamma=2$ if the right-hand side is multiplied by an additional factor $G=L/\mu$, which is the condition number of $h$. We will prove this fact in Section 2.2.

Bregman distance based on polynomial kernels. Reference functions of the form $h(x)=(1/p)\|x\|^{p}$ for some $p\geq 2$ recently attracted lots of attention following Nesterov’s work on tensor methods in convex optimization . In general, the global TSEs for the induced Bregman divergence can be less than $1$ for $p>2$. However, the modified reference function $h(x)=(1/2)\|x\|^{2}+(1/p)\|x\|^{p}$ for $p\geq 4$ has $\gamma>1$, or $\gamma=2$ with an additional factor on the right-hand side of (10), over a bounded domain. We will give detailed analysis for the case $p=4$ in Section 2.2, after introducing a relaxed version of TSE.

We observe that the largest uniform TSEs are quite different for the Bregman distances listed above. An important question is: Are these differences essential such that they lead to different convergence rates if different Bregman distances are used in an accelerated algorithm? It would be ideal to derive an intrinsic characterization that is common for most Bregman distances and essential for convergence analysis of accelerated algorithms.

Plugging the last equality into (15) and after some simple algebra, we arrive at (14). ∎

2 Bounding the triangle-scaling gain

where $w=x+t(y-x)$ for some $t\in$, which we denote as $w\in[x,y]$. Consequently, if we define

Next we consider the polynomial reference function $h(x)=(1/4)\|x\|^{4}$, which does not have bounded Hessian. In this case, we have $\nabla h(x)=\|x\|^{2}x$ and $\nabla^{2}h(x)=\|x\|^{2}\cdot I+2xx^{T}$, where $I$ is the identity matrix. Clearly $\|x\|^{2}\cdot I\preceq\nabla^{2}h(x)\preceq 3\|x\|^{2}\cdot I$. According to (16), we have

As a simple fix, we consider $h(x)=(1/2)\|x\|^{2}+(1/4)\|x\|^{4}$, whose Hessian is $\nabla^{2}h(x)=(1+\|x\|^{2})I+2xx^{T}$ and it satisfies $(1+\|x\|^{2})I\preceq\nabla^{2}h(x)\preceq(1+3\|x\|^{2})I$. Therefore, according to (16),

Accelerated Bregman proximal gradient method

In this section, we present an accelerated Bregman proximal gradient (ABPG) method for solving problem (1), and analyze its convergence rate under the uniform triangle-scaling property. Adaptive variants based on the intrinsic TSE are developed in Section 4.

To simplify notation, we define a lower approximation of $F(x)=f(x)+\Psi(x)$ by linearizing $f$ at a given point $y$:

If $f$ is $L$-smooth relative to $h$ (Definition 1), then we have both a lower and an upper approximation:

When $\gamma=2$ and $\Psi\equiv 0$, Algorithm 1 reduces to the IGA (improved interior gradient algorithm) method in , which is an extension of Nesterov’s accelerated gradient method in to the Bregman proximal setting. It was shown in that the IGA method attains $O(k^{-2})$ rate of convergence under the uniform Lipschitz condition (3). In this paper, we consider the general case $\gamma\in$ under the much weaker relatively smooth condition.

We show that the ABPG method converges with a sublinear rate of $O(k^{-\gamma})$. First, we state a basic property of optimization with Bregman distance [13, Lemma 3.2].

and $h$ is differentiable at $z_{+}$, then

The following lemma establishes a relationship between the two consecutive steps of Algorithm 1. It is an extension of Proposition 1 in , which uses $\gamma=2$ under the assumption (3).

First, using the upper approximation in (17) and line 1 and line 1 in Algorithm 1, we have

where in the last inequality we used the lower bound in (17). Subtracting $F(x)$ from both sides of the inequality above, we obtain

Dividing both sides by $\theta_{k}^{\gamma}$ and rearranging terms yield

Finally applying the condition (18) gives the desired result. ∎

The sequence $\theta_{k}=\frac{\gamma}{k+\gamma}$ for $k=0,1,2,\ldots$ satisfies the condition (18).

With $\theta_{k}=\frac{\gamma}{k+\gamma}$, we have

Recall the weighted arithmetic mean and geometric mean inequality (see, e.g., [18, Section 2.5].), i.e., for any positive real numbers $a$, $b$, $\alpha$ and $\beta$, it holds that

Setting $a=k+1$, $b=k+1+\gamma$, $\alpha=1$ and $\beta=\gamma-1$, we arrive at

which, together with (24) and (25), implies the inequality (18). ∎

A slightly faster converging sequence $\theta_{k}$ can be obtained by solving the equality in (18). Since there is no closed-form solution in general, we can find $\theta_{k+1}$ as the root of

numerically, say, using Newton’s method with $\theta_{k}$ as the starting point.

Let $\theta_{0}=1$ and $\theta_{k+1}$ be the solution to (27) for all $k\geq 0$. Then $\theta_{k}\leq\frac{\gamma}{k+\gamma}$ for all $k\geq 0$.

Let $\vartheta_{k}=\frac{\gamma}{k+\gamma}$ and define another sequence $\xi_{k}$ such that $\xi_{0}=1$ and

is monotone decreasing in $\theta$. Since $\omega(\vartheta_{k+1})\leq 1/\vartheta_{k}^{\gamma}$ by Lemma 3 and $\omega(\xi_{k+1})=1/\vartheta_{k}^{\gamma}$ by (28), we have $\xi_{k+1}\leq\vartheta_{k+1}$ for all $k\geq 0$.

Next we prove $\theta_{k}\leq\vartheta_{k}$ for all $k\geq 0$ by mathematical induction. This obviously holds for $k=0$ since $\theta_{0}=\vartheta_{0}=1$. Suppose $\theta_{k}\leq\vartheta_{k}$ holds for some $k\geq 0$. Then using the facts $\omega(\theta_{k+1})=1/\theta_{k}^{\gamma}$ and $\omega(\xi_{k+1})=1/\vartheta_{k}^{\gamma}$, we obtain $\omega(\theta_{k+1})\geq\omega(\xi_{k+1})$. Since $\omega$ is monotone decreasing, we conclude that $\theta_{k+1}\leq\xi_{k+1}$. Combining with $\xi_{k+1}\leq\vartheta_{k+1}$ obtained above, we have $\theta_{k+1}\leq\vartheta_{k+1}$. This completes the induction. ∎

Using $D_{h}(x,z_{k+1})\geq 0$ and the initializations $\theta_{0}=1$ and $z_{0}=x_{0}$, we obtain

It remains to apply the condition $\theta_{k}\leq\frac{\gamma}{k+\gamma}$. ∎

2 ABPG method with exponent adaptation

The best convergence rate of the ABPG method is obtained with the largest uniform TSE for the Bregman distance. Since it is often hard to determine the largest TSE, we present in Algorithm 2 a variant of the ABPG method with automatic exponent adaptation, called the ABPG-e method.

This method starts with a large $\gamma_{0}\geq 2$. During each iteration $k$, it reduces $\gamma_{k}$ by a small amount $\delta>0$ until some stopping criterion is satisfied. An obvious choice for the stopping criterion is the local triangle-scaling property

where $x_{k+1}=(1-\theta_{k})x_{k}+\theta_{k}z_{k+1}$ and $y_{k}=(1-\theta_{k})x_{k}+\theta_{k}z_{k}$. According to the proof of Lemma 2, we can also use the inequality (21) as stopping criterion, which is implied by (29) and the relatively smooth assumption. For convergence analysis, we only need (21) to hold, which can be less conservative than (29). In Algorithm 2, we use the following inequality as the stopping criterion

which is equivalent to (21) (by subtracting $\Psi(x_{k+1})$ from both sides of the inequality). In practice, this condition often leads to much faster convergence than using (29). Computationally, it is slightly more expensive since it needs to evaluate $f(x_{k+1})$ in addition to $\nabla f(y_{k})$ during each inner loop, while (29) does not.

By replacing inequality (18) with the one above and repeating the analysis in Section 3.1, we obtain the following result.

ABPG methods with gain adaptation

In this section, we present and analyze an adaptive ABPG method based on the concept of intrinsic TSE developed in Section 2.1. Instead of searching for the largest uniform TSE as in Algorithm 2, we can replace line 1 in Algorithm 1 by

Algorithm 3 is such a method with gain adaptation. During each iteration, the algorithm uses an inner loop to search for an value of $G_{k}$ that satisfies

which is true if the following local triangle-scaling property holds:

For any $\alpha,\beta>0$ and $\gamma\geq 1$, the following inequality holds:

The case of $\gamma=1$ is obvious. Assume $\gamma>1$. The desired inequality is equivalent to

Applying the weighted arithmetic and geometric mean inequality (26), we have

where $\overline{G}_{k}$ is a weighted geometric mean of the gains at each step:

We follow the same steps as in Section 3.1. In light of (31), the inequality (21) becomes

Then the same arguments in the proof of Theorem 1 lead to

Next we derive an upper bound for $G_{k}\theta_{k}^{\gamma}$. For convenience, let’s define for $k=0,1,2,\ldots$,

Then (32) implies $a_{k+1}=A_{k+1}-A_{k}$. Moreover, we have

Applying Lemma 5 with $\alpha=A_{k+1}^{1/\gamma}$ and $\beta=A_{k}^{1/\gamma}$, we obtain

We can eliminate the common factor $A_{k+1}$ on both sides of the above inequality to obtain

Summing the above inequality from step $0$ to $k-1$ and using $A_{0}=1/G_{0}$, we have

Using the weighted arithmetic and geometric mean inequality (e.g., [18, Section 2.5]) gives

Combining the last two inequalities above, we arrive at

Finally, substituting the inequality above into (38) gives the desired result. ∎

We note that the geometric mean $\overline{G}_{k}$ in (34) can be much smaller than the average (arithmetic mean) of $\{G_{0},G_{1},\ldots,G_{k}\}$. Under the assumption of uniform Lipschitz smoothness (3), Nesterov proposed an accelerated gradient method with non-monotone line search. However, the complexity obtained there still depends on the global Lipschitz constant $L$, more specifically, replacing $\overline{G}_{k}L$ in (33) with $\rho L$ when $\gamma=2$. Our result in (33) can be tighter if the local Lipschitz constants are smaller than $L$ (equivalently with $G_{k}<1$).

In order to estimate the overhead of the gain-adaptation procedure, we follow the approach of [27, Lemma 4]. Notice that each inner loop needs to call a gradient oracle to compute $\nabla f(y_{k})$, and also $f(x_{k+1})$ when we use (30) as the stopping criterion for gain adaptation. Let $n_{i}\geq 1$ be the number of calls of the oracle (for $\nabla f(y_{k})$) at the $i$th iteration, for $i=0,\ldots,k$. Then

Therefore, the total number of oracle calls is

Roughly speaking, on average each iteration need two oracle calls (unless $G_{k}$ becomes very large).

As an alternative to calculating $\theta_{k+1}$ by solving the equation (32), we can also use the following explicit update rule:

1 Towards the O⁡(k−2)O(k^{-2}) convergence rate

where $u_{k}\in\bigl[(1-\theta_{k})x_{k}+\theta_{k}z_{k},(1-\theta_{k})x_{k}+\theta_{k}z_{k+1}\bigr]$ and $v_{k}\in\bigl[z_{k},z_{k+1}\bigr]$. Suppose the sequence $\{x_{k}\}$ converges to the optimal solution $x_{\star}$ and $\theta_{k}\to 0$, then have $u_{k}\to x_{\star}$. If $x_{\star}$ is an interior point of the positive orthant or the simplex, meaning $x_{\star}^{(i)}>0$ for all coordinates $i$, Then we see from (40) that the bound on $G_{\theta_{k}}(x_{k},z_{k},z_{k+1})$ depends on how close $x_{\star}$ is close to the boundary (assuming $v_{k}$ is bounded).

The most interesting case is when the optimal solution $x_{\star}$ is on the boundary, i.e., when $x_{\star}^{(i)}=0$ for some coordinates $i$. In fact, in our numerical examples on the D-optimal design problem and Poisson linear inverse problem, most of the solutions are on the boundary. However, we emphasize that the iterates generated by the ABPG algorithm is never on the boundary, but may only converge to the boundary; see Assumption A, especially A.5. Our analysis applies to this case as well. In particular, we can show that if both sequences $\{x_{k}\}$ and $\{z_{k}\}$ converge, then $x_{k}^{(i)}\to 0$ implies $z_{k}^{(i)}\to 0$ and the convergence rate of $x_{k}^{(i)}$ is no faster than that of $z_{k}^{(i)}$. (Here $x_{k}^{(i)}\to 0$ means $\lim_{k\to\infty}x_{k}^{(i)}=0$.) More precisely, we have the following lemma.

Suppose an algorithm generates two sequences $\{x_{k}\}$ and $\{z_{k}\}$ in the strictly positive orthant, satisfying $x_{0}=z_{0}$ and $x_{k+1}=(1-\theta_{k})x_{k}+\theta_{k}z_{k+1}$ for all $k\geq 0$. Then

If $\{x_{k}\}$ converges and $x_{k}^{(i)}\to 0$ for some coordinate $i$, then there must exists an subsequence of $\{z_{k}^{(i)}\}$ that converges to $0$.

Suppose both sequences $\{x_{k}\}$ and $\{z_{k}\}$ converge. If $x_{k}^{(i)}\to 0$ for some $i$, then it converges at a rate that is no faster than $z_{k}^{(i)}$ in the following sense: For any monotone decreasing sequence $\{r_{k}\}$ that converges to $0$ and satisfies $z_{k}^{(i)}\geq r_{k}$ for all $k\geq 0$, we have $x_{k}^{(i)}\geq r_{k}$ for all $k\geq 0$. In particular, we can choose $r_{k}$ to be the monotone lower envelop of $z_{k}^{(i)}$, i.e., $r_{k}=\min\{z_{0}^{(i)},z_{1}^{(i)},\ldots,z_{k}^{(i)}\}$.

By the update rule $x_{k+1}=(1-\theta_{k})+\theta_{k}z_{k+1}$, we know that each $x_{k}$ is a convex combination of the points $\{z_{0}=x_{0},z_{1},\ldots,z_{k}\}$, which all lie in the strictly positive orthant. Therefore,

Part (a). Suppose $x_{k}^{(i)}\to 0$ but there is no subsequence of $\{z_{k}^{(i)}\}$ converging to $0$. Then there must exist an $\epsilon>0$ such that $z_{k}^{(i)}>\epsilon$ for all $k\geq 0$. Since $x_{k}$ is a convex combination of $\{z_{0},z_{1},\ldots,z_{k}\}$, this implies

which contradicts with the assumption that $x_{k}^{(i)}\to 0$. Therefore, there must exists an subsequence of $\{z_{k}^{(i)}\}$ that converges to $0$.

Part (b). Suppose $\{r_{k}\}$ is monotone decreasing and converges to $0$. If $z_{k}^{(i)}\geq r_{k}$ for all $k\geq 0$, then

In this sense, $x_{k}^{(i)}\to 0$ at a rate that is no faster than $z_{k}^{(i)}$. ∎

According to the bound in (40) and Lemma 6, if both sequences $\{x_{k}\}$ and $\{z_{k}\}$ converge, then

which can be bounded by a constant, since $x_{k}^{(i)}\to 0$ at a rate that is no faster than $z_{k}^{(i)}$. In this case, we have $G_{k}$ in Algorithm 3 bounded by a constant asymptotically, thus the convergence rate is $O(k^{2})$.

However, we are not able to prove the convergence of the sequences $\{x_{k}\}$ and $\{z_{k}\}$ without additional assumptions (such as relative strong convexity). Indeed, to the best of our knowledge, convergence of these sequences have not been established even under the classical uniform Lipschitz condition. Therefore, an a priori theoretical guarantee of the $O(k^{-2})$ rate seems to be out of reach in general, which seems to coroborate the recent result in that the $O(k^{-1})$ rate cannot be improved in in general for the class of relatively smooth functions.

Nevertheless, we would like to reiterate the remarks at the end of Sections 1.1. In particular, the class of relatively smooth functions is very large, and the lower bound in is established with a worst-case function with pathological nonsmooth behavior. In practical applications, we always work with one particular reference function which may possess structural properties that allow fast convergence. In Algorithm 3, the sequence $\{G_{k}\}$ is readily available as part of the computation and we can easily check the magnitude of $\overline{G}_{k}$. Whenever it is small, we obtain a numerical certificate that the algorithm did converge with the $O(k^{-2})$ rate. This is exactly what we observe in the numerical experiments in Section 6.

Accelerated Bregman dual averaging method

In this section, we present an accelerated Bregman dual averaging (ABDA) method under the relative smoothness assumption. This method extends Nesterov’s accelerated dual averaging method ( and [33, Algorithm 3]) to the relatively smooth setting. Here we focus on a simple variant in Algorithm 4 based on the uniform triangle-scaling property, although it is also possible to develop more sophisticated variants with automatic exponent or gain adaptation.

In other words, $\psi_{k+1}$ is a weighted sum of the lower approximations in (17) constructed at $y_{0},\ldots,y_{k}$:

When implementing Algorithm 4, we only need to keep track of $g_{k}$ and $\vartheta_{k}$, and there is no need to maintain the abstract form of $\psi_{k}(x)$. Here our assumption of $C$ and $\Psi$ being simple means that the minimization problem in (43) can be solved efficiently. This requirement is equivalent to that for the BPG method (8) and all variants of the ABPG methods in this paper.

To see this, we use induction. Clearly it holds for $k=0$ if we choose $\theta_{0}=1$. Suppose it holds for some $k\geq 0$, then in light of (45) and (44),

Therefore the inequality (45) holds for all $k\geq 0$.

To analyze the convergence of Algorithm 4, we need the following simple variant of Lemma 1.

Notice that for $k\geq 1$, $z_{k}$ is the minimizer of $\psi_{k}(z)+Lh(z)$ over $C$. We use Lemma 7 to obtain

Combining the inequalities (47) and (48), we obtain

where in the last equality we used recursive definition of $\psi_{k+1}$ in (41). Dividing both sides of the above inequality by $\theta_{k}^{\gamma}$, we have

Using (44) and rearranging terms gives the desired result (46), which holds for $k\geq 1$. ∎

Suppose Assumption A holds, $f$ is $L$-smooth relative to $h$ on $C$, and $\gamma$ is a uniform TSE of $D_{h}$. The sequences generated by Algorithm 4 satisfy:

In this case, we can extend the result of Lemma 8 to hold for all $k\geq 0$. Applying the inequality (46) for iterations $0,1,\ldots,k$, we obtain

where we used $\theta_{0}=1$ and $\psi_{0}\equiv 0$. Next using (44) and rearranging terms, we have

According to Lemma 4, we have $\theta_{k}\leq\frac{\gamma}{k+\gamma}$ if (44) holds, which gives (49).

where the last inequality repeats the arguments from (51) to (52). Rearranging terms leads to

and further applying Lemma 4 gives the desired result (50). ∎

where the second inequality used the upper bound in (17), and the last inequality used the lower bound in (17). Therefore, for any $x$ such that $F(x)<F(z_{1})+LD_{h}(x,z_{1})$, we have

Numerical experiments

We consider three applications of relatively smooth convex optimization: D-optimal experiment design, Poisson linear inverse problem, and relative-entropy nonnegative regression. For each application, we compare the algorithms developed in this paper with the BPG method (8) and demonstrate significant performance improvement. Our implementations and experiments are shared through an open-source repository at https://github.com/linxiaolx/accbpg.

Figure 2(b) shows that for $\gamma=1.0$ and $1.5$, $\widehat{G}_{k}$ is mostly much smaller than $1$. For $\gamma=2$, $\widehat{G}_{k}$ is much closer to $1$ but always less than $1$. This gives a numerical certificate that the ABPG method converged with $O(k^{-2})$ rate. For $\gamma=2.2$, $\widehat{G}_{k}$ stayed close to $1$ for the first 700 iterations and then jumped to $3$ and stayed around. The method diverges with larger value of $\gamma$. We didn’t plot the ABDA method (Algorithm 4) because it overlaps with ABPG for the same value of $\gamma$ when the initial point is taken as the center of the simplex, see part (a) of Theorem 4.

We also show the comparison in terms of CPU time in Figures 3(c) and 3(d). As remarked at the end of Section 3.2, the ABPG-e method only take a constant number more iterations than ABPG, thus its their comparison is very similar to the case with number of iterations. For ABPG-g, the analysis in Section 4 on page 4 shows that the number of gradient calls and proximal computations is roughly twice of the ABPG method with the same number of iterations. This is exactly what we observe in Figures 3(c) and 3(d). Given such predictable scaling between number of iterations and CPU time, we only show comparisons in the number of iterations in the rest numerical experiments.

Figure 4 shows the comparison of different methods on another random problem instance with $m=80$ and $n=120$. All methods converge much faster and reach very high precision. In particular, BPG and BPG-LS look to have linear convergence. This indicates that this problem instance is much better conditioned and the objective function may be strongly convex relative to Burg’s entropy. In this case, it is shown in that the BPG method attains linear convergence. The ABPG and ABPG-g methods demonstrate periodic non-monotone behavior. A well-known technique to avoid such oscillations and attain fast linear convergence is to restart the algorithm whenever the function value starts to increase . We applied restart (RS) to both ABPG and ABPG-g, which resulted in a much faster convergence as shown in Figure 4.

1.2 Experiment on real data

In our second experiment, we construct D-optimal design instances from LibSVM data . In particular, we consider several regression datasets – the goal is to find the most relevant data points where one shall run the experiment to evaluate the corresponding label.

Figure 5 shows the results on four different datasets: abalone ($n=4177,m=8$), bodyfat ($n=252,m=14$), mpg ($n=392,m=7$) and housing ($n=506,m=13$). The left column indicates that in each case, the best performance of ABPG is achieved with large TSE $\gamma=2$ and $\gamma=2.2$. Furthermore, ABPG with $\gamma>1$ always compared favorably over plain BPG.

Next, the second column of Figure 5 shows that both ABPG-g and ABPG (with $\gamma=2$) always significantly outperform BPG and BPG-LS. We have chosen log-log scale of the plot to contrast the $O(k^{-1})$ convergence rate of BPG (with line search) with the $O(k^{-2})$ convergence rate of ABPG and ABPG-g. In the third column, we plot the local triangle-scaling gains. They serve as numerical certificates of the empirical $O(k^{-2})$ convergence rate of ABPG and its variants. In particular, we see that $G_{k}$ for the ABPG-g algorithm is mostly flat and less than one.

2 Poisson linear inverse problem

Figure 6 shows our computational results for a randomly generated instance with $m=200$ and $n=100$ and $\Psi\equiv 0$ (no regularization). The entries of $A$ and $b$ are generated following independent uniform distribution over the interval $$.

Figure 7 shows the results for a randomly generated instance with $m=100$ and $n=1000$. In this case, since $m<n$, we added a regularization $\Psi(x)=(\lambda/2)\|x\|^{2}$ with $\lambda=0.001$. ABPG-g has the best performance. Again we observe that $G_{k}\ll 1$ most of the time, which gives a numerical certificate that the ABPG methods do converge with $O(k^{-2})$ rate.

3 Relative-entropy nonnegative regression

Acknowledgments

We thank Haihao Lu, Robert Freund and Yurii Nesterov for helpful conversations. Peter Richtárik acknowledges the support of the KAUST Baseline Research Funding Scheme.

References