A Generalized Accelerated Composite Gradient Method: Uniting Nesterov's Fast Gradient Method and FISTA

TL;DR

This paper introduces a unified accelerated gradient method that combines the strengths of FGM and FISTA, providing improved robustness, flexibility, and convergence analysis for large-scale convex optimization problems.

Contribution

The authors develop a generalized accelerated composite gradient method uniting FGM and FISTA, with a Lyapunov-based analysis and enhanced monotonicity and line-search features.

Findings

The new method outperforms FGM and FISTA in robustness and usability.

Monotonicity stabilizes convergence in composite optimization.

Simulation results validate the method's effectiveness on benchmark problems.

Abstract

The most popular first-order accelerated black-box methods for solving large-scale convex optimization problems are the Fast Gradient Method (FGM) and the Fast Iterative Shrinkage Thresholding Algorithm (FISTA). FGM requires that the objective be finite and differentiable with known gradient Lipschitz constant. FISTA is applicable to the more broad class of composite objectives and is equipped with a line-search procedure for estimating the Lipschitz constant. Nonetheless, FISTA cannot increase the step size and is unable to take advantage of strong convexity. FGM and FISTA are very similar in form. Despite this, they appear to have vastly differing convergence analyses. In this work we generalize the previously introduced augmented estimate sequence framework as well as the related notion of the gap sequence. We showcase the flexibility of our tools by constructing a Generalized Accelerated Composite Gradient Method, that unites FGM and FISTA, along with their most popular variants. The Lyapunov property of the generalized gap sequence used in deriving our method implies that both FGM and FISTA are amenable to a Lyapunov analysis, common among optimization algorithms. We further showcase the flexibility of our tools by endowing our method with monotonically decreasing objective function values alongside a versatile line-search procedure. By simultaneously incorporating the strengths of FGM and FISTA, our method is able to surpass both in terms of robustness and usability. We support our findings with simulation results on an extensive benchmark of composite problems. Our experiments show that monotonicity has a stabilizing effect on convergence and challenge the notion present in the literature that for strongly convex objectives, accelerated proximal schemes can be reduced to fixed momentum methods.

Figures

Tables

Table 1. Table 1: Algorithm stall time in WTU based on the outcome of LSSC and MC

	MC passed	MC failed
LSSC passed	$0$	$\max \{t_f,t_{\Psi }\}$
LSSC failed	$t_f+t_g+t_p$	N / A

Table 2. Table 2: Per-iteration cost of FISTA, AMGS, and generalized ACGM in the non-monotone setting

	FISTA	AMGS	ACGM
Base cost	$t_g+t_p$	$2t_g+2t_p$	$t_g+t_p$
LSSC $t_d$	$t_f$	$t_g$	$t_f$
LSSC $t_c$	$t_p$	$t_g+t_p$	$t_g+t_p$
Backtrack cost	$t_f+t_p$	$2t_g+t_p$	$t_f+t_g+t_p$

Table 3. Table 3: Oracle functions of the five test problems

	$f\left(𝒙\right)$	$\Psi \left(𝒙\right)$	$\nabla f\left(𝒙\right)$	$pro{x}_{\tau \Psi }\left(𝒙\right)$
LASSO	$\frac{1}{2}{\Vert 𝑨𝒙-𝒃\Vert }_2^2$	${\lambda }_1{\Vert 𝒙\Vert }_1$	${𝑨}^T\left(𝑨𝒙-𝒃\right)$	${𝒯}_{\tau {\lambda }_1}\left(𝒙\right)$
NNLS	$\frac{1}{2}{\Vert 𝑨𝒙-𝒃\Vert }_2^2$	${\sigma }_{{ℝ}_{+}^n}\left(𝒙\right)$	${𝑨}^T\left(𝑨𝒙-𝒃\right)$	${\left(𝒙\right)}_{+}$
L1LR	$ℐ\left(𝑨𝒙\right)-{𝒚}^T𝑨𝒙$	${\lambda }_1{\Vert 𝒙\Vert }_1$	${𝑨}^T\left(ℒ\left(𝑨𝒙\right)-𝒚\right)$	${𝒯}_{\tau {\lambda }_1}\left(𝒙\right)$
RR	$\frac{1}{2}{\Vert 𝑨𝒙-𝒃\Vert }_2^2$	$\frac{{\lambda }_2}{2}{\Vert 𝒙\Vert }_2^2$	${𝑨}^T\left(𝑨𝒙-𝒃\right)$	$\frac{1}{1+\tau {\lambda }_2}𝒙$
EN	$\frac{1}{2}{\Vert 𝑨𝒙-𝒃\Vert }_2^2$	${\lambda }_1{\Vert 𝒙\Vert }_1+\frac{{\lambda }_2}{2}{\Vert 𝒙\Vert }_2^2$	${𝑨}^T\left(𝑨𝒙-𝒃\right)$	$\frac{1}{1+\tau {\lambda }_2}{𝒯}_{\tau {\lambda }_1}\left(𝒙\right)$

Table 4. Table 6: FGM and FISTA, along with their common variants, can be considered instances of generalized ACGM with certain restrictions applied.

Algorithm	Restriction
	Smooth objective	$\mu =0$	$\mu >0$	$A_0=0$	$A_0>0$	Fixed step size	Non-monotone
FGM [4]	yes	no	no	no	yes	yes	yes
FGM [22]	yes	yes	no	unclear	unclear	no	yes
FISTA [11]	no	yes	no	yes	no	partial	yes
FISTA [23, 24]	no	yes	no	yes	no	no	yes
MFISTA [26]	no	yes	no	yes	no	yes	no
scAPG [15]	no	no	yes	no	yes	no	yes
FISTA-CP [7]	no	no	no	no	no	yes	no

Equations

$$\underset{\mathbf{x}\in {\mathbb{R}}^n}{\min }F(\mathbf{x})\stackrel{\mathrm{def}}{=}f(\mathbf{x})+\Psi (\mathbf{x}),$$

$$\text{mbox}{prox}_{\tau \Psi }(\mathbf{x})\stackrel{\mathrm{def}}{=}\underset{\mathbf{z}\in {\mathbb{R}}^n}{\mathrm{argmin}}\left(\Psi (\mathbf{z})+\frac{1}{2\tau }\Vert \mathbf{z}-\mathbf{x}{\Vert }_2^2\right),$$

$$\psi (\mathbf{x})\stackrel{\mathrm{def}}{=}{\psi }^{\ast }+\frac{\gamma }{2}\Vert \mathbf{x}-\mathbf{v}{\Vert }_2^2,\mathbf{x}\in {\mathbb{R}}^n,$$

$$Q_{\lambda ,\mathbf{y}}(\mathbf{x})\stackrel{\mathrm{def}}{=}f(\mathbf{y})+⟨\nabla f(\mathbf{y}),\mathbf{x}-\mathbf{y}⟩+\frac{\lambda }{2}\Vert \mathbf{x}-\mathbf{y}{\Vert }_2^2.$$

$$T_L(\mathbf{y})\stackrel{\mathrm{def}}{=}\underset{\mathbf{x}\in {\mathbb{R}}^n}{\mathrm{argmin}}\left(Q_{L,\mathbf{y}}(\mathbf{x})+\Psi (\mathbf{x})\right),\mathbf{y}\in {\mathbb{R}}^n,$$

$$T_L(\mathbf{y})=\text{mbox}{prox}_{\frac{1}{L}\Psi }\left(\mathbf{y}-\frac{1}{L}\nabla f(\mathbf{y})\right),\mathbf{y}\in {\mathbb{R}}^n,L>0.$$

$$A_k(F({\mathbf{x}}_k)-F({\mathbf{x}}^{\ast }))\le A_0(F({\mathbf{x}}_0)-F({\mathbf{x}}^{\ast }))+\frac{{\gamma }_0}{2}\Vert {\mathbf{x}}_0-{\mathbf{x}}^{\ast }{\Vert }_2^2.$$

$$A_{k}F({\mathbf{x}}_k)\le H_k,$$

$$H_k\stackrel{\mathrm{def}}{=}(A_k-A_0)F({\mathbf{x}}^{\ast })+A_{0}F({\mathbf{x}}_0)+\frac{{\gamma }_0}{2}\Vert {\mathbf{x}}_0-{\mathbf{x}}^{\ast }{\Vert }_2^2$$

$$F(\mathbf{x})\ge W_k(\mathbf{x})\text{mbox},\mathbf{x}\in {\mathbb{R}}^n,k\ge 1.$$

$${\mathcal{H}}_k\stackrel{\mathrm{def}}{=}(A_k-A_0)W_k({\mathbf{x}}^{\ast })+A_{0}F({\mathbf{x}}_0)+\frac{{\gamma }_0}{2}\Vert {\mathbf{x}}^{\ast }-{\mathbf{x}}_0{\Vert }_2^2,$$

$${\psi }_k(\mathbf{x})\stackrel{\mathrm{def}}{=}(A_k-A_0)W_k(\mathbf{x})+A_{0}F({\mathbf{x}}_0)+\frac{{\gamma }_0}{2}\Vert \mathbf{x}-{\mathbf{x}}_0{\Vert }_2^2,$$

