Advances in Low-Memory Subgradient Optimization

TL;DR

This paper reviews recent advances in low-memory subgradient algorithms for non-smooth convex optimization, highlighting techniques that improve speed and efficiency for large-scale problems with minimal storage requirements.

Contribution

It introduces modern methods like Nesterov smoothing, Universal Mirror Prox, and adaptive Mirror Descent, emphasizing their theoretical complexity bounds and practical applications.

Findings

Universal Mirror Prox algorithm for variational inequalities.

Primal-dual Mirror Descent method optimal for Lipschitz problems.

Application to sparse Truss Topology Design problem.

Abstract

This chapter is devoted to the black-box subgradient algorithms with the minimal requirements for the storage of auxiliary results, which are necessary to execute these algorithms. It starts with the original result of N.Z. Shor which open this field with the application to the classical transportation problem. To discuss the fundamentals of non-smooth optimization the theoretical complexity bounds for smooth and non-smooth convex and quasi-convex optimization problems are briefly exposed with the special attention given to adaptive step-size policy. Than this chapter contains descriptions of different modern techniques that allow to solve non-smooth convex optimization problems faster then lower complexity bounds: Netserov smoothing technique, Netserov Universal approach, Legendre (saddle point) representation approach. We also describe recent Universal Mirror Prox algorithm for variational inequalities. To demonstrate application of non-smooth convex optimization algorithms for solution of huge-scale extremal problems we consider convex optimization problems with non-smooth functional constraints and describe two adaptive Mirror Descent methods. The first method is of primal-dual nature and proved to be optimal in terms of lower oracle bounds for the class of Lipschitz-continuous convex objective and constraints. The advantages of application of this method to sparse Truss Topology Design problem are discussed in some details. The second method can be applied for solution of convex and quasi-convex optimization problems and is optimal in a sense of complexity bounds. The conclusion part of the survey contains the important references that characterize recent developments of non-smooth convex optimization.

Tables

Table 1. Table 1: Optimal estimates for the number of oracle calls

$N\le n$	$\left|f\left(\overrightarrow{y}\right)-f\left(\overrightarrow{x}\right)\right|\le M\Vert \overrightarrow{y}-\overrightarrow{x}\Vert$	${\Vert \nabla f\left(\overrightarrow{y}\right)-\nabla f\left(\overrightarrow{x}\right)\Vert }_{\ast }\le L\Vert \overrightarrow{y}-\overrightarrow{x}\Vert$
$f\left(\overrightarrow{x}\right)$ convex	$O\left(\frac{M^2{R}^2}{{\epsilon }^2}\right)$	$O\left(\sqrt{\frac{L{R}^2}{\epsilon }}\right)$
$f\left(\overrightarrow{x}\right)\mu -$strongly convex in $\parallel \cdot \parallel$-norm	$\tilde{O}\left(\frac{M^2}{\mu \epsilon }\right)$	$\tilde{O}\left(\sqrt{\frac{L}{\mu }}\lceil \ln \left(\frac{\mu R^2}{\epsilon }\right)\rceil \right)\left(\forall N\right)$

Equations

$$\underset{x\in E}{\min }f(x)=f_{\star }=f({\vec{x}}^{\star }),{\vec{x}}^{\star }\in X_{\star },$$

$${\begin{array}[]{c}\min\ \sum_{i=1}^{m}\sum_{j=1}^{n}c_{ij}x_{ij}\\ \sum_{i=1}^{m}x_{ij}=a_{j},\ j=1,2,\dots,n;\\ \sum_{j=1}^{n}x_{ij}=b_{i},\ i=1,2,\dots,m\\ x_{ij}\geq 0,i=1,2,\dots,m;j=1,2,\dots,n\end{array}}$$

$$\max \Phi (\vec{u},\vec{v})$$

$$\Phi (\vec{u},\vec{v})=\underset{\vec{x}\ge 0}{\inf }L(\vec{x},\vec{u},\vec{v})$$

$$\begin{gathered} L(\vec{x},\vec{u},\vec{v})=\sum _{i=1}^m\sum _{j=1}^n{c}_{ij}{x}_{ij}+ \\ \sum _{j=1}^n{u}_j(\sum _{i=1}^m{x}_{ij}-a_j)+\sum _{i=1}^m{v}_i(\sum _{j=1}^n{x}_{ij}-b_i). \end{gathered}$$

$${\begin{array}[]{c}\Phi({\vec{u}},{\vec{v}})=-m\sum_{j=1}^{n}u_{j}a_{j}-n\sum_{i=1}^{m}v_{i}b_{i}+\sum_{i=1}^{m}\sum_{j=1}^{n}\inf_{{\vec{x}}\geq 0}x_{ij}\{c_{ij}+u_{j}+v_{i}\}=\\ -m\sum_{j=1}^{n}u_{j}a_{j}-n\sum_{i=1}^{m}v_{i}b_{i}-\operatorname{Ind}_{D}({\vec{u}},{\vec{v}}),\end{array}}$$

