Dual subgradient algorithms for large-scale nonsmooth learning problems

Introduction

The problem of interest in this paper is a convex optimization problem in the form

where $X$ is a nonempty closed and bounded subset of Euclidean space $E_{x}$, and $f_{*}$ is concave and Lipschitz continuous function on $X$. We are interested in the situation where the sizes of the problem put it beyond the “practical grasp” of polynomial time interior point methods with their rather computationally expensive in the large scale case iterations. In this case the methods of choice are the First Order (FO) optimization techniques. The state of the art of these techniques can be briefly summarized as follows: $\bullet$ The most standard FO approach to ($P$) requires to provide $X,\;E_{x}$ with proximal setup $\|\cdot\|,\omega(\cdot)$, that is to equip the space $E_{x}$ with a norm $\|\cdot\|$, and the domain $X$ of the problem – with a strongly convex, modulus 1, w.r.t. $\|\cdot\|$ distance-generating function (d.-g.f.) $\omega(\cdot)$ with variation over $X$ bounded by some $\Omega^{2}_{X}$. After such a setup is fixed, generating an $\epsilon$-solution to the problem (i.e., a point $x_{\epsilon}\in X$ satisfying $\hbox{\rm Opt}(P)-f_{*}(x_{\epsilon})\leq\epsilon$) costs at most $N(\epsilon)$ steps, where

$N(\epsilon)=O(1){\Omega_{X}^{2}L^{2}\over\epsilon^{2}}$ in the nonsmooth case, where $f_{*}$ is Lipschitz continuous, with constant $L$ w.r.t. $\|\cdot\|$ (Mirror Descent (MD) algorithm, see, e.g., [11, Chapter 5]), and

$N(\epsilon)=O(1){\Omega_{X}\sqrt{{\cal L}\over\epsilon}}$ in the smooth case, where $f_{*}$ possesses Lipschitz continuous, with constant ${\cal L}$, gradient: $\|f_{*}^{\prime}(x)-f_{*}^{\prime}(x^{\prime})\|_{*}\leq{\cal L}\|x-x^{\prime}\|$, where $\|\cdot\|_{*}$ is the norm conjugate to $\|\cdot\|$ (Nesterov’s algorithm for smooth convex optimization, see, e.g., ). Here and in the sequel $O(1)$’s stand for positive absolute constants. Note that in the large scale case, these convergence rates are the best allowed under circumstances by Information-Based Complexity Theory (for details, see, e.g., ).

A step of a FO method essentially reduces to computing $f_{*}$, $f_{*}^{\prime}$ at a point and computing prox-mapping (“Bregman projection”) $(x\in X,\xi\in\ E_{x})\mapsto\mathop{\hbox{\rm argmin}}_{u\in X}\left\{\omega(u)+\langle\xi-\omega^{\prime}(x).u\rangle\right\}$, construction originating from J.-J. Moreau and L. Bregman . $\bullet$ A different way of processing ($P$) by FO methods, originating in the breakthrough paper of Nesterov , is to use Fenchel-type representation of $f_{*}$:

where $Y$ is a closed and bounded subset of Euclidean space $E_{y}$, $A:\;E_{y}\mapsto E_{x}$ is a linear mapping, and $\psi(y)$ is a convex function. Representations of this type are readily available for a wide family of “well-structured” nonsmooth objectives $f_{*}$; moreover, usually we can make $\psi$ to possess Lipschitz continuous gradient or even to be linear (for instructive examples, see, e.g., or [11, Chapter 6]). Whenever this is the case, and given proximal setups $(\|\cdot\|_{x},\omega_{x}(\cdot))$ and $(\|\cdot\|_{y},\omega_{y}(\cdot))$ for $(E_{x},X)$ and for $(E_{y},Y)$, proximal type algorithms like Nesterov’s smoothing and Dual Extrapolation , or Mirror Prox allow to find an $\epsilon$-solution to ($P$) in $O(1/\epsilon)$ proximal type iterations, the factors hidden in $O(\cdot)$ being explicitly given functions of the variations $\Omega_{X}^{2}$ of $\omega_{x}$ and $\Omega_{Y}^{2}$ of $\omega_{y}$ and of the partial Lipschitz constants of $\nabla F$ w.r.t. $x$ and to $y$.

Clearly, to be practical, methods of the outlined type should rely on “good” proximal setups – those resulting in “moderate” values of $\Omega_{X}$ and $\Omega_{Y}$ and not too difficult to compute prox-mappings, associated with $\omega_{X}$ and $\omega_{Y}$. This is indeed the case for domains $X$ arising in numerous applications (for instructive examples, see, e.g., [11, Chapter 5]). The question addressed in this paper is what to do when one of the domains, namely, $X$ does not admit a “good” proximal setup. Here are two instructive examples:

$X$ is the unit ball of the nuclear norm $\|\sigma(\cdot)\|_{1}$ in the space ${\mathbf{R}}^{p\times q}$ of $p\times q$ matrices (from now on, for a $p\times q$ matrix $x$, $\sigma(x)=[\sigma_{1}(x);...;\sigma_{\min[p,q]}(x)]$ denotes the vector comprised by singular values of $x$ taken in the non-ascending order). This domain arises in various low-rank-oriented problems of matrix recovery. In this case, $X$ does admit a proximal setup with $\Omega_{X}=O(1)\sqrt{\ln(pq)}$. However, computing prox-mapping involves full singular value decomposition (SVD) of a $p\times q$ matrix and becomes prohibitively time-consuming when $p,q$ are in the range of tens of thousand. Note that this hardly is a shortcoming of the existing proximal setups, since already computing the nuclear norm (that is, checking the inclusion $x\in X$) requires computing the SVD of $x$.