$${\psi }_k^{\ast }\stackrel{\mathrm{def}}{=}\underset{\mathbf{x}\in {\mathbb{R}}^n}{\min }{\psi }_k(\mathbf{x}),k\ge 0,$$

$${\psi }_k^{\ast }=\underset{\mathbf{x}\in {\mathbb{R}}^n}{\min }{\psi }_k(\mathbf{x})\le {\psi }_k({\mathbf{x}}^{\ast })={\mathcal{H}}_k,k\ge 0.$$

$$A_{k}F({\mathbf{x}}_k)\le {\psi }_k^{\ast },k\ge 0,$$

$$A_{k}F({\mathbf{x}}_k)\stackrel{\text{eqref}labelaas}{\le }{\psi }_k^{\ast }\stackrel{\text{eqref}labelaar}{\le }{\psi }_k({\mathbf{x}}^{\ast })={\mathcal{H}}_k\stackrel{\text{eqref}labelaao}{\le }{H}_k,k\ge 0.$$

$${\psi }_k^{\prime }(\mathbf{x})\stackrel{\mathrm{def}}{=}{\psi }_k(\mathbf{x})+H_k-{\mathcal{H}}_k,k\ge 0,$$

$${\psi }_k^{\prime }(\mathbf{x})={\psi }_k(\mathbf{x})+(A_k-A_0)(F({\mathbf{x}}^{\ast })-W_k({\mathbf{x}}^{\ast })),k\ge 0.$$

$$A_{k}F({\mathbf{x}}_k)\le {\psi }_k^{\prime \ast },k\ge 0,$$

$${\psi }_k^{\prime \ast }\stackrel{\mathrm{def}}{=}\underset{\mathbf{x}\in {\mathbb{R}}^n}{\min }{\psi }_k^{\prime }(\mathbf{x}),k\ge 0.$$

$$f({\mathbf{z}}_{k+1})\le Q_{L_{k+1},{\mathbf{y}}_{k+1}}({\mathbf{z}}_{k+1}),k\ge 0,$$

$${\mathbf{z}}_{k+1}\stackrel{\mathrm{def}}{=}{T}_{L_{k+1}}({\mathbf{y}}_{k+1}),$$

$$\begin{array}{c}F(\mathbf{x})\ge w_{k+1}(\mathbf{x})\stackrel{\mathrm{def}}{=}F({\mathbf{z}}_{k+1})+\frac{L_{k+1}+{\mu }_{\Psi }}{2}\Vert {\mathbf{z}}_{k+1}-{\mathbf{y}}_{k+1}{\Vert }_2^2\\ +(L_{k+1}+{\mu }_{\Psi })⟨{\mathbf{y}}_{k+1}-{\mathbf{z}}_{k+1},\mathbf{x}-{\mathbf{y}}_{k+1}⟩+\frac{\mu }{2}\Vert \mathbf{x}-{\mathbf{y}}_{k+1}{\Vert }_2^2.\end{array}$$

$${\mathbf{z}}_{k+1}=\underset{\mathbf{x}\in {\mathbb{R}}^n}{\mathrm{argmin}}\left(f({\mathbf{y}}_{k+1})+⟨\nabla f({\mathbf{y}}_{k+1}),\mathbf{x}-{\mathbf{y}}_{k+1}⟩+\frac{L_{k+1}}{2}\Vert \mathbf{x}-{\mathbf{y}}_{k+1}{\Vert }_2^2+\Psi (\mathbf{x})\right).$$

$$\mathbf{\xi }=L_{k+1}({\mathbf{y}}_{k+1}-{\mathbf{z}}_{k+1})-\nabla f({\mathbf{y}}_{k+1}),$$

$$f(\mathbf{x})$$

$$\Psi (\mathbf{x})$$

$$F(\mathbf{x})\ge$$

$$+\frac{{\mu }_{\Psi }}{2}\Vert \mathbf{x}-{\mathbf{z}}_{k+1}{\Vert }_2^2-\frac{L_{k+1}}{2}\Vert {\mathbf{y}}_{k+1}-{\mathbf{z}}_{k+1}{\Vert }_2^2.$$

$$W_{k+1}(\mathbf{x})={\left(\sum _{i=1}^{k+1}{a}_i\right)}^{-1}\left(\sum _{i=1}^{k+1}{a}_i{w}_i(\mathbf{x})\right),k\ge 0,$$

$$A_{k+1}=A_k+a_{k+1},k\ge 0.$$

Full text

\headers

Generalized Accelerated Composite Gradient MethodMihai I. Florea and Sergiy A. Vorobyov

{AMS}

90C06, 68Q25, 90C25

1 Introduction

Numerous large-scale convex optimization problems have recently emerged in a variety of fields, including signal and image processing, statistical inference, computer vision, and machine learning. Often, little is known about the actual structure of the objective function. Therefore, optimization algorithms used in solving such problems can only rely (e.g., by means of callback functions) on specific black-box methods, called oracle functions [1]. The term “large-scale” refers to the tractability of certain computational primitives (see also [2]). In the black-box setting, it means that the oracle functions of large-scale problems usually only include scalar functions and operations that resemble first-order derivatives.

Many large-scale applications were rendered practical to address with the advent of Nesterov’s Fast Gradient Method (FGM) [3]. The canonical form of FGM, formulated in [4], requires that the objective be differentiable with Lipschitz gradient. Many optimization problems, particularly inverse problems in fields such as sparse signal processing, linear algebra, matrix and tensor completion, and digital imaging (see [5, 6, 7, 8, 9] and references therein), have a composite structure. In these composite problems, the objective $F$ is the sum of a function $f$ with Lipschitz gradient (Lipschitz constant $L_{f}$) and a simple but possibly non-differentiable regularizer $\Psi$. The regularizer $\Psi$ embeds constraints by being infinite outside the feasible set. Often, $L_{f}$ is not known in advance. The composite problem oracle functions are the scalar $f(\bm{x})$ and $\Psi(\bm{x})$, as well as the gradient $\nabla f(\bm{x})$ and the proximal operator $prox_{\tau\Psi}(\bm{x})$.

To address composite problems, Nesterov has devised an Accelerated Multistep Gradient Scheme (AMGS) [10]. This method updates a Lipschitz constant estimate (LCE) at every iteration using a subprocess commonly referred to in the literature as “line-search” [8, 11]. However, the generation of a new iterate (advancement phase of an iteration) and line-search are interdependent and cannot be executed in parallel. Moreover, AMGS utilizes only the gradient-type oracle functions $\nabla f(\bm{x})$ and $prox_{\tau\Psi}(\bm{x})$. In many applications, including compressed sensing (e.g., LASSO [12]) and many classification tasks (e.g., $l_{1}$-regularized logistic regression), the evaluation of $\nabla f(\bm{x})$ is more computationally expensive than $f(\bm{x})$. An alternative to AMGS that uses $f(\bm{x})$ calls in line-search has been proposed by Beck and Teboulle in the form of the Fast Iterative Shrinkage-Thresholding Algorithm (FISTA) [11].111FISTA has very similar functionality and convergence guarantees to those of the first variant of FGM introduced in [3]. For clarity, we denote the latter method also as FISTA, whereas we use FGM to designate the canonical form developed in [4]. FISTA also benefits from having line-search decoupled from advancement. However, FISTA is unable to decrease the LCE at run-time and does not take into account strong convexity.

A strongly convex extension of FISTA, which we designate as FISTA Chambolle-Pock (FISTA-CP), has been recently introduced in [7], but without line-search. Two works [13, 14] propose a similar method that accommodates only feasible start, assumes $\mu_{\Psi}=0$, and does not feature line-search nor monotonically decreasing objective function values (monotonicity). In the context of the composite problem class, this method is merely a restricted version of FISTA-CP. Consequently, we also treat it as FISTA-CP. Another study [15] proposes a variant of FISTA-CP equipped with a line-search procedure capable of decreasing the LCE (a fully adaptive line-search strategy). However, this strongly convex Accelerated Proximal Gradient (scAPG) method is guaranteed to converge only if $f$ is strongly convex and if the starting point is feasible.

Apart from the above methods, which are applicable to a wide range of problems and feature convergence guarantees, the literature abounds with algorithms that either employ heuristics or are highly specialized. For instance, restart techniques have been developed to address strong convexity. However, the adaptive restart proposed in [16] lacks rigorous convergence guarantees for non-quadratic objectives and the periodic restart introduced in [3] requires prior knowledge of the Lipschitz constant. A very interesting study [17] proposes an over-relaxed variant of FISTA, a first-order method that appears to define a distinct category of algorithms. However, the study fails to argue whether the convergence guarantees are superior to those of FISTA. Despite being tested only on instances of the LASSO problem [12], the performance benefit is unclear even in this restricted scenario. Numerous other methods successfully exploit additional problem structure, such as sparsity of optimal points (e.g. [18, 19, 20, 21]), however they cannot be applied to the entire composite problem class.