$${\operatorname{Ind}_{D}({\vec{u}},{\vec{v}})=\left\{\begin{array}[]{ll}0&\mbox{ when }c_{ij}+u_{i}+v_{j}\geq 0;\\ \infty&\mbox{ otherwise.}\end{array}\right.}$$

$$\Phi (\vec{u},\vec{v})=m\sum _{j=1}^n{u}_j{a}_j+n\sum _{i=1}^m{v}_i{b}_i+$$

$$\sum _{i=1}^m\sum _{j=1}^n\underset{\vec{x}\ge 0}{\sup }{x}_{ij}\{c_{ij}+u_j+v_i\}=$$

$$=m\sum _{j=1}^n{u}_j{a}_j+n\sum _{i=1}^m{v}_i{b}_i+{\mathrm{Ind}}_D(\vec{u},\vec{v}),$$

$${\vec{u}}^{k+1}={\vec{u}}^k-\lambda g_L^u={\vec{u}}^k-\lambda \vec{a};{\vec{v}}^{k+1}={\vec{v}}^k-\lambda g_L^v={\vec{v}}^k-\lambda \vec{b};k=0,1,\dots$$

$${\Phi }_V(\vec{u},\vec{v})=m\sum _{j=1}^n{u}_j{a}_j+n\sum _{i=1}^m{v}_i{b}_i-$$

$$\sum _{i=1}^m\sum _{j=1}^n\underset{0\le \vec{x}\le V}{\max }{x}_{ij}\{c_{ij}+u_j+v_i\}=$$

$$=m\sum _{j=1}^n{u}_j{a}_j+n\sum _{i=1}^m{v}_i{b}_i+P_V(\vec{u},\vec{v})$$

$$P_V(\vec{u},\vec{v})=\sum _{i=1}^m\sum _{j=1}^n\underset{x_{ij}\in [0,V]}{\max }{x}_{ij}\{c_{ij}+u_j+v_i\}=$$

$$\sum _{i=1}^m\sum _{j=1}^{n}V\{c_{ij}+u_j+v_i{\}}_{+}$$

$${\vec{g}}_{ij}=\left\{\begin{array}[]{ll}V&\mbox{ if }c_{ij}+u_{j}+v_{i}>0\\ 0&\mbox{ if }c_{ij}+u_{j}+v_{i}<0\\ \operatorname{cone}(0,V)&\mbox{ if }c_{ij}+u_{j}+v_{i}=0\end{array}\right.$$

$${\vec{x}}^{k+1}={\vec{x}}^k-\lambda {\overline{\stackrel{⃗}{g}}}^k,k=0,1,\dots$$

$$x^{k+1}={\vec{x}}^k-\lambda g_{\nu }^k,k=0,1,\dots$$

$$f(\widetilde{{\stackrel{⃗}{x}}^k})=f({\vec{x}}^k)\text{mbox}and\mathrm{dist}(\widetilde{{\stackrel{⃗}{x}}^k},X_{\star })\le \lambda (1+ϵ)/2.$$

$$\Vert {\vec{x}}^{k+1}-{\vec{x}}^{\star }{\Vert }^2=\Vert {\vec{x}}^k-{\vec{x}}^{\star }-\lambda {\vec{g}}_{\nu }^k{\Vert }^2=\Vert {\vec{x}}^k-{\vec{x}}^{\star }{\Vert }^2+{\lambda }^2-2\lambda {\overline{\stackrel{⃗}{g}}}^k({\vec{x}}^k-{\vec{x}}^{\star }).$$

$$\begin{array}[t]{c}\min\\ {\vec{z}}\in H_{k}\end{array}\|{{\vec{x}}^{\star}}-{\vec{z}}\|^{2}=\|{{\vec{x}}^{\star}}-{\vec{z}}^{k}\|^{2}=\delta_{k},$$

$$\Vert {\vec{x}}^{k+1}-{\vec{x}}^{\star }{\Vert }^2=\Vert {\vec{x}}^k-{\vec{x}}^{\star }{\Vert }^2+{\lambda }^2-2\lambda {\delta }_k,k=0,1,2,\dots$$

$$\Vert {\vec{x}}^{k+1}-{\vec{x}}^{\star }{\Vert }^2\le \Vert {\vec{x}}^k-{\vec{x}}^{\star }{\Vert }^2-{\lambda }^2ϵ,k=0,1,2,\dots$$

$$0\le \Vert {\vec{x}}^{k+1}-{\vec{x}}^{\star }{\Vert }^2\le \Vert {\vec{x}}^0-{\vec{x}}^{\star }{\Vert }^2\le -k{\lambda }^2ϵ\to -\infty$$

$$\underset{z:f(z)=f({\vec{x}}^{k_0})}{\min }\Vert {\vec{x}}^{\star }-z{\Vert }^2=\Vert {\vec{x}}^{\star }-{\bar{z}}^{k_0}{\Vert }^2=\underset{z:f(z)\ge f(x^{k_0})}{\min }\Vert {\vec{x}}^{\star }-z{\Vert }^2\le \Vert {\vec{x}}^{\star }-z^{k_0}{\Vert }^2={\delta }_{k_0}.$$

$$N\ge n=\dim \vec{x},$$

$$f\left(\vec{x}\right)\to \underset{\vec{x}\in Q}{\min },$$

$$f\left({\vec{x}}^N\right)-f_{\ast }\le \epsilon ,$$

$$N\sim n\ln \left(\Delta f/\epsilon \right),$$

$$\Lambda =\sum _{i,j=1,1}^{m,n}{\log }_2\left|a_{ij}\right|+\sum _{i=1}^m{\log }_2\left|{\vec{b}}_i\right|+{\log }_2\left(mn\right)+1.$$

$$\underset{\vec{x}}{\min }{\Vert (A\vec{x}-\vec{b}{)}_{+}\Vert }_{\infty }\ge 2^{-\left(\Lambda -1\right)}.$$

$${\Vert (A\vec{x}-\vec{b}{)}_{+}\Vert }_{\infty }\to \underset{{\Vert {\vec{x}}_{\ast }\Vert }_{\infty }\le 2^{\Lambda }}{\min }.$$

Full text

Introduction

We consider a finite-dimensional non-differentiable convex optimization problem (COP)

[TABLE]

where $E$ denotes a finite-dimensional space of primal variables and $f:E\to\mathbb{R}$ is a finite convex function, not necessarily differentiable. For a given point ${\vec{x}}$ the subgradient oracul returns value of objective function at that point $f({\vec{x}})$ and subgradient $g\in\partial f({\vec{x}})$. We do not make any assumption about the choice of $g$ from $\partial f({\vec{x}})$. As we are interested in computational issues related to solving (1) mainly we assume that this problem is solvable and has nonempty and bounded set of solutions $X_{\star}$.

This problem enjoys a considerable popularity due to its important theoretical properties and numerous applications in large-scale structured optimization, discrete optimization, exact penalization in constrained optimization, and others. Non-smooth optimization theory made it possible to solve in an efficient way classical discrete min-max problems ddm2002 , $l_{1}$-approximation and others, at the same time opening new approaches in bi-level, monotropic programming, two-stage stochastic optimization, to name a few.

As a major steps in:the development of different algorithmic ideas we can start with the subgradient algorithm due to Shor (see Shor2012 for the overview and references to earliest works).

1 Example Application: Transportation Problem and The First Subgradient Algorithm

From utilitarian point of view the development of non-smooth (convex) optimization started with the classical transportation problem

[TABLE]

which is widely used in many applications.

By dualizing this problem with respect to balancing constrains we can convert (2) into dual problem of the kind

[TABLE]

where ${\vec{u}}=(u_{i},i=1,2,\dots,m);{\vec{v}}=(v_{j},j=1,2,\dots,n)$ are dual variables associated with the balancing constraints in (2) and $\Phi({\vec{u}},{\vec{v}})$ is the dual function defined as

[TABLE]

and $L({\vec{x}},{\vec{u}},{\vec{v}})$ is the Lagrange function of the problem:

[TABLE]

By rearranging terms in this expression we can obtain the following expression for the dual function

[TABLE]

where

[TABLE]

is the indicator function of the set $D=\{{\vec{u}},{\vec{v}}:c_{ij}+u_{j}+v_{i}\geq 0,i=1,2,\dots,m;j=1,2,\dots,n\}$ which is the feasible set of the dual problem.

Of course, by explicitely writing feasibility constraints for (3) we obtain the linear dual transportation problem with a fewer variables but with much higher number of constraints. This is bad news for textbook simplex method so many specialized algorithms were developed, one of them was simple-minded method of generalized gradient which started the development of non-smooth optimization.

This method relies on subgradient of concave function $\Phi({\vec{u}},{\vec{v}})$ which we can transform into convex just by changing signs and replacing $\inf$ with $\sup$

[TABLE]

and ask for its minimization.

According to convex analysis TR1970 the subdifferential ${\partial_{c}}\Phi({\vec{u}},{\vec{v}})$ exists for any ${\vec{v}},{\vec{u}}\in\operatorname{int}\operatorname{dom}(\operatorname{Ind}_{D})$, and in this case just equals to the (constant) vector $g_{L}=({\vec{g}}_{\vec{u}},{\vec{g}}_{\vec{v}})=({\vec{a}},{\vec{b}})$ of a linear objective in the interior of $D$. The situation becomes more complicated when ${\vec{u}},bv$ happens to be at the boundary of $D$, the subdifferential set ceases to be a singleton and becomes even unbounded, roughly speaking certain linear manifolds are added to $g_{L}$ but we will not go into details here. The difficulty is that if we mimic gradient method of the kind

[TABLE]

with a certain step-size $\lambda>0$, we inevitably violate the dual feasibility constraints as ${\vec{a}},{\vec{b}}>0.$ Than the dual function (1) becomes undefined and correspondently the subdifferential set becomes undefined as well.

There are at least two simple ways to overcome this difficulty. One is to incorporate in the gradient method certain operations which restore feasibility and the appropriate candidate for it is the orthogonal projection operation where one can make use of the special structure of constraints and sparsity. However it will still require computing projection operator of the kind $B^{T}(BBT)^{-1}B$ for basis matrices $B$ with rather uncertain number of iteration and of matrices of the size around $(n+m)\times(n+m)$. Neither computers speed nor memory sizes at that time where not up to demands to solve problems of $n+m\approx 10^{4}$ which was required by GOSPLAN!

The second ingenious way was to take into account that if $\sum_{j=1}^{n}a_{j}=\sum_{i=1}^{m}b_{i}=V$, which is required anyway for solvability of transportation problem in a closed form. The flow variables may be uniformally bounded by $V$ and the dual function can be redefined as

[TABLE]

where the penalty function $P_{V}({\vec{u}},{\vec{v}})$ is easily computed by component-wise maximization:

[TABLE]

where $\{\cdot\}_{+}=\max\{0,\cdot\}$. Than the dual objective function becomes finite, the optimization problem — unconstrained and we can use simple subgradient method with very low requirements for memory and computations.

Actually even tighter bounds ${\vec{x}}_{ij}\leq\min(a_{i},b_{j})$ can be imposed on the flow variables which may be advantageous for computational reasons.

In both cases there is a fundamental problem of recovering optimal primal $n\times m$ primal solution from $n+m$ dual. This problem was studied by many authors and recent advances in this area can be studied from the excellent paper by A. Nedic and A. Ozdoglar nedoz09 . Theoretically speaking, nonzero positive values of $c_{ij}+u_{j}^{\star}+v_{i}^{\star}$, where ${\vec{u}}^{\star},{\vec{v}}^{\star}$ are the exact optimal solutions of the dual problem (3) signal that the corresponding optimal primal flow $x_{ij}^{\star}$ is equal to zero. Hopefully after excluding these variables we obtain nondegenerate basis and can compute the remaining variables by simple and efficient linear algebra, especially taking into account the uni-modularity of basis.

However the theoretical gap between zeros and non-zeros is exponentially small even for modest length integer data therefore we need an accuracy unattainable by modern 64-128 bits hardware. Also the real life computations are always accompanied by numerical noise and we face the hard choice in fact guessing which dual constraints are active and which are not.

To connect the transportation problem with non-smooth optimization notice that the penalty function $P_{V}({\vec{u}},{\vec{v}})$ is finite with the subdifferential ${\partial_{c}}P_{V}({\vec{u}},{\vec{v}})$ which can be represented as a set of $n\times m$ matrices

[TABLE]

so the subdifferential set is a convex hull of up to $2^{(n+m)}$ extreme points — enormous number even for a modest size transportation problem. Nevertheless it is easy to get at least single member of subdifferential and consider the simplest version of subgradient method:

[TABLE]

where ${\vec{x}}^{0}$ is a given starting point, $\lambda>0$ — fixed step-size and $\bar{\vec{g}}^{k}={\vec{g}}^{k}/\|{\vec{g}}^{k}\|$ is a normalized subgradient ${\vec{g}}^{k}\in\partial f({\vec{x}}^{k})$. Of course we assume that ${\vec{g}}^{k}\neq 0$ otherwise ${\vec{x}}^{k}$ is already a solution.

Of course, there is no hope of classical convergence result such that ${\vec{x}}^{k}\to{{\vec{x}}^{\star}}\in{X_{\star}}$, but the remarkable theorem of Shor Shor79 establishes that this simplest algorithm determines at least the approximate solution. As a major step in the development of different algorithmic ideas we can start with the subgradient algorithm due to Shor (see Shor2012 for the overview and references to earliest works). Of course, there is no hope of classical convergence result such that ${\vec{x}}^{k}\to{{\vec{x}}^{\star}}\in{X_{\star}}$, but the remarkable theorem of Shor Shor79 establishes that this very simple algorithm provides an approximate solution of (1) at least theoretically.

Theorem 1.1

Let $f$ is a finite convex function with a subdifferential $\partial f$ and the sequence $\{{{\vec{x}}^{k}}\}$ is obtained by the recursive rule

[TABLE]

with $\lambda>0$ and $g^{k}_{\nu}=g^{k}/\|g^{k}\|,g^{k}\in\partial f({{\vec{x}}^{k}}),\ g^{k}\neq 0$ is a normalized subgradient at the point $x^{k}$. Then for any $\epsilon>0$ there is an infinite set $Z_{\epsilon}\subset Z$ such that for any $k\in Z_{\epsilon}$

[TABLE]

The statement of the theorem is illustrated on Fig. 1 together with the idea of the proof.

The detailed proof of the theorem goes like following: Let ${{\vec{x}}^{\star}}\in{X_{\star}}$ and estimate

[TABLE]

The last term in fact equals

[TABLE]

where $H_{k}=\{{\vec{z}}:{{\vec{z}}g^{k}_{\nu}}={{{\vec{x}}^{k}}g^{k}_{\nu}}$ is a hyperplane, orthogonal to $g^{k}_{\nu}$ and passing through the point ${{\vec{x}}^{k}}$, so

[TABLE]

If $\lambda^{2}-2\lambda\delta_{k}\leq-\lambda^{2}\epsilon$ for any $k\in Z$ then

[TABLE]

therefore

[TABLE]

when $k\to\infty$. This contradiction proves that there is $k_{0}$ such that $\lambda^{2}-2\lambda\delta_{k_{0}}>-\lambda^{2}\epsilon$ or $\delta_{k_{0}}<\lambda(1+\epsilon)/2$.

To complete the proof notice that by convexity $f(z^{k_{0}})\geq f({\vec{x}}_{k_{0}})$ and therefore

[TABLE]

By setting $\tilde{\vec{x}}^{0}=z^{k_{0}}$ we obtain $\|{{\vec{x}}^{\star}}-\tilde{\vec{x}}^{0}\|^{2}<\lambda(1+\epsilon)/2$.

By replacing $x^{0}$ in (11) by $\tilde{\vec{x}}^{0}$ and repeating the reasoning above we obtain $\tilde{x}^{1}$ such that $\|{{\vec{x}}^{\star}}-\tilde{\vec{x}}^{1}\|^{2}<\lambda(1+\epsilon)/2$, then in the same manner $\tilde{\vec{x}}^{2},\tilde{\vec{x}}^{3}$ and so on with $\|{{\vec{x}}^{\star}}-\tilde{\vec{x}}^{k}\|^{2}<\lambda(1+\epsilon)/2,k=2,3,\dots$ which complete the proof.

2 Complexity Results for Convex Optimization

At this section we describe the complexity results for non-smooth convex optimization problems. Most of the results mentioned below can be found in books nemirovsky1983problem ; polyak1987introduction ; nesterov2018lectures ; bubeck2015 ; ben-tal2015lectures . We start with the ‘small dimensional problems’, when

[TABLE]

where $N$ is a number of oracle calls (number of subgradient calculations or/and calculations of separation hyperplane to some simple set at a given point).

Let’s consider convex optimization problem

[TABLE]

where $Q$ – is a compact and simple set (it’s significant here!). Based on at least $N$ subgradient calculations (in general, oracle calls) we would like to find such a point ${\vec{x}}^{N}$ that

[TABLE]

where $f_{\ast}=f\left({{\vec{x}}_{\ast}}\right)$ is an optimal value of function in (13), ${\vec{x}}_{\ast}$ – the solution of (13). The lower and the upper bounds for the oracle complexity is (up to a multiplier, which has logarithmic dependence on some characteristic of the set $Q$)

[TABLE]

where \Delta f=\mathop{\sup}\limits_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\vec{x}},{\vec{y}}}\in Q}\left\{{f\left({\vec{y}}\right)-f\left({\vec{x}}\right)}\right\}. The center of gravity method Levin1965 ; Newman1965 converges according to this estimate. The center of gravity method in $n=1$ is a simple binary search method brent1973algorithms . But in $n>1$ this method is hard to implement. The complexity of iteration is too high, because we required center of gravity oracle bubeck2015 . Well known ellipsoid method Shor85 ; nemirovsky1983problem requires111Here and below for all (large) $n$: $\tilde{O}(g(n))\leq C\cdot(\ln n)^{r}g(n)$ with some constants $C>0$ and $r\geq 0$. Typically, $r=1$. If $r=0$, then $\tilde{O}(\cdot)=O(\cdot)$. $N=\tilde{O}\left({n^{2}\ln\left({{\Delta f}\mathord{\left/{\vphantom{{\Delta f}\varepsilon}}\right.\kern-1.2pt}\varepsilon}\right)}\right)$ oracle calls and $O\left({n^{2}}\right)$ iteration complexity. In Vaidya1989 ; bubeck2015 a special version of cutting plane method was proposed. This method (Vayda’s method) requires $N=\tilde{O}\left({n\ln\left({{\Delta f}\mathord{\left/{\vphantom{{\Delta f}\varepsilon}}\right.\kern-1.2pt}\varepsilon}\right)}\right)$ oracle calls and has iteration complexity $\tilde{{\rm O}}\left({n^{2.37}}\right)$. In the work lee2015faster there proposed a method with $N=\tilde{O}\left({n\ln\left({{\Delta f}\mathord{\left/{\vphantom{{\Delta f}\varepsilon}}\right.\kern-1.2pt}\varepsilon}\right)}\right)$ oracle calls and iteration complexity $\tilde{O}\left({n^{2}}\right)$. Unfortunately, for the moment it’s not obvious that this method is very practical one due to the large log-factors in $\tilde{O}\left(\right)$.