Note that whenever a prox-mapping associated with $X$ is “easy to compute,” it is equally easy to maximize over $X$ a linear form (since ${\hbox{}\rm Prox}_{x}(-t\xi)$ converges to the maximizer of $\langle\xi,x\rangle$ over $X$ as $t\to\infty$). In such a case, we have at our disposal an efficient Linear Optimization (LO) oracle – a routine which, given on input a linear form $\xi$, returns a point $x_{X}(\xi)\in\mathop{\hbox{\rm Argmax}}_{x\in X}\langle\xi,x\rangle$. This conclusion, however, cannot be reversed – our abilities to maximize, at a reasonable cost, linear functionals over $X$ does not imply the possibility to compute a prox-mapping at a comparable cost. For example, when $X\in{\mathbf{R}}^{p\times q}$ is the unit ball of the nuclear norm, maximizing a linear function $\langle\xi,x\rangle=\hbox{\rm Tr}(\xi x^{T})$ over $X$ requires finding the largest singular value of a $p\times q$ matrix and associated left singular vector. For large $p$ and $q$, solving the latter problem is by orders of magnitude cheaper than computing full SVD of a $p\times q$ matrix. This and similar examples motivate the interest, especially in Machine Learning community, in optimization techniques solving ($P$) via an LO oracle for $X$. In particular, the only “classical” technique of this type – the Conditional Gradient (CG) algorithm going back to Frank and Wolfe – has attracted much attention recently. In the setting of CG method it is assumed that $f$ is smooth (with Hölder continuous gradient), and the standard result here (which is widely known, see, e.g., ) is the following.

Let $X$ be a closed and bounded convex set in a Euclidean space $E_{x}$ such that $X$ linearly spans $E_{x}$. Assume that we are given an LO oracle for $X$, and let $f_{*}$ be a concave continuously differentiable function on $X$ such that for some ${\cal L}<\infty$ and $q\in(1,2]$ one has

where $\|\cdot\|_{X}$ is the norm on $E_{x}$ with the unit ball ${1\over 2}[X-X]$. Consider a recurrence of the form

where $x_{X}(\xi)\in\mathop{\hbox{\rm Argmax}}_{x\in X}\langle\xi,x\rangle$ and $x_{1}\in X$. Then for all $t=2,3,...$ one has

Assuming an LO oracle for $X$ available, the major limitation in solving ($P$) by the Conditional Gradient method is the requirement for problem objective $f_{*}$ to be smooth (otherwise, there are no particular requirements to the problem geometry). What to do if this requirement is not satisfied? In this paper, we investigate two simple options for processing this case, based on Fenchel-type representation (1) of $f_{*}$ which we assume to be available. Fenchel-type representations are utilized in various ways by a wide spectrum of duality-based convex optimization algorithms (see, e.g., and references therein). Here we primarily focus on “nonsmooth” case, where the representation involves a Lipschitz continuous convex function $\psi$ given by a First Order oracle. Besides this, we assume that $Y$ (but not $X$!) does admit a proximal setup $(\|\cdot\|_{y},\omega_{y}(\cdot))$. In this case, we can pass from the problem of interest ($P$) to its dual

Clearly, the LO oracle for $X$ along with the FO oracle for $\psi$ provide a FO oracle for $(D)$:

Since $Y$ admits a proximal setup, this is enough to allow to get an $\epsilon$-solution to $(D)$ in $N(\epsilon)=O(1){L^{2}\Omega^{2}_{Y}\over\epsilon^{2}}$ steps, $L$ being the Lipschitz constant of $f$ w.r.t. $\|\cdot\|_{y}$. Whatever slow the resulting rate of convergence could look, we shall see in the mean time that there are important applications where this rate seems to be the best known so far. When implementing the outlined scheme, the only nontrivial question is how to recover a good optimal solution to the problem $(P)$ of actual interest from a good approximate solution to its dual problem $(D)$. The proposed answer to this question stems from the pretty simple at the first glance machinery of accuracy certificates proposed recently in , and closely related to the work . The summary of our approach is as follows. When solving $(D)$ by a FO method, we generate search points $y_{\tau}\in Y$ where the subgradients $f^{\prime}(y_{\tau})$ of $f$ are computed; as a byproduct of the latter computation, we have at our disposal the points $x_{\tau}=x(y_{\tau})$. As a result, after $t$ steps we have at our disposal execution protocol $y^{t}=\{y_{\tau},f^{\prime}(y_{\tau})\}_{\tau=1}^{t}$. An accuracy certificate associated with this protocol is, by definition, a collection $\lambda^{t}=\{\lambda^{t}_{\tau}\}_{\tau=1}^{t}$ of nonnegative weights $\lambda^{t}_{\tau}$ summing up to 1: $\sum_{\tau=1}^{t}\lambda^{t}_{\tau}=1$. The resolution of the certificate is, by definition, the quantity

An immediate observation is (see section 2) that setting $\widehat{y}^{t}=\sum_{\tau=1}^{t}\lambda^{t}_{\tau}y_{\tau}$, $\widehat{x}^{t}=\sum_{\tau=1}^{t}\lambda^{t}_{\tau}x_{\tau}$, we get a pair of feasible solutions to $(D)$ and to $(P)$ such that

Thus, assuming that the FO method in question produces, in addition to search points, accuracy certificates for the resulting execution protocols and that the resolution of these certificates goes to 0 as $t\to\infty$ at some rate, we can use the certificates to build feasible approximate solutions to $(D)$ and to $(P)$ with nonoptimalities, in terms of the objectives of the respective problems, going to 0, at the same rate, as $t\to\infty$.

The scope of the outlined approach depends on whether we are able to equip known methods of nonsmooth convex minimization with computationally cheap mechanisms for building “good” accuracy certificates. The meaning of “good” in this context is exactly that the rate of convergence of the corresponding resolution to 0 is identical to the standard efficiency estimates of the methods (e.g., for MD this would mean that $\epsilon(y^{t},\lambda^{t})\leq O(1)L\Omega_{Y}t^{-1/2}$). provides a positive answer to this question for the most attractive academically polynomial time oracle-oriented algorithms for convex optimization, like the Ellipsoid method. These methods, however, usually are poorly suited for large-scale applications. In this paper, we provide a positive answer to the above question for the three most attractive oracle-oriented FO methods for large-scale nonsmooth convex optimization known to us. Specifically, we consider $\bullet$ MD (where accuracy certificates are easy to obtain, see also ), $\bullet$ Full Memory Mirror Descent Level (MDL) method (a Mirror Descent extension of the Bundle-Level method ; to the best of our knowledge, this extension was not yet described in the literature), and $\bullet$ Non-Euclidean Restricted Memory Level method (NERML) originating from , which we believe is the most attractive tool for large-scale nonsmooth oracle-based convex optimization. To the best of our knowledge, equipping NERML with accuracy certificates is a novel development.

