Linearization Algorithms for Fully Composite Optimization

TL;DR

This paper introduces new first-order algorithms for fully composite optimization that efficiently handle differentiable and non-differentiable components, extending classical methods with improved convergence and practical implementation.

Contribution

It generalizes classical Frank-Wolfe and Conditional Gradient methods for non-differentiable problems using a stronger linear minimization oracle and provides convergence analysis and acceleration.

Findings

Global convergence rates for convex and non-convex objectives

An affine-invariant analysis of the basic method

An accelerated convex optimization algorithm with improved complexity

Abstract

This paper studies first-order algorithms for solving fully composite optimization problems over convex and compact sets. We leverage the structure of the objective by handling its differentiable and non-differentiable components separately, linearizing only the smooth parts. This provides us with new generalizations of the classical Frank-Wolfe method and the Conditional Gradient Sliding algorithm, that cater to a subclass of non-differentiable problems. Our algorithms rely on a stronger version of the linear minimization oracle, which can be efficiently implemented in several practical applications. We provide the basic version of our method with an affine-invariant analysis and prove global convergence rates for both convex and non-convex objectives. Furthermore, in the convex case, we propose an accelerated method with correspondingly improved complexity. Finally, we provide illustrative experiments to support our theoretical results.

Tables

Table 1. Table 1: Summary of convergence complexities for solving non-smooth problems. Note (1) (1) \color[rgb]{0.0,0.47,0.44}{{}^{\text{(1)}}} marks complexities reaching an ε 𝜀 \varepsilon functional residual. Note (2) (2) \color[rgb]{0.0,0.47,0.44}{{}^{\text{(2)}}} marks complexities for reaching Clarke-stationary points. Note (3) (3) \color[rgb]{0.0,0.47,0.44}{{}^{\text{(3)}}} marks complexities for obtaining d − limit-from 𝑑 d- stationary points. Note (4) (4) \color[rgb]{0.0,0.47,0.44}{{}^{\text{(4)}}} marks the complexity for reaching a small norm of the gradient mapping. Finally, note (5) (5) \color[rgb]{0.0,0.47,0.44}{{}^{\text{(5)}}} marks the complexity of minimizing the positive quantity ( 14 ).

Reference	$\phi$ subclass	Use structure?	# FO	# PO/LMO	Observations
Shor et al. (1985)	cvx, L-cont	no	$𝒪\Leftarrow {\epsilon }^{\nwarrow \in }\Rightarrow$${}^{\text{(1)}}$	$𝒪\Leftarrow {\epsilon }^{\nwarrow \in }\Rightarrow$${}^{\text{(1)}}$	projection
Thekumparampil et al. (2020)	cvx, L-cont	no	$𝒪\Leftarrow {\epsilon }^{\nwarrow \in }\Rightarrow$${}^{\text{(1)}}$	$𝒪\Leftarrow {\epsilon }^{\nwarrow \in }\Rightarrow$${}^{\text{(1)}}$	smoothing, vanilla LMO
Doikov and Nesterov (2022)	cvx, fully-comp	yes	$𝒪\Leftarrow {\epsilon }^{\nwarrow \infty }\Rightarrow$${}^{\text{(1)}}$	$𝒪\Leftarrow {\epsilon }^{\nwarrow \infty }\Rightarrow$${}^{\text{(1)}}$	modif. LMO
(this work) Alg. 2	cvx, fully-comp.	yes	$𝒪\Leftarrow {\epsilon }^{\nwarrow \infty \propto \in }\Rightarrow$${}^{\text{(1)}}$	$𝒪\Leftarrow {\epsilon }^{\nwarrow \infty }\Rightarrow$${}^{\text{(1)}}$	modif. LMO
De Oliveira (2023)	non-cvx, upper-$C^{1,\alpha }$	no	$𝒪\Leftarrow {\epsilon }^{\nwarrow \in }\Rightarrow$ ${}^{\text{(2)}}$	$𝒪\Leftarrow {\epsilon }^{\nwarrow \in }\Rightarrow$ ${}^{\text{(2)}}$	vanilla LMO
Kreimeier et al. (2023)	non-cvx, abs-smooth	no	$𝒪\Leftarrow {\epsilon }^{\nwarrow \in }\Rightarrow$ ${}^{\text{(3)}}$	$𝒪\Leftarrow {\epsilon }^{\nwarrow \in }\Rightarrow$ ${}^{\text{(3)}}$	modif. LMO
Drusvyatskiy and Paquette (2019)	non-cvx, comp	yes	$𝒪\Leftarrow {\epsilon }^{\nwarrow \in }\Rightarrow$ ${}^{\text{(4)}}$	$𝒪\Leftarrow {\epsilon }^{\nwarrow \in }\Rightarrow$ ${}^{\text{(4)}}$	prox. steps
(this work) Alg. 1	non-cvx, fully-comp	yes	$\widetilde{𝒪}\left({\epsilon }^{-2}\right)$ ${}^{\text{(5)}}$	$\widetilde{𝒪}\left({\epsilon }^{-2}\right)$ ${}^{\text{(5)}}$	modif. LMO

Equations

$$\underset{\mathbf{x}\in \text{mathcalorigin}X}{\min }\left[\phi (\mathbf{x})\stackrel{\mathrm{def}}{=}F(\mathbf{f}(\mathbf{x}),\mathbf{x})\right],$$

$$\begin{array}[]{rcl}\mathbf{y}_{k+1}&=&\pi_{\mathcalorigin{X}}\bigl{(}\mathbf{y}_{k}-\alpha_{k}\nabla\varphi(\mathbf{y}_{k})\bigr{)},\qquad\alpha_{k}>0,\end{array}$$

$$\begin{array}[]{rcl}\mathbf{y}_{k+1}&\in&\mathop{\rm Argmin}\limits_{\mathbf{x}}\bigl{\{}\,\langle\nabla\varphi(\mathbf{y}_{k}),\mathbf{x}\rangle\;:\;\mathbf{x}\in\mathbf{y}_{k}+\gamma_{k}(\mathcalorigin{X}-\mathbf{y}_{k})\,\bigr{\}},\qquad\gamma_{k}\in(0,1].\end{array}$$

$$\begin{array}[]{rcl}\varphi^{\star}&=&\min\limits_{\mathbf{x}\in\mathcalorigin{X}}\Bigl{[}\,\varphi(\mathbf{x})\;\;\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\;\;F(\mathbf{f}(\mathbf{x}),\mathbf{x})\,\Bigr{]},\qquad\mathcalorigin{X}\subset\mathbb{R}^{d},\end{array}$$

$$\begin{array}[]{rcl}\mathop{\rm Argmin}\limits_{\mathbf{x}\in\mathcalorigin{X}}F\bigl{(}\mathbf{A}\mathbf{x}+\mathbf{b},\mathbf{x}\bigr{)}+\langle\mathbf{u},\mathbf{x}\rangle\end{array}$$

$$\begin{array}[]{rcl}F(\gamma\mathbf{u},\mathbf{x})&\leq&\gamma F(\mathbf{u},\mathbf{x}),\qquad\forall\mathbf{u}\in\mathbb{R}^{n},\;\mathbf{x}\in\mathcalorigin{X},\;\gamma\geq 1.\end{array}$$

$$\begin{array}[]{rcl}\mathcal{S}&\!\!=\!\!&\mathcal{S}_{\mathbf{f},F,\mathcalorigin{X}}\;\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\par\sup\limits_{\begin{subarray}{c}\mathbf{x},\mathbf{y}\in\mathcalorigin{X},\,\gamma\in(0,1]\\ \mathbf{y}_{\gamma}=\mathbf{x}+\gamma(\mathbf{y}-\mathbf{x})\end{subarray}}F\bigl{(}\frac{2}{\gamma^{2}}\bigl{[}\mathbf{f}(\mathbf{y}_{\gamma})-\mathbf{f}(\mathbf{x})-\nabla\mathbf{f}(\mathbf{x})(\mathbf{y}_{\gamma}-\mathbf{x})\bigr{]},\mathbf{y}_{\gamma}\bigr{)}\;<\;+\infty.\end{array}$$

$$\Vert \nabla f_i(\mathbf{x})-\nabla f_i(\mathbf{y})\Vert \le L_i\Vert \mathbf{x}-\mathbf{y}\Vert \forall \mathbf{x},\mathbf{y}\in \text{mathcalorigin}X,\forall i\in [n].$$

$$\begin{array}[]{rcl}F(\mathbf{u},\mathbf{x})&\leq&F(\mathbf{v},\mathbf{x}).\end{array}$$