Based on ellipsoid method in the late 70-th Leonid Khachyan showed khachiyan1979polynomial that LP is in P in byte complexity. Let us shortly explain the idea. The main question is whether $A{\vec{x}}\leq{\vec{b}}$ is solvable or not, where $n=\dim{\vec{x}}$, $m=\dim{\vec{b}}$ and all elements of $A$ and ${\vec{b}}$ are integers. We would like also to find one of the exact solutions ${\vec{x}}_{\ast}$. This problem up to a logarithmic factor in complexity is equivalent to the problem to find the exact solution of LP problem $\left\langle{{\vec{c}},{\vec{x}}}\right\rangle\to\mathop{\min}\limits_{A{\vec{x}}\leq{\vec{b}}}$ with integer $A$, ${\vec{b}}$ and ${\vec{c}}$. We consider only inequality constraints as it is known that to find the exact solution of $A{\vec{x}}={\vec{b}}$ one can use polynomial Gauss elimination algorithm with $O\left({n^{3}}\right)$ arithmetic operations (a.o.) complexity.

Let us introduce

[TABLE]

If $A{\vec{x}}\leq{\vec{b}}$ is compatible, then there exists such ${\vec{x}}_{\ast}$ that $\left\|{{\vec{x}}_{\ast}}\right\|_{\infty}\leq 2^{\Lambda}$, $A{\vec{x}}_{\ast}\leq{\vec{b}}$ otherwise