We also consider a different approach to non-smooth convex optimization over a domain given by LO oracle, approach mimicking Nesterov’s smoothing . Specifically, assuming, as above, that $f_{*}$ is given by Fenchel-type representation (1) with $Y$ admitting a proximal setup, we use this setup, exactly in the same way as in , to approximate $f_{*}$ by a smooth function which then is minimized by the CG algorithm. Therefore, the only difference with is in replacing Nesterov’s optimal algorithm for smooth convex optimization (which requires a good proximal point setup for $X$) with although slower, but less demanding (just LO oracle for $X$ is enough) CG method. We shall see in the mean time that, unsurprisingly, the theoretical complexity of the two outlined approaches – “nonsmooth” and “smoothing” ones – are essentially the same.

The main body of the paper is organized as follows. In section 2, we develop the components of the approach related to duality and show how an accuracy certificate with small resolution yields a pair of good approximate solutions to $(P)$ and $(D)$. In section 3, we show how to equip the MD, MDL and NERML algorithms with accuracy certificates. In section 4, we investigate the outlined smoothing approach. In section 5, we consider examples, primarily of Machine Learning origin, where we prone the usage of the proposed algorithms. Section 6 reports some preliminary numerical results. Some technical proofs are relegated to the appendix.

Duality and accuracy certificates

Let $E_{x}$ be a Euclidean space, $X\subset E_{x}$ be a nonempty closed and bounded convex set equipped with LO oracle – a procedure which, given on input $\xi\in E_{x}$, returns a maximizer $x_{X}(\xi)$ of the linear form $\langle\xi,x\rangle$ over $x\in X$. Let $f_{*}(x)$ be a concave function given by Fenchel-type representation:

where $Y$ is a convex compact subset of a Euclidean space $E_{y}$ and $\psi$ is a Lipschitz continuous convex function on $Y$ given by a First Order oracle.

By the standard saddle point argument, we have $\hbox{\rm Opt}(P)=\hbox{\rm Opt}(D)$.

2 Main observation

Observe that the First Order oracle for $\psi$ along with the LO oracle for $X$ provide a First Order oracle for $(D)$; specifically, the vector field

where $\psi^{\prime}(y)\in\partial\psi(y)$ is a subgradient field of $f$.

Consider a collection $y^{t}=\{y_{\tau}\in Y,f^{\prime}(y_{\tau})\}_{\tau=1}^{t}$ along with a collection $\lambda^{t}=\{\lambda_{\tau}\geq 0\}_{\tau=1}^{t}$ such that $\sum_{\tau=1}^{t}\lambda_{\tau}=1$, and let us set

In the sequel, the components $y_{\tau}$ of $y^{t}$ will be the search points generated by a First Order minimization method as applied to $(D)$ at the steps $1,...,t$. We call $y^{t}$ the associated execution protocol, call a collection $\lambda^{t}$ of $t$ nonnegative weights summing up to 1 an accuracy certificate for this protocol, and refer to the quantity $\epsilon(y^{t},\lambda^{t})$ as to the resolution of the certificate $\lambda^{t}$ at the protocol $y^{t}$.

Our main observation (cf. ) is as follows:

Let $y^{t}$, $\lambda^{t}$ be as above. Then $\widehat{x}:=x(y^{t},\lambda^{t})$, $\widehat{y}:=y(y^{t},\lambda^{t})$ are feasible solutions to problems $(P)$, $(D)$, respectively, and

Proof. Let $F(x,y)=\langle x,Ay+a\rangle+\psi(y)$ and $x(y)=x_{X}(Ay+a)$, so that $f(y)=F(x(y),y)$. Observe that $f^{\prime}(y)=F^{\prime}_{y}(x(y),y)$, where $F^{\prime}_{y}(x,y)$ is a selection of the subdifferential of $F$ w.r.t. $y$, that is, $F^{\prime}_{y}(x,y)\in\partial_{y}F(x,y)$ for all $x\in X$, $y\in Y$. Setting $x_{\tau}=x(y_{\tau})$, we have for all $y\in Y$:

The inclusions $\widehat{x}\in X$, $\widehat{y}\in Y$ are evident. o

In the proof of Proposition 2.1, the linearity of $F$ w.r.t. $x$ was never used, so that in fact we have proved a more general statement: Given a concave in $x\in X$ and convex in $y\in Y$ Lipschitz continuous function $F(x,y)$, let us associate with it a convex function $f(y)=\max_{x\in X}F(x,y)$, a concave function $f_{*}(x)=\min_{y\in Y}F(x,y)$ and problems $(P)$ and $(D)$. Let $F^{\prime}_{y}(x,y)$ be a vector field with $F^{\prime}_{y}(x,y)\in\partial_{y}F(x,y)$, so that with $x(y)\in\mathop{\hbox{\rm Argmax}}_{x\in X}F(x,y)$, the vector $f^{\prime}(y)=F^{\prime}_{y}(x(y),y)$ is a subgradient of $f$ at $y$. Assume that problem $(D)$ associated with $F$ is solved by a FO method using $f^{\prime}(y)=F^{\prime}_{y}(x(y),y)$ which produced execution protocol $y^{t}$ and accuracy certificate $\lambda^{t}$. Then setting

Moreover, let $\delta\geq 0$, and let $x_{\delta}(y)$ be a $\delta$-maximizer of $F(x,y)$ in $x\in X$: for all $y\in Y$,

Proposition 2.1 says that whenever we can equip the subsequent execution protocols generated by a FO method, as applied to the dual problem $(D)$, with accuracy certificates, we can generate solutions to the primal problem $(P)$ of inaccuracy going to 0 at the same rate as the certificate resolution. In the sequel, we shall point out some “good” accuracy certificates for several most attractive FO algorithms for nonsmooth convex minimization.

Accuracy certificates in oracle-oriented methods for large-scale nonsmooth convex optimization

As it was mentioned in the introduction, the Mirror Descent (MD) algorithm solving $(D)$ is given by a norm $\|\cdot\|$ on $E_{y}$ and a distance-generating function (d.-g.f.) $\omega(y):Y\to{\mathbf{R}}$ which should be continuous and convex on $Y$, should admit a continuous in $y\in Y^{o}=\{y\in Y:\partial\omega(y)\neq\emptyset\}$ selection of subdifferentials $\omega^{\prime}(y)$, and should be strongly convex, modulus 1, w.r.t. $\|\cdot\|$, that is,