$$\begin{array}[]{rcl}\mathcal{S}&\leq&\sup\limits_{\begin{subarray}{c}\mathbf{x},\mathbf{y}\in\mathcalorigin{X},\,\gamma\in[0,1]\\ \mathbf{y}_{\gamma}=\mathbf{x}+\gamma(\mathbf{y}-\mathbf{x})\end{subarray}}F(\nabla^{2}\mathbf{f}(\mathbf{y}_{\gamma})[\mathbf{y}-\mathbf{x}]^{2},\mathbf{x}).\end{array}$$

$$\begin{array}[]{rcl}\mathcal{S}&\leq&F(\mathbf{L}\mathcalorigin{D}_{\mathcalorigin{X}}^{2})\;\;\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\;\;\sup\limits_{\mathbf{x}\in\mathcalorigin{X}}F(\mathbf{L}\mathcalorigin{D}_{\mathcalorigin{X}}^{2},\mathbf{x}).\end{array}$$

$$\begin{array}[]{c}\min\limits_{\mathbf{x}\in\mathcalorigin{X}}\max\limits_{1\leq i\leq n}f_{i}(\mathbf{x}),\end{array}$$

$$\begin{array}[]{rcl}\min\limits_{\mathbf{x}\in\mathcalorigin{X}}\max\limits_{1\leq i\leq n}\langle\mathbf{a}_{i},\mathbf{x}\rangle+b_{i}&\quad\Leftrightarrow&\min\limits_{\mathbf{x}\in\mathcalorigin{X},t\in\mathbb{R}}\bigl{\{}t\;:\;\langle\mathbf{a}_{i},\mathbf{x}\rangle+b_{i}\leq t,\,1\leq i\leq n\bigr{\}}.\end{array}$$

$$\begin{array}[]{rcl}\min\limits_{\mathbf{x}\in\mathcalorigin{X}}\max\limits_{1\leq i\leq n}\langle\mathbf{a}_{i},\mathbf{x}\rangle+b_{i}&=&\min\limits_{\mathbf{x}\in\mathcalorigin{X}}\max\limits_{\bm{\lambda}\in\Delta_{n}}\sum\limits_{i=1}^{n}\lambda^{(i)}\bigl{[}\langle\mathbf{a}_{i},\mathbf{x}\rangle+b_{i}\bigr{]}=\max\limits_{\bm{\lambda}\in\Delta_{n}}g(\bm{\lambda}),\end{array}$$

$$\begin{array}[]{rcl}\min\limits_{\mathbf{x}\in\mathcalorigin{X}}\|\mathbf{f}(\mathbf{x})\|,\end{array}$$

$${\mathbf{x}}_{k+1}\in \underset{\mathbf{x}\in \text{mathcalorigin}X}{\mathrm{Argmin}}\Vert \mathbf{f}({\mathbf{y}}_k)+\nabla \mathbf{f}({\mathbf{y}}_k)(\mathbf{x}-{\mathbf{y}}_k)\Vert ,\text{ and }{\mathbf{y}}_{k+1}=(1-{\gamma }_k){\mathbf{y}}_k+{\gamma }_k{\mathbf{x}}_{k+1}.$$

$$\begin{array}[]{rcl}\mathbf{x}_{k+1}&\in&\mathop{\rm Argmin}\limits_{\mathbf{x}\in\mathcalorigin{X}}F\bigl{(}\mathbf{f}(\mathbf{y}_{k})+\nabla\mathbf{f}(\mathbf{y}_{k})(\mathbf{x}-\mathbf{y}_{k}),\,\mathbf{x}\bigr{)}\end{array}$$

$$\Delta_{k}\;\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\;\varphi(\mathbf{y}_{k})-F\big{(}\mathbf{f}(\mathbf{y}_{k})+\nabla\mathbf{f}(\mathbf{y}_{k})(\mathbf{x}_{k+1}-\mathbf{y}_{k}),\mathbf{x}_{k+1}\big{)}.$$

$$\begin{array}[]{rcl}\Delta_{k}&\geq&\max\limits_{\mathbf{x}\in\mathcalorigin{X}}\Bigl{[}\varphi(\mathbf{y}_{k})-F\big{(}\mathbf{f}(\mathbf{y}_{k})+\nabla\mathbf{f}(\mathbf{y}_{k})(\mathbf{x}-\mathbf{y}_{k}),\mathbf{x}\big{)}\Bigr{]}\\[10.0pt] &\geq&\max\limits_{\mathbf{x}\in\mathcalorigin{X}}\Bigl{[}\varphi(\mathbf{y}_{k})-F(\mathbf{f}(\mathbf{x}),\mathbf{x})\Bigr{]}\;\;=\;\;\varphi(\mathbf{y}_{k})-\varphi^{\star}.\end{array}$$

$$\phi ({\mathbf{y}}_k)-{\phi }^{\star }\le \frac{2\mathcal{S}}{1+k}\text{ and }\underset{1\le i\le k}{\min }{\Delta }_i\le \frac{6\mathcal{S}}{k}.$$

$$\begin{array}[]{rcl}\min\limits_{0\leq i\leq k}\Delta_{i}&\leq&\frac{\varphi(\mathbf{y}_{0})-\varphi^{\star}+0.5\mathcal{S}(1+\ln(k+1))}{\sqrt{k+1}}.\end{array}$$

$$\begin{array}[]{rcl}\mathbf{y}_{k+1}&\in&\mathop{\rm Argmin}\limits_{\mathbf{y}\in\mathbf{y}_{k}+\gamma_{k}(\mathcalorigin{X}-\mathbf{y}_{k})}\|\mathbf{f}(\mathbf{y}_{k})+\nabla\mathbf{f}(\mathbf{y}_{k})(\mathbf{y}-\mathbf{y}_{k})\|.\end{array}$$

$$\begin{array}[]{rcl}\min\limits_{0\leq i\leq k}\max\limits_{\mathbf{y}\in\mathcalorigin{X}}\langle\nabla\Phi(\mathbf{y}_{i}),\mathbf{y}_{i}-\mathbf{y}\rangle&\leq&\mathcal{O}\bigl{(}\frac{\ln(k)}{\sqrt{k}}\bigr{)}.\end{array}$$

$$\begin{array}[]{rcl}\mathop{\rm Argmin}\limits_{\mathbf{u}\in\mathcalorigin{X}}\Bigl{\{}P(\mathbf{u})\stackrel{{\scriptstyle\mathrm{def}}}{{=}}F(\mathbf{f}(\mathbf{z})+\nabla\mathbf{f}(\mathbf{z})(\mathbf{u}-\mathbf{z}),\mathbf{u})+\frac{\beta}{2}\|\mathbf{u}-\mathbf{x}\|_{2}^{2},\;\;\beta>0\Bigr{\}}.\end{array}$$

$$F(\mathbf{f}(\mathbf{z})+\nabla f(\mathbf{z})({\mathbf{u}}^{+}-\mathbf{z}),{\mathbf{u}}^{+})+\beta ⟨{\mathbf{u}}^{+}-\mathbf{x},{\mathbf{u}}^{+}⟩$$

$$\le F(\mathbf{f}(\mathbf{z})+\nabla f(\mathbf{z})(\mathbf{u}-\mathbf{z}),\mathbf{u})+\beta ⟨{\mathbf{u}}^{+}-\mathbf{x},\mathbf{u}⟩+\eta ,\forall \mathbf{u}\in \text{mathcalorigin}X.$$

$$\begin{array}[]{rcl}\varphi(\mathbf{y}_{k})-\varphi^{\star}&\leq&\frac{\delta+8cF(\mathbf{L})\mathcalorigin{D}_{\mathcalorigin{X}}^{2}}{(k+2)(k+3)}+\frac{2\max\{0,1-c\}F(\mathbf{L})\mathcalorigin{D}_{\mathcalorigin{X}}^{2}}{k+3}.\end{array}$$

$$\begin{array}[]{c}P(\mathbf{u}_{t})-P^{\star}\;\;\leq\;\;\frac{2\beta\mathcalorigin{D}_{\mathcalorigin{X}}^{2}}{t+1}\qquad\text{ and }\qquad\min\limits_{1\leq i\leq t}\Delta_{t}\;\;\leq\;\;\frac{6\beta\mathcalorigin{D}_{\mathcalorigin{X}}^{2}}{t}.\end{array}$$

$$\underset{\mathbf{x}\in \text{mathcalorigin}X}{\min }\left\{\underset{i=1,n}{\max }{\mathbf{x}}^{\top }{\mathbf{A}}_i\mathbf{x}-{\mathbf{b}}_i^{\top }\mathbf{x}\right\},\text{ for }\text{mathcalorigin}X={\Delta }_d,\subseteq {\mathbb{R}}^d,$$

$$\underset{\mathbf{X}\in \text{mathcalorigin}X}{\min }\left\{\underset{i=1,n}{\max }\sum _{(k,l)\in {\Omega }_i}{\left({\mathbf{X}}_{k,l}-{\mathbf{A}}_{k,l}^{(i)}\right)}^2\right\},\text{ for }\text{mathcalorigin}X:=\{\mathbf{X}\in {\mathbb{R}}^{d\times m},{\Vert \mathbf{X}\Vert }_{\ast }\le r\}$$

