An Optimal High-Order Tensor Method for Convex Optimization

TL;DR

This paper introduces an optimal high-order tensor algorithm for convex optimization that achieves the best possible iteration complexity, extending previous methods to the composite convex case with both smooth and non-smooth functions.

Contribution

It proposes a new high-order tensor algorithm with optimal iteration complexity for composite convex functions, matching established lower bounds.

Findings

Achieves iteration complexity of O(1/k^{(3d+1)/2}) for d-th order methods.

Matches the lower bounds for optimal convex optimization algorithms.

Uses an accelerated hybrid proximal extragradient framework with bisection steps.

Abstract

This paper is concerned with finding an optimal algorithm for minimizing a composite convex objective function. The basic setting is that the objective is the sum of two convex functions: the first function is smooth with up to the d-th order derivative information available, and the second function is possibly non-smooth, but its proximal tensor mappings can be computed approximately in an efficient manner. The problem is to find -- in that setting -- the best possible (optimal) iteration complexity for convex optimization. Along that line, for the smooth case (without the second non-smooth part in the objective), Nesterov (1983) proposed an optimal algorithm for the first-order methods (d=1) with iteration complexity O( 1 / k^2 ). A high-order tensor algorithm with iteration complexity of O( 1 / k^{d+1} ) was proposed by Baes (2009) and Nesterov (2018). In this paper, we propose a new high-order tensor algorithm for the general composite case, with the iteration complexity of O( 1 / k^{(3d+1)/2} ), which matches the lower bound for the d-th order methods as established in Nesterov (2018), and Shamir et al. (2018), and hence is optimal. Our approach is based on the Accelerated Hybrid Proximal Extragradient (A-HPE) framework proposed in Monteiro and Svaiter (2013), where a bisection procedure is installed for each A-HPE iteration. At each bisection step a proximal tensor subproblem is approximately solved, and the total number of bisection steps per A-HPE iteration is bounded by a logarithmic factor in the precision required.

Equations

$$\underset{x\in {\mathbb{R}}^n}{\min }F(x):=f(x)+h(x),$$

$${\nabla }^{d}f(x{)}_{i_1,\dots ,i_d}=\frac{{\partial }^{d}f}{\partial x_{i_1}\cdots \partial x_{i_d}}(x),\forall 1\le i_1,...,i_d\le n.$$

$$\mathcal{T}[z^1,\dots ,z^d]:=\sum _{i_1,\dots ,i_d=1}^n{\mathcal{T}}_{i_1,\dots ,i_d}{z}_{i_1}^1\dots z_{i_d}^d.$$

$$\Vert \mathcal{T}\Vert :=\underset{\Vert z^i\Vert =1,i=1,...,d}{\max }\mathcal{T}[z^1,...,z^d].$$

$$\mathcal{T}[z^{k+1},\dots ,z^d{]}_{i_1,\cdots ,i_k}:=\sum _{i_{k+1},\dots ,i_d=1}^n{\mathcal{T}}_{i_1,\dots ,i_k,i_{k+1},\dots ,i_d}{z}_{i_{k+1}}^{k+1}\dots z_{i_d}^d$$

$$(z_{\ast }^1,\dots ,z_{\ast }^k):=\underset{\Vert y^i\Vert =1,i=1,...,k}{\mathrm{argmax}}\left(\mathcal{T}[z^{k+1},...,z^d]\right)[y^1,...,y^k].$$

$$\Vert \mathcal{T}[z^{k+1},...,z^d]\Vert =\mathcal{T}[z_{\ast }^1,\dots ,z_{\ast }^k,z^{k+1},...,z^d]\le \Vert \mathcal{T}\Vert \Vert z^{k+1}\Vert \cdots \Vert z^d\Vert .$$

$$\Vert {\nabla }^{d}f(x)-{\nabla }^{d}f(y)\Vert \le L_d\Vert x-y\Vert \forall x,y\in {\mathbb{R}}^n,$$

$$\underset{y\in {\mathbb{R}}^n}{\min }f(y)+h(y)+\frac{1}{2\lambda }\Vert y-x{\Vert }^2.$$

$$f_x(y):=f(x)+\nabla f(x{)}^{\top }(y-x)+\frac{1}{2}{\nabla }^{2}f(x)[y-x{]}^2+\cdots +\frac{1}{d!}{\nabla }^{d}f(x)[y-x{]}^d+\frac{M}{(d+1)!}\Vert y-x{\Vert }^{d+1},$$

$$\Vert \nabla f(y)-\nabla f_x(y)\Vert \le \frac{L_d+M}{d!}\Vert y-x{\Vert }^d.$$

$$\underset{y\in {\mathbb{R}}^n}{\min }{f}_x(y)+h(y)+\frac{1}{2\lambda }\Vert y-x{\Vert }^2.$$

$$u\in (\nabla f_x+\partial h)(y),\lambda u+y-x=0.$$

$${\partial }_{ϵ}h(x):=\{u|h(y)\ge h(x)+⟨y-x,u⟩-ϵ,\forall y\in {\mathbb{R}}^n\}.$$

$$u\in (\nabla f_x+{\partial }_{ϵ}h)(y)\text{mbox}and\Vert \lambda u+y-x{\Vert }^2+2\lambda ϵ\le {\hat{\sigma }}^2\Vert y-x{\Vert }^2.$$

$$v_{k+1}\in {\partial }_{{ϵ}_{k+1}}F({\tilde{y}}_{k+1}),$$

$$\Vert {\lambda }_{k+1}{v}_{k+1}+{\tilde{y}}_{k+1}-{\tilde{x}}_k{\Vert }^2+2{\lambda }_{k+1}{ϵ}_{k+1}\le {\sigma }^2\Vert {\tilde{y}}_{k+1}-{\tilde{x}}_k{\Vert }^2$$

$${\tilde{x}}_k$$

$$a_{k+1}$$

$$A_{k+1}$$

$$x_{k+1}$$

$$\frac{1}{2}\Vert x_{\ast }-x_k{\Vert }^2+A_k(F(y_k)-F_{\ast })+\frac{1-{\sigma }^2}{2}\sum _{j=1}^k\frac{A_j}{{\lambda }_j}\Vert {\tilde{y}}_j-{\tilde{x}}_{j-1}{\Vert }^2\le \frac{1}{2}{D}^2.$$

$$\sum _{j=1}^k\frac{A_j}{{\lambda }_j}\Vert {\tilde{y}}_j-{\tilde{x}}_{j-1}{\Vert }^2\le \frac{D^2}{1-{\sigma }^2}.$$

$$A_k\ge \frac{1}{4}{\left(\sum _{j=1}^k\sqrt{{\lambda }_j}\right)}^2.$$

$$\Vert y_k-x_{\ast }\Vert \le \left(\frac{2}{\sqrt{1-{\sigma }^2}}+1\right)D.$$

$$\frac{d!{\sigma }_l}{L_d+M}\le {\lambda }_{k+1}\Vert y_{k+1}-{\tilde{x}}_k{\Vert }^{d-1}\le \frac{d!{\sigma }_u}{L_d+M}$$

$${\tilde{x}}_k=\frac{A_k}{A_k+a_{k+1}}{y}_k+\frac{a_{k+1}}{A_k+a_{k+1}}{x}_k$$

$$a_{k+1}=\frac{{\lambda }_{k+1}+\sqrt{{\lambda }_{k+1}^2+4{\lambda }_{k+1}{A}_k}}{2}.$$

$$v_{k+1}$$

$$A_{k+1}$$

$$x_{k+1}$$

$$v_{k+1}$$

Full text

Keywords: convex optimization; tensor method; acceleration; iteration complexity.