A proximal setup $(\|\cdot\|,\omega(\cdot))$ for $Y,E_{y}$ gives rise to several entities, namely,

Bregman distance $V_{y}(z)=\omega(z)-\omega(y)-\langle\omega^{\prime}(y),z-y\rangle$ ($y\in Y^{o},z\in Y$). Due to strong convexity of $\omega$, we have

$\omega$-center $y_{\omega}=\mathop{\hbox{\rm argmin}}_{y\in Y}\omega(y)$ of $Y$ and $\omega$-diameter

which combines with the inequality $V_{y}(z)\geq{1\over 2}\|z-y\|^{2}$ to yield the relation

where $\xi\in E_{y}$ and $y\in Y^{o}$. This mapping takes its values in $Y^{o}$ and satisfies the relation

1.2 Mirror Descent algorithm

which is oracle represented, meaning that we have access to an oracle which, given on input $y\in Y$, returns $g(y)$. From now on we assume that this field is bounded:

where $\|\cdot\|_{*}$ is the norm conjugate to $\|\cdot\|$. The algorithm is the recurrence

where $\gamma_{\tau}>0$ are stepsizes. Let us equip this recurrence with accuracy certificates, setting

of $\lambda^{t}$ on the execution protocol $y^{t}=\{y_{\tau},g(y_{\tau})\}_{\tau=1}^{t}$ satisfies the standard MD efficiency estimate

In particular, if $\gamma_{\tau}={\gamma(t)\over\|g(y_{\tau})\|_{*}}$, $\gamma(t):={\Omega\over\sqrt{t}}$ for $1\leq\tau\leq t$ We assume here that $g_{\tau}\neq 0$ for all $\tau\leq t$. In the opposite case, the situation is trivial: when $g(y_{\tau_{*}})=0$, for some $\tau*\leq t$, setting $\lambda^{t}_{\tau}=0$ for $\tau\neq\tau_{*}$ and $\lambda^{t}_{\tau_{*}}=1$, we ensure that $\epsilon(y^{t},\lambda^{t})=0$.,

The proof of the proposition follows the lines of the “classical” proof in the case when $g(\cdot)$ is the (sub)gradient field of the objective of a convex minimization problem (see, e.g., Proposition 5.1 of ), and is omitted.

In order to solve $(D)$, we apply MD to the vector field $g=f^{\prime}$. Assuming that

we can set $L[g]=L_{f}$. For this setup, Proposition 2.1 implies that the MD accuracy certificate $\lambda^{t}$, as defined in (15), taken together with the MD execution protocol $y^{t}=\{y_{\tau},g(y_{\tau})=f^{\prime}(y_{\tau}):=A^{*}\underbrace{x(y_{\tau})}_{x_{\tau}}+\psi^{\prime}(y_{\tau})\}_{\tau=1}^{t}$, yield the primal-dual feasible approximate solutions

Combining this conclusion with Proposition 3.1, we arrive at the following result:

In the case of (18), for every $t=1,2,...$ the $t$-step MD with the stepsize policyWe assume that $f^{\prime}(y_{\tau})\neq 0$, for $\tau\leq t$; otherwise, as we remember, the situation is trivial.

as applied to $(D)$ yields feasible approximate solutions $\widehat{x}^{t}$, $\widehat{y}^{t}$ to $(P)$, $(D)$ such that

In particular, given $\epsilon>0$, it takes at most

2 Convex minimization with certificates, II: Mirror Descent with full memory

Algorithm MDL – Mirror Descent Level method – is a non-Euclidean version of (a variant of) the Bundle-Level method ; to the best of our knowledge, this extension was not presented in the literature. Another novelty in what follows is equipping the method with accuracy certificates.

MDL is a version of MD with “full memory”, meaning that the first order information on the objective being minimized is preserved and utilized at subsequent steps, rather than being “summarized” in the current iterate, as it is the case for MD. While the guaranteed accuracy bounds for MDL are similar to those for MD, the typical practical behavior of the algorithm is better than that of the “memoryless” MD.

MDL with certificates which we are about to describe is aimed at processing an oracle-represented vector field (13) satisfying (14), with the same assumptions on $Y$ and the same proximal setup as in the case of MD.

We associate with $y\in Y$ the affine function

and with a finite set $S\subset Y$ the family ${\cal F}_{S}$ of affine functions on $E_{y}$ which are convex combinations of the functions $h_{y}(z)$, $y\in S$. In the sequel, the words “we have at our disposal a function $h(\cdot)\in{\cal F}_{S}$” mean that we know the functions $h_{y}(\cdot)$, $y\in S$, and nonnegative weights $\lambda_{y}$, $y\in S$, summing up to 1, such that $h(z)=\sum_{y\in S}\lambda_{y}h_{y}(z)$.

of the algorithm is, given a tolerance $\epsilon>0$, to find a finite set $S\subseteq Y$ and $h\in{\cal F}_{S}$ such that

Note that our target $S$ and $h$ are of the form $S=\{y_{1},...,y_{t}\}$, $h(y)=\sum_{\tau=1}^{t}\lambda_{\tau}\langle g(y_{\tau}),y_{\tau}-y\rangle$ with nonnegative $\lambda_{\tau}$ summing up to 1. In other words, our target is to build an execution protocol $y^{t}=\{y_{\tau},g(y_{\tau})\}_{\tau=1}^{t}$ and an associated accuracy certificate $\lambda^{t}$ such that $\epsilon(y^{t},\lambda^{t})\leq\epsilon$.

2.2 Construction

As applied to (13), MDL at a step $t=1,2,...$ generates search point $y_{t}\in Y$ where the value $g(y_{t})$ of $g$ is computed; it provides us with the affine function $h_{t}(z)=\langle g(y_{t}),y_{t}-z\rangle$. Besides this, the method generates finite sets $I_{t}\subset\{1,...,t\}$ and

Steps of the method are split into subsequent phases numbered $s=1,2,...$, and every phase is associated with optimality gap $\Delta_{s}\geq 0$.