$$\begin{array}[]{rcl}\varphi(\mathbf{y}_{\gamma})&\leq&F\bigl{(}\mathbf{f}(x)+\nabla\mathbf{f}(\mathbf{x})(\mathbf{y}_{\gamma}-\mathbf{x}),\,\mathbf{y}_{\gamma}\bigr{)}\;+\;\frac{\gamma^{2}}{2}\mathcal{S}.\end{array}$$

Full text

1 Introduction

In this paper we consider fully composite optimization problems of the form

[TABLE]

where $\mathcalorigin{X}$ is a convex and compact set, $F:\mathbb{R}^{n}\times\mathcalorigin{X}\to\mathbb{R}$ is a simple but possibly non-differentiable convex function and $\mathbf{f}:\mathcalorigin{X}\to\mathbb{R}^{n}$ is a smooth mapping, which is the main source of computational burden.

Problems of this type cover and generalize many classical use-cases of composite optimization and are often encountered in applications. In this work, we develop efficient algorithms for solving (1) by leveraging the structure of the objective and using the linearization principle. Our method generalizes the well-known Frank-Wolfe algorithm (Frank and Wolfe, 1956) and ensures asymptotically faster convergence rates compared to methods treating $\varphi$ in a black-box fashion.

A classical algorithm for solving smooth optimization problems is the Gradient Descent method (GD), proposed by Cauchy in 1847 (see historical note by Lemaréchal (2012)). It rests on the idea of linearizing the function around the current iterate, taking a step in the negative gradient direction and projecting the result onto the feasible set $\mathcalorigin{X}$ for $k\geq 0$:

[TABLE]

where $\pi_{\mathcalorigin{X}}$ is the projection operator onto $\mathcalorigin{X}$. Surprisingly, the same kind of iterations can minimize general non-smooth convex functions by substituting $\nabla\varphi(\mathbf{y}_{k})$ with any subgradient in the subdifferential $\partial\varphi(\mathbf{y}_{k})$. The resulting Subgradient method was proposed by Shor et al. (1985).

Another notable example for smooth optimization over a convex and bounded constraint set $\mathcalorigin{X}$ is the Frank-Wolfe (FW) method (Frank and Wolfe, 1956). Again, a linearization of the objective around the current iterate is used to query the so-called linear minimization oracle (LMO) associated with $\mathcalorigin{X}$, for every $k\geq 0$:

[TABLE]

Steps of type (3) can be significantly cheaper than those involving projections (2) for a few notable domains such as nuclear norm balls and spectrahedrons (Combettes and Pokutta, 2021), making FW the algorithm of choice in such scenarios. Moreover, the solutions found by FW methods can benefit from additional properties such as sparsity (Jaggi, 2013). These desirable features make FW methods suitable for large scale optimization, a fact which led to an increased interest in recent years (we point the reader to the monograph of Braun et al. (2022) for a detailed presentation). Unfortunately, the vanilla FW algorithm does not extend to non-differentiable problems in the same straightforward manner as GD – a counterexample is given by Nesterov (2018a). The question of developing non-smooth versions of the FW algorithm therefore remains open, and is the main focus of this article.

Finally, we touch on the issue of convergence rates – a principal means of theoretically characterizing optimization algorithms. The classical monograph of Nemirovski and Yudin (1983) establishes that the $\mathcal{O}(1/\sqrt{k})$ rate of the Subgradient method is optimal for general non-differentiable convex problems, while the $\mathcal{O}(1/k)$ rate of its counterpart GD is far from the lower bound of $\Omega(1/k^{2})$ for $L-$smooth convex functions. Similar results are established by Lan (2013) for LMO-based algorithms, although in this case the $\mathcal{O}(1/k)$ rate is matched by a lower bound for smooth convex minimization. This relatively slow convergence of FW algorithms is a result of their affine-invariant oracle, which is independent on the choice of norm. In light of these lower bounds, one can only hope to improve convergence rates by imposing additional structure on the problem to be solved.

The present work leverages this observation and studies a subclass of (possibly) non-smooth and non-convex problems with the specific structure of (1). Our methods require only linearizations of the differentiable component $\mathbf{f}$, while the non-differentiable function $F$ is kept as a part of the subproblem solved within oracle calls. We show that this approach is a viable way of generalizing FW methods to address problem (1), with the possibility of acceleration. Our contributions can be summarized as follows.

• •

We propose a basic method for problem (1), which is affine-invariant and equipped with accuracy certificates. We prove the global convergence rate of $\mathcal{O}(1/k)$ in the convex setting, and of $\tilde{\mathcal{O}}(1/\sqrt{k})$ in the non-convex case.

• •

We propose an accelerated method with inexact proximal steps which attains a convergence rate of $\mathcal{O}(1/k^{2})$ for convex problems. Our algorithm achieves the optimal $\mathcalorigin{O}\left(\varepsilon^{-1/2}\right)$ oracle complexity for smooth convex problems in terms of the number of computations of $\nabla\mathbf{f}$.

• •

We provide proof-of-concept numerical experiments, that demonstrate the efficiency of our approach for solving composite problems.

Related Work.

The present work lies at the intersection of two broad lines of study: general methods for composite optimization and FW algorithms. The former category encompasses many approaches that single out non-differentiable components in the objective’s structure, and leverage this knowledge in the design of efficient optimization algorithms. This approach originated in the works of Burke (1985, 1987); Nesterov (1989); Nemirovski (1995); Pennanen (1999); Boţ et al. (2007, 2008). A popular class of additive composite optimization problems was proposed by Beck and Teboulle (2009); Nesterov (2013) and the modern algorithms for general composite formulations were developed by Cui et al. (2018); Drusvyatskiy and Lewis (2018); Drusvyatskiy and Paquette (2019); Bolte et al. (2020); Burke et al. (2021); Doikov and Nesterov (2022).

The primitive on which most of the aforementioned methods rely is a proximal-type step – a generalization of (2). Depending on the geometry of the set $\mathcalorigin{X}$, such steps may pose a significant computational burden. Doikov and Nesterov (2022) propose an alternative contracting-type method for fully composite problems, which generalizes the vanilla FW algorithm. Their method relies on a simpler primitive built on the linearization principle, which can be much cheaper in practice. We study the same problem structure as Doikov and Nesterov (2022) and devise methods with several advantages over the aforementioned approach, including an affine-invariant analysis, accuracy certificates, convergence guarantees for non-convex problems and, in the convex case, an accelerated convergence. Moreover, we decouple stepsize selection from the computational primitive, to enable efficient line search procedures.

Our methods are also intimately related to FW algorithms, which they generalize. For smooth and convex problems, vanilla FW converges at the cost of $\mathcalorigin{O}\left(\varepsilon^{-1}\right)$ LMO and first order oracle (FO) calls in terms the Frank-Wolfe gap – an accuracy measure bounding functional suboptimality (Jaggi, 2013). For smooth non-convex problems, a gap value of at most $\varepsilon$ is attained after $\mathcalorigin{O}\left(\varepsilon^{-2}\right)$ LMO and FO calls (Lacoste-Julien, 2016). Due to the relatively slow convergence of LMO-based methods, recent efforts have gone into devising variants with improved guarantees. The number of FO calls was reduced to the lower bound for smooth convex optimization by Lan and Zhou (2016), local acceleration was achieved following a burn-in phase by Diakonikolas et al. (2020); Carderera et al. (2021); Chen and Sun (2022), and empirical performance was enhanced by adjusting the update direction with gradient information by Combettes and Pokutta (2020). Of the aforementioned works, closest to ours is the Conditional Gradient Sliding (CGS) algorithm proposed by Lan and Zhou (2016) and further studied by Yurtsever et al. (2019); Qu et al. (2018). CGS uses the acceleration framework of Nesterov (1983) and solves the projection subproblem inexactly via the FW method, achieving the optimal complexity of $\mathcalorigin{O}\left(\varepsilon^{-1/2}\right)$ FO calls for smooth convex problems. We rely on a similar scheme for improving FO complexity in the convex case.

In the context of generic non-smooth convex objectives, the FW algorithm was studied by Lan (2013), who proposes a smoothing-based approach matching the lower bound of $\Omega(\varepsilon^{-2})$ LMO calls. The method however requires $\mathcalorigin{O}\left(\varepsilon^{-4}\right)$ FO calls, a complexity which is later improved to $\mathcalorigin{O}\left(\varepsilon^{-2}\right)$ by Garber and Hazan (2016) through a modified LMO for polytopes, by Ravi et al. (2019) with a (differently) modified LMO, and finally by Thekumparampil et al. (2020) through a combination of smoothing and the CGS algorithm. Our algorithm, instead, leverages the structure of problem (1) and a modified LMO to achieve improved rates, with the added benefit of an affine invariant method and analysis. We also mention FW methods for additive composite optimization (Argyriou et al., 2014; Yurtsever et al., 2018, 2019; Zhao and Freund, 2022), with the former three relying on proximal steps and the latter assuming a very restricted class of objectives.