Mathematics Subject Classification: 90C60, 90C25, 47H05

d

1 Introduction

In this paper, we consider the following composite unconstrained convex optimization:

[TABLE]

where $f$ is differentiable and convex, and $h$ is convex but possibly non-smooth. In this context, we assume that convex tensor (polynomial) proximal mappings regarding $h$ can be approximately computed efficiently. Given that structure, a fundamental quest is to find an optimal algorithm that solves the above problem, using the available derivative information of the smooth part $f$.

In case $F(x)=f(x)$, and only the gradient information of $f$ is available, Nesterov [25] proposed a gradient-type algorithm, which achieves the overall iteration complexity of $O(1/k^{2})$, matching the lower bound on the iteration complexity of this class of solution methods, hence is known to be an optimal algorithm among all the first-order methods. Since Nesterov’s seminal work [25], especially in the recent years when the large scale machine learning applications have come under the spotlight, there has been a surge of research effort to extend Nesterov’s approach to more general settings; see e.g. [4, 12, 19, 14, 30], and/or to incorporate certain adaptive strategies to enhance the practical performances of the acceleration; see e.g. [20, 29, 13]. At the same time, there has also been a considerable research effort to fully understand the underpinning mechanism of the first-order acceleration phenomenon; see e.g. [7, 32, 33, 34].

When the Hessian information is available, Nesterov [26] proposed an acceleration scheme for cubic regularized Newton’s method, and he showed that the iteration complexity bound improves from $O\left(1/k^{2}\right)$ to $O\left(1/k^{3}\right)$. A few years later, Monteiro and Svaiter [24] proposed a totally different acceleration scheme, which they termed as Accelerated Hybrid Proximal Extragradient Method (A-HPE) framework, and they proved that if the second-order information is incorporated into the A-HPE framework then the corresponding accelerated Newton proximal extragradient method has a superior iteration complexity bound of $O\left(1/k^{7/2}\right)$ over $O\left(1/k^{3}\right)$. In 2018, Arjevani, Shamir and Shiff [31] showed that $O\left(1/k^{7/2}\right)$ is actually a lower bound for the oracle complexity of the second-order methods for convex smooth optimization. This shows that the accelerated Newton proximal extragradient method is an optimal second-order method.

As evidenced by the special cases $d=1$ and $d=2$, there is a clear tradeoff between the level of derivation information required and the overall iteration complexity improved. Therefore, a natural and important question arises:

What is the exact tradeoff relationship between $d$ and the worst-case iteration complexity?

Such question has been in fact raised and addressed in some way in recent works [5, 10, 11, 22] in the context of nonconvex optimization. For convex optimization, the accelerated cubic regularized Newton method was generalized to the general high-order case [3, 27] with the iteration complexity being $O\left(1/k^{d+1}\right)$, where $d$ is the order of derivative information used in the algorithm. Jiang, Lin and Zhang [18] extended Nesterov’s approach to accommodate the composite optimization (1.1) and relaxed the requirement on the knowledge of problem parameters such as the Lipschitz constants and the requirement on the exact solutions of the subproblems while maintaining the same iteration bound as in [3, 27]. Along the line of bounding the worst case iteration complexity using up to the $d$-th order derivative information, there have also been significant progresses as well. Arjevani, Shamir and Shiff [31] showed that the worst case iteration complexity of any algorithm in that setting cannot be better than $O\left(1/k^{(3d+1)/2}\right)$. A simplified analysis of the bound can be found in Nesterov [27]. So, there was a gap between the achieved iteration bound $O\left(1/k^{d+1}\right)$ and the best possible bound of $O\left(1/k^{(3d+1)/2}\right)$. Clearly at least one of the two bounds is improvable. In this paper, we aim to settle the above theoretical quest by providing a new implementable algorithm whose iteration complexity is precisely $O\left(1/k^{(3d+1)/2}\right)$. As a result, the tradeoff relationship discussed above is pinned down to be exactly $O\left(1/k^{(3d+1)/2}\right)$.

Our algorithm is based on the A-HPE framework of Monteiro and Svaiter [24], which is presented as Algorithm 1 in this paper. In fact, our algorithm specifies a way to generate an approximate solution through the use of high order derivative information by Taylor expansion. In each iteration, such approximate solution is computed by means of a bisection process. At each bisection step, a regulated convex tensor (polynomial) optimization subproblem is approximately solved. Moreover, we show that, to implement one A-HPE iteration, the number of bisection steps – each calling to solve a convex tensor subproblems – is upper bounded by a logarithmic factor in the inverse of the required precision. Our bisection procedure is similar to the one proposed in [24] for the case $d=2$; however, a key modification is applied which enables the removal of the so-called “bracketing stage” used in [24].

After submitting the first version of the paper, we became aware of two other independent works [15, 8] establishing similar iteration bounds as ours, with the main difference being that the focus of [15, 8] is on the smooth case: $F(x)=f(x)$, while our method accommodates a composite objective function. The common theoretical development by the three groups was subsequently jointly announced in the form of abstract at Conference on Learning Theory (COLT) [17]. It is also worth mentioning that other than the afore-mentioned three papers there are some other related works on high-order optimization methods [6, 1, 2] based on large-step A-HPE framework.

The rest of the paper is organized as follows. In Section 2, we introduce some preliminaries including the assumptions and the high-order oracle model used throughout this paper. Then we present our optimal tensor method and its iteration complexity analysis in Section 3. The line search subroutine being used in the main procedure of our optimal tensor method is presented and analyzed in Section 4. Finally, some technical proofs and lemmas are provided in the appendix.

2 Preliminaries

2.1 Notations

We denote $\nabla^{d}f(x)$ to be the $d$-th order derivative tensor at point $x$ of function $f$ with the $(i_{1},...,i_{d})$ component given as:

[TABLE]

Given a $d$-th order tensor $\mathcal{T}$ and vectors $z^{1},\dots,z^{d}\in\mathbb{R}^{n}$, we denote

[TABLE]

The operator norm associated with $\mathcal{T}$ is defined as:

[TABLE]

For given $z^{k+1},\dots,z^{d}$, $\mathcal{T}[z^{k+1},\dots,z^{d}]$ is a $k$-th order tensor with the associated $(i_{1},\cdots,i_{k})$ component defined as:

[TABLE]

for $1\leq i_{1},...,i_{k}\leq n$. Denote

[TABLE]

One has

[TABLE]

As a matter of convention, for quantities $x$ and $y$, we use the notation $y=\Theta(x)$ to indicate the relation that there are positive constants $a$ and $b$ such that $ax\leq y\leq bx$. If $a$ is absent, then we shall indicate the relation as $y=O(x)$.

2.2 High-Order Oracle Model and Regularized Tensor Approximation

In this paper, we consider the following high-order oracle model and the algorithm we are going to propose is such oracle model.

$d$-th Order Oracle Model

•

$f$ is $d$ times Lipschitz-continuous and differentiable with Lipschitz constant $L_{d}$ for $d$-th order derivative tensor; i.e.

$\|\nabla^{d}f(x)-\nabla^{d}f(y)\|\leq L_{d}\|x-y\|\ \ \ \forall x,y\in\mathbb{R}^{n},$

(2.2)

where the left side is the $d$-th order tensor operator norm.

•

Given any $x$, the oracle returns $f(x),\nabla f(x),\nabla^{2}f(x),...,\nabla^{d}f(x)$.

•

At iteration $k$, $x_{k}$ is generated from a deterministic function $h$ and the oracle’s responses at any linear combination of $x_{1},x_{2},...,x_{k-1}$ and $\nabla^{i}f(x_{1}),\nabla^{i}f(x_{2}),...,\nabla^{i}f(x_{k-1})$, where $1\leq i\leq d$.