To initialize the method, we set $y_{1}=y_{\omega}$, $I_{0}=\emptyset$ (whence $S_{0}=\emptyset$ as well), $\Delta_{0}=+\infty$.

given $y_{t}$, we compute $g(y_{t})$, thus getting $h_{t}(\cdot)$, and set $I_{t}^{+}=I_{t-1}\cup\{t\}$, $S_{t}^{+}=S_{t-1}\cup\{y_{t}\}$;

By the von Neumann lemma, an optimal solution to this (auxiliary) problem is associated with nonnegative and summing up to 1 weights $\lambda_{\tau}^{t}$, $\tau\in I_{t}^{+}$ such that

and we assume that as a result of solving (24), both $\epsilon_{t}$ and $\lambda^{t}_{\tau}$ become known. We set $\lambda^{t}_{\tau}=0$ for all $\tau\leq t$ which are not in $I_{t}^{+}$, thus getting an accuracy certificate $\lambda^{t}=[\lambda^{t}_{1};...;\lambda^{t}_{t}]$ for the execution protocol $y^{t}=\{y_{\tau},g(y_{\tau})\}_{\tau=1}^{t}$ along with $h^{t}(\cdot)=\sum_{\tau=1}^{t}\lambda^{t}_{\tau}h_{\tau}(\cdot)$. Note that by construction

If $\epsilon_{t}\leq\epsilon$, we terminate – $h(\cdot)=h^{t}(\cdot)$ satisfies (23). Otherwise we proceed as follows:

If (case A) $\epsilon_{t}\leq\gamma\Delta_{s-1}$, $\gamma\in(0,1)$ being method’s control parameter, we say that step $t$ starts phase $s$ (e.g., step $t=1$ starts phase 1), set

2.3 Efficiency estimate

Given on input a target tolerance $\epsilon>0$, the MDL algorithm terminates after finitely many steps, with the output $y^{t}=\{Y\ni y_{\tau},g(y_{\tau})\}_{\tau=1}^{t}$, $\lambda^{t}=\{\lambda^{t}_{\tau}\geq 0\}_{\tau=1}^{t}$, $\sum_{\tau}\lambda^{t}_{\tau}=1$ such that

The number of steps of the algorithm does not exceed

Assume that $\omega(\cdot)$ is continuously differentiable on the entire $Y$, so that the quantity

is finite. From the proof of Proposition 3.2 it follows immediately that one can substitute the rule “$\widehat{y}_{t}=y_{\omega}$ when $t$ starts a phase and $\widehat{y}_{t}=y_{t}$ otherwise” with a simpler one “$\widehat{y}_{t}=y_{t}$ for all $t$,” at the price of replacing $\Omega$ in (28) with $\Omega^{+}$.

is completely similar to the case of MD: given a desired tolerance $\epsilon>0$, one applies MDL to the vector field $g(y)=f^{\prime}(y)$ until the target (23) is satisfied. Assuming (18), we can set $L[g]=L_{f}$, so that by Proposition 3.2 our target will be achieved in

steps, with $L_{f}$ given by (18). Assuming that the target is attained at a step $t$, we have at our disposal the execution protocol $y^{t}=\{y_{\tau},f^{\prime}(y_{\tau})\}_{\tau=1}^{t}$ along with the accuracy certificate $\lambda^{t}=\{\lambda^{t}_{\tau}\}$ such that $\epsilon(y^{t},\lambda^{t})\leq\epsilon$ (by the same Proposition 3.2). Therefore, specifying $\widehat{x}^{t}$, $\widehat{y}^{t}$ according to (19) and invoking Proposition 2.1, we ensure (22). Note that the complexity $t=t(\epsilon)$ of finding these solutions, as given by (29), is completely similar to the complexity bound (21) of MD.

3 Convex minimization with certificates, III: restricted memory Mirror Descent

The fact that the number of linear functions $h_{\tau}(\cdot)$ involved into the auxiliary problems (24), (26) (and thus computational complexity of these problems) grows as the algorithm proceeds is a serious shortcoming of MDL from the computational viewpoint. NERML (Non-Euclidean Restricted Memory Level Method) algorithm, which originates from is a version of MD “with restricted memory”. In this algorithm the number of affine functions $h_{\tau}(\cdot)$ participating in (24), (26) never exceeds $m+1$, where $m$ is a control parameter which can be set to any desired positive integer value. The original NERML algorithm, however, was not equipped with accuracy certificates, and our goal here is to correct this omission.

Same as MDL, NERML processes an oracle-represented vector field (13) satisfying the boundedness condition (14), with the ultimate goal to ensure (23). The setup for the algorithm is identical to that for MDL.

The algorithm builds search sequence $y_{1}\in Y,y_{2}\in Y,...$ along with the sets $S_{\tau}=\{y_{1},...,y_{\tau}\}$, according to the following rules: A. Initialization. We set $y_{1}=y_{\omega}:=\mathop{\hbox{\rm argmin}}_{y\in Y}\omega(y)$, compute $g(y_{1})$ and set $f_{1}=\max\limits_{y\in Y}h_{y_{1}}(y)$. We clearly have $f_{1}\geq 0$.

In the case of $f_{1}=0$, we terminate and output $h(\cdot)=h_{y_{1}}(\cdot)\in{\cal F}_{S_{1}}$, thus ensuring (23) with $\epsilon=0$.

When $f_{1}>0$, we proceed. Our subsequent actions are split into phases indexed with $s=1,2,...$.

B. Phase $s=1,2,...$ At the beginning of phase $s$, we have at our disposal

the set $S^{s}=\{y_{1},...,y_{t_{s}}\}\subset Y$ of already built search points, and

an affine function $h^{s}(\cdot)\in{\cal F}_{S^{s}}$ along with the real $f_{s}:=\max\limits_{y\in Y}h^{s}(y)\in(0,f_{1}]$.

To save notation, we denote the search points generated at phase $s$ as $u_{1},u_{2},...$, so that $y_{t_{s}+\tau}=u_{\tau}$, $\tau=1,2,...$. B.1. Initializing phase $s$. We somehow choose collection of $m$ functions $h_{0,j}^{s}(\cdot)\in{\cal F}_{S^{s}}$, $1\leq j\leq m$, such that the set

B.2. Step $\tau=1,2,...$ of phase $s$: B.2.1. At the beginning of step $\tau$, we have at our disposal