Finally, two concurrent works study FW methods for some restricted classes of non-smooth and non-convex problems. De Oliveira (2023) shows that vanilla FW with line-search can be applied to the special class of upper$-C^{1,\alpha}$ functions, when one replaces gradients with an arbitrary element in the Clarke subdifferential. A rate of $\mathcalorigin{O}\left(\varepsilon^{-2}\right)$ is shown for reaching a Clarke-stationary point in a setting comparable to ours. A similar rate is shown by Kreimeier et al. (2023) for reaching a $d$-stationary point of abs-smooth functions through the use of a modified LMO. Both these algorithms are structure-agnostic. A summary of method complexities for solving non-smooth problems is provided in Table 1.

Notation.

We denote by $[n]$ the set $\{1,\ldots n\}$ and by $\left\lVert\,\cdot\,\right\rVert$ the standard Euclidean norm, unless explicitly stated otherwise. We define the diameter of a bounded set $\mathcalorigin{X}$ as $\mathcalorigin{D}_{\mathcalorigin{X}}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\max_{\mathbf{z},\mathbf{y}\in\mathcalorigin{X}}\{\left\lVert\,\mathbf{z}-\mathbf{y}\,\right\rVert\}$. We use the notation $\Delta_{n}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\{\bm{\lambda}\in\mathbb{R}^{n}_{+}\,:\,\langle\bm{\lambda},\mathbf{e}\rangle=1\}$ to denote the standard $n$-dimensional simplex, where $\mathbf{e}$ is the vector of all ones. For a differentiable, scalar-valued function $f:\mathbb{R}^{d}\to\mathbb{R}$ we use $\nabla f(\mathbf{x})\in\mathbb{R}^{d}$ to denote its gradient vector and $\nabla^{2}f(\mathbf{x})\in\mathbb{R}^{d\times d}$ to denote its Hessian matrix. For a differentiable vector-valued function $\mathbf{f}:\mathbb{R}^{d}\to\mathbb{R}^{n}$ defined as $\mathbf{f}=(f_{1},f_{2},\ldots f_{n})$ we denote by $\nabla\mathbf{f}(\mathbf{x})$ its Jacobian matrix defined as $\nabla\mathbf{f}(\mathbf{x})=\sum_{i=1}^{n}\mathbf{e}_{i}\nabla f_{i}(\mathbf{x})^{\top}\in\mathbb{R}^{n\times d},$ where $\mathbf{e}_{i}$ are the standard basis vectors in $\mathbb{R}^{n}$. We represent the second directional derivatives applied to the same direction $\mathbf{h}\in\mathbb{R}^{d}$ as $\nabla^{2}f(\mathbf{x})[\mathbf{h}]^{2}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\langle\nabla^{2}f(\mathbf{x})\mathbf{h},\mathbf{h}\rangle\in\mathbb{R}$, and $\nabla^{2}\mathbf{f}(\mathbf{x})[\mathbf{h}]^{2}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\sum_{i=1}^{n}\mathbf{e}_{i}\nabla^{2}f_{i}(\mathbf{x})[\mathbf{h}]^{2}\in\mathbb{R}^{n}$.

2 Problem Setup, Assumptions and Examples

The problems addressed by this work are represented by the following structured objective

[TABLE]

where $\mathcalorigin{X}$ is a convex and compact set and the inner mapping $\mathbf{f}:\mathcalorigin{X}\to\mathbb{R}^{n}$ is differentiable and defined as $\mathbf{f}(\mathbf{x})=(f_{1}(\mathbf{x}),\ldots,f_{n}(\mathbf{x}))\,\in\,\mathbb{R}^{n},$ where each $f_{i}:\mathcalorigin{X}\to\mathbb{R}$ is differentiable. We assume access to a first-order oracle $\nabla\mathbf{f}$, which is the main source of computational burden. The outer component $F:\mathbb{R}^{n}\times\mathcalorigin{X}\to\mathbb{R}$, on the other hand, is directly accessible to the algorithm designer and is simple (see assumptions). However, $F$ is possibly non-differentiable.

In this work, we propose two algorithmic solutions addressing problem (4), which we call a fully composite problem. Our methods importantly assume that subproblems of the form

[TABLE]

are efficiently solvable, where $\mathbf{A}\in\mathbb{R}^{n\times d}$ and $\mathbf{b}\in\mathbb{R}^{n},\mathbf{u}\in\mathbb{R}^{d}$. Oracles of type  (5) are sequentially called during the optimization procedure and take as arguments linearizations of the difficult nonlinear components of (4). Naturally, solving (5) cheaply is possible only when $F$ is simple and $\mathcalorigin{X}$ has an amenable structure.

In particular, template (4) encompasses to some standard problem formulations. For example, the classical Frank-Wolfe setting is recovered when $F(\mathbf{u},\mathbf{x})\equiv u^{(1)}$, in which case problem (4) becomes $\min_{\mathbf{x}\in\mathcalorigin{X}}f_{1}(\mathbf{x})$ and subproblem (5) reduces to a simple LMO: $\mathop{\rm argmin}_{\mathbf{x}\in\mathcalorigin{X}}\langle\mathbf{u},\mathbf{x}\rangle$. The setting of proximal-gradient methods is similarly covered, by letting $F(\mathbf{u},\mathbf{x})\equiv u^{(1)}+\psi(\mathbf{x})$ for a given convex function $\psi$ (e.g., a regularizer). Then, problem (4) reduces to additive composite optimization \min_{\mathbf{x}\in\mathcalorigin{X}}\Bigl{\{}f_{1}(\mathbf{x})+\psi(\mathbf{x})\Bigr{\}}, and subproblem (5) becomes \mathop{\rm argmin}_{\mathbf{x}\in\mathcalorigin{X}}\Bigl{\{}\langle\mathbf{u},\mathbf{x}\rangle+\psi(\mathbf{x})\Bigr{\}}.

We now formally state the assumptions on the fully composite problem (4), along with commentary and examples.

Assumption 1

The outer function $F:\mathbb{R}^{n}\times\mathcalorigin{X}\to\mathbb{R}$ is jointly convex in its arguments. Additionally, $F(\mathbf{u},\mathbf{x})$ is subhomogeneous in $\mathbf{u}$:

[TABLE]

Assumption 2a

The inner mapping $\mathbf{f}:\mathcalorigin{X}\to\mathbb{R}^{n}$ is differentiable and the following affine-invariant quantity is bounded:

[TABLE]

Assumption 2b

Each component $f_{i}(\cdot)$ has a Lipschitz continuous gradient on $\mathcalorigin{X}$ with constant $L_{i}$:

[TABLE]

We denote the vector of Lipschitz constants by $\mathbf{L}=(L_{1},\ldots,L_{n})\in\mathbb{R}^{n}$.

Assumption 3

Each component $f_{i}:\mathcalorigin{X}\to\mathbb{R}$ is convex. Moreover, $F(\cdot,\mathbf{x})$ is monotone $\forall\mathbf{x}\in\mathcalorigin{X}$. Thus, for any two vectors $\mathbf{u},\mathbf{v}\in\mathbb{R}^{n}$ such that $\mathbf{u}\leq\mathbf{v}$ (component-wise), it holds that

[TABLE]

A few comments are in order. Assumption 1, which is also required by Doikov and Nesterov (2022), represents the formal manner in which we ask that $F$ be simple – through convexity and bounded growth in $\mathbf{u}$. This assumption ensures convexity of subproblem (5), irrespective of the nature of $\mathbf{f}$.

Assumption 2a is a generalization of the standard bounded curvature premise typical for Frank-Wolfe settings (Jaggi, 2013). Requirement (7) is mild, as it only asks that the curvature of $\mathbf{f}$ remains bounded under $F$ over $\mathcalorigin{X}$. Importantly, the quantity $\mathcal{S}$ is affine-invariant (remains unchanged under affine reparametrizations of $\mathcalorigin{X}$), which enables us to obtain convergence rates with the same property. Further discussion on the importance of affine-invariant analysis for FW algorithms is provided by Jaggi (2013). For mappings $\mathbf{f}$ that are twice differentiable, we can bound the quantity $\mathcal{S}$ from Assumption 2a using Taylor’s formula and the second derivatives, as follows

[TABLE]

This quantity is reminiscent of the quadratic upper-bound used to analyze smooth optimization methods. In particular, for monotone non-decreasing $F$, a compact $\mathcalorigin{X}$ and Lipschitz continuous $\nabla f_{i}$ with respect to a fixed norm $\|\cdot\|$, the assumption is satisfied with