[TABLE]

Thus, the question of compatibility of $A{\vec{x}}\leq{\vec{b}}$ is equivalent to the problem of finding minimum of the following non-smooth convex optimization problem

[TABLE]

The approach of khachiyan1979polynomial is to apply ellipsoid method for this problem with $\varepsilon=2^{-\Lambda}$. From the complexity of this method, it follows that in $O\left({n\Lambda}\right)$-bit arithmetic with $\tilde{O}\left({mn+n^{2}}\right)$ cost of PC memory one can find ${\vec{x}}_{\ast}$ (if it exists) in $\tilde{O}\left({n^{3}\left({n^{2}+m}\right)\Lambda}\right)$ a.o.

Note, that in the ideal arithmetic with real numbers it is still an open question blum2012complexity whether it is possible to find the exact solution of LP an problem (with the data given by real numbers) in polynomial time in the ideal arithmetic ($\pi\cdot e$ – costs $O\left(1\right))$.

Now let us consider ‘large dimensional problems’

[TABLE]

Table 1 describes (for more details see ben-tal2015lectures ; bubeck2015 ; nesterov2018lectures ) optimal estimates for the number of oracle calls for convex optimization problem (13) in the case when $N\leq n.$ Now $Q$ is not necessarily compact set.

Here $R$ is a “distance” (up to a $\ln n$-factor) between starting point and the nearest solution

[TABLE]

Let’s describe optimal method in the most simple case: $Q={\mathbb{R}}^{n}$, $\left\|\cdot\right\|=\left\|\cdot\right\|_{2}$ polyak1987introduction ; nesterov2009primal-dual . Define

[TABLE]

The main iterative process is (for simplicity we’ll denote arbitrary element of $\partial f\left({\vec{x}}\right)$ as $\nabla f\left({\vec{x}}\right))$

[TABLE]

Assume that under ${\vec{x}}\in B_{2}^{n}\left({{\vec{x}}_{\ast},\sqrt{2}R}\right)$

[TABLE]

where $R=\left\|{{\vec{x}}^{0}-{\vec{x}}_{\ast}}\right\|_{2}$.

Hence, from (14), (15) we have

[TABLE]

[TABLE]

[TABLE]

Here we choose ${\vec{x}}={\vec{x}}_{\ast}$ (if {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\vec{x}}}_{\ast} isn’t unique, we choose the nearest ${\vec{x}}_{\ast}$ to ${\vec{x}}^{0})$

[TABLE]

[TABLE]

[TABLE]

If

[TABLE]

then

[TABLE]

Note that the precise lower bound for fixed steps first-order methods for the class of convex optimization problems with (15) Drori-Teboulle2016

[TABLE]

Inequality (18) means that (see also Table 1)

[TABLE]

So, one can mentioned that if we will use in (14)

[TABLE]

the result (18) holds with nesterov2009primal-dual

[TABLE]

If we put in (19),

[TABLE]

like in (17), the result similar to (18) also holds

[TABLE]

not only for the convex functions, but also for quasi-convex functions Polyak1969 ; nesterov1989 :

[TABLE]

Note that

[TABLE]

Hence, for all $k=0,...,N$,

[TABLE]

therefore

[TABLE]

Inequality (20) justifies that we need assumption (15) holds only with ${\vec{x}}\in B_{2}^{n}\left({{\vec{x}}_{\ast},\sqrt{2}R}\right)$.

For the general (constrained) case (13) we introduce a norm $\left\|\cdot\right\|$ and some prox-function $d\left({\vec{x}}\right)\geq 0$, which is continuous and 1-strongly convex with respect to $\left\|\cdot\right\|$, i.e. $d(y)-d(x)-\langle d(x),y-x\rangle\geq\frac{1}{2}\|x-y\|^{2}$, for all $x,y\in Q$. We also introduce Bregman’s divergence ben-tal2015lectures

[TABLE]

We set $R^{2}=V[{\vec{x}}^{0}]({\vec{x}}_{\ast})$, where ${\vec{x}}_{\ast}$ – is solution of (13) (if ${\vec{x}}_{\ast}$ isn’t unique then we assume that ${\vec{x}}_{\ast}$ is minimized $V[{\vec{x}}^{0}]({\vec{x}}_{\ast})$. The natural generalization of iteration process (14) is Mirror Descent algorithm nemirovskii1979efficient ; ben-tal2015lectures which iterates as

[TABLE]

For this iteration process instead of (16) we have

[TABLE]

where $\left\|{\nabla f\left({\vec{x}}\right)}\right\|_{\ast}\leq M$ for all x:{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}V[{\vec{x}}]({\vec{x}}_{*})\leq 2V[{\vec{x}}^{0}]({\vec{x}}_{*})}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}=2R^{2}} , see also Section 4.

Analogues of formulas (17), (18), (20) are also valid

[TABLE]

where

[TABLE]

and

[TABLE]

In ben-tal2015lectures authors discus how to choose $d({\vec{x}})$ for different simple convex sets $Q$. One of these examples (unit simplex) will considered below. Note, that in all these examples one can guarantees that ben-tal2015lectures :

[TABLE]

Note, that if $Q={\mathbb{R}}^{n}$, $\left\|\cdot\right\|=\left\|\cdot\right\|_{2}$ then $d({\vec{x}})=\frac{1}{2}\|{\vec{x}}\|_{2}^{2}$, $V[{\vec{x}}]({\vec{y}})=\frac{1}{2}\|{\vec{y}}-{\vec{x}}\|_{2}^{2}$,

[TABLE]

[TABLE]

that corresponds to the standard gradient-type iteration process (14).

Example (unit simplex). We have

[TABLE]

[TABLE]

For $k=0,...,N-1$, $i=1,...,n$

[TABLE]

The main result here is

[TABLE]

Note, that if we use $\left\|\cdot\right\|_{2}$-norm and $d\left({\vec{x}}\right)=\frac{1}{2}\left\|{{\vec{x}}-{\vec{x}}^{0}}\right\|_{2}^{2}$ here, we will have higher iteration complexity (2-norm projections on unit simplex) and

[TABLE]

Since typically $M_{2}={\rm O}\left({\sqrt{n}M_{\infty}}\right)$, it is worth to use $\left\|\cdot\right\|_{1}$-norm.

Assume now that $f(x)$ in (13) is additionally $\mu$-strongly convex in $\left\|\cdot\right\|_{2}$ norm:

[TABLE]

Let

[TABLE]

where

[TABLE]

Then Lacost-Julien2012

[TABLE]

Hence (see also Table 1),

[TABLE]

This bound is also un-improvable up to a constant factor nemirovsky1983problem ; nesterov2018lectures .

3 Looking into the Black-Box

In this section we consider how problem special structure can be used to solve non-smooth optimization problems with the convergence rate $O\left(\frac{1}{k}\right)$, which is faster than the lover bound $O\left(\frac{1}{\sqrt{k}}\right)$ for general class of non-smooth convex problems nemirovsky1983problem . Nevertheless, there is no contradiction as additional structure is used and we are looking inside the black-box.

3.1 Nesterov’s smoothing

In this subsection, following nesterov2005smooth , we consider the problem

[TABLE]

where $A:E_{1}\to E_{2}^{*}$ is a linear operator, $\phi({\vec{u}})$ is a continuous convex function on $Q_{2}$, $Q_{1},Q_{2}$ are convex compacts, $h$ is convex function with $L_{h}$-Lipschitz-continuous gradient.