the set $S^{s}_{\tau-1}$ of all previous search points;

a collection of functions $\{h^{s}_{\tau-1,j}(\cdot)\in{\cal F}_{S^{s}_{\tau-1}}\}_{j=1}^{m}$ such that the set

current search point $u_{\tau}\in Y^{s}_{\tau-1}$ such that

Note that this relation is trivially true when $\tau=1$.

B.2.2. Our actions at step $\tau$ are as follows. B.2.2.1. We compute $g(u_{\tau})$ and set

$\lambda_{j}^{\tau}\geq 0$ and $\sum_{j=1}^{m+1}\lambda^{\tau}_{j}=1$. We assume that when solving the auxiliary problem, we compute the above weights $\lambda^{\tau}_{j}$, and thus have at our disposal the function

due to $Y_{\tau}\subset Y^{s}_{\tau-1}$. B.2.2.5. By optimality conditions for (31) (see Lemma A.1), for certain nonnegative $\mu_{j}$, $1\leq j\leq m+1$, such that

$\bullet$ In the case of $\mu=\sum_{j}\mu_{j}>0$, we set

We then discard from the collection $\{h^{s}_{\tau-1,j}(\cdot)\}_{j=1}^{m+1}$ two (arbitrarily chosen) elements and add to $h^{s}_{\tau,1}$ the remaining $m-1$ elements of the collection, thus getting an $m$-element collection $\{h^{s}_{\tau,j}\}_{j=1}^{m}$ of elements of ${\cal F}_{S^{s}_{\tau}}$.

In both cases (those of $\mu>0$ and of $\mu=0$), we have built the data required to start step $\tau+1$ of phase $s$, and we proceed to this step. The description of the algorithm is completed.

Same as MDL, the outlined algorithm requires solving at every step two nontrivial auxiliary optimization problems – (30) and (31). It is explained in that these problems are relatively easy, provided that $m$ is moderate (note that this parameter is under our full control) and $Y$ and $\omega$ are “simple and fit each other,” meaning that we can easily solve problems of the form

(that is, our proximal setup for $Y$ results in easy-to-compute prox-mapping).

By construction, the presented algorithm produces upon termination (if any)

an execution protocol $y^{t}=\{y_{\tau},g(y_{\tau})\}_{\tau=1}^{t}$, where $t$ is the step where the algorithm terminates, and $y_{\tau}$, $1\leq\tau\leq t$, are the search points generated in course of the run; by construction, all these search points belong to $Y$;

an accuracy certificate $\lambda^{t}$ – a collection of nonnegative weights $\lambda_{1},...,\lambda_{t}$ summing up to 1 – such that the affine function $h(y)=\sum_{\tau=1}^{t}\lambda_{\tau}\langle g(y_{\tau}),y_{\tau}-y\rangle$ satisfies the relation $\epsilon(y^{t},\lambda^{t}):=\max\limits_{x\in Y}h(x)\leq\epsilon,$ where $\epsilon$ is the target tolerance, exactly as required in (23).

3.2 Efficiency estimate

Given on input a target tolerance $\epsilon>0$, the NERML algorithm terminates after finitely many steps, with execution protocol $y^{t}$ and accuracy certificate $\lambda^{t}$, described in Remark 3.4. The number of steps of the algorithm does not exceed

Inspecting the proof of Proposition 3.3, it is immediately seen that when $\omega(\cdot)$ is continuously differentiable on the entire $Y$, one can replace the rule (31) with

where $y^{s}$ is an arbitrary point of $Y^{o}=Y$. The cost of this modification is that of replacing $\Omega$ in the efficiency estimate with $\Omega^{+}$, see Remark 3.1. Computational experience shows that a good choice of $y^{s}$ is the best, in terms of the objective, search point generated before the beginning of phase $s$.

is completely similar to the case of MDL, with the bound

Observe that Propositions 3.1-3.3 do not impose restrictions of the vector field $g(\cdot)$ processed by the respective algorithms aside from the boundedness assumption (14). Invoking Remark 2.1, we arrive at the following conclusion:

in the situation of section 2.1 and given $\delta\geq 0$, let instead of exact maximizers $x(y)\in\mathop{\hbox{\rm Argmax}}_{x\in X}\langle x,Ay+a\rangle$, approximate maximizers $x_{\delta}(y)\in X$ such that $\langle x_{\delta}(y),Ay+a\rangle\geq\langle x(y),Ay+a\rangle-\delta$ for all $y\in Y$ be available. Let also

be the associated approximate subgradients of the objective $f$ of $(D)$. Assuming

let MD/MDL/NERML be applied to the vector field $g(\cdot)=f^{\prime}_{\delta}(\cdot)$. Then the number of steps of each method before termination remains bounded by the respective bound (17), (28) or (36), with $L_{f,\delta}$ in the role of $L_{f}$. Besides this, defining the approximate solutions to $(P)$, $(D)$ according to (19), with $x_{\tau}=x_{\delta}(y_{\tau})$, we ensure the validity of $\delta$-relaxed version of the accuracy guarantee (22), specifically, the relation

An alternative: Smoothing

An alternative to the approach we have presented so far is based on the use of the proximal setup for $Y$ to smoothen $f_{*}$ and then to maximize the resulting smooth approximation of $f_{*}$ by the Conditional Gradient (CG) algorithm. This approach is completely similar to the one used by Nesterov in his breakthrough paper , with the only difference that since in our situation domain $X$ admits LO oracle rather than a good proximal setup, we are bounded to replace the $O(1/t^{2})$-converging Nesterov’s method for smooth convex minimization with $O(1/t)$-converging CG.

Let us describe the CG implementation in our setting. Suppose that we are given a norm $\|\cdot\|_{x}$ on $E_{x}$, a representation (1) of $f_{*}$, a proximal point setup $(\|\cdot\|_{y},\;\omega(\cdot))$ for $Y$ and a desired tolerance $\epsilon>0$. We assume w.l.o.g. that $\min\limits_{y\in Y}\omega(y)=0$ and set, following Nesterov ,

From (1), the definition of $\Omega$ and the relation $\min_{Y}\omega=0$ it immediately follows that