[TABLE]

Assumption 2b is standard and considered separately from Assumption 2a to allow for different levels of generality in our results. The restriction to $\mathcalorigin{X}$ makes this a locally-Lipschitz gradient assumption on $f_{i}$.

Finally, Assumption 3, which is also made by Doikov and Nesterov (2022), is required whenever we need to ensure the overall convexity of $\varphi(\mathbf{x})$. The monotonicity of $F$ is necessary in addition to convexity of each $f_{i}$, since the composition of convex functions is not necessarily convex (Boyd and Vandenberghe, 2004). We rely on this assumption for deriving asymptotically faster convergence rates in the convex setting (Section 4).

To conclude this section, we provide the main application examples that fall under our fully composite template and which satisfy to our assumptions. Further examples can be found in Appendix C.

Example 2.1**.**

Let $F(\mathbf{u},\mathbf{x})\equiv\max\limits_{1\leq i\leq n}u^{(i)}$. Function $F$ satisfies Assumptions 1 and 3 and problem (4) becomes

[TABLE]

while oracle (5) becomes

[TABLE]

Max-type minimization problems of this kind result from scalarization approaches in multi-objective optimization, and their solutions were shown to be (weakly) Pareto optimal (Chapter 3.1 in Miettinen, 1999). As such, problem (9) is relevant to a wide variety of applications requiring optimal trade-offs amongst several objective functions, and appears in areas such as machine learning, science and engineering (see the introductory sections of, e.g., Daulton et al., 2022; Zhang and Golovin, 2020). Problem (9) also covers some instances of constrained $\ell_{\infty}$ regression.

When $\mathcalorigin{X}$ is a polyhedron, subproblem (10) can be solved via Linear Programming, while for general $\mathcalorigin{X}$ one can resort to Interior-Point Methods (Nesterov and Nemirovski, 1994). Another option for solving (10) is to note that under strong duality (Rockafellar, 1970) we have

[TABLE]

where g(\bm{\lambda})\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\min\limits_{\mathbf{x}\in\mathcalorigin{X}}\sum_{i=1}^{n}\lambda^{(i)}\bigl{[}\langle\mathbf{a}_{i},\mathbf{x}\rangle+b_{i}\bigr{]}. The maximization of $g$ in (11) can be done very efficiently for small values of $n$ (with, e.g., the Ellipsoid Method or the Mirror Descent algorithm), since evaluating $g(\bm{\lambda})$ and $\partial g(\bm{\lambda})$ reduces to a vanilla LMO call over $\mathcalorigin{X}$. An interesting case is $n=2$, for which (11) becomes a univariate maximization problem and one may use binary search to solve it at the expense of a logarithmic number of LMOs.

Example 2.2**.**

Let $F(\mathbf{u},\mathbf{x})\equiv\|\mathbf{u}\|$ for an arbitrary fixed norm $\|\cdot\|$. Function $F$ satisfies Assumption 1 and problem (4) can be interpreted as solving a system of non-linear equations over $\mathcalorigin{X}$

[TABLE]

while oracle (5) amounts to solving the (constrained) linear system $\min\limits_{\mathbf{x}\in\mathcalorigin{X}}\|\mathbf{A}\mathbf{x}+\mathbf{b}\|$. Problems of this kind can be encountered in applications such as robust phase retrieval (Duchi and Ruan, 2019) with phase constraints.

The iterations of Algorithm 1 can be interpreted as a variant of the Gauss-Newton method (Burke and Ferris, 1995; Nesterov, 2007; Tran-Dinh et al., 2020), solving the (constrained) linear systems:

[TABLE]

In the particular case of solving systems of non-linear equations over compact convex sets, our algorithms can be seen as modified Gauss-Newton methods with global convergence guarantees.

3 The Basic Method

We present the first new method for solving problem (4) in Algorithm 1. The central idea is to linearize the differentiable components of the objective and then to minimize this new model over the constraint $\mathcalorigin{X}$, via calls to an oracle of type (5). The next iterate is defined as a convex combination with coefficient (or stepsize) $\gamma$ between the computed minimizer and the preceding iterate.

A similar method for tackling problems of type (4) in the convex setting was proposed by Doikov and Nesterov (2022). Different from theirs, our method decouples the parameter $\gamma_{k}$ from the minimization subproblem. This change is crucial since it allows us to choose the parameter $\gamma_{k}$ after minimizing the model, thus enabling us to use efficient line search rules. Moreover, we provide Algorithm 1 with a more advanced affine-invariant analysis and establish its convergence in the non-convex setup.

We also mention that for solving problems of type (9), oracle (5) reduces to the minimization of a piecewise linear function over $\mathcalorigin{X}$. Therefore, it has the same complexity as the modified LMOs of Kreimeier et al. (2023) and, moreover, subproblem (5) is convex irrespective of the nature of $\mathbf{f}$.

Accuracy Certificates.

The standard progress metric of FW algorithms, which Algorithm 1 generalizes, is the ‘Frank-Wolfe gap’ or Hearn gap (Hearn, 1982). For smooth objectives, it is defined as $\Delta_{k}=\max_{\mathbf{y}\in\mathcalorigin{X}}\langle\nabla\varphi(\mathbf{y}_{k}),\mathbf{y}_{k}-\mathbf{y}\rangle$, for each iterate $\mathbf{y}_{k}$. This quantity is computed cost-free during the algorithm’s iterations and has the desirable property of upper-bounding the suboptimality of the current iterate: $\Delta_{k}\geq\varphi(\mathbf{y}_{k})-\varphi^{\star}$. Notably, its semantics straightforwardly extend to the non-convex setting (Lacoste-Julien, 2016). Additionally, convergence guarantees on the gap are desirable due to its affine invariance, which aligns with the affine invariance of classical FW algorithm.

Our setting does not permit a direct generalization of the FW gap with all of the above properties. Rather, we introduce the following accuracy certificate, which is readily available in each iteration:

[TABLE]

For minimization of a smooth (not necessarily convex) function, quantity (14) indeed reduces to the standard FW gap. Moreover, for convex $\varphi(\mathbf{x})$ (Assumption 3) we can conclude that

[TABLE]

Hence, for a tolerance $\varepsilon>0$, the criterion $\Delta_{k}\leq\varepsilon$ can be used as the stopping condition for our method in convex scenarios. Moreover, the value of $\Delta_{k}$ can be used for computing the parameter $\gamma_{k}$ through line search.

Convergence on Convex Problems.

In the following, we prove the global convergence of Algorithm 1 in case when $\varphi(\mathbf{x})$ is convex.

Theorem 3.1

Let Assumptions 1, 2a, and 3 be satisfied. Let $\gamma_{k}\coloneqq\min\{1,\frac{\Delta_{k}}{\mathcal{S}}\}$ or $\gamma_{k}\coloneqq\frac{2}{2+k}$. Then, for $k\geq 1$ it holds that

[TABLE]

Our method recovers the rate of classical FW methods for smooth problems, while being applicable to the wider class of fully composite problems (4). Thus, our $\mathcal{O}(1/k)$ rate improves upon the $\mathcal{O}(1/\sqrt{k})$ of black-box non-smooth optimization. Clearly, the improvement is achievable by leveraging the structure of the objective within the algorithm.

Convergence on Non-convex Problems.

In this case, $\Delta_{k}$ has different semantics and no longer provides an accuracy certificate for the functional residual. This quantity is nevertheless important, since it enables us to quantify the algorithm’s progress in the non-convex setting, while maintaining an affine-invariant analysis. The following theorem states the convergence guarantee on $\Delta_{k}$ for non-convex problems.

Theorem 3.2

Let Assumptions 1 and 2a be satisfied. Let $\gamma_{k}\coloneqq\min\{1,\frac{\Delta_{k}}{\mathcal{S}}\}$ or $\gamma_{k}\coloneqq\frac{1}{\sqrt{1+k}}$. Then, for all $k\geq 1$ it holds that

[TABLE]

Theorem 3.2 recovers a similar rate to the classical FW methods (Lacoste-Julien, 2016). The line search rule for parameter $\gamma_{k}$ makes our method universal, thereby allowing us to attain practically faster rates automatically when the iterates lie within a convex region of the objective.

As previously mentioned, the progress measure (14) does not upper-bound functional suboptimality in the general non-convex setting. However, in some cases, we may still be able to establish convergence of meaningful quantities for non-convex fully composite problems with the linearization method. Namely, let us consider problem (12) in Example 2.2 for the Euclidean norm, i.e., $F(\mathbf{u},\mathbf{x})=\left\lVert\,\mathbf{u}\,\right\rVert$, and the following simple iterations:

[TABLE]