Let us consider an example of $f(x)=\|A{\vec{x}}-{\vec{b}}\|_{\infty}$ with $A\in\mathbb{R}^{m\times n}$. Then,

[TABLE]

$h=0$, $E_{2}=\mathbb{R}^{m}$, $\phi({\vec{u}})=\langle{\vec{u}},{\vec{b}}\rangle$ and $Q_{2}$ is the ball in 1-norm.

The main idea of Nesterov is to add regularization inside the definition of $f$ in (21). More precisely, a prox-function $d_{2}({\vec{u}})$ (see definition in Section 2) is introduced for the set $Q_{2}$ and a smoothed counterpart $f_{\mu}({\vec{x}})$ for $f$ is defined as

[TABLE]

and ${\vec{u}}_{\mu}({\vec{x}})$ is the optimal solution of this maximization problem.

Theorem 3.1 (nesterov2005smooth )

The function $f_{\mu}({\vec{x}})$ is well defined, convex and continuously differentiable at any ${\vec{x}}\in E_{1}$ with $\nabla f_{\mu}({\vec{x}})=\nabla h({\vec{x}})+A^{*}{\vec{u}}_{\mu}({\vec{x}})$. Moreover, $\nabla f_{\mu}({\vec{x}})$ is Lipschitz continuous with constant $L_{\mu}=L_{h}+\frac{\|A\|_{1,2}^{2}}{\mu}$.

Here the adjoint operator $A^{*}$ is defined by equality $\langle A{\vec{x}},{\vec{u}}\rangle=\langle A^{*}{\vec{u}},{\vec{x}}\rangle$, ${\vec{x}}\in E_{1},{\vec{u}}\in E_{2}$ and the norm of the operator $\|A\|_{1,2}$ is defined by $\|A\|_{1,2}=\max_{{\vec{x}},{\vec{u}}}\{\langle A{\vec{x}},{\vec{u}}\rangle:\|{\vec{x}}\|_{E_{1}}=1,\|{\vec{u}}\|_{E_{2}}=1\}$.

Since $Q_{2}$ is bounded, $f_{\mu}({\vec{x}})$ is a uniform approximation for the function $f$, namely, for all ${\vec{x}}\in Q_{1}$,

[TABLE]

where $D_{2}=\max\{d_{2}({\vec{u}}):{\vec{u}}\in Q_{2}\}$.

Then, the idea is to choose $\mu$ sufficiently small and apply accelerated gradient method to minimize $f_{\mu}({\vec{x}})$ on $Q_{1}$. We use accelerated gradient method from dvurechensky2017adaptive ; dvurechensky2018computational which is different from the original method of nesterov2005smooth .

Theorem 3.2 (dvurechensky2017adaptive ; dvurechensky2018computational )

Let the sequences $\{{\vec{x}}^{k},{\vec{y}}^{k},{\vec{z}}^{k},\alpha_{k},C_{k}\}$, $k\geq 0$ be generated by Algorithm 1. Then, for all $k\geq 0$, it holds that

[TABLE]

Following the same steps as in the proof of Theorem 3 in nesterov2005smooth , we obtain

Theorem 3.3

Let Algorithm 1 be applied to minimize $f_{\mu}({\vec{x}})$ on $Q_{1}$ with $\mu=\frac{2\|A\|_{1,2}}{N+1}\sqrt{\frac{D_{1}}{D_{2}}}$, where $D_{1}=\max\{d_{1}({\vec{x}}):{\vec{x}}\in Q_{1}\}$. Then, after $N$ iterations, we have

[TABLE]

Proof

Applying Theorem 3.2 to $f_{\mu}$, and using (22), we obtain

[TABLE]

Substituting the value of $\mu$ from the theorem statement, we finish the proof.

A generalization of the smoothing technique for the case of non-compact sets $Q_{1},Q_{2}$, which is especially interesting when dealing with problems dual to problems with linear constraints, can be found in tran-dinh2015smooth . Ubiquitous entropic regularization of optimal transport cuturi2013sinkhorn can be seen as a particular case of the application of smoothing technique, especially in the context of Wasserstein barycenters cuturi2014fast ; uribe2018distributed ; dvurechensky2018decentralize .

3.2 Nemirovski’s Mirror Prox

In his paper nemirovski2004prox , Nemirovski considers problem (21) in the following form

[TABLE]

pointing to the fact that this problem is as general as (21). Indeed, the change of variables ${\vec{u}}\leftarrow({\vec{u}},t)$ and the feasible set $Q_{2}\leftarrow\{({\vec{u}},t):\min_{{\vec{u}}^{\prime}\in Q_{2}}\phi({\vec{u}}^{\prime})\leq t\leq\phi({\vec{u}})\}$ allows to make $\phi$ linear. His idea is to consider problem (29) directly as a convex-concave saddle point problem and associated weak variational inequality (VI).

[TABLE]

where the operator

[TABLE]

is monotone, i.e. $\langle\Phi({\vec{z}}_{1})-\Phi({\vec{z}}_{2}),{\vec{z}}_{1}-{\vec{z}}_{2}\rangle\geq 0$, and Lipschitz-continuous, i.e. $\|\Phi({\vec{z}}_{1})-\Phi({\vec{z}}_{2})\|_{*}\leq L\|{\vec{z}}_{1}-{\vec{z}}_{2}\|$. With the appropriate choice of norm on $E_{1}\times E_{2}$ and prox-function for $Q_{1}\times Q_{2}$, see Section 5 in nemirovski2004prox , the Lipschitz constant for $\Phi$ can be estimated as $L=2\|A\|_{1,2}\sqrt{D_{1}D_{2}}+L_{h}D_{1}$.

Theorem 3.4 (nemirovski2004prox )

Assume that $\Phi({\vec{z}})$ is monotone and $L$-Lipschitz-continuous. Then, for any $k\geq 1$ and any ${\vec{u}}\in Q$,

[TABLE]

Moreover, if the VI is associated with a convex-concave saddle point problem, i.e.

• •

$E=E_{1}\times E_{2}$,

• •

$Q=Q_{1}\times Q_{2}$* with convex compact sets $Q_{1}\subset E_{1}$, $Q_{2}\subset E_{2}$*

• •

\Phi({\vec{z}})=\Phi({\vec{x}},{\vec{u}})=\left(\begin{array}[]{c}\nabla_{{\vec{x}}}f({\vec{x}},{\vec{u}})\\ -\nabla_{{\vec{u}}}f({\vec{x}},{\vec{u}})\end{array}\right)* for a continuously differentiable function $f({\vec{x}},{\vec{u}})$ which is convex in ${\vec{x}}\in Q_{1}$ and concave in ${\vec{u}}\in Q_{2}$,*

then

[TABLE]

Choosing appropriately the norm in the space $E_{1}\times E_{2}$ and applying Mirror Prox algorithm to solve problem (29) as a saddle point problem, we obtain that the saddle point error in the l.h.s. of (35) decays as $\frac{2\|A\|_{1,2}\sqrt{D_{1}D_{2}}+L_{h}D_{1}}{k}$. This is slightly worse than the rate in (27) since the accelerated gradient method allows the faster decay for the smooth part $h({\vec{x}})$. An accelerated Mirror Prox method with the same rate as in (27) can be found in chen2017accelerated .

4 Non-Smooth Optimization in Large Dimensions

The optimization of non-smooth functionals with constraints attracts widespread interest in large-scale optimization and its applications ben-tal1997robust ; nesterov2014primal-dual . Subgradient methods for non-smooth optimization have a long history starting with the method for deterministic unconstrained problems and Euclidean setting in shor1967generalized and the generalization for constrained problems in polyak1967general , where the idea of steps switching between the direction of subgradient of the objective and the direction of subgradient of the constraint was suggested. Non-Euclidean extension, usually referred to as Mirror Descent, originated in nemirovskii1979efficient ; nemirovsky1983problem and was later analyzed in beck2003mirror . An extension for constrained problems was proposed in nemirovsky1983problem , see also recent version in beck2010comirror . To prove faster convergence rate of Mirror Descent for strongly convex objective in an unconstrained case, the restart technique nemirovskii1985optimal ; nemirovsky1983problem ; nesterov1983method was used in juditsky2012first-order . Usually, the step-size and stopping rule for Mirror Descent requires to know the Lipschitz constant of the objective function and constraint, if any. Adaptive step-sizes, which do not require this information, are considered in nemirovskii1979efficient for problems without inequality constraints, and in beck2010comirror for constrained problems.