Recall that the exact proximal minimization at point $x$ with stepsize $\lambda>0$ is defined as

[TABLE]

To utilize all the derivative information, we consider the regularized tensor approximation of $f(y)$ at point $x$:

[TABLE]

where $M>0$ is the parameter of the high-order regularization term $\|y-x\|^{d+1}$. Then, by (2.2) and the Taylor expansion, we can bound the gap between $f_{x}(\cdot)$ and $f(\cdot)$ for any $x$ (see Nesterov [27]):

Lemma 2.1

For every $x,y\in\mathbb{R}^{n}$,

[TABLE]

Therefore, it is natural to consider the tensor approximation of (2.3):

[TABLE]

In fact, (2.5) is the subproblem to be solved in the Optimal Tensor Method that will be introduced later. Note that similar subproblems have appeared in [24] and [27]. Specifically, the one used in [24] corresponds to $d=2$ in (2.5) without the term involving $\|y-x\|^{d+1}$ (i.e., $M=0$), while [27] uses the subproblem that only minimizes $f_{x}(y)$ (i.e., without the nonsmooth term $h(y)$ and the quadratic regularization term $\frac{1}{2\lambda}\|y-x\|^{2}$). In contrast, our above subproblem installs both the high-order and quadratic regularization terms.

Note that the unique solution $y$ of (2.5) is characterized by the following optimality condition:

[TABLE]

For a scalar $\epsilon\geq 0$, the $\epsilon$-subdifferential of a proper closed convex function $h$ is defined as:

[TABLE]

With the above notion in mind, let us consider the following approximate solution for (2.6) (hence (2.5)).

Definition 2.1

Given $(\lambda,x)\in\mathbb{R}_{++}\times\mathbb{R}^{n}$ and $\hat{\sigma}\geq 0$, the triplet $(y,u,\epsilon)\in\mathbb{R}^{n}\times\mathbb{R}^{n}\times\mathbb{R}_{+}$ is called a $\hat{\sigma}$-approximate solution of (2.5) at $(\lambda,x)$ if

[TABLE]

Obviously, if $(y,u)$ is the solution pair of (2.6), then $(y,u,0)$ is a $\hat{\sigma}$-approximate solution of (2.5) at $(\lambda,x)$ for any $\hat{\sigma}\geq 0$. In the rest of our analysis, we assume the availability of a subroutine which, for given $(\lambda,x)$ and $\hat{\sigma}>0$, returns a $\hat{\sigma}$-approximate solution $(y,u,\epsilon)$. Let us call this subroutine ATS (Approximate Tensor Subroutine). Different from [27], where a similar subproblem as (2.5) without a possible nonsmooth function $h(\cdot)$ and regularization term $\frac{1}{2\lambda}\|y-x\|^{2}$ is exactly solved, we only assume an approximate solution in the form of (2.7) is available and no further assumption on $h(\cdot)$ is required. Note that the possibly nonsmooth function $h(\cdot)$ can be viewed as a fixed parameter in ATS. Once $h(\cdot)$ is given, ATS could be called in each step of the bisection search, which itself is a subroutine in the main procedure of our algorithm.

3 The Optimal Tensor Method

3.1 The tensor algorithm and its iteration complexity

Our bid to the optimal tensor algorithm is based on the so-called Accelerated Hybrid Proximal Extragradient (A-HPE) framework proposed by Monteiro and Svaiter [24] for problem (1.1), whose main steps can be schematically sketched below:

Note that in STEP 2, the stopping condition is $0\in\partial F(y_{k})$. However, in practice, the condition is replaced by an approximate version of it (see Algorithm 2 below).

In the following, we quote some technical results derived in [24] for A-HPE. Since our proposed algorithm is within that framework, the results in Lemma 3.1 hold true for our method as well, and they will be used in the subsequent analysis.

Lemma 3.1

Suppose the sequence $\{x_{k},y_{k},\tilde{x}_{k},\tilde{y}_{k}\}$ is genernated from Algorithm 1. Let $x_{*}$ be the projection of $x_{0}$ onto the set of optimal value points $X_{*}$, $F_{*}$ be the optimal value, and $D$ be the distance from $x_{0}$ to $X_{*}$. Then for any integer $k\geq 1$, it holds that (Theorem 3.6 in [24]),

[TABLE]

Therefore,

[TABLE]

Furthermore, $A_{k}$ and $\lambda_{k}$ has the following relation (Lemma 3.7 in [24])

[TABLE]

If $y_{k}$ is chosen as $\tilde{y}_{k}$ for all k, the distance between $y_{k}$ and $x_{*}$ can be bounded as follows (Theorem 3.10 in [24]),

[TABLE]

Now we are ready to propose our optimal tensor method in Algorithm 2.

At this point, neither Algorithm 1 nor Algorithm 2 has been shown to be implementable. In fact, STEP 3 in both algorithms presented above remain unspecified. Since $\lambda_{k+1}$ appears in both (3.7) and (3.9), it is even unclear why such solutions as required by STEP 3 exist at all. Actually, the double roles played by $\lambda_{k+1}$ in (3.7) and (3.9) are crucial for the overall $O\left(1/k^{(3d+1)/2}\right)$ convergence rate. As a tradeoff, such $\lambda_{k+1}$ is not easy to find. In Section 4, we shall discuss a practical method to find a proper $\lambda_{k+1}$ (and thus establish a practical implementation of STEP 3 in Algorithm 2) via the Approximate Tensor Subroutine (ATS) in combination with a line-search subroutine.

First, let us remark that Algorithm 2 is indeed a specialization of A-HPE. For simplicity, we let $y_{k+1}=\tilde{y}_{k+1}$ in STEP 4 of Algorithm 1. Because $(y_{k+1},u_{k+1},\epsilon_{k+1})$ is a $\hat{\sigma}$-approximate solution at $(\lambda_{k+1},\tilde{x}_{k})$, one has that $u_{k+1}\in(\nabla f_{\tilde{x}_{k}}+\partial_{\epsilon_{k+1}}h)(y_{k+1})$, and so we have

[TABLE]

which satisfies (3.1). To establish (3.2), we need the following proposition.

Proposition 3.2

Let $(y,u,\epsilon)$ be a $\hat{\sigma}$-approximate solution of (2.5) at $(\lambda,{\tilde{x}})$ such that (3.7) holds. Define $v:=\nabla f(y)+u-\nabla f_{\tilde{x}}(y)$. Then,

[TABLE]

Consequently,

[TABLE]

where $\sigma_{u}$ is a input paramter in Algorithm 2 and also appears in (3.7).

Proof. First of all, according to Lemma 2.1, it follows that

[TABLE]

Combining the above inequality with (2.7), one has that

[TABLE]

proving the first inequality. Then, by the right hand side of (3.7), $\lambda\frac{L_{d}+M}{d!}\|y-\tilde{x}\|^{d-1}\leq\sigma_{u}$, and so the second inequality follows. $\Box$

We summarize the above discussion in the theorem below.

Theorem 3.3

Algorithm 2 is a manifestation of the A-HPE framework, and thus the results of Lemma 3.1 hold for the sequence generated by Algorithm 2.

Before addressing the implementation of STEP 3 in Algorithm 2, let us first present the overall iteration complexity of Algorithm 2, assuming STEP 3 could be implemented.

Theorem 3.4

*Let $D$ be the distance of $x_{0}$ to $X_{*}$. Then, for any integer $k\geq 1$, the iterate $y_{k}$ generated by Algorithm 2 satisfies:

[TABLE]