Note that in (18), differently from (13), the value of $\gamma_{k}$ is selected prior to the oracle call. Denoting the squared objective as \Phi(\mathbf{x})\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\frac{1}{2}\bigl{[}\varphi(\mathbf{x})\bigr{]}^{2}=\frac{1}{2}\|\mathbf{f}(\mathbf{x})\|^{2} and following our analysis, we can state the convergence of process (18) in terms of the classical FW gap with respect to $\Phi$. The proof is deferred to Appendix A.5.

Proposition 3.1

Let $\gamma_{k}\coloneqq\frac{1}{\sqrt{1+k}}$. Then, for the iterations (18), under Assumption 2b and for all $k\geq 1$, it holds that

[TABLE]

We further show in Appendix B that $\Delta_{k}$ can be related to the classical FW gap, when our iterates lie in a smooth region of $F$. Whether we can provide a meaningful interpretation of $\Delta_{k}$ in the general non-convex case, however, remains an interesting open question.

4 The Accelerated Method

We now move away from the affine-invariant formulation of Algorithm 1 to a setting in which, by considering regularized minimization subproblems along with convexity and Lipschitz continuity of gradients, we can accelerate the Basic Method. We achieve acceleration by resorting to the well-known three-point scheme of Nesterov (1983), in which the proximal subproblem is solved inexactly via calls to oracles of type (5). This approach was first analyzed in the context of FW methods by Lan and Zhou (2016).

We propose Algorithm 2 which consists of a two-level scheme: an outer-loop computing the values of three iterates $\mathbf{y}$, $\mathbf{x}$ and $\mathbf{z}$ in $\mathcalorigin{X}$, and a subsolver computing inexact solutions to the proximal subproblem

[TABLE]

Note that the minimization in (19) does not conform to our oracle model (5) due to the quadratic regularizer. However, we can approximate its solution by iteratively solving subproblems in which we linearize the squared norm to match the template of (5). This procedure, denoted as InexactProx in Algorithm 2, returns a point $\mathbf{u}^{+}$ satisfying the optimality condition $\eta$-inexactly for some $\eta>0$:

[TABLE]

Note that condition (4) implies $P(\mathbf{u}^{+})\leq P(\mathbf{u})+\eta,\;\forall\mathbf{u}\in\mathcalorigin{X}$. Formally, the main convergence result characterizing Algorithm 2 is the following.

Theorem 4.1

Let Assumptions 1,  2b, and 3 be satisfied. We choose $\gamma_{k}\coloneqq\frac{3}{k+3}$, $\beta_{k}\coloneqq cF(\mathbf{L})\gamma_{k}$ and $\eta_{k}\coloneqq\frac{\delta}{3(k+1)(k+2)}$ where $\delta>0$ and $c\geq 0$ are chosen constants, and $F(\mathbf{L})\coloneqq\sup_{\mathbf{x}\in\mathcalorigin{X}}F(\mathbf{L},\mathbf{x})$. Then, for all $k\geq 1$ it holds that

[TABLE]

The proof of Theorem 4.1 comes from a natural sequence of steps involving the properties of the operators and the approximate optimality of $\mathbf{x}_{k+1}$. The crucial step in attaining the improved convergence is the choice of parameters $\gamma_{k}$, $\beta_{k}$ and $\eta_{k}$. Notably, the decay speed required of $\eta_{k}$ is quadratic, meaning that the subproblems are solved with fast-increasing accuracy and at the cost of additional time spent in the subsolver. The constant $\delta$ allows us to fine-tune the accuracy required for the first several iterations of the algorithm, where we can demand a lower accuracy. In practice, we can always choose $\delta=1$ as a universal rule, and the optimal choice is $\delta=F(\mathbf{L})\mathcalorigin{D}_{\mathcalorigin{X}}^{2}$ when these parameters are known. The factor $cF(\mathbf{L})$ in the definition of $\beta_{k}$ can be interpreted as the quality of the approximation of the Lipschitz constant for our problem. Namely, it is exactly computed for $c=1$, and over or underestimated for $c>1$ and $c\in(0,1)$ respectively.

We describe each of the bounding terms independently: the first is highly reminiscent of the usual bounds accompanying FW-type algorithms in terms of constants, albeit now with quadratic decay speed. The second term indicates the behavior of the algorithm as a function of $c$: overestimation of $F(\mathbf{L})$ ensures quadratic rates of convergence, since the second term becomes negative. Conversely, underestimation of $F(\mathbf{L})$ brings us back into the familiar FW convergence regime of $\mathcalorigin{O}\left(1/k\right)$ as the second term becomes positive. The extreme case $c=0$ (and hence $\beta_{k}=0$) essentially reduces Algorithm 2 to Algorithm 1, since the projection subproblem reduces to problem (5) which we assume to be easily solvable. We therefore have robustness in terms of choosing the parameter $c$ and the exact knowledge of $F(\mathbf{L})$ is not needed, even though it may come at the cost of a slower convergence. In contrast, classical Fast Gradient Methods are usually very sensitive to such parameter choices (Devolder, 2013).

Theorem 4.1 provides an accelerated rate on the iterates $\mathbf{y}_{k}$ – an analogous result to that of Lan and Zhou (2016) albeit under a different oracle. This convergence rate is conditioned on the subsolver returning an $\eta_{k}$-inexact solution to the projection subproblem and therefore any subsolver satisfying the condition can achieve this rate. As with any optimization algorithm, convergence guarantees may also be stated in terms of the oracle complexity required to reach $\varepsilon$ accuracy. For Algorithm 2 all the oracle calls are deferred to the subsolver InexactProx, which we describe and analyze in the next section.

5 Solving the Proximal Subproblem

We now provide an instance of the InexactProx subsolver which fully determines the oracle complexity of the Accelerated Method (Algorithm 2). It relies on a specific adaptation of Algorithm 1 to the structure of (19). The quadratic regularizer is linearized and oracles of type (5) are called once per inner iteration, while the Jacobian $\nabla\mathbf{f}(\mathbf{z}_{k})$ is computed once per subsolver call. The main challenge here is to find a readily available quantity defining the exit condition of the subsolver, which we denote by $\Delta_{t}$.

The parameters of Algorithm 3 are fully specified, and the stopping condition depends on $\Delta_{t}\geq P(\mathbf{u}_{t})-P^{\star}$, which is a meaningful progress measure. The algorithm selects its stepsize via closed-form line search to improve practical performance. When $F(\mathbf{u})\equiv\mathbf{u}^{(1)}$, this procedure recovers the classical FW algorithm with line search applied to problem (19).

We prove two results in relation to Algorithm 3: its convergence rate and the total oracle complexity of Algorithm 2 when using Algorithm 3 as the subsolver. The rate and analysis are similar to the ones for the Basic Method, utilizing additionally the form of the proximal subproblem.

Theorem 5.1

Let Assumptions 1,  2b, and 3 be satisfied. Then, for all $t\geq 1$ it holds that

[TABLE]

Consequently, Algorithm 3 returns an $\eta$-approximate solution according to condition (4) after at most \mathcal{O}\Bigl{(}\frac{\beta\mathcalorigin{D}_{\mathcalorigin{X}}^{2}}{\eta}\Bigr{)} iterations.

We note that oracle (5) is called once per inner iteration and the Jacobian $\nabla\mathbf{f}(\mathbf{z}_{k})$ is computed once per subsolver call. In particular, when using Algorithm 3 as a subsolver, our Accelerated Method achieves the optimal number of $\mathcalorigin{O}\left(\varepsilon^{-1/2}\right)$ Jacobian computations typical of smooth and convex optimization, while maintaining a $\mathcalorigin{O}\left(\varepsilon^{-1}\right)$ complexity for the number of calls to oracle (5). The results are stated in the following corollary.

Corollary 5.1

Consider the optimal choice of parameters for Algorithm 2, that is $c\coloneqq 1$ and $\delta\coloneqq F(\mathbf{L})\mathcalorigin{D}_{\mathcalorigin{X}}^{2}$. Then, solving problem (4) with $\varepsilon$ accuracy $\varphi(\mathbf{y}_{k})-\varphi^{\star}\leq\varepsilon$, requires \mathcal{O}\bigl{(}\sqrt{\frac{F(\mathbf{L})\mathcalorigin{D}_{\mathcalorigin{X}}^{2}}{\varepsilon}}\bigr{)} computations of $\nabla\mathbf{f}$. In addition, the total number of calls to oracle (5) is \mathcal{O}\bigl{(}\frac{F(\mathbf{L})\mathcalorigin{D}_{\mathcalorigin{X}}^{2}}{\varepsilon}\bigr{)}.

Finally, we note that for smaller values of parameter $c\in[0,1]$ in Algorithm 2 (underestimating the Lipschitz constant), the complexity of InexactProx procedure improves. Thus, for $c=0$ we have $\beta=0$ (no regularization) and Algorithm 3 finishes after just one step.

6 Experiments