Formally speaking, we consider the following convex constrained minimization problem

[TABLE]

where $X$ is a convex closed subset of a finite-dimensional real vector space $E$, $f:X\to\mathbb{R}$, $g:E\to\mathbb{R}$ are convex functions.

We assume $g$ to be a non-smooth Lipschitz-continuous function and the problem (3) to be regular. The last means that there exists a point $\bar{{\vec{x}}}$ in relative interior of the set $X$, such that $g(\bar{{\vec{x}}})<0$.

Note that, despite problem (36) contains only one inequality constraint, considered algorithms allow to solve more general problems with a number of constraints given as $\{g_{i}({\vec{x}})\leq 0,i=1,...,m\}$. The reason is that these constraints can be aggregated and represented as an equivalent constraint given by $\{g({\vec{x}})\leq 0\}$, where $g({\vec{x}})=\max_{i=1,...,m}g_{i}({\vec{x}})$.

We consider two adaptive Mirror Descent methods bayandina2018mirror for the problem (36). Both considered methods have complexity $O\left(\frac{1}{\varepsilon^{2}}\right)$ and optimal.

We consider algorithms, which are based on Mirror Descent method. Thus, we start with the description of proximal setup and basic properties of Mirror Descent step. Let $E$ be a finite-dimensional real vector space and $E^{*}$ be its dual. We denote the value of a linear function $g\in E^{*}$ at ${\vec{x}}\in E$ by $\langle{\vec{g}},{\vec{x}}\rangle$. Let $\|\cdot\|_{E}$ be some norm on $E$, $\|\cdot\|_{E,*}$ be its dual, defined by \|{\vec{g}}\|_{E,*}=\max\limits_{{\vec{x}}}\big{\{}\langle{\vec{g}},{\vec{x}}\rangle,\|{\vec{x}}\|_{E}\leq 1\big{\}}. We use $\nabla f({\vec{x}})$ to denote any subgradient of a function $f$ at a point ${\vec{x}}\in{\rm dom}f$.

Given a vector ${\vec{x}}\in X^{0}$, and a vector ${\vec{p}}\in E^{*}$, the Mirror Descent step is defined as

[TABLE]

We make the simplicity assumption, which means that $\mathrm{Mirr}[{\vec{x}}]({\vec{p}})$ is easily computable.

The following lemma ben-tal2015lectures describes the main property of the Mirror Descent step.

Lemma 1

Let $f$ be some convex function over a set $X$, $h>0$ be a step-size, ${\vec{x}}\in X^{0}$. Let the point ${\vec{x}}^{+}$ be defined by ${\vec{x}}^{+}=\mathrm{Mirr}[{\vec{x}}](h(\nabla f({\vec{x}})))$. Then, for any ${\vec{z}}\in X$,

[TABLE]

The following analog of Lemma 1 for $\delta$-subgradient $\nabla_{\delta}f$ holds.

Lemma 2

Let $f$ be some convex function over a set $X$, $h>0$ be a step-size, ${\vec{x}}\in X^{0}$. Let the point ${\vec{x}}^{+}$ be defined by ${\vec{x}}^{+}=\mathrm{Mirr}[{\vec{x}}](h\cdot(\nabla_{\delta}f({\vec{x}})))$. Then, for any ${\vec{z}}\in X$,

[TABLE]

We consider problem (36) in two different settings, namely, non-smooth Lipschitz-continuous objective function $f$ and general objective function $f$, which is not necessarily Lipschitz-continuous, e.g. a quadratic function. In both cases, we assume that $g$ is non-smooth and is Lipschitz-continuous

[TABLE]

Let ${\vec{x}}_{*}$ be a solution to (36). We say that a point $\tilde{{\vec{x}}}\in X$ is an $\varepsilon$-solution to (36) if

[TABLE]

All considered in this section methods (Algorithms 3 and 4) are applicable in the case of using $\delta$-subgradient instead of usual subgradient. For this case we can get an $\varepsilon$-solution $\tilde{{\vec{x}}}\in X$:

[TABLE]

The methods we describe are based on the of Polyak’s switching subgradient method polyak1967general for constrained convex problems, also analyzed in nesterov2010introduction , and Mirror Descent method originated in nemirovsky1983problem ; see also nemirovskii1979efficient .

4.1 Convex Non-Smooth Objective Function

In this subsection, we assume that $f$ is a non-smooth Lipschitz-continuous function

[TABLE]

Let ${\vec{x}}_{*}$ be a solution to (36) and assume that we know a constant $\Theta_{0}>0$ such that

[TABLE]

For example, if $X$ is a compact set, one can choose $\Theta_{0}^{2}=\max_{{\vec{x}}\in X}d({\vec{x}})$.

Theorem 4.1

Assume that inequalities (40) and (43) hold and a known constant $\Theta_{0}>0$ is such that $d({\vec{x}}_{*})\leq\Theta_{0}^{2}$. Then, Algorithm 3 stops after not more than

[TABLE]

iterations and $\bar{{\vec{x}}}^{k}$ is an $\varepsilon$-solution to (36) in the sense of (41).

Let us now show that Algorithm 3 allows to reconstruct an approximate solution to the problem, which is dual to (36). We consider a special type of problem (36) with $g$ given by

[TABLE]

Then, the dual problem to (36) is

[TABLE]

where $\lambda_{i}\geq 0,i=1,...,m$ are Lagrange multipliers.

We slightly modify the assumption (44) and assume that the set $X$ is bounded and that we know a constant $\Theta_{0}>0$ such that

[TABLE]

As before, denote $[k]=\{j\in\{0,...,k-1\}\}$, $J=[k]\setminus I$. Let $j\in J$. Then a subgradient of $g({\vec{x}})$ is used to make the $j$-th step of Algorithm 3. To find this subgradient, it is natural to find an active constraint $i\in 1,...,m$ such that $g({\vec{x}}^{j})=g_{i}({\vec{x}}^{j})$ and use $\nabla g({\vec{x}}^{j})=\nabla g_{i}({\vec{x}}^{j})$ to make a step. Denote $i(j)\in 1,...,m$ the number of active constraint, whose subgradient is used to make a non-productive step at iteration $j\in J$. In other words, $g({\vec{x}}^{j})=g_{i(j)}({\vec{x}}^{j})$ and $\nabla g({\vec{x}}^{j})=\nabla g_{i(j)}({\vec{x}}^{j})$. We define an approximate dual solution on a step $k\geq 0$ as

[TABLE]

and modify Algorithm 3 to return a pair $(\bar{{\vec{x}}}^{k},\bar{\lambda}^{k})$.

Theorem 4.2

Assume that the set $X$ is bounded, the inequalities (40) and (43) hold and a known constant $\Theta_{0}>0$ is such that $d({\vec{x}}_{*})\leq\Theta_{0}^{2}$. Then, modified Algorithm 3 stops after not more than

[TABLE]

iterations and the pair $(\bar{{\vec{x}}}^{k},\bar{\lambda}^{k})$ returned by this algorithm satisfies

[TABLE]

Now we consider an interesting example of huge-scale problem nesterov2014subgradient ; nesterov2014primal-dual with a sparse structure. We would like to illustrate two important ideas. Firstly, consideration of the dual problem can simplify the solution, if it is possible to reconstruct the solution of the primal problem by solving the dual problem. Secondly, for a special sparse non-smooth piece-wise linear functions we suggest a very efficient implementation of one subgradient iteration nesterov2014subgradient . In such cases simple subgradient methods (for example, Algorithm 3) can be useful due to the relatively inexpensive cost of iterations.

Recall (see e.g. nesterov2014primal-dual ) that Truss Topology Design problem consists in finding the best mechanical structure resisting to an external force with an upper bound for the total weight of construction. Its mathematical formulation looks as follows:

[TABLE]

where $\overline{{\vec{f}}}$ is a vector of external forces, ${\vec{z}}\in R^{2n}$ is a vector of virtual displacements of $n$ nodes in $R^{2}$, ${\vec{w}}$ is a vector of $m$ bars, and $T$ is the total weight of construction. The compliance matrix $A({\vec{w}})$ has the following form:

[TABLE]

where ${\vec{a}}_{i}\in R^{2n}$ are the vectors describing the interactions of two nodes connected by an arc. These vectors are very sparse: for 2D-model they have at most 4 nonzero elements.