FISTA (when the objective is non-strongly convex), FISTA-CP, and scAPG (when $f$ is strongly-convex) are identical in form to certain variants of FGM. Therefore, FISTA, FISTA-CP, and scAPG can be viewed as extensions of FGM for composite objectives. However, whereas FGM was derived using the estimate sequence [4], FISTA, FISTA-CP (including the variant found in [13, 14]), and scAPG were proposed without derivation. No details have been provided on how the update rules of each method were constructed, apart from their obvious validation in the corresponding convergence analysis. Consequently, new features of FGM cannot be directly incorporated into FISTA, FISTA-CP, and scAPG. For instance, Nesterov has proposed in [22] a line-search variant of FGM with an implicit derivation based on the estimate sequence, albeit only for non-strongly convex objectives. Had FISTA been derived in this way, a fully adaptive line-search variant of FISTA could simply be obtained from [22]. Instead, a sophisticated fully adaptive line-search extension was proposed in [23], with a technical derivation based on the mathematical constructs of [11]. On the other hand, through partial adoption of the estimate sequence, i.e., relating FISTA to the “constant step scheme I” variant of FGM in [4] and AMGS, we have arrived at a similar but simpler fully adaptive line-search scheme for FISTA [24], but again only in the non-strongly convex scenario.

In [25], we have introduced the augmented estimate sequence framework and used it to derive the Accelerated Composite Gradient Method (ACGM), which incorporates by design a fully adaptive line-search procedure. ACGM has the convergence guarantees of FGM, the best among primal first-order methods, while being as broadly applicable as AMGS. In addition, FISTA-CP and FISTA, along with the fully adaptive line-search extensions in [23] and [24], are particular cases of ACGM [25]. However, to accommodate infeasible start, we have imposed restrictions on the input parameters. Algorithms such as the “constant step scheme III” variant of FGM [4] and scAPG [15] are guaranteed to converge only when the starting point is feasible, and thus do not correspond to any instance of ACGM in [25].

1.1 Contributions

• •

In this paper, by lifting the restrictions on the input parameters imposed by infeasible start, we generalize the augmented estimate sequence framework and derive a generalization of ACGM that encompasses FGM and FISTA, along with their variants.

• •

We further showcase the flexibility and power of the augmented estimate sequence framework by endowing ACGM with monotonicity alongside its adaptive line-search procedure. Monotonicity is a desirable property, particularly when dealing with proximal operators that lack a closed form expression or other kinds of inexact oracles [7, 26]. Even when dealing with exact oracles, monotonicity leads to a more stable and predictable convergence rate.

• •

Our generalized ACGM, with its additional features, including the fully adaptive estimation of the Lipschitz constant, is therefore superior to FGM and FISTA in terms of flexibility and usability.

• •

We further generalize the definition of the wall-clock time unit (previously introduced in [25]), to account for arbitrary oracle function costs. We analyze the per-iteration complexity of our method in this generalized context, taking advantage of parallelization when possible.

• •

We support our theoretical findings with simulation results on an extended benchmark of signal processing and machine learning problems.

1.2 Assumptions and notation

We consider the composite optimization problems of the form

[TABLE]

where $\bm{x}$ is a vector of $n$ optimization variables (bold type indicates a vector in $\mathbb{R}^{n}$), and $F$ is the objective function. The constituents of the objective $F$ are the convex differentiable function $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ and the convex lower semicontinuous regularizer function $\Psi:\mathbb{R}^{n}\rightarrow\mathbb{R}\cup\{\infty\}$. Function $f$ has Lipschitz gradient (Lipschitz constant $L_{f}>0$) and a strong convexity parameter $\mu_{f}\geq 0$. Note that the existence of a strong convexity parameter does not imply the strong convexity property since any convex function has a strong convexity parameter of zero. Regularizer $\Psi$ has a strong convexity parameter $\mu_{\Psi}\geq 0$, entailing that objective $F$ has a strong convexity parameter $\mu=\mu_{f}+\mu_{\Psi}$. Constraints are enforced by making $\Psi$ infinite outside the feasible set, which is closed and convex.

Apart from the above properties, nothing is assumed known about functions $f$ and $\Psi$, which can only be accessed in a black-box fashion [1] by querying oracle functions $f(\bm{x})$, $\nabla f(\bm{x})$, $\Psi(\bm{x})$, and $\mbox{prox}_{\tau\Psi}(\bm{x})$, with arguments $\bm{x}\in\mathbb{R}^{n}$ and $\tau>0$. The proximal operator $\mbox{prox}_{\tau\Psi}(\bm{x})$ is given by

[TABLE]

for all $\bm{x}\in\mathbb{R}^{n}$ and $\tau>0$, with $\|.\|_{2}$ designating the standard Euclidean norm in $\mathbb{R}^{n}$.

Central to our derivation are generalized parabolae, quadratic functions whose Hessians are multiples of the identify matrix. We refer to the strongly convex ones simply as parabolae. Any parabola $\psi:\mathbb{R}^{n}\rightarrow\mathbb{R}$ can be written in canonical form as

[TABLE]

where $\gamma>0$ denotes the curvature, $\bm{v}$ is the vertex, and $\psi^{*}$ is the optimum value.

For conciseness, we introduce the generalized parabola expression $Q_{\lambda,\bm{y}}(\bm{x})$ for all $\bm{x},\bm{y}\in\mathbb{R}^{n}$ and $\lambda\geq 0$ as

[TABLE]

The proximal gradient operator $T_{L}(\bm{y})$ [10] can be expressed succinctly using (4) as

[TABLE]

where $L>0$ is a parameter corresponding to the inverse of the step size. Operator $T_{L}(\bm{y})$ can be evaluated in terms of oracle functions as

[TABLE]

2 Generalizing ACGM

In our previous formulation of the augmented estimate sequence framework and subsequent derivation of the Accelerated Composite Gradient Method (ACGM) in [25], we have chosen to accommodate infeasible start by default. In the sequel, we argue why this assumption was necessary and provide a generalized framework, wherein this assumption is only a particular case.

2.1 FGM estimate sequence

We begin the generalization of the estimate sequence by noting that a rate of convergence can only be obtained on function values and not iterates. It may be argued that, when the objective function is strongly convex, many first-order schemes, including the non-accelerated fixed-point methods, guarantee the linear convergence of the iterates to the optimal point. However, when the problem is non-strongly convex, the optimization landscape may contain a high-dimensional subspace of very low curvature in the vicinity of the set of optimal points. In this case, a convergence rate of iterates remains a difficult open problem [27]. For instance, Nesterov has provided in [4] an ill-conditioned quadratic problem where a convergence rate for iterates cannot be formulated for all first-order schemes of a certain structure.

Hence, we choose to measure convergence using the image space distance (ISD), which is the distance between the objective values at iterates and the optimal value. The decrease rate of an upper bound on the ISD gives the convergence guarantee (provable convergence rate). The generalization of the estimate sequence framework follows naturally from the formulation of such guarantees. Specifically, we interpret the image space distance upper bound (ISDUB), provided by Nesterov for FGM in [4], for all $k\geq 0$ as

[TABLE]

Point $\bm{x}^{*}$ can be any optimal point. However, we consider it fixed throughout this work. The convergence guarantee is given by the sequence $\{A_{k}\}_{k\geq 0}$ with $A_{0}\geq 0$. Each iteration may not be able to improve the objective function value, but at least it must improve the convergence guarantee, meaning that $\{A_{k}\}_{k\geq 0}$ must be an increasing sequence. From $A_{0}\geq 0$ it follows that $A_{k}>0$ for all $k\geq 1$.

The right-hand side of (7) is a weighted sum between the initial ISD and the corresponding domain space term (DST), with weights given by $A_{0}$ and $\gamma_{0}$, respectively. In the derivation of FGM, the weights are constrained as $A_{0}>0$ and $\gamma_{0}\geq A_{0}\mu$ whereas in the original ACGM [25] we have enforced $A_{0}=0$ and $\gamma_{0}\leq 1$. Given that our current aim is to provide a generic framework, we impose no restrictions on the weights, apart from $A_{0}\geq 0$ and $\gamma_{0}>0$. The former restriction follows from the convexity of $F$ while the latter is required by the estimate sequence, along with its augmented variant, as we shall demonstrate in the sequel.

The ISDUB expression can be rearranged to take the form

[TABLE]

where

[TABLE]

is the highest upper bound that can be placed on weighted objective values $A_{k}F(\bm{x}_{k})$ to satisfy (7).

The value of $H_{k}$ depends on the optimal value $F(\bm{x}^{*})$, which is an unknown quantity. The estimate sequence provides a computable, albeit more stringent, replacement for $H_{k}$. It is obtained as follows. The convexity of the objective implies the existence of a sequence $\{W_{k}\}_{k\geq 1}$ of convex global lower bounds on $F$, written as

[TABLE]

For instance, when querying $\nabla f(\bm{x})$ at a test point $\bm{y}_{k}$, we obtain a global lower bound on $f$ as $Q_{\mu_{f},\bm{y}_{k}}(\bm{x})$ for all $\bm{x}\in\mathbb{R}^{n}$. Consequently $Q_{\mu_{f},\bm{y}_{k}}+\Psi$ is a global lower bound of $F$. The lower bounds obtained from test points up to and including $\bm{y}_{k}$ can be combined to produce $W_{k}$. The drawback of incorporating $\Psi$ in the lower bounds lies in the introduction of additional $\mbox{prox}_{\tau\Psi}(\bm{x})$ calls, as in [10]. In Subsection 2.3 we will formulate more computationally efficient bounds. Notwithstanding, the framework we introduce here does not assume any structure on $W_{k}$ apart from convexity.