The experiments are implemented in Python 3.9 and run on a MacBook Pro M1 with 16 GB RAM. For both experiments we use the Projected Subgradient Method as a baseline (Shor et al., 1985), with a stepsize of $\frac{p}{\sqrt{k}}$ where $p$ is tuned for each experiment. The CVXPY library (Diamond and Boyd, 2016) is used to solve subproblems of type (5). The random seed for our experiments is always set to 666013, and we set $c=1$ since we can analytically compute the Lipschitz constants or their upper bounds.

6.1 Minimization Over the Simplex

We consider the following optimization problem

[TABLE]

where $\mathbf{A}_{i}\in\mathbb{R}^{d\times d}$ are random PSD matrices and $\mathbf{b}\in\mathbb{R}^{d}$. The problem conforms to Example 2.1, and we use $d=500$ and $n=10$. We generate $\mathbf{A}_{i}=\mathbf{Q}_{i}\mathbf{D}\mathbf{Q}_{i}^{\top}$, where $\mathbf{D}$ is a diagonal matrix of eigenvalues decaying linearly from $1$ to $10^{-6}$, and $\mathbf{Q}_{i}$ is a randomly generated orthogonal matrix using the scipy.stats.ortho_group method (Mezzadri, 2006). The vectors $\mathbf{b}_{i}$, which determine the position of the quadratics in space, are set as follows: $\mathbf{b}_{i}=\mathbf{e}_{i}\cdot 10$, $\mathbf{b}_{9}=\mathbf{0}$ (the origin), $\mathbf{b}_{10}=\mathbf{1}\cdot 10$ (the all ones vector multiplied by $10$). We set $\delta=0.2$ in the Accelerated Method (see Theorem 4.1) and settle for $p=1.42$ following tuning of the Subgradient Method. Finally, we set $\mathbf{x}_{0}=\mathbf{e}_{3}\in\Delta_{d}$ for all methods.

The convergence results in terms of FO oracles and oracles of type (5) are shown in Figure 1 and 1, respectively. The figures highlight the improvement in terms of the number of FO calls, while showing comparable performance in terms of subproblem oracle calls, as predicted by our theory.

6.2 Minimization Over the Nuclear Norm Ball

We consider the following optimization problem

[TABLE]

Formulation (22) models a matrix completion scenario where we wish to recover an $\mathbf{X}^{\star}$ that minimizes the largest error within a given set of matrices $\mathbf{A}^{(i)}$. The matrices $\mathbf{A}^{(i)}$ are only partially revealed through a set of corresponding indices $\Omega_{i}$. This problem conforms to Example 2.1 and we use $d=30$, $m=10$, $r=7$, where $r$ is the rank of matrices $\mathbf{A}^{(i)}$. The data is generated in identical fashion as in Section 5.2 of (Lan and Zhou, 2016) on Matrix Completion. We set $\delta=100$ in the Accelerated Method (see Theorem 4.1) and settle for $p=0.2$ following tuning of the Subgradient Method. Finally, we set $\mathbf{x}_{0}=\mathbf{0}^{d\times m}\in\mathcalorigin{X}$ for all methods.

The convergence results in terms of FO oracles and oracles of type (5) are shown in Figure 2 and 2, respectively. The figures highlight the improvement in terms of the number of FO calls, while showing comparable performance in terms of subproblem oracle calls, as predicted by our theory.

7 Conclusion

Our work illustrates how improved convergence rates may be attained by assuming precise structure within a class of objectives (e.g., non-differentiable ones). Moreover, it shows how a simple principle such as linearizing the differentiable components of a function composition can be used to create more benign subproblems that are efficiently solved. Interesting future work may address relaxing the assumptions on $F$, extending this framework to stochastic settings, and meaningfully interpreting the quantity $\Delta_{k}$ in the non-convex setting.

\acks

The authors thank Aditya Varre and the anonymous reviewers for their helpful feedback on the manuscript. This work was partly supported by the Swiss State Secretariat for Education, Research and Innovation (SERI) under contract number 22.00133.

Appendix A Proofs

Assumption 2a implies global progress bounds on our fully composite objective with an inner linearization of $\mathbf{f}$, as stated in the following Lemma A.1. This lemma provides a basis for all our convergence results.

Lemma A.1

Let $\mathbf{x},\mathbf{y}\in\mathcalorigin{X}$ and $\gamma\in[0,1]$. Denote $\mathbf{y}_{\gamma}=\mathbf{x}+\gamma(\mathbf{y}-\mathbf{x})$. Then, it holds

[TABLE]

Proof A.1**.**

Note that the subhomogenity assumption (6) is equivalent to the following useful inequality for the outer component of the objective see (Theorem 7.1 in Doikov and Nesterov (2022)):

[TABLE]

Then, we have

[TABLE]

which is the desired bound.

A.1 Proof of Theorem 3.1

See 3.1

Proof A.2**.**

Indeed, for one iteration of the method, we have

[TABLE]

where we used in $(*)$ that $F(\cdot,\cdot)$ is jointly convex. Hence, we obtain the following inequality for the progress of one step, for $k\geq 0$:

[TABLE]

Now, let us choose use an auxiliary sequence $A_{k}:=k\cdot(k+1)$ and $a_{k+1}:=A_{k+1}-A_{k}=2(k+1)$. Then,

[TABLE]

which is one of the possible choices for $\gamma_{k}$. Note that for the other choice, we set \gamma_{k}=\min\bigl{\{}1,\frac{\Delta_{k}}{\mathcal{S}}\bigr{\}}, which maximizes the right hand side of (25) with respect to $\gamma_{k}\in[0,1]$. Hence, in both cases we have that

[TABLE]

or, rearranging the terms,

[TABLE]

Telescoping this bound for the first $k\geq 1$ iterations, we get

[TABLE]

It remains to prove the convergence in terms of the accuracy measure $\Delta_{k}$. For that, we telescope the bound (26), which is

[TABLE]

for the $k\geq 1$ iterations, and use the convergence for the functional residual (27):

[TABLE]

To finish the proof, we need to divide both sides by $\sum\limits_{i=1}^{k}a_{i+1}=A_{k+1}-a_{1}=k(3+k)$.

A.2 Proof of Theorem 3.2

See 3.2

Proof A.3**.**

As in the proof of the previous theorem, our main inequality (25) on the progress of one step is:

[TABLE]

where we can substitute $\gamma_{k}=\frac{1}{\sqrt{k+1}}$ for the both strategies of choosing this parameter.

Summing up this bound for the first $k+1$ iterations, we obtain

[TABLE]

Using the bound $\varphi(\mathbf{y}_{k+1})\geq\varphi^{\star}$ and our value of $\gamma_{i}$, we get

[TABLE]

which is (17).

A.3 Proof of Theorem 4.1

See 4.1

Proof A.4**.**

Let us consider one iteration of the method, for some $k\geq 0$.

Since all the components of $\mathbf{f}$ have the Lipschitz continuous gradients, it hold that

[TABLE]

where the vector inequality is component-wise. Then, using the properties of $F$, we have

[TABLE]

where the second equality comes from the update rule of $\mathbf{y}_{k+1}$, the second inequality comes from the joint convexity in Assumption 1, the third inequality comes from convexity of the components of $\mathbf{f}$ and monotonicity of $F$.

Since we are introducing a norm-regularized minimization subproblem for the purpose of acceleration, the term $\frac{\gamma_{k}^{2}F(\mathbf{L})}{2}\|\mathbf{x}_{k+1}-\mathbf{x}_{k}\|^{2}$ can be further upper bounded using $\eta_{k}$-approximate guarantee (4), as follows, for any $\mathbf{x}\in\mathcalorigin{X}$:

[TABLE]

where we used our choice $\beta_{k}=cF(\mathbf{L})\gamma_{k}$.

Therefore, combining these two bounds together, we obtain

[TABLE]

for all $\mathbf{x}\in\mathcalorigin{X}$, where we used convexity of $\mathbf{f}$ and monotonicity of $F$.

We now subtract $\varphi(\mathbf{x})$ from both sides, let $\mathbf{x}=\mathbf{x}^{\star}$ and denote $\varepsilon_{k}:=\varphi(\mathbf{y}_{k})-\varphi^{\star}$, which gives

[TABLE]

We now move on to choosing the parameters $\gamma_{k}$, $\eta_{k}$ and $\beta_{k}$ in a way that allows us to accelerate. For more flexibility, we let $\gamma_{k}:=\frac{a_{k+1}}{A_{k+1}}$, for some sequences $\{a_{k}\}_{k\geq 0}$ and $\{A_{k}\}_{k\geq 0}$ that will be defined later. Then (30) becomes:

[TABLE]

and therefore we have

[TABLE]

We wish to choose sequences $A_{k}$, $a_{k}$, $\beta_{k}$ and $\eta_{k}$ such that we obtain a $\mathcalorigin{O}\left(1/k^{2}\right)$ rate of convergence on the functional residual of $\varphi(\cdot)$. The constraint we require on the sequences is $\gamma_{k}F(\mathbf{L})\leq\beta_{k}$. The following choices