Let us rewrite the problem (50) as a Linear Programming problem.

[TABLE]

Note that for the inequality in the third line we do not need any assumption.

Denote by ${\vec{y}}^{*}$ the optimal solution of the optimization problem in the brackets. Then there exist multipliers ${\vec{x}}^{*}\in R^{m}_{+}$ such that

[TABLE]

where $J_{+}=\{i:\,\langle{\vec{a}}_{i},{\vec{y}}^{*}\rangle=1\}$, and $J_{-}=\{i:\,\langle{\vec{a}}_{i},{\vec{y}}^{*}\rangle=-1\}$. Multiplying the first equation in (52) by ${\vec{y}}^{*}$, we get

[TABLE]

Note that the first equation in (52) can be written as

[TABLE]

Let us reconstruct now the solution of the primal problem. Denote

[TABLE]

Then, in view of (54) we have $\overline{{\vec{f}}}=A({\vec{w}}^{*}){\vec{z}}^{*}$, and ${\vec{w}}^{*}\in\triangle(T)$. Thus, the pair (55) is feasible for the primal problem. On the other hand,

[TABLE]

Thus, the duality gap in the chain (51) is zero, and the pair $({\vec{w}}^{*},{\vec{z}}^{*})$, defined by (55) is the optimal solution of the primal problem.

The above discussion allows us to concentrate on the following (dual) Linear Programming problem:

[TABLE]

which we can solve by the primal-dual Algorithm 3.

Assume that we have local truss: each node is connected only with few neighbors. It allows to apply the property of sparsity for vectors ${\vec{a}}_{i}$ ($1\leq i\leq m$). In this case the computational cost of each iteration grows as $O(\log_{2}m)$ nesterov2014subgradient ; nesterov2014primal-dual .

In nesterov2014subgradient a special class of huge-scale problems with sparse subgradient was considered. According to nesterov2014subgradient for smooth functions this is a very rare feature. For example, for quadratic function $f({\vec{y}})=\frac{1}{2}\langle A{\vec{y}},{\vec{y}}\rangle$ the gradient $\nabla f({\vec{y}})=A{\vec{y}}$ usually is dense even for a sparse matrix $A$.

However, the subgradient of non-smooth function $f({\vec{y}})=\max_{1\leq i\leq m}\langle{\vec{a}}_{i},{\vec{y}}\rangle$ (see (56) above) are sparse provided that all vectors ${\vec{a}}_{i}$ share this property. This fact is based on the following observation. For the function $f({\vec{y}})=\max_{1\leq i\leq m}\langle{\vec{a}}_{i},{\vec{y}}\rangle$ with sparse matrix $A=({\vec{a}}_{1},{\vec{a}}_{2},...,{\vec{a}}_{m})$ the vector $\nabla f({\vec{y}})={\vec{a}}_{i({\vec{y}})}$ is a subgradient at point ${\vec{y}}$. Then the standard subgradient step

[TABLE]

changes only a few entries of vector ${\vec{y}}$ and the vector ${\vec{z}}_{+}=A^{T}{\vec{y}}_{+}$ differs from ${\vec{z}}=A^{T}{\vec{y}}$ also in a few positions only. Thus, the function value $f({\vec{y}}_{+})$ can be easily updated provided that we have an efficient procedure for recomputing the maximum of $m$ values.

Note the objective functional in (56) is linear and the costs of iteration of Algorithm 3 and considered in nesterov2014subgradient switching simple subgradient scheme is comparable. At the same time, the step productivity condition is simpler for Algorithm 3 as considered in nesterov2014subgradient switching subgradient scheme. Therefore main observations for nesterov2014subgradient are correct for Algorithm 3.

4.2 General Convex and Quasi-Convex Objective Functions

In this subsection, we assume that the objective function $f$ in (36) might not satisfy (43) and, hence, its subgradient could be unbounded. One of the examples is a quadratic function. We also assume that inequality (44) holds.

We further consider ideas in nesterov2010introduction ; nesterov2015subgradient and adapt them for problem (36), in a way that our algorithm allows to use non-Euclidean proximal setup, as does Mirror Descent, and does not require to know the constant $M_{g}$. Following nesterov2010introduction , given a function $f$ for each subgradient $\nabla f({\vec{x}})$ at a point ${\vec{y}}\in X$, we define

[TABLE]

The following result gives complexity estimate for Algorithm 4 in terms of $v_{f}[{\vec{x}}_{*}]({\vec{x}})$. Below we use this theorem to establish complexity result for smooth objective $f$.

Theorem 4.3

Assume that inequality (40) holds and a known constant $\Theta_{0}>0$ is such that $d({\vec{x}}_{*})\leq\Theta_{0}^{2}$. Then, Algorithm 4 stops after not more than

[TABLE]

iterations and it holds that $\min_{i\in I}v_{f}[{\vec{x}}_{*}]({\vec{x}}^{i})\leq\varepsilon$ and $g(\bar{{\vec{x}}}^{k})\leq\varepsilon$.

To obtain the complexity of our algorithm in terms of the values of the objective function $f$, we define non-decreasing function

[TABLE]

and use the following lemma from nesterov2010introduction .

Lemma 3

Assume that $f$ is a convex function. Then, for any ${\vec{x}}\in X$,

[TABLE]

Corollary 1

Assume that the objective function $f$ in (36) is given as $f({\vec{x}})=\max_{i\in\{1,...,m\}}f_{i}({\vec{x}})$, where $f_{i}({\vec{x}})$, $i=1,...,m$ are differentiable with Lipschitz-continuous gradient

[TABLE]

Then $\bar{{\vec{x}}}^{k}$ is $\widetilde{\varepsilon}$-solution to (36) in the sense of (41), where

[TABLE]

Remark 1

According to nesterov1989 ; nesterov2018lectures main lemma 3 holds for quasi-convex objective functions Polyak1969 too:

[TABLE]

This means that results of this subsection are valid for quasi-convex objectives.

Remark 2

In view of the Lipschitzness and, generally speaking, non-smoothness of functional limitations, the obtained estimate for the number of iterations means that the proposed method is optimal from the point of view of oracle evaluations: $O\left(\frac{1}{\varepsilon^{2}}\right)$ iterations are sufficient for achieving the required accuracy $\varepsilon$ of solving the problem for the class of target functionals considered in this section of the article. Note also that the considered algorithm 3 applies to the considered classes of problems with constraints with convex objective functionals of different smoothness levels. However, the non-fulfillment, generally speaking, of the Lipschitz condition for the objective functional $f$ does not allow one to substantiate the optimality of the algorithms 3 in the general situation (for example, with a Lipschitz-continuous gradient). More precisely, situations are possible when the productive steps of the norm (sub)gradients of the objective functional $\|\nabla f({\vec{x}}^{k})\|_{*}$ are large enough and this will interfere with the speedy achievement of the stopping criterion of the 3.

5 Universal Methods

In this section we consider problem

[TABLE]

where $Q$ is a convex set and $f$ is a convex function with Hölder-continuous subgradient

[TABLE]

with $\nu\in[0,1]$. The case $\nu=0$ corresponds to non-smooth optimization and the case $\nu=1$ corresponds to smooth optimization. The goal of this section is to present the Universal Accelerated Gradient method first proposed by Nesterov nesterov2015universal . This method is a black-box method which does not require the knowledge of constants $\nu,L_{\nu}$ and works in accordance with the lower complexity bound $O\left(\left(\frac{L_{\nu}R^{1+\nu}}{{\epsilon}}\right)^{\frac{2}{1+3\nu}}\right)$ obtained in nemirovsky1983problem .

The main idea is based on the observation that a non-smooth convex function can be upper bounded by a quadratic objective function slightly shifted above. More precisely, for any ${\vec{x}},{\vec{y}}\in Q$,

[TABLE]

where

[TABLE]

The next idea is to apply an accelerated gradient method with backtracking procedure to adapt for the unknown $L(\delta)$ with appropriately chosen $\delta$. The method we present is based on accelerated gradient method from dvurechensky2017adaptive ; dvurechensky2018computational and, thus is different from the original method of nesterov2015universal .

Inequality (64) guarantees that the backtracking procedure in the inner cycle is finite.

Theorem 5.1 (nesterov2015universal )

Let $f$ satisfy (63). Then,

[TABLE]

Moreover, the number of oracle calls is bounded by

[TABLE]