By substituting the optimal value terms $F(\bm{x}^{*})$ in (8) with $W_{k}(\bm{x}^{*})$, we obtain $\mathcal{H}_{k}$, a lower bound on $H_{k}$, given by

[TABLE]

for all $k\geq 0$. This still depends on $\bm{x}^{*}$. However, $\mathcal{H}_{k}$ can be viewed as the value of an estimate function, taken at an optimal point $\bm{x}^{*}$. The estimate functions $\psi_{k}$ are defined as functional extensions of $\mathcal{H}_{k}$, namely

[TABLE]

for all $\bm{x}\in\mathbb{R}^{n}$ and $k\geq 0$. Note that the first estimate function $\psi_{0}$ does not contain a lower bound term. Therefore, it is not necessary to define $W_{0}$. The collection of estimate functions $\{\psi_{k}\}_{k\geq 0}$, is referred to as the estimate sequence.

The estimate function optimum value, given by

[TABLE]

is guaranteed to be lower than $\mathcal{H}_{k}$, since

[TABLE]

As such, $\psi_{k}^{*}$ provides the sought after computable replacement of $H_{k}$. Note that if the lower bounds $W_{k}$ are linear, the estimate functions are generalized parabolae with the curvature given by $\gamma_{0}$. In this case, the existence of $\psi_{k}^{*}$ is conditioned by $\gamma_{0}>0$, explaining the assumption made in (7).

Thus, it suffices to maintain the estimate sequence property, given by

[TABLE]

to satisfy the ISDUB expression (7). The proof follows from the above definitions as

[TABLE]

2.2 Generalizing the augmented estimate sequence

The interval between the maintained upper bound $\psi_{k}^{*}$ and the highest allowable bound $H_{k}$ contains $\mathcal{H}_{k}$. This allows us to produce a relaxation of the estimate sequence by forcibly closing the gap between $\mathcal{H}_{k}$ and $H_{k}$. Namely, we define the augmented estimate functions as

[TABLE]

with $\{\psi^{\prime}_{k}(\bm{x})\}_{k\geq 0}$ being the augmented estimate sequence. We expand definition (17) as

[TABLE]

The augmented estimate sequence property is formulated as

[TABLE]

where the augmented estimate function optimum value $\psi^{\prime*}_{k}$ is given by

[TABLE]

2.3 Quadratic lower bounds

As we have seen in Subsection 2.1, the structure of the estimate sequence and, subsequently, the augmented estimate sequence, depends on the global lower bounds $W_{k}$. These can be obtained by combining simple lower bounds, each produced at each iteration $k$ using test point $\bm{y}_{k}$ and denoted as $w_{k}$. In the following, we propose simple lower bounds $w_{k}$ that can be obtained in a computationally efficient manner, each requiring only one call to $T_{L}(\bm{x})$.

Lemma 2.1**.**

If at every iteration $k\geq 0$, the auxiliary points and LCEs obey the descent condition, stated as

[TABLE]

where

[TABLE]

then the objective is lower bounded for all $\bm{x}\in\mathbb{R}^{n}$ as

[TABLE]

Proof 2.2**.**

By expanding (22) we obtain

[TABLE]

This subproblem is unconstrained and has a strongly convex and subdifferentiable objective. The first-order optimality condition is therefore equivalent to

[TABLE]

where $\bm{\xi}$ is a specific subgradient of function $\Psi$ at point $\bm{z}_{k+1}$.

Moreover, from the definition of the strong convexity parameters of $f$ and $\Psi$ we have that

[TABLE]

By adding together (21), (26), and (27), canceling and rearranging terms we obtain

[TABLE]

*Further applying (25) to (28) and rearranging terms gives the desired result. *

Once we have obtained the simple lower bounds, we need to determine a means of combining them. We propose the convex combination given by

[TABLE]

where

[TABLE]

We have assumed that $\{A_{k}\}_{k\geq 0}$ is increasing which implies that the weights are strictly positive, namely

[TABLE]

From definition (12), we have that the initial estimate function is a parabola, given by

[TABLE]

Applying (29) and (30) to (12) we get simple recursion rule in the form of

[TABLE]

Lemma 2.3**.**

When the simple lower bounds are given by Lemma 2.1 and are combined as in (29), then the estimate functions and the augmented estimate functions at every iteration are parabolic, and thus can be written in canonical form for all $\bm{x}\in\mathbb{R}^{n}$ and $k\geq 0$ as

[TABLE]

with $\psi_{k}^{*}$ and $\psi_{k}^{\prime*}$ given by (13) and (20), respectively, and

[TABLE]

Moreover, the curvatures and vertices obey the following recursion rules

[TABLE]

Proof 2.4**.**

The Hessian of $\psi_{0}(\bm{x})$ is given by $\gamma_{0}I_{n}$ (where $I_{n}$ is the identity matrix of size $n$) with $\gamma_{0}>0$ while the Hessian of $w_{k+1}(\bm{x})$ is $\mu I_{n}$. From the recursion rule in (33), it follows by simple induction that all estimate functions are parabolae, and as such, can be written as in (34) with (13) and (36). The relation in (18) means that the augmented estimate functions differ from the estimate function by a constant term. It follows that (35) with (20) holds.

Expanding the recursion rule in (33) using (34) and taking the gradient we obtain

[TABLE]

*Matching the first-order coefficients and constant terms on both sides of (39) gives (37) and (38), respectively. *

By iterating (37), we obtain a closed form expression of the estimate function curvature as

[TABLE]

2.4 Generalizing the gap sequence

We introduce the augmented estimate sequence gaps $\{\Gamma_{k}\}_{k\geq 0}$ and the gap sequence $\{\Delta_{k}\}_{k\geq 0}$, respectively, as

[TABLE]

Proposition 2.5**.**

When the estimate functions are parabolic, as in Lemma 2.3, the augmented estimate sequence gaps can be expressed more succinctly as

[TABLE]

Proof 2.6**.**

From (18), using (34) and (35) in Lemma 2.3, it follows that

[TABLE]

Therefore, the gap $\Gamma_{k}$ can be expressed as

[TABLE]

Note that (18), (32), and Lemma 2.3 imply that $\psi_{0}^{\prime*}=\psi_{0}^{\prime}(\bm{x}_{0})=\psi_{0}(\bm{x}_{0})=A_{0}F(\bm{x}_{0})$ which, together with (41), guarantees that the initial augmented estimate sequence gap $\Gamma_{0}$ is zero, regardless of the structure of the lower bounds. Therefore, a sufficient condition for the preservation of the augmented estimate sequence property (19) as the algorithm progresses is that the augmented estimate sequence gaps do not increase. Proposition 2.5 states that when the lower bounds are generalized parabolae, the sufficient condition for (19) becomes

[TABLE]

Therefore, (19) and the subsequent the convergence of the algorithm are implied by the Lyapunov (non-increasing) property of the gap sequence. Such Lyapunov functions have been widely employed in the theoretical analysis of optimization schemes (e.g., [28, 29, 30]). Moreover, the simple structure and generic nature of the gap sequence further justifies the augmentation of the estimate sequence.

2.5 A design pattern for first-order accelerated methods

The augmented estimate sequence property places an upper bound on $F(\bm{x}_{k+1})$, which must be satisfied before the generation of $\bm{x}_{k+1}$. The black-box nature of the objective constitutes an obstacle in this respect and it is therefore advantageous to employ a surrogate upper bound on the objective. Numerous algorithms employ this technique, and are often referred to as majorization minimization (MM) algorithms (see, e.g., [31, 32] and references therein). New iterates are generated as

[TABLE]

where $u_{k+1}$ is a local upper bound on the objective at the new iterate, namely

[TABLE]

Combining (30), (33), and (46) yields the algorithm design pattern outlined in Pattern 1. Note that Nesterov’s FGM and AMGS adhere to Pattern 1. This pattern will form the scaffolding of our generalized ACGM.

Pattern 1 takes as input the starting point $\bm{x}_{0}$, the curvature $\gamma_{0}$, the initial guarantee $A_{0}$ and, if the Lipschitz constant is not known in advance, an initial LCE $L_{0}>0$. In line 7 of Pattern 1, the main iterate is given by the minimum of the local upper bound $u_{k+1}$. Alongside the main iterate, the algorithm maintains an estimate function $\psi_{k+1}$, obtained from the previous one by adding a simple global lower bound $w_{k+1}$ weighted by $a_{k+1}>0$ (line 8 of Pattern 1). Function $w_{k+1}$ need not be given by Lemma 2.1, but we assume that it is uniquely determined by an auxiliary point $\bm{y}_{k+1}$. The current LCE $L_{k+1}$, weight $a_{k+1}$, and auxiliary point $\bm{y}_{k+1}$ are computed using algorithm specific methods $\mathcal{S}$, $\mathcal{F}_{a}$, and $\mathcal{F}_{y}$, respectively (lines 3, 4, and 5 of Pattern 1). These methods take as parameters the state of the algorithm, given by current values of the main iterate, LCE, weight, and estimate function.