The above theorem establishes the $O(1/k^{\frac{3d+1}{2}})$ iteration complexity for Algorithm 2. Since Algorithm 2 falls into the category of the High-Order Oracle Model, whose iteration complexity has a lower bound of $O(1/k^{\frac{3d+1}{2}})$; see Arjevanim, Shamir and Shiff [31] and Nesterov [27]. The worst-case iteration complexity of Algorithm 2 matches this lower bound and it is therefore an optimal method.

3.2 Proof of Theorem 3.4

We first provide a recursive bound on $A_{k}$ as an intermediate step.

Proposition 3.5

Let $D$ be the distance of $x_{0}$ to $X_{*}$. Suppose $\{A_{k}\}_{k=1}^{\infty}$ is generated from Algorithm 2, then

[TABLE]

where $q=\frac{3d+1}{d-1}$ , $p=\frac{3d+1}{2d+2}$ and $C=\frac{D^{2}}{(1-({\hat{\sigma}}+\sigma_{u})^{2})}\left(\frac{d!\sigma_{l}}{L_{d}+M}\right)^{-\frac{2}{d-1}}$.

Proof. Suppose $\{x_{k},y_{k},\tilde{x}_{k}\}$ is the sequence generated by Algorithm 2. Then, according to (3.4) and Proposition 3.2, it holds that

[TABLE]

which together with the left hand side of (3.7) implies

[TABLE]

By the definition of $p$ and $q$, we have $\frac{1}{p}+\frac{1}{q}=1$. Using Hölder’s inequality, together with (3.14), we have

[TABLE]

Finally, by (3.5) we obtain

[TABLE]

$\Box$

Proof of Theorem 3.4. Let $p$, $q$ and $C$ be defined as in Proposition 3.5. Construct $\{B_{k}\}$ such that $B_{1}=A_{1}$ and $B_{i}=T^{\frac{1-(2p/q)^{i-1}}{1-2p/q}}(A_{1})^{(2p/q)^{i-1}}$ for $i\geq 2$, where $T:=\frac{1}{4}(\frac{1}{C})^{\frac{2p}{q}}(\frac{2}{d+1})^{2p}$. Next, we shall apply induction to show that for any $k\geq 1$,

[TABLE]

where $r_{i}=\frac{3d+1}{2}\left[1-(2p/q)^{i-1}\right]$. When $i=1$, this is obvious because $A_{k}\geq A_{1}=B_{1}k^{r_{1}}$. Now suppose that for any $k\geq 1$, $A_{k}\geq B_{i}k^{r_{i}}$ for some $i$. Then, by the induction hypothesis and (3.13) it holds that

[TABLE]

where the last inequality follows from

[TABLE]

Let us further simplify the expression. First of all, from the definition of $T$ and $B_{i}$, one observes that

[TABLE]

Then, the construction of $q$ and $r_{i}$ implies that

[TABLE]

where the second last equality holds true due to the fact that $2p/q=(d-1)/(d+1)$. Now the desired inequality (3.15) follows by combining (3.16), (3.17) and (3.18). Observe that $2p/q=(d-1)/(d+1)<1$, and so $\lim\limits_{i\to\infty}B_{i}=T^{\frac{1}{1-2p/q}}=T^{\frac{d+1}{2}}$ and $\lim\limits_{i\to\infty}r_{i}=\frac{3d+1}{2}$. Finally, by letting $i\rightarrow\infty$ in (3.15) and using the definition of $C$ in (3.14), we have

[TABLE]

Combining it with (3.3), we have

[TABLE]

$\Box$

3.3 Comparison with Nesterov’s Accelerated Tensor Method

In Nesterov’s accelerated tensor method [27], an auxiliary function

[TABLE]

with $l_{k}$ being some linear function, is constructed to satisfy

[TABLE]

where $A_{k}=\Theta(k^{d+1})$. In fact, the function $\psi_{k}(x)$ serves as a bridge to guarantee the following relation:

[TABLE]

As a result, $F(y_{k})-F_{*}\leq\frac{M}{A_{k}}\|x_{*}-x_{0}\|^{d+1}$ yielding the iteration complexity of $O(1/k^{d+1})$.

In the implementation of high-order A-HPE framework, it is crucial to ensure that condition (3.7) is satisfied. In the remainder of the paper, we shall focus on how to satisfy (3.7) in STEP 3 of Algorithm 2. Our bid is to use bisection on a parameter $\lambda$ (to be introduced later), while calling an Approximate Tensor Subroutine (ATS). Observe that $\tilde{x}$, which is the point to define $f_{\tilde{x}}(y)$ in (2.4) to approximate the smooth function $f(y)$, is indeed heavily dependent on $\lambda$. In other words, we need to search for the point where the Taylor expansion (2.4) is to be computed. This is a key difference between the A-HPE framework and Nesterov’s approach [27]. Once condition (3.7) is satisfied, then inequality (3.3) would follow, which leads to the following tighter estimation than (3.20):

[TABLE]

as $\beta_{k}=\frac{1-\sigma^{2}}{2}\sum^{k}_{j=1}\frac{A_{j}}{\lambda_{j}}\|\tilde{y}_{j}-\tilde{x}_{j-1}\|^{2}\geq 0$ is totally missing in (3.20). The above inequality also gives an upper bound on $\beta_{k}$. Together with the lower bound (3.6) this gives a better lower bound on $A_{k}$, namely $A_{k}\geq O(k^{\frac{3d+1}{2}})$, which leads to the optimal iteration complexity presented in Theorem 3.4.

4 A Line Search Subroutine and Its Iteration Complexity

After establishing the overall iteration complexity for Algorithm 2, it remains to find a way to implement STEP 3 of the algorithm. In this section we discuss how this can be done, from a special case to the general one. The idea is better illustrated by considering the special case. Finally, for the general composite objective function, assuming the tensor proximal mapping regarding $h(x)$ is possible, our approach is based on a line-search procedure for the point on which the Taylor expansion is computed.

4.1 The Non-Composite Case

Let us first consider a special case for Algorithm 2 where $F(x)=f(x)$ in the objective function and $y_{k+1}$ is the exact solution of the following convex tensor proximal point problem:

[TABLE]

We shall discuss how to find $\lambda_{k+1}$ to satisfy the alternative condition in STEP 3 of Algorithm 2.

Note that for fixed $x_{k}$ and $y_{k}$, $\tilde{x}_{k}$ and $y_{k+1}$ are uniquely determined by $\lambda_{k+1}$. Therefore the functions $\tilde{x}_{k}(\lambda)$ and $y_{k+1}(\lambda)$ are continuous with respect to $\lambda$ (where we denote $\lambda_{k+1}$ to be $\lambda$). Next, we show that:

(i)

$\lambda\|y_{k+1}(\lambda)-\tilde{x}_{k}(\lambda)\|^{d-1}\rightarrow 0$, as $\lambda\rightarrow 0$;

(ii)

Either there exists an increasing sub-sequence $\lambda_{j}\uparrow\infty$, such that $\lambda_{j}\|y_{k+1}(\lambda_{j})-\tilde{x}_{k}(\lambda_{j})\|^{d-1}\rightarrow\infty$ as $j\rightarrow\infty$, or there exists $\hat{\lambda}$ such that $\|{\nabla f(y_{k+1}(\lambda))}\|\leq\bar{\rho}$ for any $\lambda\geq\hat{\lambda}$.

Observe that

[TABLE]