Translating this rate of convergence to the language of complexity, we obtain that to obtain a solution with an accuracy ${\epsilon}$ the number of iterations is no more than

[TABLE]

i.e. is optimal.

In his paper, Nesterov considers a more general composite optimization problem

[TABLE]

where $h$ is a simple convex function, and obtains the same complexity guarantees. Universal methods were extended for the case of strongly convex problems by a restart technique in roulet2017sharpness , for non-convex optimization in ghadimi2015generalized and for the case of non-convex optimization with inexact oracle in dvurechensky2017gradient . As we can see from (64), universal accelerated gradient method is connected to smooth problems with inexact oracle. The study of accelerated gradient methods with inexact oracle was first proposed in aspremont2008smooth and was very well developed in devolder2014first ; dvurechensky2016stochastic ; bogolubsky2016learning ; dvurechensky2017gradient including stochastic optimization problems and strongly convex problems. A universal method with inexact oracle can be found in dvurechensky2017universal . Experiments show nesterov2015universal that universal method accelerates to $O\left(\frac{1}{k}\right)$ rate for non-smooth problems with a special ”smoothing friendly” (see Section 3) structure. This is especially interesting for traffic flow modeling problems, which possess such structure baimurzina2017universal .

Now we consider universal analog of A.S. Nemirovski’s proximal mirror method for variational inequalities with a Holder-continuous operator. More precisely, we consider universal extension of Algorithm 2 which allows to solve smooth and non-smooth variational inequalities without the prior knowledge of the smoothness. Main idea of the this method is the adaptive choice of constants and level of smoothness in minimized prox-mappings at each iteration. These constants are related to the Hölder constant of the operator and this method allows to find a suitable constant at each iteration.

Theorem 5.2 (dvurechensky2018generalized )

For any $k\geq 1$ and any ${\vec{z}}\in Q$,

[TABLE]

Note that if $\max_{{\vec{z}}\in Q}V[{\vec{z}}^{0}]({\vec{z}})\leq D$, we can construct the following adaptive stopping criterion for our algorithm

[TABLE]

Next, we consider the case of Hölder-continuous operator $\Phi$ and show that Algorithm 6 is universal. Assume for some $\nu\in[0,1]$ and $L_{\nu}\geq 0$

[TABLE]

holds. The following inequality is a generalization of (64) for VI. For any ${\vec{x}},{\vec{y}},{\vec{z}}\in Q$ and $\delta>0$,

[TABLE]

[TABLE]

where

[TABLE]

So, we have

[TABLE]

Let us consider estimates of the necessary number of iterations are obtained to achieve a given quality of the variational inequality solution.

Corollary 2 (Universal Method for VI)

Assume that the operator $\Phi$ is Hölder continuous with constant $L_{\nu}$ for some $\nu\in[0,1]$ and $M_{-1}\leq\left(\frac{2}{\varepsilon}\right)^{\frac{1-\nu}{1+\nu}}L_{\nu}^{\frac{2}{1+\nu}}$. Also assume that the set $Q$ is bounded. Then, for all $k\geq 0$, we have

[TABLE]

As it follows from (77), if $M_{k}\geq L(\frac{\varepsilon}{2})$, (74) holds. Thus, for all $i=0,...,k-1$, we have $M_{i}\leq 2\cdot L(\frac{\varepsilon}{2})$ and

[TABLE]

(78) holds. Here $L(\cdot)$ is defined in (76). ∎

Let us add some remarks.

Remark 3

Since the algorithm does not use the values of parameters $\nu$ and $L_{\nu}$, we obtain the following iteration complexity bound

[TABLE]

to achieve

[TABLE]

Using the same reasoning as in nesterov2015universal , we estimate the number of oracle calls for Algorithm 6. The number of oracle calls on each iteration $k$ is equal to $2i_{k}$. At the same time, $M_{k}=2^{i_{k}-2}M_{k-1}$ and, hence, $i_{k}=2+\log_{2}\frac{M_{k}}{M_{k-1}}$. Thus, the total number of oracle calls is

[TABLE]

where we used that $M_{k}\leq 2L(\frac{\varepsilon}{2})$.

Thus, the number of oracle calls of the Algorithm 6 does not exceed:

[TABLE]

Remark 4

We can apply this method to convex-concave saddle problems of the form

[TABLE]

where $Q_{1,2}$ are convex compacts in $\mathbb{R}^{n}$, $f$ is convex in ${\vec{x}}$ and concave in ${\vec{y}}$, there is $\nu\in[0,1]$ and constants $L_{11,\nu},L_{12,\nu},L_{21,\nu},L_{22,\nu}\textless+\infty$:

[TABLE]

[TABLE]

for all ${\vec{x}},{\vec{x}}+\Delta{\vec{x}}\in Q_{1},{\vec{y}},{\vec{y}}+\Delta{\vec{y}}\in Q_{2}$.

It is possible to achieve an acceptable approximation $(\widehat{{\vec{x}}},\widehat{{\vec{y}}})\in Q_{1}\times Q_{2}$:

[TABLE]

for the saddle point $({\vec{x}}_{*},{\vec{y}}_{*})\in Q_{1}\times Q_{2}$ of the (80) problem in no more than

[TABLE]

iterations, which indicates the optimality of the proposed method, at least for $\nu=0$ and $\nu=1$. However, in practice experiments show that (81) can be achieved much faster due to the adaptability of the method.

6 Concluding remarks

Modern numerical methods for non-smooth convex optimization problems are typically based on the structure of the problem. We start with one of the most powerful example of such type. For geometric median search problem there exists efficient method that significantly outperform described above lower complexity bounds cohen2016geometric . In Machine Learning we typically meet the problems with hidden affine structure and small effective dimension (SVM) that allow us to use different smoothing techniques allen2016optimal . Description of one of these techniques (Nesterov’s smoothing technique) one can find in this survey. The other popular technique is based on averaging of the function around the small ball with the center at the point in consideration duchi2012randomized . A huge amount of data since applications lead to composite optimization problems with non smooth composite (LASSO). For this class of problems accelerated (fast) gradient methods are typically applied beck2009fast , nesterov2013gradient , lan2016gradient . This approach (composite optimization) have been recently expanded for more general class of problems tyurin2017fast . In different Image Processing applications one can find a lot of non-smooth problems formulations with saddle-point structure. That is the goal function has Legendre representation. In this case one can apply special versions of accelerated (primal-dual) methods chambolle2011first-order , chen2014optimal , lan2016accelerated . Universal Mirror Prox method described above demonstrates the alternative approach which can be applied in rather general context. Unfortunately, the most of these tricks have proven to be beyond the scope of this survey. But we include in the survey the description of the Universal Accelerated Gradient Descent algorithm tyurin2017fast which in the general case can also be applied to a wide variety of problems.

Another important direction in Non-smooth Convex Optimization is huge-scale optimization for sparse problems nesterov2014subgradient . The basic idea that reduce huge dimension to non-smoothness is as follows:

[TABLE]

is equivalent to the single non-smooth constraint:

[TABLE]

We demonstrated this idea above on Truss Topology Design example.

One should note that we concentrate in this survey only on deterministic convex optimization problems, but the most beautiful things in non smooth optimization is that stochasticity nemirovsky1983problem , duchi2016introductory , juditsky2011firstI , juditsky2011firstII and online context hazan2016introduction in general doesn’t change (up to a logarithmic factor in the strongly convex case) anything in complexity estimates. As an example, of stochastic (randomized) approach one can mentioned the work anikin2015efficient where one can find reformulation of Google problem as non smooth convex optimization problem. Special randomized Mirror Descent algorithm allows to solve this problem almost independently on the number of vertexes.

Finally, let’s note that in the decentralized distributed non smooth (stochastic) convex optimization for the last few years there appear optimal methods lan2017communication , uribe2017optimal , bubeckoptimal .

Acknowledgements.

The article was supported in its major parts by the grant 18-29-03071 mk from Russian Foundation for Basic Research. E. Nurminski acknowledges the partial support from the project 1.7658.2017/6.7 of Ministry of Science and Higher Professional Education in Section 1. The work of A. Gasnikov, P. Dvurechensky and F. Stonyakin in Section 3 was partially supported by Russian Foundation for Basic Research grant 18-31-20005 mol_a_ved. The work of F. Stonyakin in Subsection 4.1 was supported by Russian Science Foundation grant 18-71-00048.