2.6 Upper bounds

Interestingly, all of the above results do not rely on a specific form of the local upper bound $u_{k+1}$, as long as assumption (21) holds for all $k\geq 0$. The enforced descent condition (21) can be equivalently expressed in terms of composite objective values (by adding $\Psi(\bm{\bm{z}}_{k+1})$ to both sides) as

[TABLE]

We want our algorithm to converge as fast as possible while maintaining the monotonicity property, expressed as

[TABLE]

Then, without further knowledge of the objective function, (48) and (49) suggest a simple expression of the upper bound for all $\bm{x}\in\mathbb{R}^{n}$ and $k\geq 0$ in the form of

[TABLE]

where $\sigma_{X}$ is the indicator function [33] of set $X$, given by

[TABLE]

2.7 Towards an algorithm

With the lower and upper bounds fully defined, we select functions $\mathcal{S}$, $\mathcal{F}_{a}$, and $\mathcal{F}_{y}$ so as to preserve the Lyapunov property of the gap sequence (45).

Theorem 2.7**.**

If at iteration $k\geq 0$ we have

[TABLE]

then

[TABLE]

where subexpressions $\mathcal{A}_{k+1}$, $\mathcal{B}_{k+1}$, $\bm{g}_{k+1}$, $\bm{s}_{k+1}$, and $Y_{k+1}$ are, respectively, defined as

[TABLE]

Proof 2.8**.**

We assume $k\geq 0$ throughout this proof. The descent condition assumption implies that the lower bound property in Lemma 2.1 holds. We define

[TABLE]

which simplifies the following subexpression of $w_{k+1}$ in Lemma 2.1:

[TABLE]

The tightness of lower bound $w_{k+1}$ in Lemma 2.1 is given by the residual $R_{k+1}$ as

[TABLE]

From Lemma 2.1 it follows that

[TABLE]

By expanding terms in (55) using the identity in (53) we obtain

[TABLE]

where

[TABLE]

The vertex update in (38) implies that

[TABLE]

Using (58) and (59) in (57) and rearranging terms yields

[TABLE]

with $S_{k+1}$ and $V_{k+1}$ defined as

[TABLE]

By further applying (37) and (58), $S_{k+1}$ and $V_{k+1}$, respectively, become

[TABLE]

*Putting together (52), (56), (60), (63), and (64) with $F(\bm{x}_{k+1})\leq F(\bm{z}_{k+1})$ and the fact that the second term in the right-hand side of (63) is always non-negative, we obtain the desired result. More details can be found in proofs of Theorems 3, 4, and 5 in [34]. *

Theorem 2.7 provides a simple sufficient condition for the monotonicity of the gap sequence.

Proposition 2.9**.**

The monotonic decrease of the gap sequence in (45) is guaranteed regardless of the algorithmic state if, for all $k\geq 0$, the following hold:

[TABLE]

where

[TABLE]

Proof 2.10**.**

Herein, we consider all $k\geq 0$. According to Theorem 2.7, an obvious sufficient condition for (45) is $\mathcal{A}_{k+1}\geq 0$ combined with $\mathcal{B}_{k+1}\geq 0$. Due to the black-box nature of the oracle functions and the assumption of arbitrary algorithmic state, we allow for $\bm{g}_{k+1}$ to be any vector in $\mathbb{R}^{n}$. Quantity $\mathcal{B}_{k+1}$ is a scalar product between two vectors, the one containing $\bm{g}_{k+1}$ being arbitrary. Therefore, a sufficient condition for $\mathcal{B}_{k+1}\geq 0$ is that $\bm{s}_{k+1}=0$, which we can always guarantee by setting (65). Moreover, the value of $\|\bm{g}_{k+1}\|_{2}^{2}$ could be very large, but is always non-negative. Therefore, a sufficient condition for $\mathcal{A}_{k+1}\geq 0$ is that

[TABLE]

By expanding (68) using (30) and (37) and rearranging terms, we obtain that $a_{k+1}$ needs to satisfy the inequality

[TABLE]

*Because $L_{k+1}>\mu_{f}$, the solutions of (69) lie between the two roots of the corresponding equation, only one of which is positive and given by $\mathcal{E}(\gamma_{k},A_{k},L_{k+1})$. *

Proposition 2.9 allows us to select $\mathcal{F}_{a}$ and $\mathcal{F}_{y}$ as to preserve the Lyapunov property in (45). For $\mathcal{F}_{y}$, the obvious choice is (65). Out of the multitude of potential $\mathcal{F}_{a}$, we choose the one that yields the largest convergence guarantees $A_{k+1}$, given by

[TABLE]

which corresponds to equality in (66), namely

[TABLE]

For determining the LCE, we select the backtracking line-search method $\mathcal{S}_{A}$ employed by AMGS  [10] and the original ACGM [25]. The search parameters comprise the LCE increase rate $r_{u}>1$ and the LCE decrease rate $0<r_{d}<1$. The search terminates when the line-search stopping criterion (LSSC) in (21) is satisfied.

2.8 Putting it all together

We have thus determined a search strategy $\mathcal{S}_{A}$, the initial estimate function $\psi_{0}$ in (32), the upper bounds $u_{k+1}$ in (50), the lower bounds $w_{k+1}$ in Lemma 2.1, function $\mathcal{F}_{a}$ in (70), and function $\mathcal{F}_{y}$ in (65). Substituting these expressions in the design model outlined in Pattern 1, we can write down a generalization of ACGM in estimate sequence form, as listed in Algorithm 2.

Non-monotone generalized ACGM can be obtained by enforcing $\bm{x}_{k+1}=\bm{z}_{k+1}$ for all $k\geq 0$, accomplished by replacing line 16 of Algorithm 2 with

[TABLE]

3 Complexity analysis

3.1 Worst-case convergence guarantees

Algorithm 2 maintains the convergence guarantee in (7) explicitly at run-time as state variable $A_{k}$. Moreover, if sufficient knowledge of the problem is available, it is possible to formulate a worst-case convergence guarantee before running the algorithm.

For our analysis, we will need to define a number of curvature-related quantities, namely the local inverse condition number $q_{k+1}$ for all $k\geq 0$, the worst-case LCE $L_{u}$, and the worst-case inverse condition number $q_{u}$, respectively, given by

[TABLE]

The worst-case convergence rate for generalized ACGM is stated in the following theorem.

Theorem 3.1**.**

If $\gamma_{0}\geq A_{0}\mu$, the generalized ACGM algorithm generates a sequence $\{\bm{x}_{k}\}_{k\geq 1}$ that satisfies

[TABLE]

where

[TABLE]

Proof 3.2**.**

The non-negativity of the weights (31) implies that $\gamma_{k}\geq\gamma_{0}$ for all $k\geq 0$. Combined with (71) and $A_{k}\mu\geq 0$, we have for all $k\geq 0$ that

[TABLE]

As we can see from Algorithm 2, scaling $A_{0}$ and $\gamma_{0}$ by a fixed factor does not alter the behavior of generalized ACGM. Additionally, $\gamma_{0}$ is guaranteed to be non-zero. To simplify calculations, we introduce the normalized convergence guarantees $\bar{A}_{k}\operatorname{\ \overset{def}{=}\ }A_{k}/\gamma_{0}$ for all $k\geq 0$.

Regardless of the outcome of individual line-search calls, the growth of the normalized accumulated weights obeys

[TABLE]

Taking into account that $A_{0}\geq 0$, we obtain by induction that

[TABLE]

From assumption $\gamma_{0}\geq A_{0}\mu$, (40) implies $\gamma_{k}\geq A_{k}\mu$ for all $k\geq 0$. Hence

[TABLE]

Since $A_{0}\geq 0$, we have that $\bar{A}_{1}\geq\frac{1}{L_{u}-\mu_{f}}$. By induction, it follows that

[TABLE]

*Substituting (78) and (80) in (7) completes the proof. *

Theorem 3.1 shows that generalized ACGM has the best convergence guarantees among its class of algorithms (see [25] for an in-depth discussion).

Note that the assumption $\gamma_{0}\geq A_{0}\mu$ always holds for non-strongly convex objectives and that ACGM is guaranteed converge for strongly convex objectives also when $\gamma_{0}<A_{0}\mu$. However, in the latter case, it is more difficult to obtain simple lower bounds on the convergence guarantees. We leave such an endeavor to future research.

3.2 Wall-clock time units

So far, we have measured the theoretical performance of algorithms in terms of convergence guarantees (including the worst-case ones) indexed in iterations. This does not account for the complexity of individual iterations. In [25], we have introduced a new measure of complexity, the wall-clock time unit (WTU), to compare optimization algorithms more reliably. We thus distinguish between two types of convergence guarantees. One is the previously used iteration convergence guarantee, indexed in iterations and a new computational convergence guarantee, indexed in WTU.

The WTU is a measure of running time in a shared memory parallel scenario. The computing environment consists of a small number of parallel processing units (PPU). Each PPU may be a parallel machine itself. The number of parallel units is considered sufficient to compute any number of independent oracle functions simultaneously. The shared-memory system does not impose constraints on parallelization, namely, it is uniform memory access (UMA) [35] and it is large enough to store the arguments and results of oracle calls for as long as they are needed.