where $f(\tilde{x}_{k}(\lambda))$ is bounded, since $\tilde{x}_{k}(\lambda)$ is a convex combination of $x_{k}$ and $y_{k}$. Letting $\lambda\rightarrow 0$ in the above inequality leads to $\|y_{k+1}(\lambda)-\tilde{x}_{k}(\lambda)\|^{2}\rightarrow 0$, which implies $\lambda\|y_{k+1}(\lambda)-\tilde{x}_{k}(\lambda)\|^{d-1}\rightarrow 0$ as $\lambda\rightarrow 0$, proving (i).

To prove (ii), it suffices to show that if the “either” part does not hold, then the “or” part must hold. In this case, there must exist $C_{1}>0$ such that when $\lambda\rightarrow\infty$, $\lambda\|y_{k+1}(\lambda)-\tilde{x}_{k}(\lambda)\|^{d-1}\leq C_{1}$, and thus $\|y_{k+1}(\lambda)-\tilde{x}_{k}(\lambda)\|\rightarrow 0$. Moreover, for any $\lambda>0$ the optimality condition is

[TABLE]

Letting $\lambda\rightarrow\infty$ in the above identity yields that $\nabla f_{\tilde{x}_{k}(\lambda)}(y_{k+1}(\lambda))\rightarrow 0$. Recall that in this case we have $\|y_{k+1}(\lambda)-\tilde{x}_{k}(\lambda)\|\rightarrow 0$, thus $\nabla f(y_{k+1}(\lambda))\rightarrow 0$ proving the “or” part.

To summarize, either we have $\lambda\|y_{k+1}(\lambda)-\tilde{x}_{k}(\lambda)\|^{d-1}\rightarrow 0$ as $\lambda\rightarrow 0$ and $\lambda_{j}\|y_{k+1}(\lambda_{j})-\tilde{x}_{k}(\lambda_{j})\|^{d-1}\rightarrow\infty$ as $j\rightarrow\infty$, which guarantees the existence of $\lambda$ to satisfy (3.7) due to the continuity of $\lambda\|y_{k+1}(\lambda)-\tilde{x}_{k}(\lambda)\|^{d-1}$ on $\lambda$. Or we have a $\lambda_{k+1}$ such that $\|{\nabla f(y_{k+1}(\lambda))}\|\leq\bar{\rho}$. In this case, since $h(x)$ is not present, $u_{k+1}=\nabla f_{\tilde{x}_{k}}(y_{k+1})$ and $\|\nabla f(y_{k+1})+u_{k+1}-\nabla f_{\tilde{x}_{k}}(y_{k+1})\|=\|\nabla f(y_{k+1})\|\leq\bar{\rho}$. Therefore, we have shown that the alternative condition in STEP 3 is actually satisfied.

4.2 A Bisection Subroutine

To present the algorithm that computes $\lambda$ satisfing the conditions in STEP 3, we first construct $\beta_{k+1}=\frac{a_{k+1}}{A_{k}+a_{k+1}}$. From (3.9), we can see that $\lambda_{k+1}=\frac{a_{k+1}^{2}}{A_{k}+a_{k+1}}$. Therefore, we are able to represent $\lambda_{k+1}$ and $\tilde{x}_{k}$ by means of $\beta_{k+1}$:

[TABLE]

In the $k$-th iteration, we denote

[TABLE]

Its inverse on the domain $\lambda>0$ is

[TABLE]

which is monotonically increasing.

We shall perform bisection on $\beta$ instead of $\lambda$ in STEP 3 of Algorithm 2 to search for $\lambda_{k+1}$. In that way, the initial interval for the bisection is $[0,1]$. (Monteiro and Svaiter [24] presented a bisection process for their A-HPE algorithm too. However, we can skip what they called the bracketing stage in [24]).

We remark that the conditions on $\bar{\rho}$ and $\bar{\epsilon}$ are only used in the final stage of the algorithm to decide the point that is close to optimum. In the implementation, it is reasonable to set a lower precision at the beginning stage of the algorithm. Now an upper bound for the overall number of iterations required by Algorithm 3 is presented in the following theorem, whose proof will be postponed to the subsequent section.

Theorem 4.1

Algorithm 3 needs to perform no more than

[TABLE]

bisection steps before reaching $\lambda_{k+1}>0$ and a $\hat{\sigma}$-approximate solution $(y_{k+1},u_{k+1},\epsilon_{k+1})$ at $(\lambda_{k+1},\tilde{x}_{k}(\lambda_{k+1}))$ satisfying

[TABLE]

or to return $v_{k+1}$ and $\epsilon_{k+1}$ such that $\|v_{k+1}\|\leq\bar{\rho}$ and $\epsilon_{k+1}\leq\bar{\epsilon}$.

4.3 The Iteration Complexity Analysis

In this subsection, we establish the iteration bound of Algorithm 3 and give a proof for Theorem 4.1. First, we review some facts for maximal monotone operator. For a point-to-set operator $T:\mathbb{R}^{n}\rightrightarrows\mathbb{R}^{n}$, its graph is defined as:

[TABLE]

and the operator $T$ is called monotone if

[TABLE]

and $T$ is maximal monotone if it is monotone and maximal in the family of monotone operators with respect to the partial order of inclusion. Given a maximal monotone operator $T:\mathbb{R}^{n}\rightrightarrows\mathbb{R}^{n}$ and a scalar $\epsilon$, the associated $\epsilon$-enlargement $T^{\epsilon}:\mathbb{R}^{n}\rightrightarrows\mathbb{R}^{n}$ is defined as:

[TABLE]

For a convex function $f$, its subdifferential $\partial f$ is monotone if $f$ is a proper function. If $f$ is a proper lower semicontinuous convex function, then $\partial f$ is maximal monotone [28].

Recall that the optimality condition of subproblem (2.5) is characterized by (2.6), which is:

[TABLE]

Furthermore, $x$ is optimal to (1.1) if and only if $y=x$. Therefore, it is natural to consider the residual

[TABLE]

for any $\lambda>0,x\in\mathbb{R}^{n}$. The above residual was adopted in [24] for the quadratic subproblem. In this paper, to accommodate the high-order information, we consider the following modified residual:

[TABLE]

We have an immediate property regarding $\psi(\cdot)$.

Proposition 4.2

Let $x\in\mathbb{R}^{n},\lambda>0$ and $\hat{\sigma}\geq 0$. If $(y,u,\epsilon)$ is a $\hat{\sigma}$-approximate solution of (2.5) at $(\lambda,x)$, then

[TABLE]

Proof. From proposition 7.3 in [24], it holds that

[TABLE]

Notice $\varphi^{d-1}(\lambda;x)=\lambda^{d-2}\psi(\lambda;x)$, and so (4.4) readily follows by combining the above inequalities and identity. $\Box$

Lemma 4.3

Let scalars $\bar{\rho}>0$, $\bar{\epsilon}>0$, $\hat{\sigma}\geq 0$ and $\alpha>0$ be given and satisfy $\hat{\sigma}+\frac{L_{d}+M}{d!}\alpha:=\sigma<1$. Suppose

[TABLE]

and $(y,u,\epsilon)$ is a $\hat{\sigma}$-approximate solution of (2.5) at $(\lambda,x)$ for some vector $x\in\mathbb{R}^{n}$. Then, one of the following holds: either (a) $\lambda\|y-x\|^{d-1}>\alpha$; or (b) the vector $v:=\nabla f(y)-\nabla f_{x}(y)+u$ satisfies

[TABLE]

Proof. Suppose that $\lambda$ satisfies (4.6) but not (a), namely

[TABLE]