and $f_{*}^{\beta}$ clearly is concave. It is well known (the proof goes back to J.-J. Moreau ) that strong convexity, modulus 1 w.r.t. $\|\cdot\|_{y}$, of $\omega(y)$ implies smoothness of $f_{*}^{\beta}$, specifically,

Observe also that under the assumption that an optimal solution $y(x)$ of the right hand side minimization problem in (38) is available at a moderate computational cost, In typical applications, $\psi$ is just linear, so that computing $y(x)$ is as easy as computing the value of the prox-mapping associated with $Y$, $\omega(\cdot)$. we have at our disposal a FO oracle for $f_{*}^{\beta}$:

We can now use this oracle, along with the LO oracle for $X$, to solve $(P)$ by CG If our only goal were to approximate $f_{*}$ by a smooth concave function, we could use other options, most notably the famous Moreau-Yosida regularization . The advantage of the smoothing based on Fenchel-type representation (1) of $f_{*}$ and proximal setup for $Y$ is that in many cases it indeed allows for computationally cheap approximations, which is usually not the case for Moreau-Yosida regularization.. In the sequel, we refer to the outlined algorithm as to SCG (Smoothed Conditional Gradient).

for SCG is readily given by Proposition 1.1. Indeed, assume from now on that $X$ is contained in $\|\cdot\|_{x}$-ball of radius $R$ of $E_{x}$. It is immediately seen that under this assumption, (40) implies the validity of the condition (cf. (2) with $q=2$)

In order to find an $\epsilon$-maximizer of $f_{*}$, it suffices, by (39), to find an $\epsilon/2$-maximizer of $f_{*}^{\beta}$; by (4) (where one should set $q=2$), what takes

Let us assume, as above, that $X$ is contained in the centered at the origin $\|\cdot\|_{x}$-ball of radius $R$, and let us compare the essentially identical to each otherprovided the parameters $\gamma\in(0,1)$, $\theta\in(0,1)$ in (29), (37) are treated as absolute constants. complexity bounds (21), (29), (37), with the bound (42). Under the natural assumption that the subgradients of $\psi$ we use satisfy the bounds $\|\psi^{\prime}(y)\|_{y,*}\leq L_{\psi}$, where $L_{\psi}$ is the Lipschitz constant of $\psi$ w.r.t. the norm $\|\cdot\|_{y}$, (18) implies that

Thus, the first three complexity bounds reduce to

while the conditional gradients based complexity bound is

We see that assuming $L_{\psi}\leq O(1)R\|A\|_{y;x,*}$ (which indeed is the case in many applications, in particular, in the examples we are about to consider), the complexity bounds in question are essentially identical. This being said, we believe that the two approaches in question seem to have their own advantages and disadvantages. Let us name just a few:

Formally, the SCG has a more restricted area of applications than MD/MDL/NERML, since relative simplicity of the optimization problem in (38) is a more restrictive requirement than relative simplicity of computing prox-mapping associated with $Y,\omega(\cdot)$. At the same time, in most important applications known to us $\psi$ is just linear, and in this case the just outlined phenomenon disappears.

An argument in favor of SCG is its insensitivity to the Lipschitz constant of $\psi$. Note, however, that in the case of linear $\psi$ (which, as we have mentioned, is the case of primary interest) the nonsmooth techniques admit simple modifications (not to be considered here) which make them equally insensitive to $L_{\psi}$.

Our experience shows that the convergence pattern of nonsmooth methods utilizing memory (MDL and NERML) is, at least at the beginning of the solution process, much better than is predicted by their worst-case efficiency estimates. It should be added that in theory there exist situations where the nonsmooth approach “most probably,” or even provably, significantly outperforms the smooth one. This is the case when $E_{y}$ is of moderate dimension. A well-established experimental fact is that when solving $(D)$ by MDL, every $\dim E_{y}$ iterations of the method reduce the inaccuracy by an absolute constant factor, something like 3. It follows that if $n$ is in the range of few hundreds, a couple of thousands of MDL steps can yield a solution of accuracy which is incomparably better than the one predicted by the theoretical worst-case oriented $O(1/\epsilon^{2})$ complexity bound of the algorithm. Moreover, in principle one can solve $(D)$ by the Ellipsoid method with certificates , building accuracy certificate of resolution $\epsilon$ in polynomial time $O(1)n^{2}\ln(L_{f}\Omega[Y,\omega(\cdot)]/\epsilon)$. It follows that when $\dim E_{y}$ is in the range of few tens, the nonsmooth approach allows to solve, in moderate time, problems $(P)$ and $(D)$ to high accuracy. Note that low dimensionality of $E_{y}$ by itself does not prevent $X$ to be high-dimensional and “difficult;” how frequent are these situations in actual applications, this is another story.

We believe that the choice of one, if any, of the outlined approaches to use, is the issue which should be resolved, on the case-by-case basis, by computational practice. We believe, however, that it makes sense to keep them both in mind.

Application examples

In this section we work out some application examples, with the goal to demonstrate that the approach we are proposing possesses certain application potential.

Our first example (for its statistical motivation, see ) is as follows: given a symmetric $p\times p$ matrix $b$ and a positive real $R$, we want to find the best entrywise approximation of $b$ by a positive semidefinite matrix $x$ of given trace $R$, that is, to solve the problem

where ${\mathbf{S}}^{p}$ is the space of $p\times p$ symmetric matrices. Note that with our $X$, computing prox-mappings associated with all known proximal setups needs eigenvalue decomposition of a $p\times p$ symmetric matrix and thus becomes computationally demanding in the large scale case. On the other hand, to maximize a linear form $\langle\xi,x\rangle=\hbox{\rm Tr}(\xi x)$ over $x\in X$ requires computing the maximal eigenvalue of $\xi$ along with corresponding eigenvector. In the large scale case this task is by orders of magnitude less demanding than computing full eigenvalue decomposition. Note that our $f_{*}$ admits a simple Fenchel-type representation:

Equipping $E_{y}={\mathbf{S}}^{p}$ with the norm $\|\cdot\|_{y}=\|\cdot\|_{1}$, and $Y$ with the d.-g.f.

where $\alpha$ is an appropriately chosen constant of order of 1 (induced by the necessity to make $\omega(\cdot)$ strongly convex, modulus 1, w.r.t. $\|\cdot\|_{1}$), we get a proximal setup for $Y$ such that