In order to compare algorithms based on a unified benchmark, in [25] we have assumed that $f(\bm{x})$ and $\nabla f(\bm{x})$ require $1$ WTU each while all other operations are negligible and amount to [math] WTU. In this work, we generalize the analysis. We attribute finite non-negative costs $t_{f}$, $t_{g}$, $t_{\Psi}$, and $t_{p}$ to $f(\bm{x})$, $\nabla f(\bm{x})$, $\Psi(\bm{x})$, and $\mbox{prox}_{\tau\Psi}(\bm{x})$, respectively. However, since we are dealing with large-scale problems, we maintain the assumption that element-wise vector operations, including scalar-vector multiplications, vector additions, and inner products, have negligible complexity when compared to oracle functions and assign a cost of [math] WTU to each. Synchronization of PPUs also incurs no cost. Consequently, when computed in isolation, an objective function value call costs $t_{F}=\max\{t_{f},t_{\Psi}\}$, ascribable to separability, while a proximal gradient operation costs $t_{T}=t_{g}+t_{p}$, due to computational dependencies.

3.3 Per-iteration complexity

We measure this complexity in WTU on the shared memory system described in the previous subsection and consider a parallel implementation involving speculative execution [35].

The advancement phase of a generalized ACGM iteration consists of one proximal gradient step (line 9 of Algorithm 2). Hence, every iteration has a base cost of $t_{T}=t_{g}+t_{p}$. LSSC and the monotonicity condition (MC) in line 16 of Algorithm 2 can be evaluated in parallel with subsequent iterations. Both rely on the computation of $f(\hat{\bm{z}}_{k+1})$, which in the worst case requires $\lceil t_{f}/t_{T}\rceil$ dedicated PPUs. In addition, MC may need up to $\lceil t_{\Psi}/t_{T}\rceil$ PPUs.

Backtracks stall the algorithm in a way that cannot be alleviated by parallelization or intensity reduction. Therefore, it is desirable to make them a rare event. Assuming that the local curvature of $f$ varies around a fixed value, this would mean that $\log(r_{u})$ should be significantly larger than $-log(r_{d})$. With such a parameter choice, the algorithm can proceed from one iteration to another by speculating that backtracks do not occur at all. Let the current iteration be indexed by $k$. If the LSSC of iteration $k$ fails, then the algorithm discards all the state information pertaining to all iterations made after $k$, reverts to iteration $k$, and performs the necessary computation to correct the error. We consider that a misprediction incurs a detection cost $t_{d}$ and a correction cost $t_{c}$. LSSC requires the evaluation of $f(\hat{\bm{z}}_{k+1})$ and incurs a detection cost of $t_{d}=t_{f}$. A backtrack entails recomputing $\hat{\bm{y}}_{k+1}$, yielding an LSSC $t_{c}=t_{T}$ correction time.

Overshoots are assumed to occur even less often. Similarly, the algorithm proceeds speculating that MC always passes and defaults to (72). Hence, MC has $t_{d}=t_{F}$, due to its dependency on $\Psi(\hat{\bm{z}}_{k+1})$, but once the algorithmic state of iteration $k$ has been restored, no additional oracle calls are needed, leading to $t_{c}=0$. MC and LSSC can be fused into a single condition, giving rise to the scenarios outlined in Table 1. Note that if LSSC fails, MC is not evaluated.

For non-monotone generalized ACGM, each backtrack adds $t_{f}+t_{T}$ WTU to a base iteration cost of $t_{T}$. A comparison to other methods employing line-search is shown in Table 2.

4 Extrapolated form

4.1 Monotonicity and extrapolation

In the original ACGM [25], the auxiliary point can be obtained from two successive main iterates through extrapolation. Interestingly, this property is preserved for any step size. We show in the following how monotonicity alters this property and bring generalized monotone ACGM to a form in which the auxiliary point is an extrapolation of state variables. We begin with the following observation.

Lemma 4.1**.**

In Algorithm 2, the estimate sequence vertices can be obtained from other state variables through extrapolation as

[TABLE]

Proof 4.2**.**

In this proof we consider all $k\geq 0$. The vertex update in (38) can be rewritten using (71) as

[TABLE]

The auxiliary point update in (65) can also be written as

[TABLE]

From (37) we have that

[TABLE]

Putting together (81), (82), and eliminating $\bm{y}_{k+1}$ using (83) we obtain

[TABLE]

Lemma 4.1 enables us to write down the corresponding expression for the auxiliary point.

Proposition 4.3**.**

In Algorithm 2, the auxiliary points $\bm{y}_{k+1}$ obey the following extrapolation rule:

[TABLE]

where

[TABLE]

and $\bm{1}_{X}$ denotes the membership function of set $X$, namely

[TABLE]

Proof 4.4**.**

The results herein apply to all $k\geq 1$. Lemma 4.1 applied for $\bm{v}_{k}$, combined with the auxiliary point update in (65), leads to

[TABLE]

Depending on the outcome of the update in line 16 of Algorithm 2, we distinguish two situations.

If MC passes at iteration $k-1$ ($F(\bm{z}_{k})\leq F(\bm{x}_{k-1})$), we set $\bm{x}_{k}=\bm{z}_{k}$, hence

[TABLE]

If the algorithm overshoots ($F(\bm{z}_{k})>F(\bm{x}_{k-1})$) then, by monotonicity, we impose $\bm{x}_{k}=\bm{x}_{k-1}$, which leads to

[TABLE]

*Combining (88) and (89) completes the proof. *

Until this point we have assumed that the first iteration $k=0$ does not use the auxiliary point extrapolation rule in Proposition 4.3. To write our generalized ACGM in a form similar to monotone FISTA (MFISTA [26]) and the monotone version of FISTA-CP [7], we define the vertex extrapolation factor in Lemma 4.1 as

[TABLE]

Proposition 4.5**.**

For $k\geq 1$, the vertex extrapolation factor can be determined using a recursion rule that does not depend on weights $a_{k}$ and $A_{k}$, given by

[TABLE]

Subexpression $\omega_{k}$ and auxiliary point extrapolation factor $\beta_{k}$ can also be written for all $k\geq 1$ as

[TABLE]

Proof 4.6**.**

In this proof we take all $k\geq 1$. From (73) it follows that

[TABLE]

Moreover,

[TABLE]

Putting together (94), (95), (96), and (97) we obtain

[TABLE]

Since $\gamma_{k+1}>0$, (91) holds. Subexpression $\omega_{k}$ is similarly obtained as

[TABLE]

*Substituting (90) in (85) gives (93). *

For simplicity, we wish to extend the update rules from Propositions 4.3 and 4.5 to the first iteration $k=0$. The missing parameters follow naturally from this extension. First, $t_{0}$ can be obtained by setting $k=0$ in (91) as

[TABLE]

Next, we introduce a “phantom iteration” $k=-1$ with the main iterate as the only state parameter. We set $x_{-1}\operatorname{\overset{def}{=}}x_{0}$ so that any value of $\beta_{0}$ will satisfy Proposition 4.3. For brevity, we choose to obtain $\beta_{0}$ from expression (93) with $k=0$.

Thus, with the initialization in (98) and the recursion in (91), we have completely defined the vertex extrapolation factor sequence $\{t_{k}\}_{k\geq 0}$, and derived from it the auxiliary extrapolation factor expression in (93). Now, we no longer need to maintain weight sequences $\{a_{k}\}_{k\geq 1}$ and $\{A_{k}\}_{k\geq 0}$. We simplify our generalized ACGM further by noting that, to produce the auxiliary point, the extrapolation rule in Proposition 4.3 depends on three vector parameters. However, it is not necessary to store both $\bm{z}_{k}$ and $\bm{x}_{k-1}$ across iterations. To address applications where memory is limited, we only maintain the difference term $\bm{d}_{k}$, given by

[TABLE]

The extrapolation rule in Proposition 4.3 becomes

[TABLE]

Note that subexpression $\omega_{k}$ contains only recent information whereas $\bm{d}_{k}$ needs only to access the state of the preceding iterations.

The above modifications yield a form of generalized ACGM based on extrapolation, which we list in Algorithm 3. To obtain a non-monotone algorithm, it suffices to replace line 16 of Algorithm 3 with (72).

We stress that while Algorithms 2 and 3 carry out different computations, they are mathematically equivalent with respect to the main iterate sequence $\{\bm{x}_{k}\}_{k\geq 0}$. The oracle calls and their dependencies in Algorithm 3 are also identical to those in Algorithm 2. Therefore the per-iteration complexity is the same.

4.2 Retrieving the convergence guarantee

For Algorithm 2, the convergence guarantee in (7) is obtained directly from the state variable $A_{k}$. For Algorithm 3, we need the following result.

Lemma 4.7**.**

The vertex extrapolation factor expression in (98) generalizes to arbitrary $k\geq 0$ as

[TABLE]

Proof 4.8**.**

For $k=0$, (98) implies that Lemma 4.7 holds. For $k\geq 1$ we have that

[TABLE]

From Lemma 4.7, we distinguish two scenarios.

If $\gamma_{0}\neq A_{0}\mu$, from (40), (73), and Lemma 4.7 it follows that the convergence guarantee can be derived directly from the state parameters without alterations to Algorithm 3 as