In that case, recall that $\partial h_{\epsilon}$ is the $\epsilon$-subdifferential of $h$ and $(\partial h)^{\epsilon}$ is the $\epsilon$-enlargement of operator $\partial h$. According to Proposition 3 in [9], one has $\partial h_{\epsilon}(x)\subseteq(\partial h)^{\epsilon}(x)$ for any $\epsilon\geq 0$ and $x\in\mathbb{R}^{n}$. Therefore, the inclusion in (4.7) directly follows from Proposition 3.2. Moreover, inequality (3.11) leads to

[TABLE]

Together with (4.6) and (4.8), the above inequality yields

[TABLE]

On the other hand, inequality (3.11) also implies that

[TABLE]

Combined with (4.6) and (4.8) this leads to

[TABLE]

Hence, (b) must hold in this case. $\Box$

In the rest of this section, we simply let $\alpha=\alpha_{-}$ in Lemma 4.3 and denote

[TABLE]

Lemma 4.3 implies that if $\lambda$ is sufficiently large, then either Algorithm 3 stops because (4.7) is satisfied, or $\lambda\|y-x\|^{d-1}\geq\alpha_{-}$, which achieves half of the bisection goal. Now we are ready to prove Theorem 4.1.

Proof of Theorem 4.1. Suppose that Algorithm 3 has performed $j$ bisection steps before triggering the stopping criteria. We aim to show $j\leq\Theta\left(\max\{\log_{2}(\bar{\epsilon}^{-1}),\log_{2}(\bar{\rho}^{-1})\}\right)$. At that iteration let us denote $x_{+}=x_{\beta_{+}}$, $x_{-}=x_{\beta_{-}}$, $y_{+}=y_{\beta_{+}}$ and $y_{-}=y_{\beta_{-}}$, and we also have $\beta_{+}-\beta_{-}=\frac{1}{2^{j}}$. Denote $\bar{\beta}=\beta(\bar{\lambda})$, where $\bar{\lambda}$ is as defined in (4.9). If $\bar{\beta}\leq\frac{1}{2}$ then $\frac{1}{1-\bar{\beta}}\leq 2$; if $\bar{\beta}>\frac{1}{2}$, then (4.1) gives

[TABLE]

Therefore, in the rest of the proof we may assume $j\geq\log_{2}(2/(1-\bar{\beta}))$, for otherwise $j<\log_{2}(2/(1-\bar{\beta}))\leq\Theta(\max\{\log_{2}(\bar{\rho}^{-1}),\log_{2}(\bar{\epsilon}^{-1})\})$ already holds.

Note that the bisection search starts with $\beta_{+}=1$, corresponding to $\lambda_{+}=+\infty$ according to (4.1) when $\beta_{+}$ is not updated during the procedure. However, the following lemma tells us that after running Algorithm 3 for a number of iterations, $\lambda_{+}$ will be reduced and upper bounded by some constant depending on $\bar{\epsilon}$ and $\bar{\rho}$.

Lemma 4.4

Suppose that Algorithm 3 has performed $j$ bisection steps with $j\geq\log_{2}(2/(1-\bar{\beta}))$, where $\bar{\beta}=\beta(\bar{\lambda})$ and $\bar{\lambda}$ are as defined in (4.9). Then we have

[TABLE]

We shall continue our discussion without disruption here and leave the proof of Lemma 4.4 to the appendix. Since Algorithm 3 did not stop before iteration $j$, the bound on $\beta_{+}$ must have been previously updated, and so

[TABLE]

where $\lambda_{+}$ is upper bounded due to Lemma 4.4.

By Proposition 4.2, we have that

[TABLE]

Consequently,

[TABLE]

The parameters $\alpha_{+}$ and $\alpha_{-}$ are pre-specified. Therefore, it suffices to show that $\psi_{+}-\psi_{-}$ is upper bounded by $\beta_{+}-\beta_{-}$ multiplied by some constant factor and hence the number of bisection search $j$ can be bounded as well. To this end, denote

[TABLE]

Then, there exist

[TABLE]

and

[TABLE]

To proceed, we have the following bound on $\lambda_{-}^{2}\|\bar{u}_{+}-\bar{u}_{-}\|$ whose proof can be found in the appendix.

Lemma 4.5

It holds that

[TABLE]

Note that

[TABLE]

for any $a,b>0$. Now combining (4.13), (4.14), (4.15) and (4.16) we have

[TABLE]

Next, by applying (1.9), (4.9), Lemma A.2, Lemma A.5 and Lemma A.6, we have

[TABLE]

Combining the bounds above with (4.17) yields

[TABLE]

where the last inequality is due to $d\geq 2$. Because $\beta_{+}-\beta_{-}=\frac{1}{2^{j}}$, from (4.11) we have

[TABLE]

The left hand side of the above inequality is a positive constant. Therefore,

[TABLE]

as required. $\Box$

Remark 4.6

In fact, we can quantify the constants in the proof of Theorem 4.1 more explicitly, and obtain the exact form of the bound $\Theta\left(\max\{\log_{2}(\bar{\epsilon}^{-1}),\log_{2}(\bar{\rho}^{-1})\}\right)$. Recall that

[TABLE]

Introduce the following constants

[TABLE]

where

[TABLE]

and $B_{1},...,B_{d}$ is a sequence defined by

[TABLE]

Then the complexity bound in Theorem 4.4 can be explicitly expressed by

[TABLE]

It is clear that the dependence of the resulting bound depends logarithmically on the parameters $L_{d},D$, and input parameters $\alpha_{+},\alpha_{-},\sigma_{u}$, and polynomially on $d$. The derivation of (4.18) is skipped for the sake of succinctness.

Now, for a given $\epsilon>0$, we denote

[TABLE]

Combining the bounds provided in Theorem 3.4, Theorem 4.1 and (4.18), we obtain the overall iteration bound for Algorithm 2 in terms of the ATS calls as follows:

Theorem 4.7

Given $\epsilon>0$. Assume that Algorithm 2 is implemented with $M\leq 2L_{d}$. Set $\bar{\epsilon}=\epsilon/2$, $\bar{\rho}\leq\min\left\{\frac{\epsilon}{2\bar{D}_{\epsilon}},\ \epsilon\right\}$, and define

[TABLE]

where

[TABLE]

Then, a point $z\in\mathbb{R}^{n}$ satisfying

[TABLE]

can be found by Algorithm 2 with no more than $K_{\epsilon}$ calls of ATS.

Proof. We consider two cases separately. In the first case, Algorithm 2 terminates because we find a $k\leq K_{\epsilon}$ such that $\|v_{k}\|\leq\bar{\rho}$ and $\|\epsilon_{k}\|\leq\bar{\epsilon}$. As $v_{k}=\nabla f(y_{k})-\nabla f_{x_{k}}(y_{k})+u_{k}$ and $u_{k}\in\nabla f_{x_{k}}(y_{k})+\partial_{\epsilon_{k}}h(y_{k})$, we have $v_{k}\in\nabla f(y_{k})+\partial_{\epsilon_{k}}h(y_{k})$. Let $x_{*}$ be the projection of $x_{0}$ onto $X_{*}$. By the convexity of $f$ and $h$,

[TABLE]

Summing up the two inequalities above yields

[TABLE]

By the construction of $\bar{\rho}$, we have that $\|v_{k}\|\leq\bar{\rho}\leq\epsilon$. Together with the definition of $\bar{D}_{\epsilon}$, this implies that $\|y_{k}-x_{*}\|\leq\bar{D}_{\epsilon}$. Again, by evoking the construction of $\bar{\rho}$ and $\bar{\epsilon}$, it holds

[TABLE]

In the other case, condition (3.7) holds for every $k\leq K_{\epsilon}$. Then, according to Theorem 3.4, for all $k\leq K_{\epsilon}$ we have

[TABLE]

By the definition of $K_{\epsilon}$ and letting $k=K_{\epsilon}$ in the above inequality, one has