[TABLE]

give us the desired outcome, since equation (A.4) becomes:

[TABLE]

since $a_{i+1}\beta_{i}-a_{i}\beta_{i-1}=\frac{9cF(\mathbf{L})\left(i^{2}+5i+4\right)}{i^{2}+5i+6}<9cF(\mathbf{L})$ and $\left\lVert\,\mathbf{x}_{i}-\mathbf{x}^{\star}\,\right\rVert^{2}\leq\mathcalorigin{D}_{\mathcalorigin{X}}^{2}$.

A.4 Proof of Theorem 5.1

See 5.1

Proof A.5**.**

Let us introduce our subproblem, in a general form, that is

[TABLE]

where $r(\cdot)$ is a differentiable convex function, whose gradient is Lipschitz continuous with constant $\beta>0$, and $h(\mathbf{u})$ is a general proper closed convex function, not necessarily differentiable.

In our case, for computing the inexact proximal step (19), we set

[TABLE]

for a fixed $\mathbf{x}$ and $\mathbf{z}$.

Then, in each iteration of Algorithm 3, we compute, for $t\geq 0$:

[TABLE]

The optimality condition for this operation is (see, e.g. Theorem 3.1.23 in Nesterov (2018b))

[TABLE]

Therefore, employing the Lipschitz continuity of the gradient of $r(\cdot)$, we have

[TABLE]

where the last equality comes from the definition of $\Delta_{t}$ in Algorithm 3.

Note that $\alpha_{t}$ is defined as the minimizer of $\frac{\beta\alpha_{t}^{2}}{2}\|\mathbf{v}_{t+1}-\mathbf{u}_{t}\|^{2}-\alpha_{t}\Delta_{t}$ and hence, for any other $\rho_{t}\in[0,1]$ it will hold that:

[TABLE]

At the same time,

[TABLE]

where the last line follows from the convexity of $r(\mathbf{u})$. Letting $\mathbf{u}:=\mathbf{u}^{\star}$ (solution to (31)) in (36) and further substituting it into (35) and subtracting $\mathbf{s}^{\star}$ from both sides, we obtain

[TABLE]

Now, let us choose $\rho_{t}:=\frac{a_{t+1}}{A_{t+1}}$ for sequences $A_{t}:=t\cdot(t+1)$, and $a_{t+1}:=A_{t+1}-A_{t}=2(t+1)$. Then, \begin{array}[]{rcl}\rho_{t}&:=&\frac{2}{2+t},\;\;t\geq 0.\end{array} Using this choice, inequality (37) can be rewritten as

[TABLE]

Telescoping this inequality for the first iterations, we obtain, for $t\geq 1$:

[TABLE]

This is the global convergence in terms of the functional residual. It remains to justify the convergence for the accuracy certificates $\Delta_{t}$. Multiplying (35) by $A_{t+1}$, we obtain

[TABLE]

Telescoping this bound, we get, for $t\geq 1$:

[TABLE]

Dividing both sides by $\sum_{i=1}^{t}a_{i+1}=A_{t+1}-A_{1}=t(3+t)$ completes the proof we finally get:

[TABLE]

A.5 Proof of Proposition 3.1

See 3.1

Proof A.6**.**

In our case, we have $\varphi(\mathbf{x})\equiv\|\mathbf{f}(\mathbf{x})\|_{2}$. Using Lemma A.1, we obtain

[TABLE]

where $\mathbf{x}_{k+1}\in\mathcalorigin{X}$ is the point such that $\mathbf{y}_{k+1}=\mathbf{y}_{k}+\gamma_{k}(\mathbf{x}_{k+1}-\mathbf{y}_{k})$. Using convexity of the function $g(\mathbf{x})\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\|\mathbf{f}(\mathbf{y}_{k})+\gamma_{k}\nabla\mathbf{f}(\mathbf{y}_{k})(\mathbf{x}-\mathbf{y}_{k})\|_{2}$, we get that

[TABLE]

where the subgradient $g^{\prime}(\mathbf{x}_{k+1})=\gamma_{k}\nabla\mathbf{f}(\mathbf{y}_{k})^{\top}\frac{\mathbf{f}_{k+1}}{\|\mathbf{f}_{k+1}\|_{2}}$ with $\mathbf{f}_{k+1}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\mathbf{f}(\mathbf{y}_{k})+\gamma_{k}\nabla\mathbf{f}(\mathbf{y}_{k})(\mathbf{x}_{k+1}-\mathbf{y}_{k})$, satisfies the stationary condition for the method step:

[TABLE]

A few comments are in order now about the use of the subgradient above. Note that we wish to impose an assumption on $\mathbf{f}$ which can ensure that $\mathbf{f}(\mathbf{y}_{k})+\gamma_{k}\nabla\mathbf{f}(\mathbf{y}_{k})(\mathbf{x}-\mathbf{y}_{k})\neq\mathbf{0}\in\mathbb{R}^{n}$. First, some preliminaries. Under Assumption 2b on $\mathbf{f}$, it holds that:

[TABLE]

From here, we can bound the products between Jacobians and iterates as follows:

[TABLE]

Thus, without loss of generality, we can shift $\mathbf{f}$ by a constant vector of identical values depending on $\mathcal{G}\mathcalorigin{D}_{\mathcalorigin{X}}$ such that we ensure, for example, $\mathbf{f}(\mathbf{y}_{k})+\gamma_{k}\nabla\mathbf{f}(\mathbf{y}_{k})(\mathbf{x}-\mathbf{y}_{k})>\mathbf{0}$ component-wise. Hence, combining these observations with (40), we have

[TABLE]

Then, by lower bounding appropriately using (42) and (43), we get:

[TABLE]

Substituting $\gamma_{k}:=\frac{1}{\sqrt{1+k}}$ and telescoping this bound would lead to the desired global convergence (for the details, see the end of the proof of Theorem 3.2).

Appendix B Interpretation of $\Delta_{k}$ in the non-convex setting

While we cannot make any strong claims about the meaning of $\Delta_{k}$ in general, we can provide an additional observation for this quantity when the outer component $F$ is smooth inside a ball included in $\mathcalorigin{X}$.

Thus, consider a ball of radius $\varepsilon$ centered at $\mathbf{y}_{k}$ denoted by $B(\mathbf{y}_{k},\varepsilon)=\{\mathbf{x}\in\mathbb{R}^{d}\;:\;\|\mathbf{x}-\mathbf{y}_{k}\|\leq\varepsilon\}$, and set $\mathcal{B}=B(\mathbf{y}_{k},\varepsilon)\cap\mathcalorigin{X}$. Assuming that $F(\mathbf{u},\mathbf{x})$ is differentiable at all points from $\mathbb{R}^{n}\times\mathcal{B}$, and that its gradient is Lipschitz continuous with constant $L_{F}$, we have for any $\mathbf{x}\in\mathcal{B}\subseteq\mathcal{X}$:

[TABLE]

Hence, for a small enough ball, $\Delta_{k}$ is an $\mathcalorigin{O}\left(\varepsilon^{2}\right)$-approximation of the original FW gap restricted to the considered neighborhood. If, in addition, the composite function $\varphi$ is convex in $\mathcal{B}$ and there is a local optimum $\mathbf{x}^{\star}\in\mathcal{B}$, then $\Delta_{k}$ is an $\mathcalorigin{O}\left(\varepsilon^{2}\right)$-approximation of functional suboptimality.

Appendix C Additional Application Examples

Example C.1**.**

We define a generalized nonlinear model as,

[TABLE]

where $\phi:\mathbb{R}\to\mathbb{R}$ is a fixed convex loss function, and $n$ is the number of data points. Problem (4) then reduces to training a (non-convex) model, for example a neural network, with respect to the constraint set $\mathcalorigin{X}$: \begin{array}[]{rcl}\min\limits_{\mathbf{x}\in\mathcalorigin{X}}\sum\limits_{i=1}^{m}\phi(f_{i}(\mathbf{x})).\end{array}

Solving this problem then involves training a linear model within the basic subroutine (5) \displaystyle\begin{array}[]{rcl}\min\limits_{\mathbf{x}\in\mathcalorigin{X}}\sum\limits_{i=1}^{m}\phi(\langle\mathbf{a}_{i},\mathbf{x}\rangle+b_{i}),\end{array} which is a convex problem. Amongst the loss functions relevant to Machine Learning, the following are convex and subhomogeneous thus making $F$ in (45) satisfy Assumption 1:

• •

$\ell_{1}$-regression: $\phi(t)=|t|$

• •

Hinge loss (SVM): $\phi(t)=\max\{0,t\}$

• •

Logistic loss: $\phi(t)=\log(1+e^{t})$