[TABLE]

4.3 Border-case

However, if $\gamma_{0}=A_{0}\mu$, then (40), (73), and Lemma 4.7 imply that

[TABLE]

Therefore, the state parameters of Algorithm 3 no longer contain information on the convergence guarantee but can be brought to a simpler form. It follows from (92), (93), and (102) that the auxiliary point extrapolation factor is given by

[TABLE]

The sequence $\{t_{k}\}_{k\geq 0}$ does not store any relevant information and can be left out. This means that the convergence guarantee $A_{k}$ requires a dedicated update. A simple recursion rule follows from (30), (40), and (71) as

[TABLE]

Due to scaling invariance, we can select any pair $(A_{0},\gamma_{0})$ that is a positive multiple of $(1,\mu)$. For simplicity, we choose $A_{0}=1$ and $\gamma_{0}=\mu$.

To reduce computational intensity, we modify subexpressions $\bm{d}_{k}$ and $\omega_{k}$ as

[TABLE]

The local inverse condition number sequence $\{q_{k}\}_{k\geq 0}$ does not appear in updates (104), (105), and (106). Hence, it can also be abstracted away. The form taken by generalized ACGM in this border-case, after simplifications, is listed in Algorithm 4.

A non-monotone variant can be obtained by replacing line 13 of Algorithm 4 with (72). The per-iteration complexity of Algorithm 4 matches the one of Algorithms 2 and 3.

5 Simulation results

5.1 Benchmark setup

We have tested the variants of generalized ACGM introduced in this work against the methods considered at the time of writing to be the state-of-the-art on the problem class outlined in Subsection 1.2. The proposed methods included in the benchmark are non-monotone ACGM (denoted as plain ACGM), monotone ACGM (MACGM), and, for strongly-convex problems, border-case non-monotone ACGM (BACGM) as well as border-case monotone ACGM (BMACGM). The state-of-the-art methods are FISTA-CP, monotone FISTA-CP (MFISTA-CP) [7], AMGS [10], and FISTA with backtracking line-search (FISTA-BT) [11].

We have selected as test cases five synthetic instances of composite problems in the areas of statistics, inverse problems, and machine learning. Three are non-strongly convex: least absolute shrinkage and selection operator (LASSO) [12], non-negative least squares (NNLS), and $l_{1}$-regularized logistic regression (L1LR). The other two are strongly-convex: ridge regression (RR) and elastic net (EN) [36]. Table 3 lists the oracle functions of all above mentioned problems.

Here, the sum softplus function $\mathcal{I}(\bm{x})$, the element-wise logistic function $\mathcal{L}(\bm{x})$, and the shrinkage operator $\mathcal{T}_{\tau}(\bm{x})$ are, respectively, given for all $\bm{x}\in\mathbb{R}^{n}$ and $\tau>0$ by

[TABLE]

To attain the best convergence guarantees for AMGS, Nesterov suggests in [10] that all known global strong convexity be transferred to the simple function $\Psi$. When line-search is enabled, generalized ACGM also benefits slightly from this arrangement when $r_{u}>1$ (Theorem 3.1). Without line-search, the convergence guarantees of generalized ACGM do not change as a result of strong convexity transfer, in either direction. Thus, for fair comparison, we have incorporated in $\Psi$ the strongly-convex quadratic regularization term for the RR and EN problems. In the following, we describe in detail each of the five problem instances. All random variables are independent and identically distributed, unless stated otherwise.

LASSO. Real-valued matrix $\bm{A}$ is of size $m=500$ by $n=500$, with entries drawn from $\mathcal{N}(0,1)$. Vector $\bm{b}\in\mathbb{R}^{m}$ has entries sampled from $\mathcal{N}(0,9)$. Regularization parameter $\lambda_{1}$ is $4$. The starting point $\bm{x}_{0}\in\mathbb{R}^{n}$ has entries drawn from $\mathcal{N}(0,1)$.

NNLS. Sparse $m=1000\times n=10000$ matrix $\bm{A}$ has approximately $10\%$ of entries, at random locations, non-zero. The non-zero entries are drawn from $\mathcal{N}(0,1)$ after which each column $j\in\{1,...,n\}$ is scaled independently to have an $l_{2}$ norm of $1$. Starting point $\bm{x}_{0}$ has $10$ entries at random locations all equal to $4$ and the remainder zero. Vector $\bm{b}$ is obtained from $\bm{b}=\bm{A}\bm{x}_{0}+\bm{z}$, where $\bm{z}$ is standard Gaussian noise.

L1LR. Matrix $\bm{A}$ has $m=200\times n=1000$ entries sampled from $\mathcal{N}(0,1)$, $\bm{x}_{0}$ has exactly $10$ non-zero entries at random locations, each entry value drawn from $\mathcal{N}(0,225)$, and $\lambda_{1}=5$. Labels $\bm{y}_{i}\in\{0,1\},i\in\{1,...,m\}$ are selected with probability $\mathbb{P}(\bm{Y}_{i}=1)=\mathcal{L}(\bm{A}\bm{x})_{i}$.

RR. Dimensions are $m=500\times n=500$. The entries of matrix $\bm{A}$, vector $\bm{b}$, and starting point $\bm{x}_{0}$ are drawn from $\mathcal{N}(0,1)$, $\mathcal{N}(0,25)$, and $\mathcal{N}(0,1)$, respectively. Regularizer $\lambda_{2}$ is given by $10^{-3}(\sigma_{max}(\bm{A}))^{2}$, where $\sigma_{max}(\bm{A})$ is the largest singular value of $\bm{A}$.

EN. Matrix $\bm{A}$ has $m=1000\times n=500$ entries sampled from $\mathcal{N}(0,1)$. Starting point $\bm{x}_{0}$ has $20$ non-zero entries at random locations, each entry value drawn from $\mathcal{N}(0,1)$. Regularization parameter $\lambda_{1}$ is obtained according to [37] as $\lambda_{1}=1.5\sqrt{2\log(n)}$ and $\lambda_{2}$ is the same as in RR, $\lambda_{2}=10^{-3}(\sigma_{max}(\bm{A}))^{2}$.

The Lipschitz constant $L_{f}$ is given by $(\sigma_{max}(\bm{A}))^{2}$ for all problems except for L1LR where it is $\frac{1}{4}(\sigma_{max}(\bm{A}))^{2}$. Strongly convex problems RR and EN have $\mu=\mu_{\Psi}=\lambda_{2}$ and inverse condition number $q=\mu/(L_{f}+\mu_{\Psi})=1/1001$.

To be able to benchmark against FISTA-CP and FISTA-BT, which lack fully adaptive line-search, we have set $L_{0}=L_{f}$ for all tested algorithms, thus giving FISTA-CP and FISTA-BT an advantage over the proposed methods. To highlight the differences between ACGM and BACGM, we ran ACGM and MACGM with parameters $A_{0}=0$ and $\gamma_{0}=1$.

Despite the problems differing in structure, the oracle functions have the same computational costs. We consider one matrix-vector multiplication to cost 1 WTU. Consequently, for all problems, we have $t_{f}=1$ WTU, $t_{g}=2$ WTU, and $t_{\Psi}=t_{p}=0$ WTU.

The line-search parameters were selected according to the recommendation given in [5]. For AMGS and FISTA-BT we have $r_{u}^{\rm{AMGS}}=r_{u}^{\rm{FISTA}}=2.0$ and $r_{d}^{\rm{AMGS}}=0.9$. The variants of generalized ACGM and AMGS are the only methods included in the benchmark that are equipped with fully adaptive line-search. We have decided to select $r_{d}^{\rm{ACGM}}$ to ensure that ACGM and AMGS have the same overhead. We formally define the line-search overhead of method $\mathcal{M}$, denoted by $\Omega^{\rm{\mathcal{M}}}$, as the average computational cost attributable to backtracks per WTU of advancement. Assuming that the LCEs hover around a fixed value (Subsection 3.3), we thus have that

[TABLE]

From (110) and (111) we have that $r_{d}^{\rm{ACGM}}=(r_{d}^{\rm{AMGS}})^{\frac{2}{3}}$, with no difference for border-case or monotone variants.

For measuring ISDs, we have computed beforehand an optimal point estimate $\hat{\bm{x}}^{*}$ for each problem instance. Each $\hat{\bm{x}}^{*}$ was obtained as the main iterate after running MACGM for $5000$ iterations with parameters $A_{0}=0$, $\gamma_{0}=1$, $L_{0}=L_{f}$, and aggressive search parameters $r_{d}=0.9$ and $r_{u}=2.0$.

5.2 Non-strongly convex problems

The convergence results for LASSO, NNLS, and L1LR are shown in Figure 1. The LCE variation during the first 200 WTU is shown in Figure 4. For NNLS, floating point precision was exhausted after 100 WTU and the LCE variation was only plotted to this point (Figure 4(b)). In addition, the average LCEs are listed in Table 4.