[TABLE]

$\Box$

Note that the implementation of our framework is based on the assumption that the ATS can be efficiently computed. By the construction, the objective in (2.5) has a $\frac{1}{\lambda}$-strongly convex smooth part, which can be solved by many start-of-the-art optimization algorithms. However, there is few efficient algorithms customized for problem (2.5). Without further knowledge of the problem structure, the proposed approach is not necessarily more efficient as compared to, e.g., a direct application of a general-purpose convex optimization algorithm on (1.1). At the end of the paper, we shall briefly discuss a method for the subroutine in the special case $d=3$ introduced by Nesterov [27], which opens the door for this line of research. Of course, how to efficiently solve ATS in general remains a further research topic.

5 Concluding Remark

To conclude this paper, we shall discuss how to compute ATS efficiently with $d=3$. Note that in STEP 2.a of Algorithm 3, an Approximate Tensor Subroutine (ATS)

is required, which can be implemented in polynomial time in the case of convex optimization. In some applications, ATS may be implemented efficiently if some additional structures on the tensor (Taylor) expansion and/or the $h$ function exist. In this subsection, we show how ATS (i.e., solve problem (2.5)) may be computed efficiently in the absence of the non-smooth part, i.e. $F(x)=f(x)$, when $d=3$.

Note that since $h=0$, the $\epsilon_{\beta}$ in the bisection subroutine may be simply set to 0.

In this case, the objective function in (2.5) becomes: $f_{x}(y)+\frac{1}{\lambda}\|y-x\|^{2}=f(x)+\Omega(y-x)$ where

[TABLE]

Therefore, the subproblem (2.5) is equivalent to $\min_{z\in\mathbb{R}^{n}}\Omega(z)$. Let $M=3\kappa^{2}L_{3}$ with $\kappa>1$. Then, a similar argument as in Lemma 4 of [27] implies that function $\Omega(z)$ satisfies the strong relative smoothness condition

[TABLE]

with respect to function

[TABLE]

Such condition allows to minimize $\Omega(z)$ efficiently by a gradient method described in [21, 27], where we need to solve the following problem in every iteration:

[TABLE]

which was considered at the end of Section 5 in [27]. According to a min-max argument in [27], the above problem is shown to be equivalent to

[TABLE]

which is actually a univariate optimization problem with a strongly convex and analytic objective function, hence is easily solvable in practice. Note that a matrix inverse operation is required in the univariate optimization above. Since in all iterations of the gradient method described in [21, 27], the matrix $A$ is exactly $\nabla^{2}f(x)$ throughout and only the vector $a$ varies, the matrix inverse operation needs to be performed only once. As the gradient method in [21, 27] is linearly convergent, the total computational cost of an ATS call is in the order of $O(n^{3}+n^{2}\log(\bar{\rho}^{-1}))$.

Acknowledgments.

We would like to thank Tianyi Lin at UC Berkeley for the helpful discussions at the early stages of the project.

Appendix A Proofs of the lemmas in Section 4

We first establish an uniform lower bound as well as an upper bound for the sequence $\{A_{k}\}$.

Lemma A.1

Let $D$ be the distance of $x_{0}$ to $X_{*}$. Suppose $\{A_{k}\}_{k=1}^{\ell}$ is generated from Algorithm 2, and the algorithm has not stopped at iteration $\ell$. Then for any integer $1\leq k\leq\ell$, it holds that

[TABLE]

and

[TABLE]

Proof. We first establish the lower bound. Since $\{A_{k}\}$ is monotonically increasing, it suffices to lower bound $A_{1}$. Recall that $A_{0}=0$ and $A_{1}=A_{0}+a_{1}=\lambda_{1}$, and the choice of large-step (3.7) in Algorithm 2 leads to

[TABLE]

Moreover, Lemma 3.1 implies that

[TABLE]

where $x_{*}$ is the projection of $x_{0}$ onto the optimal solution set $X_{*}$. Combining the above two inequalities with the fact that $\tilde{x}_{0}=x_{0}$, it follows that

[TABLE]

Therefore,

[TABLE]

which is a uniform lower bound of the sequence $\{A_{k}\}$.

Next, we provide the upper bound. By invoking (3.12) to $(y_{k},v_{k},\epsilon_{k},\lambda_{k},\tilde{x}_{k-1})$, it holds that

[TABLE]

Then, combining (1.4) with (3.4) leads to

[TABLE]

Moreover, it follows from (3.14) that

[TABLE]

where $\alpha_{-}=\frac{d!\sigma_{l}}{L_{d}+M}$. Combining the above two inequalities yields

[TABLE]

or equivalently,

[TABLE]

On the other hand, (1.5) together with (3.4) implies that

[TABLE]

Consequently,

[TABLE]

Now, since the algorithm has not been terminated, we have either $\|v_{k}\|\geq\bar{\rho}$ or $\epsilon_{k}\geq\bar{\epsilon}$, which combined with (1.7) and (1.8) yields that

[TABLE]

and the conclusion follows. $\Box$

Below we prove Lemma 4.4 and Lemma 4.5 respectively.

Proof of Lemma 4.4. We first demonstrate that

[TABLE]

Otherwise, we have $\beta_{-}\geq\bar{\beta}$ and $\lambda_{-}>\bar{\lambda}$ as the function $\lambda(\beta)$ is strictly increasing in $\beta$. This together with Lemma 4.3 and the fact that $\lambda_{-}\|y_{-}-x_{-}\|\leq\alpha_{-}$, implies that the vector $v_{\beta_{-}}=\nabla f(y_{-})-\nabla f_{x_{-}}(y_{-})+u_{\beta_{-}}$ satisfies

[TABLE]

and thus the algorithm would have been terminated, yielding a contradiction. So we have $\beta_{-}\leq\bar{\beta}$. Since $j\geq\log_{2}(2/(1-\bar{\beta}))$, we have

[TABLE]

Together with the monotonicity of the function $\lambda(\beta)$ this implies that

[TABLE]

Therefore, $\lambda_{+}\leq A_{k}\frac{\left(1+\bar{\beta}\right)^{2}}{2(1-\bar{\beta})}\leq\left(\frac{3}{2}\right)^{2}A_{k}$ when $\bar{\beta}\leq\frac{1}{2}$ and $\lambda_{+}\leq\bar{\lambda}\frac{\left(1+\bar{\beta}\right)^{2}}{\bar{\beta}^{2}}\leq\left(\frac{2}{1/2}\right)^{2}\bar{\lambda}$ when $\bar{\beta}\geq\frac{1}{2}$. Combining the bounds in both cases, we have

[TABLE]

where the second inequality is due to the upper bounds of $\bar{\lambda}$ and $A_{k}$ in (4.9) and (1.2) respectively. $\Box$

**Proof of Lemma 4.5. ** Let $\bar{v}:=\bar{u}_{+}-\nabla f_{x_{+}}(\bar{y}_{+})+\nabla f_{x_{-}}(\bar{y}_{+})$. Then $\bar{v}\in(\nabla f_{x_{-}}+\partial h)(\bar{y}_{+})$. By (4.13) and (4.14), it holds that

[TABLE]

Recall that $\bar{v}\in(\nabla f_{x_{-}}+\partial h)(\bar{y}_{+})$ and $\bar{u}_{-}\in(\nabla f_{x_{-}}+\partial h)(\bar{y}_{-})$, and from the convexity of $f_{x_{-}}+h$ it holds that

[TABLE]

Therefore,

[TABLE]

Using the previous identity and the triangle inequality of the norms implies that

[TABLE]

$\Box$

Next we shall present the lemmas with proofs that were used in Section 4.