We see that our problem of interest fits well the setup of methods developed in this paper. Invoking the bounds (44), (45), we conclude that (46) can be solved within accuracy $\epsilon$ in at most

steps by any of methods MD, MDL or NERML, and in at most $O(1){R^{2}\ln(p)\over\epsilon^{2}}$ steps by SCG.

It is worth to mention that in the case in question, the algorithms yielded by the nonsmooth approach admit a “sparsification” as follows. We are in the case of $\langle x,Ay+a\rangle\equiv\hbox{\rm Tr}(xy)$, and $X=\{x:x\succeq 0,\hbox{\rm Tr}(x)=R\}$, so that $x(y)=Re_{y}e_{y}^{T}$, where $e_{y}$ is the leading eigenvector of a matrix $y$ normalized to have $\|e_{y}\|_{2}=1$. Given a desired accuracy $\epsilon>0$ and a unit vector $e_{y,\epsilon}$ such that $\hbox{\rm Tr}(y[e_{y,\epsilon}e_{y,\epsilon}^{T}])\geq\hbox{\rm Tr}(y[e_{y}e_{y}^{T}])-R^{-1}\epsilon$, and setting $x_{\epsilon}(y)=Re_{y,\epsilon}e_{y,\epsilon}^{T}$, we ensure that $x_{\epsilon}(y)\in X$ and that $x_{\epsilon}(y)$ is an $\epsilon$-maximizer of $\langle x,Ay+a\rangle$ over $x\in X$. Invoking Remark 3.6, we conclude that when utilizing $x_{\epsilon}(\cdot)$ in the role of $x(\cdot)$, we get $2\epsilon$-accurate solutions to $(P)$, $(D)$ in no more than $t(\epsilon)$ steps. Now, we can take as $e_{y,\epsilon}$ the normalized leading eigenvector of an arbitrary matrix $\widehat{y}(y)$ such that $\|\sigma(\widehat{y}-y)\|_{\infty}\leq\epsilon$. Assuming $R/\epsilon>1$ and given $y\in Y$, let us sort the magnitudes of entries in $y$ and build $y_{\epsilon}$ by “thresholding” – by zeroing out as many smallest in magnitude entries as possible under the restriction that the remaining part of the matrix $y$ is symmetric, and the sum of squares of the entries we have replaced with zeros does not exceed $R^{-2}\epsilon^{2}$. Since $\|y\|_{1}\leq 1$, the number $N_{\epsilon}$ of nonzero entries in $y_{\epsilon}$ is at most $O(1)R^{2}/\epsilon^{2}$. On the other hand, by construction, the Frobenius norm $\|\sigma(y-y_{\epsilon})\|_{2}$ of $y-y_{\epsilon}$ is $\leq R^{-1}\epsilon$, thus $\|\sigma(y-y_{\epsilon})\|_{\infty}\leq R^{-1}\epsilon$, and we can take as $e_{y,\epsilon}$ the normalized leading eigenvector of $y_{\epsilon}$. When the size $p$ of $y$ is $\gg R^{2}/\epsilon^{2}$ (otherwise the outlined sparsification does not make sense), this approach reduces the problem of computing the leading eigenvector to the case when the matrix is question is relatively sparse, thus reducing its computational cost.

2 Nuclear norm SVM

Our next example is as follows: we are given an $N$-element sample of $p\times q$ matrices $z_{j}$ (“images”) equipped with labels $\epsilon_{j}\in\{-1,1\}$. We assume the images to be normalized by the restriction

We want to find a linear classifier of the form

which predicts well labels of images. We assume that there exist such right and left orthogonal transformations of the image, that the label can be predicted using only a small number of diagonal elements of the transformed image. This implies that the classifier we are looking for is sparse in the corresponding basis, or that the matrix $x$ is of low rank. We arrive at the “low-rank-oriented” SVM-based reformulation of this problem:

where $\|\sigma(\cdot)\|_{1}$ is the nuclear norm, $[a]_{+}=\max[a,0]$, and $R\geq 1$ is a parameter.The restriction $R\geq 1$ is quite natural. Indeed, with the optimal choice of $x$, we want most of the terms $\left[1-\epsilon_{j}\left[\langle x,z_{j}\rangle+b\right]\right]_{+}$ to be $\ll 1$; assuming that the number of examples with $\epsilon_{j}=-1$ and $\epsilon_{j}=1$ are of order of $N$, this condition can be met only when $|\langle x,z_{j}\rangle|$ are at least of order of 1 for most of $j$’s. The latter, in view of (48), implies that $\|\sigma(x)\|_{1}$ should be at least $O(1)$.

In this case the domain $X$ of problem ($P$) is the ball of the nuclear norm in the space ${\mathbf{R}}^{p\times q}$ of $p\times q$ matrices and $p,q$ are large. As we have explained in the introduction, same as in the example of the previous section, in this case the computational complexity of LO oracle is typically much smaller than the complexity of computing prox-mapping. Thus, from practical viewpoint, in a meaningful range of values of $p,q$ the LO oracle is “affordable,” while the prox-mapping is not.

Observing that $[a]_{+}=\max_{0\leq y\leq 1}ya$, and denoting ${\mathbf{1}}=[1;...;1]\in{\mathbf{R}}^{q}$, we get

from now on we assume that $Y\neq\emptyset$. When setting

and passing from minimizing $h(x)$ to maximizing $f_{*}(x)\equiv-h(x)$, problem (49) becomes

Let us equip $E_{y}={\mathbf{R}}^{N}$ with the standard Euclidean norm $\|\cdot\|_{2}$, and $Y$ - with the Euclidean d.-g.f. $\omega(y)={1\over 2}y^{T}y$. Observe that

(we are in the case of $\|\cdot\|_{*}=\|\cdot\|_{2}$), and, besides,

We conclude that for every $\epsilon>0$, the number $t$ of MD steps needed to ensure (22) does not exceed

(see (21)), and similarly for MDL, NERML, and SCG.

3 Multi-class classification under ∞|2conditional2\infty|2 norm constraint

Our last example illustrates the potential of the proposed approach in the case when the domain $X$ of $(P)$ does not admit a proximal setup with “moderate” $\Omega_{X}$. Namely, let $(P)$ be the problem