Both variants of ACGM outperform in iterations and especially in WTU the competing methods in each of these problem instances. Even though for LASSO and NNLS, the iteration convergence rate of AMGS is slightly better in the beginning (Figures 1(a) and 1(c)), AMGS lags behind afterwards and, when measured in terms of computational convergence rate, has the poorest performance among the methods tested (Figures 1(b), 1(d), and 1(f)). FISTA-BT produces the same iterates as FISTA-CP, as theoretically guaranteed in the non-strongly convex case for $L_{0}=L_{f}$.

The overall superiority of ACGM and MACGM can be attributed to the effectiveness of line-search. Interestingly, ACGM manages to surpass FISTA-CP and MFISTA-CP even when the latter are supplied with the exact value of the global Lipschitz constant. This is because ACGM is able to accurately estimate the local curvature, which is often well below $L_{f}$. For the L1LR problem, where the smooth part $f$ is not the square of a linear function, the local curvature is substantially lower than the global Lipschitz constant with LCEs hovering around one fifth of $L_{f}$ (Figure 4(c)). One would expect AMGS to be able to estimate local curvature as accurately as ACGM. This is does not happen due to AMGS’s reliance on a “damped relaxation condition” [10] line-search stopping criterion. For LASSO and NNLS, the average LCE of AMGS is actually above $L_{f}$. ACGM has an average LCE that is roughly two thirds that of AMGS on these problems whereas for L1LR the average is more than three times lower than AMGS. The difference between the LCE averages of ACGM and MACGM is negligible.

Indeed, monotonicity, as predicted, does not alter the overall iteration convergence rate and has a stabilizing effect. MACGM overshoots do have a negative but limited impact on the computational convergence rate. We have noticed in our simulations that overshoots occur less often for larger problems, such as the tested instance of NNLS.

5.3 Strongly convex problems

The convergence results for RR and EN are shown in Figures 3(a), 3(b), 3(c), and 3(d). The LCE variation is shown in Figure 5 with LCE averages listed in Table 5.

Strongly convex problems have a unique optimum point and accelerated first-order schemes are guaranteed to find an accurate estimate of it in domain space (see [4] for a detailed analysis). Along with Theorem 3.1, it follows that

[TABLE]

Thus, we can display accurate estimates of ISDUBs in (7), of the form $U_{k}\operatorname{\overset{def}{=}}\Delta_{0}/A_{k}$, $k\geq 1$, for methods that maintain convergence guarantees at runtime. These are shown in Figures 3(e) and 3(f) as upper bounds indexed in WTU.

For the smooth RR problem, the effectiveness of each tested algorithm is roughly given by the increase rate of the accumulated weights (Figures 3(a) and 3(c)). In iterations, AMGS converges the fastest. However, in terms of WTU usage, it is the least effective of the methods designed to deal with strongly convex objectives. The reasons are the high cost of its iterations, its low asymptotic rate compared to ACGM and FISTA-CP, and the stringency of its damped relaxation criterion that results in higher LCEs (on average) than ACGM (Figure 5(a) and Table 5). The computational convergence rate of BACGM is the best, followed by ACGM and FISTA-CP. This does not, however, correspond to the upper bounds (Figure 3(e)). While BACGM produces the largest accumulated weights $A_{k}$, the high value of the ISD term in $\Delta_{0}$ causes BACGM to have poorer upper bounds than ACGM, except for the first iterations. In fact, the effectiveness of BACGM on this problem is exceptional, partly due to the regularity of the composite gradients. This regularity also ensures the monotonicity of BACGM, ACGM, and FISTA-CP. FISTA-BT does not exhibit linear convergence on this problem and it is even slower than AMGS after 500 WTU despite its lower line-search overhead and advantageous parameter choice $L_{0}=L_{f}$.

On the less regular EN problem, ACGM leads all other methods in terms of both iteration and computational convergence rates (Figures 3(b) and 3(d)). The advantage of ACGM, especially over BACGM, is accurately reflected in the upper bounds (Figure 3(f)). However, convergence is much faster than the upper bounds would imply. Even FISTA-BT has a competitive rate, due to the small number of iterations ($150$) needed for high accuracy results. The ineffectiveness of AMGS on this problem is mostly due to its high LCEs (Figure 5(b) and Table 5). The proposed ACGM and variants show comparable average LCEs. Here as well, monotonicity has a stabilizing effect and does not have a significant impact on the computational convergence rate.

6 Discussion

Our simulation results suggest that enforcing monotonicity in ACGM is generally beneficial in large-scale applications. It leads to a more predictable convergence rate and, provided that the number of overshoots per iteration is small, monotonicity has a negligible impact on the computational convergence rate as well. Our experimental results indicate that the frequency of overshoots generally decreases with problem size.

From a theoretical standpoint, the proposed method can be viewed as a unification of FGM and FISTA, in their most common forms. Specifically, the fixed-step variant ($L_{k}=L_{f}$ for all $k\geq 0$) of ACGM in extrapolated form (Algorithm 3) is equivalent to both the monotone and non-monotone variants of FISTA-CP with the theoretically optimal step size $\tau^{\rm{FISTA-CP}}=1/L_{f}$. Moreover, when $\mu=0$, the original non-monotone fixed-step ACGM coincides with the original formulation of FISTA in [11]. Adding monotonicity yields MFISTA [26]. Also for $\mu=0$, but without the fixed-step restriction, the original non-monotone ACGM in estimate sequence form reduces to the robust FISTA-like algorithm in [24], whereas in extrapolated form it is a valuable simplification of the method introduced in [23]. The original backtracking FISTA can be obtained in the same way as the variants in [24, 23] by further removing the ratio $L_{k+1}/L_{k}$ from the $t_{k+1}$ update in line 7 of Algorithm 3 (noting that $\mu_{\Psi}=0$ for $\mu=0$). Since backtracking FISTA relies on the assumption that this ratio is never less that $1$, the convergence analysis remains unaffected.

When dealing with differentiable objectives, we can assume without loss of generality that $\Psi(\bm{x})=0$ for all $\bm{x}\in\mathbb{R}^{n}$. In what follows, we consider generalized non-monotone fixed-step ACGM in estimate sequence form, unless stated otherwise. By substituting the local upper bound functions $u_{k+1}$ at every iteration $k\geq 0$ with any functions that produce iterates satisfying the descent condition, which means in this context that

[TABLE]

where $\bm{x}_{k+1}$ is given by line 7 of Pattern 1, we obtain the “general scheme of optimal method” in [4]. Both the monotone and non-monotone variants adhere to this scheme. The correspondence between Nesterov’s notation in [4] and ours is for all $k\geq 0$ given by:

[TABLE]

The remaining state parameters are identical. Note that FGM makes the assumption that $A_{0}^{\rm{ACGM}}>0$, which is incompatible with the original specification of ACGM in [25]. With the above assumption, the non-monotone variant (Algorithm 2) is in fact identical to the “constant step scheme I” in [4]. Similarly, the extrapolated form of fixed-step non-monotone ACGM (Algorithm 3) corresponds exactly to the “constant step scheme II” in [4] while fixed-step border-case non-monotone ACGM (Algorithm 4) is identical to the “constant step scheme III” in [4]. We note that for both the RR and EN problems, regardless of the actual performance of BACGM, the convergence guarantees of BACGM are poorer than those of ACGM with $A_{0}=0$. This discrepancy in guarantees is supported theoretically because, in the most common applications, the ISD term in $\Delta_{0}$ is large compared to the DST. This extends to the fixed-step setup and challenges the notion found in the literature (e.g., [38]) that for strongly-convex functions, FGM and FISTA-CP are momentum methods that take the form of the “constant step scheme III” in [4]. In fact, the border-case form may constitute the poorest choice of parameters $A_{0}$ and $\gamma_{0}$ in many applications. Indeed, the worst-case results in Theorem 3.1 favor $A_{0}=0$.

The FGM variant in [22] is a particular case of non-monotone ACGM with variable step size (Algorithm 2) when the objective is non-strongly convex ($\mu=0$) and the step size search parameters are set to $r_{u}^{\rm{ACGM}}=2$ and $r_{d}^{\rm{ACGM}}=0.5$. The notation correspondence is as follows:

[TABLE]

The remaining parameters are identical.

Thus, by relaxing the assumption that $A_{0}=0$, we have devised a generalized variant of ACGM that effectively encompasses FGM [4], with its recently introduced variant [22], as well as the original FISTA [11], including its adaptive step-size variants [24, 23], the monotone version MFISTA [26], and the strongly convex extensions FISTA-CP [7] and scAPG [15]. A summary of how the above first-order methods relate to generalized ACGM is given in Table 6.

Due to its adaptivity, generalized ACGM is not limited to the composite problem framework in Subsection 1.2. It is also guaranteed to converge on problems where the gradient of $f$ is not globally Lipschitz continuous. Constituent function $f$ needs to have Lipschitz gradient only in the area explored by the algorithm.

In terms of usability, generalized ACGM does not require a priori knowledge of Lipschitz constant $L_{f}$, or a lower estimate of it, beforehand. Thus, the proposed method can be utilized without any quantitative knowledge of the problem. Lack of information does not hinder generalized ACGM more than any other primal first-order method while additional information, such as the values of strong convexity parameters $\mu_{f}$ and $\mu_{\Psi}$ or even an accurate estimate $L_{0}$ of the curvature around $\bm{x}_{0}$, leads to a state-of-the-art performance increase unsurpassed for its class.