Lemma A.2

Suppose $\lambda_{+}$, $\lambda_{-}$, $\beta_{+}$ and $\beta_{-}$ are generated from Algorithm 3. When the nubmer of iteration $j$ in Algorithm 3 satisfying $j\geq\log_{2}(2/(1-\bar{\beta}))$ with $\bar{\beta}=\beta(\bar{\lambda})$ and $\bar{\lambda}$ defined in (4.9), we have

[TABLE]

Proof. Since $j\geq\log_{2}(2/(1-\bar{\beta}))$, inequality (1.10) holds. By the mean-value theorem and the definition of $\lambda(\beta)$, there exists $\eta\in(\beta_{-},\beta_{+})$ such that

[TABLE]

where the inequality is due to (1.10). Recall that

[TABLE]

The relation of $\beta$ and $\lambda$ in (4.1) gives

[TABLE]

Therefore, by invoking (4.9), (1.1) and (1.2), we have

[TABLE]

$\Box$

The following lemma is exactly Proposition 4.5 in [23].

Lemma A.3

Let $A:\mathbb{R}^{s}\rightrightarrows\mathbb{R}^{s}$ be a maximal monotone operator. Then for any $x,\tilde{x}\in\mathbb{R}^{s}$, we have

[TABLE]

Moreover, if $x_{*}\in A^{-1}(0)$ then

[TABLE]

Now we can bound the residual in terms of the distance between current iterate and an optimal solution.

Lemma A.4

Let $T:=\nabla f+\partial h$ and $T_{x}:=\nabla f_{x}+\partial h$. Assume that $x_{*}\in T^{-1}(0)=(\nabla f+\partial h)^{-1}(0)$ and let $\bar{x},x\in\mathbb{R}^{n}$ be given. Then,

[TABLE]

As a consequence, for every $x\in\mathbb{R}^{n}$, $x_{*}\in T^{-1}(0)$, and $\lambda>0$, it holds that

[TABLE]

Proof. Let $r$ be a constant mapping such that $r(x)=\nabla f(x_{*})-\nabla f_{\bar{x}}(x_{*})$ for any $x\in\mathbb{R}^{n}$. Then, construct $A:=T_{\bar{x}}+r$, where $A$ is also a maximal monotone operator. By Lemma A.3,

[TABLE]

Let $y=x+\lambda(\nabla f(x_{*})-\nabla f_{\bar{x}}(x_{*}))$ and $z=\left(I+\lambda r+\lambda T_{\bar{x}}\right)^{-1}(y)$. We have

[TABLE]

Canceling $\lambda(\nabla f(x_{*})-\nabla f_{\bar{x}}(x_{*}))$ on both sides leads to

[TABLE]

Combining the above inequality with (1.16) and Lemma 2.1 we have

[TABLE]

which proves (1.14), and (1.15) follows from (1.14) straightforwardly. $\Box$

Lemma A.5

Suppose $x_{+}=x_{\beta_{+}}$ and $x_{-}=x_{\beta_{-}}$ are generated from Algorithm 3, and $\bar{y}_{+}$, $\bar{y}_{-}$ are defined in (4.12). When the number of iterations $j$ in Algorithm 3 satisfies $j\geq\log_{2}(2/(1-\bar{\beta}))$ with $\bar{\beta}=\beta(\bar{\lambda})$ and $\bar{\lambda}$ defined in (4.9), we have

[TABLE]

and

[TABLE]

Proof. Let $x_{*}$ be the projection of $x_{0}$ onto the optimal solution set $X_{*}$. According to Lemma 3.1, it holds that

[TABLE]

By (4.2), we have $x_{+}=(1-\beta_{+})y_{k}+\beta_{+}x_{k}$ and $x_{-}=(1-\beta_{-})y_{k}+\beta_{-}x_{k}$. Therefore,

[TABLE]

and

[TABLE]

Recall that in the proof of Theorem 4.1, we showed that when $j\geq\log_{2}(2/(1-\bar{\beta}))$, $\beta_{-}\leq\bar{\beta}$ and inequality (4.10) holds. Applying Lemma A.4 with $x=\bar{x}=x_{+}$, inquality (4.10) and (4.5), we have

[TABLE]

Since $\lambda(\beta)$ is monotonically increasing in $\beta$, $\beta_{-}\leq\bar{\beta}$ amounts to

[TABLE]

Finally, applying Lemma A.4 again with $x=\bar{x}=x_{-}$ and using (4.5) yields

[TABLE]

$\Box$

Lemma A.6

Suppose that $x_{+}=x_{\beta_{+}}$ and $x_{-}=x_{\beta_{-}}$ are generated by Algorithm 3, and $\bar{y}_{+}$, $\bar{y}_{-}$ are defined in (4.12). If the number of iterations $j$ in Algorithm 3 satisfies $j\geq\log_{2}(2/(1-\bar{\beta}))$ with $\bar{\beta}=\beta(\bar{\lambda})$ and $\bar{\lambda}$ defined in (4.9), then we have

[TABLE]

Proof. According to the definition of function $f_{x}(\cdot)$ in (2.4), it holds that

[TABLE]

Note that (1.19) can be further bounded as follows:

[TABLE]

where the first inequality is due to (4.16), and the second last inequality is from Lemma A.5, and that $\|\bar{y}_{+}-x_{-}\|\leq\|\bar{y}_{+}-x_{+}\|+\|x_{+}-x_{-}\|$.

It remains to bound (1.19). We first show by induction that for $1\leq\ell\leq d$ and any convex combination of $x_{-},x_{+}$ and $x_{*}$ denoted by $z$,

[TABLE]

Our induction works backwardly starting from the base case: $\ell=d$. Recall that $z$ is a convex combination of $x_{-},x_{+}$ and $x_{*}$. By (1.17) we have

[TABLE]

Therefore,

[TABLE]

Moreover, by invoking (2.2) and Lemma A.5, we have

[TABLE]

Now suppose that the conclusion holds for some $\ell+1$. Consider

[TABLE]

Denote D_{1}:=\Big{(}2+\frac{2}{\sqrt{1-\sigma^{2}}}\Big{)}D. By letting $x_{t}=\frac{t_{1}}{t_{1}+t_{2}}x_{-}+\frac{t_{2}}{t_{1}+t_{2}}x_{+}$ and (1.17), we have \|x_{t}-x_{*}\|\leq\|x_{-}-x_{*}\|+\|x_{+}-x_{*}\|\leq\Big{(}2+\frac{2}{\sqrt{1-\sigma^{2}}}\Big{)}D=D_{1}. Consequently,

[TABLE]

where the second last inequality is due to (2.1) and the last inequality follows from the induction hypothesis on (1.21). Then, it follows that

[TABLE]

Now by induction on (1.21), applying Lemma A.5 and using (2.1) we have

[TABLE]

Therefore, by induction it follows that (1.21) and (1.22) hold for any $1\leq\ell\leq d$.

Now we come back to bound (1.19). For $2\leq\ell\leq d$,

[TABLE]

Applying Lemma A.5 and (1.21), the first term on the right hand side of (1.23) can be further upper bounded as follows:

[TABLE]

Moreover, applying Lemma A.5, (1.21) and (2.1) to the second term on the right hand side of (1.23) gives that

[TABLE]

Putting (1.23), (1.24) and (1.25) together yields

[TABLE]

for $\ell=2,...,d$. When $\ell=1$, (1.22) guarantees that

[TABLE]

Therefore, the quantity in (1.19) can be bounded as

[TABLE]

Finally, replacing (1.19) and (1.19) with (1.20) and (1.26) respectively leads to the desired conclusion. $\Box$