Solving Variational Inequalities with Monotone Operators on Domains Given by Linear Minimization Oracles

Introduction

The majority of First Order methods (FOM’s) for large-scale convex minimization (and all known to us FOM’s for large-scale convex-concave saddle point problems and variational inequalities with monotone operators) are of proximal type: at a step of the algorithm, one needs to compute prox-mapping – to minimize over problem’s domain the sum of a linear function and a specific for the algorithm strongly convex distance generating function (d.-g.f.), in the simplest case, just squared Euclidean norm. As a result, the practical scope of proximal algorithms is restricted to proximal-friendly domains – those allowing for d.-g.f.’s with not too expensive computationally prox-mappings. What follows is motivated by the desire to develop FOM’s for solving convex-concave saddle point problems on bounded domains with “difficult geometry” – those for which no d.-g.f.’s resulting in nonexpensive prox-mappings (and thus no “implementable” proximal methods) are known. In what follows, we relax the assumption on problem’s domain to be proximal-friendly to the weaker assumption to admit computationally nonexpensive Linear Minimization Oracle (LMO) – a routine capable to minimize a linear function over the domain. This indeed is a relaxation: to minimize within a desired, whatever high, accuracy a linear form over a bounded proximal-friendly domain is the same as to minimize over the domain the sum of large multiple of the form and the d.-g.f. Thus, proximal friendliness implies existence of a nonexpensive LMO, but not vice versa. For example, when the domain is the ball $B_{n}$ of nuclear norm in ${\mathbf{R}}^{n\times n}$, computing prox-mapping, for all known proximal setups, requires full singular value decomposition of an $n\times n$ matrix, which can be prohibitively time consuming when $n$ is large. In contrast to this, minimizing a linear form over $B_{n}$ only requires finding the leading singular vectors of an $n\times n$ matrix, which is much easier than full-fledged singular value decomposition.

Recently, there was significant interest in solving convex minimization problems on domains given by LMO’s. The emphasis in this line of research is on smooth/smooth norm-regularized convex minimization , where the main “working horse” is the classical Conditional Gradient (a.k.a. Frank-Wolfe) algorithm originating from and intensively studied in 1970’s (see and references therein). Essentially, Conditional Gradient is the only traditional convex optimization technique capable to handle convex minimization problems on LMO-represented domains. In its standard form, Conditional Gradient algorithm, to the best of our knowledge, is not applicable beyond the smooth minimization setting; we are not aware of any attempt to apply this algorithm even to the simplest – bilinear – saddle point problems. The approach proposed in this paper is different and is inspired by our recent paper , where a method for nonsmooth convex minimization over an LMO-represented convex domain was developed. The latter method utilizes Fenchel-type representations of the objective in order to pass from the problem of interest to its special dual. In many important cases the domain of the dual problem is proximal-friendly, so that the dual problem can be solved by proximal FOM’s. We then use the machinery of accuracy certificates originating from allowing to recover a good solution to the problem of interest from the information accumulated when solving the dual problem. In this paper we follow the same strategy in the context of variational inequalities (v.i.’s) with monotone operators (this covers, in particular, convex-concave saddle point problems). Specifically, we introduce the notion of a Fenchel-type representation of a monotone operator, allowing to associate with the v.i. of interest its dual, which is again a v.i. with monotone operator with the values readily given by the representation and the LMO representing the domain of the original v.i. Then we solve the dual v.i. (e.g., by a proximal-type algorithm) and use the machinery of accuracy certificates to recover a good solution to the v.i. of interest from the information gathered when solving the dual v.i.

The main body of the paper is organized as follows. Section 2 outlines the background of convex-concave saddle point problems, variational inequalities with monotone operators and accuracy certificates. In section 3, we introduce the notion of a Fenchel-type representation of a monotone operator and the induced by this notion concept of v.i. dual to a given v.i. This section also contains a simple fully algorithmic “calculus” of Fenchel-type representations of monotone operators: it turns out that basic monotonicity-preserving operations with these operators (summation, affine substitution of argument, etc.) as applied to operands given by Fenchel-type representations yield similar representation for the result of the operation. As a consequence, our abilities to operate numerically with Fenchel-type representations of monotone operators are comparable with our abilities to evaluate the operators themselves. Section 4 contains our main result – Theorem 1. It shows how information collected when solving the dual v.i. to some accuracy, can be used to build an approximate solution of the same accuracy to the primal v.i. In Section 4 we present a self-contained description of two well known proximal type algorithms for v.i.’s with monotone operators – Mirror Descent (MD) and Mirror Prox (MP) – which indeed are capable to collect the required information. Section 5 is devoted to some modifications of our approach as applied to an affine monotone operator. In the concluding Section 6, we illustrate the proposed approach by applying it to the “matrix completion problem with spectral norm fit” – to the problem

where $\|x\|_{\hbox{\scriptsize\rm nuc}}=\sum_{i}\sigma_{i}(x)$ is the nuclear norm, $\sigma(x)$ being the singular spectrum of $x$, $\|x\|_{2,2}=\max_{i}\sigma_{i}(x)$ is the spectral norm, and $u\mapsto{\cal A}u$ is a linear mapping from ${\mathbf{R}}^{n\times n}$ to ${\mathbf{R}}^{m\times m}$.

Preliminaries

Let $Y$ be a nonempty closed convex set in Euclidean space $E_{y}$ and $H(y):Y\to E_{y}$ be a monotone operator:

The variational inequality (v.i.) associated with $(H,Y)$ is

(every) $y_{*}\in Y$ satisfying the target relation in ${\hbox{\rm VI}}(H,Y)$ is called a weak solution to the v.i.; when $Y$ is convex and compact, and $H$ is monotone on $Y$, weak solutions always exist. A strong solution to v.i. is a point $y_{*}\in Y$ such that $\langle H(y_{*}),y-y_{*}\rangle\geq 0$ for all $y\in Y$; from the monotonicity of $H$ is follows that a strong solution is a weak one as well. Note that when $H$ is monotone and continuous on $Y$ (this is the only case we will be interested in), weak solutions are exactly the strong solutions.

The accuracy measure naturally quantifying the inaccuracy of a candidate solution $y\in Y$ to ${\hbox{\rm VI}}(H,Y)$ is the dual gap function

this (clearly nonnegative for $y\in Y$) quantity is zero if and only if $y$ is a weak solution to the v.i.

We will be interested also in the Special case where $Y=V\times W$ is the direct product of nonempty convex compact subsets $V\subset E_{v}$ and $W\subset E_{w}$ of Euclidean spaces $E_{v}$, $E_{w}$, and $H$ is associated with Lipschitz continuous function $f(v,w):Y=V\times W\to{\mathbf{R}}$ convex in $v\in V$ and concave in $w\in W$:

We can associate with the Special case two optimization problems

under our assumptions ($V,W$ are convex and compact, $f$ is continuous convex-concave) these problems are solvable with equal optimal values. We associate with a pair $(v,w)\in V\times W$ the saddle point inaccuracy

Given $Y,\;H$, let us call a collection ${\cal C}^{N}=\{y_{t}\in Y,\,\lambda_{t}\geq 0,\,H(y_{t})\}_{t=1}^{N}$ with $\sum_{t}\lambda_{t}=1$, an $N$-step accuracy certificate. For $Z\subset Y$, we call the quantity

the resolution of the certificate ${\cal C}^{N}$ w.r.t. $Z$.

Let us make two observations coming back to :

Let $Y$ be a closed convex set in Euclidean space $E_{y}$, $H$ be a monotone operator on $Y$, and ${\cal C}^{N}=\{y_{t}\in Y,\lambda_{i}\geq 0,H(y_{t})\}_{t=1}^{N}$ be an accuracy certificate. Setting

we have $\widehat{y}\in Y$, and for every nonempty closed convex subset $Y^{\prime}$ of $Y$ it holds

Thus, $\langle H(z),\widehat{y}-z\rangle\leq{\hbox{\rm Res}}\left({\cal C}^{N}|Y\right)$ for all $z\in Y^{\prime}$, and (1) follows. In the Special case, setting $y_{t}=[v_{t};w_{t}]$, $\widehat{y}=[\widehat{v};\widehat{w}]$, for every $y=[v;w]\in Y=V\times W$ we have

Since the resulting inequality holds true for all $v\in V$, $w\in W$, we get $\overline{f}(\widehat{v})-\underline{f}(\widehat{w})\leq{\hbox{\rm Res}}({\cal C}^{N}|Y)$, and (2) follows. $\square$

Lemma 1 can be partially inverted in the case of skew-symmetric operator $H$, that is,

with skew-symmetric $(S=-S^{*})$ From now on, for a linear mapping $x\mapsto Bx:E\to F$, where $E,F$ are Euclidean spaces, $B^{*}$ denotes the conjugate of $B$, that is, a linear mapping $y\mapsto B^{*}y:F\to E$ uniquely defined by the identity $\langle Bx,y\rangle=\langle x,B^{*}y\rangle$ for all $x\in E,$ $y\in F$. linear operator $S$. A skew-symmetric $H$ clearly satisfies the identity

Let $Y$ be a convex compact set in Euclidean space $E_{y}$, $H(y)=a+Sy$ be skew-symmetric, and let ${\cal C}^{N}=\{y_{t}\in Y,\lambda_{t}\geq 0,H(y_{t})\}_{t=1}^{N}$ be an accuracy certificate. Then for $\widehat{y}=\sum_{t}\lambda_{t}y_{t}$ it holds

Proof. We already know that $\epsilon_{\hbox{\scriptsize\rm vi}}(\widehat{y}|H,Y)\leq{\hbox{\rm Res}}\left({\cal C}^{N}|Y\right)$. To prove the inverse inequality, note that for every $y\in Y$ we have

Thus, $\sum_{t}\lambda_{t}\langle H(y_{t}),y_{t}-y\rangle\leq\epsilon_{\hbox{\scriptsize\rm vi}}(\widehat{y}|H,Y)$ for all $y\in Y$, so that ${\hbox{\rm Res}}\left({\cal C}^{N}|Y\right)\leq\epsilon_{\hbox{\scriptsize\rm vi}}(\widehat{y}|H,Y)$. $\square$

Assume we are in the Special case, so that $Y=V\times W$ is a direct product of two convex compact sets, and the monotone operator $H$ is associated with a convex-concave function $f(v,w)$. Assume also that $f$ is bilinear: $f(v,w)=\langle a,v\rangle+\langle b,w\rangle+\langle w,{\cal A}v\rangle$, so that $H$ is affine and skew-symmetric. Then for every $y\in Y$ it holds

Proof. Consider accuracy certificate ${\cal C}^{1}=\{y_{1}=y,\lambda_{1}=1,H(y_{1})\}$; for this certificate, $\widehat{y}$ as defined in Lemma 2 is just $y$. Therefore, by Lemma 2, ${\hbox{\rm Res}}({\cal C}^{1}|Y)=\epsilon_{\hbox{\scriptsize\rm vi}}(y|H,Y)$. This equality combines with Lemma 1 to imply (8). $\square$

Representations of Monotone Operators

To explain the origin of the developments to follow, let us summarize the approach to solving convex minimization problems on domains given by Linear Minimization Oracles (LMOs), developed in . The principal ingredient of this approach is a Fenchel-type representation of a convex function $f:X\to{\mathbf{R}}$ defined on a convex subset $X$ of Euclidean space $E$; by definition, such a representation is

where $Y$ is a convex subset of Euclidean space $F$ and $\psi:Y\to{\mathbf{R}}$ is convex. Assuming for the sake of simplicity that $X$, $Y$ are compact and $\psi$ is continuously differentiable on $Y$, representation (9) allows to associate with the primal problem

with the same optimal value. Observe that the first order information on the (concave) objective of $(D)$ is readily given by the first order information on $\psi$ and the information provided by an LMO for $X$. As a result, we can solve $(D)$ by, say, a proximal type First Order Method, provided that $Y$ is proximal-friendly. The crucial in this approach question of how to recover a good approximate solution to the problem of interest $(P)$ from the information collected when solving $(D)$ is addressed via the machinery of accuracy certificates .

In the sequel, we intend to apply a similar scheme to the situation where the role of $(P)$ is played by a variational inequality with monotone operator on a convex compact domain $X$ given by an LMO. Our immediate task is to outline informally what a Fenchel-type representation of a monotone operator is and how we intend to use such a representation. To this end note that $(P)$ and $(D)$ can be reduced to variational inequalities with monotone operators, specifically

the “primal” v.i. stemming from $(P)$. The domain of this v.i. is $X$, and the operator is $f^{\prime}(x)=Ay(x)+a$, where $y(x)$ is a maximizer of the function $\langle x,Ay\rangle-\psi(y)$ over $y\in Y$, or, which is the same, a (strong) solution to the v.i. given by the domain $Y$ and the monotone operator $y\mapsto G(y)-A^{*}x$, where $G(y)=\psi^{\prime}(y)$;

the “dual” v.i. stemming from $(D)$. The domain of this v.i. is $Y$, and the operator is $y\mapsto G(y)-A^{*}x(y)$, where $x(y)$ is a minimizer of $\langle x,Ay+a\rangle$ over $x\in X$.

Observe that both operators in question are described in terms of a monotone operator $G$ on $Y$ and affine mapping $y\mapsto Ay+a:F\to E$; in the above construction $G$ was the gradient field of $\psi$, but the construction of the primal and the dual v.i.’s makes sense whenever $G$ is a monotone operator on $Y$ satisfying minimal regularity assumptions. The idea of the approach we are about to develop is as follows: in order to solve a v.i. with a monotone operator $\Phi$ and domain $X$ given by an LMO,

We represent $\Phi$ in the form of $\Phi(x)=Ay(x)+a$, where $y(x)$ is a strong solution to the v.i. on $Y$ given by the operator $G(y)-A^{*}x$, $G$ being an appropriate monotone operator on $Y$. It can be shown that a desired representation always exists, but by itself existence does not help much – we need the representation to be suitable for numerical treatment, to be available in a “closed computation-friendly form.” We show that “computation-friendly” representations of monotone operators admit a kind of fully algorithmic calculus which, for all basic monotonicity-preserving operations, allows to get straightforwardly a desired representation of the result of an operation from the representations of the operands. In view of this calculus, “closed analytic form” representations, allowing to compute efficiently the values of monotone operators, automatically lead to required computation-friendly representations.

We use the representation from A to build the “dual” v.i. with domain $Y$ and the operator $\Theta(y)=G(y)-A^{*}x(y)$, with exactly the same $x(y)$ as above, that is, x(y)\in\mathop{\hbox{\rm Argmin\,}}_{x\in X}\langle x,Ay+a\rangle. We shall see that $\Theta$ is monotone, and that usually there is a significant freedom in choosing $Y$; in particular, we typically can choose $Y$ to be proximal-friendly.

We solve the dual v.i. by an algorithm, like Mirror Descent or Mirror Prox, which produce necessary accuracy certificates. We will see – and this is our main result – that such a certificate ${\cal C}^{N}$ can be converted straightforwardly into a feasible solution $x^{N}$ to the v.i. of interest such that $\epsilon_{\hbox{\scriptsize\rm vi}}(x^{N}|\Phi,X)\leq{\hbox{\rm Res}}({\cal C}^{N}|Y)$. As a result, if the certificates in question are good, meaning that the resolution of ${\cal C}^{N}$ as a function of $N$ obeys the standard efficiency estimates of the algorithm used to solve the dual v.i., we solve the v.i. of interest with the same efficiency estimate as the one for the dual v.i. It remains to note that most of the existing first order algorithms for solving v.i.’s with monotone operators (various versions of polynomial time cutting plane algorithms, like the Ellipsoid method, Subgradient/Mirror Descent, and different bundle-level versions of Mirror Descent) indeed produce good accuracy certificates, see .

2 The construction

Consider the situation where we are given

a nonempty closed convex set $Y\subset F$;

which is good w.r.t. $A,Y$, goodness meaning that the variational inequality ${\hbox{\rm VI}}(G(\cdot)-A^{*}x,Y)$ has a strong solution for every $x\in E$. Note that when $Y$ is convex compact, every continuous monotone operator on $Y$ is good, whatever be $A$;

a nonempty convex compact set $X$ in $E$.

These data give rise to two operators: “primal” $\Phi:X\to E$ which is monotone, and “dual” $\Psi:Y\to F$ which is antimonotone (that is, $-\Psi$ is monotone).

2.2 Primal monotone operator

The primal operator $\Phi:E\to E$ is defined by

Observe that required $y(x)$ do exist: these are just strong solutions to the variational inequality given by the monotone operator $G(y)-A^{*}x$ and the domain $Y$.

Now, with $x^{\prime},x^{\prime\prime}\in E$, setting $y(x^{\prime})=y^{\prime}$, $y(x^{\prime\prime})=y^{\prime\prime}$, so that $y^{\prime},y^{\prime\prime}\in Y$, we have

Thus, $\Phi(x)$ is monotone. We call (10) a representation of the monotone operator $\Phi$, and the data $F,A,a,y(\cdot),G(\cdot),Y$ – the data of the representation. We also say that these data represent $\Phi$.

Given a convex domain $X\subset E$ and a monotone operator $\bar{\Phi}$ on this domain, we say that data $F$, $A$, $a$, $y(\cdot)$, $G(\cdot)$, $Y$ of the above type represent $\bar{\Phi}$ on $X$, if the monotone operator $\Phi$ represented by these data coincides with $\bar{\Phi}$ on $X$.

2.3 Dual operator

(in words: $\Psi(y)=A^{*}x(y)-G(y)$, where $x(y)$ minimizes $\langle Ay+a,x\rangle$ over $x\in X$). This operator clearly is antimonotone, as the sum of two antimonotone operators $-G(y)$ and $A^{*}x(y)$; antimonotonicity of the latter operator stems from the fact that it is obtained from the antimonotone operator $\psi(z)$ – a section of the superdifferential \mathop{\hbox{\rm Argmin\,}}_{x\in X}\langle z,x\rangle of a concave function – by affine substitution of variables: $A^{*}x(y)=A^{*}\psi(Ay+a)$, and this substitution preserves antimonotonicity.

Note that computing the value of $\Psi$ at a point $y$ reduces to computing $G(y)$, $Ay+a$, a single call to the Linear Minimization Oracle for $X$ to get $x(y)$, and computing $A^{*}x(y)$.

3 Calculus of representations

Let $F,A,a,y(\cdot),G(\cdot),Y$ represent a monotone operator $\Phi:E\to E$:

that is, a representation of $\lambda\Phi$ is given by $F,\lambda A,\lambda a,y(\cdot),\lambda G(\cdot),Y$; note that the operator $\lambda G$ clearly is good w.r.t. $\lambda A,Y$, since $G$ is good w.r.t. $A,Y$.

3.2 Summation

Let $F_{i},A_{i},a_{i},y_{i}(\cdot),G_{i}(\cdot),Y_{i}$, $1\leq i\leq m$, represent monotone operators $\Phi_{i}(x):E\to E$:

represent $\sum_{i}\Phi_{i}(x)$. Note that the operator $G(\cdot)$ clearly is good since $G_{1},...,G_{m}$ are so.

3.3 Affine substitution of argument

Let $F,A,a,y(\cdot),G(\cdot),Y$ represent $\Phi:E\to E$, let $H$ be a Euclidean space and $h\mapsto Qh+q$ be an affine mapping from $H$ to $E$. We have

that is, $F,\widehat{A},\widehat{a},\widehat{y}(\cdot),\widehat{G}(\cdot),Y$ represent $\widehat{\Phi}$. Note that $\widehat{G}$ clearly is good since $G$ is so.

3.4 Direct sum

Let $F_{i},A_{i},a_{i},y_{i}(\cdot),G_{i}(\cdot),Y_{i}$, $1\leq i\leq m$, represent monotone operators $\Phi_{i}(x_{i}):E_{i}\to E_{i}$. Then

represent $\Phi:\;E_{1}\times...\times E_{m}\to E_{1}\times...\times E_{m}$. Note that $G$ clearly is good since $G_{1},...,G_{m}$ are so.

3.5 Representing affine monotone operators

Consider an affine monotone operator on a Euclidean space $E$:

Its Fenchel-type representation on a convex compact set $X\subset E$ is readily given by the data $F=E$, $A=S$, $G(y)=S^{*}y:F\to F$ (this operator indeed is monotone), $y(x)=x$ and $Y$ being either the entire $F$, or (any) compact convex subset of $E=F$ which contains $X$; note that $G$ clearly is good w.r.t. $A,F$, same as is good w.r.t. $A,Y$ when $Y$ is compact. To check that the just defined $F,A,a,y(\cdot),G(\cdot),Y$ indeed represent $\Phi$ on $X$, observe that when $x\in X$, $y(x)=x$ belongs to $Y\supset X$ and clearly satisfies the relation $0\leq\langle A^{*}x-G(y(x)),y(x)-y\rangle\geq 0$ for all $y\in Y$ (see (10)), since $A^{*}x-G(y(x))=S^{*}x-S^{*}x=0$. Besides this, for $x\in X$ we have

as required for a representation. The dual antimonotone operator associated with this representation of $\Phi$ on $X$ is

3.6 Representing gradient fields

Let $f(x)$ be a convex function given by Fenchel-type representation

where $Y$ is a convex compact set in Euclidean space $F$, and $\psi(\cdot):F\to{\mathbf{R}}$ is a continuously differentiable convex function. Denoting by $y(x)$ a maximizer of $\langle x,Ay\rangle-\psi(y)$ over $y$, observe that

is a subgradient field of $f$, and that this monotone operator is given by a representation with the data $F,A,a,y(\cdot),G(\cdot):=\nabla\psi(\cdot),Y$; $G$ is good since $Y$ is compact.

Main result

Consider the situation described in section 3.2. Thus, we are given Euclidean space $E$, a convex compact set $X\subset E$ and a monotone operator $\Phi:E\to E$ represented according to (10), the data being $F,A,a,y(\cdot),G(\cdot),Y$. We denote by $\Psi:Y\to F$ the dual (antimonotone) operator induced by the data $X,A,a,y(\cdot),G(\cdot)$, see (11). Our goal is to solve variational inequality given by $\Phi$, $X$, and our main observation is that a good accuracy certificate for the variational inequality given by $(-\Psi,Y)$ induces an equally good solution to the variational inequality given by $(\Phi,X)$. The exact statement is as follows.

Let $X\subset E$ be a convex compact set and $\Phi:X\to E$ be a monotone operator represented on $X$, in the above sense, by data $F,A,a,y(\cdot),G(\cdot),Y$. Let also $\Psi:Y\to F$ be the antimonotone operator as defined above by the data $X,A,a,y(\cdot),G(\cdot)$. Let, finally,

be an accuracy certificate associated with the monotone operator $[-\Psi]$ and $Y$. Setting

(these points are byproducts of computing $\Psi(y_{t})$, $1\leq t\leq N$) and

When $\Phi(x)=a+Sx$, $x\in X$, with skew-symmetric $S$, we have also

In view of Theorem 1, given a representation of the monotone operator $\Phi$ participating in the v.i. of interest ${\hbox{\rm VI}}(\Phi,X)$, we can reduce solving the v.i. to solving the dual v.i. ${\hbox{\rm VI}}(-\Psi,Y)$ by an algorithm producing good accuracy certificates. Below we discuss in details the situation when the latter algorithm is either Mirror Descent (MD) [15, Chapter 5], or Mirror Prox (MP) (, see also [15, Chapter 6] and ).

Theorem 1 may be extended to the situation where the relationships (11), defining the dual operator $\Psi$ hold only approximately. We present here the following slight extension of the main result:

Let $X\subset E$ be a convex compact set and $\Phi:X\to E$ be a monotone operator represented on $X$, in the sense of section 3, by data $F,A,a,y(\cdot),G(\cdot),Y$. Given a positive integer $N$, sequences $y_{t}\in Y$, $x_{t}\in X$, $1\leq t\leq N$, and nonnegative reals $\lambda_{t}$, $1\leq t\leq N$, summing up to 1, let us set

Proofs of Theorems 1 and 2 are given in section A.1.

Saddle Point MD and MP are algorithms for solving convex-concave saddle point problems and variational inequalities with monotone operatorsMD algorithm originates from ; its modern proximal form was developed in . MP was proposed in . For the most present exposition of the algorithms, see [15, Chapters 5,6] and .. The algorithms are of proximal type, meaning that in order to apply the algorithm to a v.i. ${\hbox{\rm VI}}(H,Y)$, where $Y$ is a nonempty closed convex set in Euclidean space $E_{y}$ and $H$ is a monotone operator on $Y$, one needs to equip $E_{y}$ with a norm $\|\cdot\|$, and $Y$ - with a continuously differentiable distance generating function (d.-g.f.) $\omega(\cdot):Y\to{\mathbf{R}}$ compatible with $\|\cdot\|$, meaning that $\omega$ is strongly convex, modulus 1, w.r.t. $\|\cdot\|$. We call $\|\cdot\|,\omega(\cdot)$ proximal setup for $Y$. This setup gives rise to

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

where the concluding inequality is due to strong convexity of $\omega$,

$\omega$-size of a nonempty subset $Y^{\prime}\subset Y$

Due to the origin of $y_{\omega}$, we have V_{y_{\omega}}(y)\leq\mbox{\small\frac{1}{2}}\Omega^{2}[Y^{\prime}] for all $y\in Y^{\prime}$, implying that $\|y-y_{\omega}\|\leq\Omega$ for all $y\in Y^{\prime}$.

Given $y\in Y$, the prox-mapping with center $u$ is defined as

Let $Y$ be a nonempty closed convex set in Euclidean space $E_{y}$, $H=\{H_{t}:Y\to E_{y}\}_{t=1}^{\infty}$ be a sequence of vector fields, and $\|\cdot\|,\omega(\cdot)$ be a proximal setup for $Y$. As applied to $(H,Y)$, MD is the recurrence

In both MD and MP, $\gamma_{t}>0$ are stepsizes. The most important to us properties of these recurrences are as follows.

For $N=1,2,...$, consider the accuracy certificate

associated with (19). Then for every $Y^{\prime}\subset Y$ one has

where $\|\cdot\|_{*}$ is the norm conjugate to $\|\cdot\|$: $\|\xi\|_{*}=\max_{\|x\|\leq 1}\langle\xi,x\rangle.$

with some finite $M\geq 0$, then, given $Y^{\prime}\subset Y$, $N$ and setting

For $N=1,2,...$, consider the accuracy certificate

where $\|\cdot\|_{*}$ is the norm conjugate to $\|\cdot\|$. In particular, if

with some finite $L\geq 0$, $M\geq 0$, then given $Y^{\prime}\subset Y$, $N$ and setting

To make the text self-contained, we provide the proofs of these known results in the appendix.

2 Intermediate summary

Theorem 1 combines with Proposition 1 to imply the following claim:

In the situation of Theorem 1, let $y_{1},...,y_{N}$ be the trajectory of $N$-step MD as applied to the stationary sequence $H=\{H_{t}(\cdot)=-\Psi(\cdot)\}_{t=1}^{\infty}$ of vector fields, and let $x_{t}=x(y_{t})$, $t=1,...,N$. Then, setting $\lambda_{t}={\gamma_{t}\over\sum_{\tau=1}^{N}\gamma_{\tau}}$, $1\leq t\leq N$, we ensure that

finite and specifying $\gamma_{t}$, $1\leq t\leq N$, according to (23) with $Y^{\prime}=Y(X)$, we ensure that

When $\Phi(x)=Sx+a$, $x\in X$, with a skew-symmetric $S$, $\epsilon_{\hbox{\scriptsize\rm vi}}(\widehat{x}|\Phi,X)$ in the latter relation can be replaced with ${\hbox{\rm Res}}(\{x_{t},\lambda_{t},\Phi(x_{t})\}_{t=1}^{N}|X)$.

In the sequel, we shall refer to the implementation of our approach presented in Corollary 2 as to our basic scheme.

Modifications in Affine case

In this section, we present some modifications of the proposed approach as applied to the case of v.i. ${\hbox{\rm VI}}(\Phi,X)$ with affine monotone operator $\Phi$ and LMO-represented convex compact domain $X$. While the worst-case complexity bounds for the modified scheme are similar to the ones stated in Corollary 2, there are reasons to believe that in practice the modified scheme could outperform the basic one.

In the rest of this section, we consider the case when the monotone operator $\Phi$ is affine:

and our goal is to solve ${\hbox{\rm VI}}(\Phi,X)$, where $X$ is a convex compact subset of $E$; w.l.o.g. we assume that $0\in X$. We suppose that $\Phi$ is given by an affine Fenchel-type representation, that is, a representation with data

$F$ is a Euclidean space, $y\mapsto Ay+a$ is an affine mapping from $F$ to $E$;

$y\mapsto Gy:F\to F$ is a linear monotone mapping, and $x\to Bx:E\to F$ is a linear mapping such that

Note that (34.$a$) implies that setting $y(x)=Bx$, $y(x)$ is a strong solution to the v.i. associated with the operator $G(y)-A^{*}x$ and $Y=F$, while (34.$b$) says that $Ay(x)+a=\Phi(x)$, $x\in E$, as required in the definition (10) of a Fenchel-type representation of a monotone operator.

2 Strategy

We intend to get an approximate solution to ${\hbox{\rm VI}}(\Phi,X)$ by applying MP to a properly built sequence $H=\{H_{t}(\cdot)\}$ of vector fields on $F$. Let us fix a proximal setup $\|\cdot\|,\omega(\cdot)$ for $Y=F$; w.l.o.g., we assume that the $\omega$-center \mathop{\hbox{\rm argmin\,}}_{F}\omega(\cdot) of $F$ is the origin, that is, $\omega^{\prime}(0)=0$. Let $L$ be the operator norm of the mapping $y\mapsto G(y):=Gy:F\to F$ from $\|\cdot\|$ to $\|\cdot\|_{*}$, so that

or, equivalently, $\langle z,Gy\rangle\leq\|z\|\|y\|$ for all $z,y\in F$. In the sequel, we set

2.2 The construction

We intend to build $H_{t}(\cdot)$ recursively, according to the recurrence

Note that independently of the choice of $x_{t}\in E$, we have

Now the relationships of the MP recurrence imply that (see (58) and (59))

The essence of the matter is how we update the vectors $x_{t}$; this is the issue we consider next.

Since $\omega(\cdot)$ is strongly convex on $F$, the function $f_{y}(\cdot)$ is well defined on $E$; $f_{y}$ is convex as the supremum of a family of affine functions of $x$. Moreover, it is well known that in fact $f_{y}(\cdot)$ possesses Lipschitz continuous gradient. Specifically, let $\|\cdot\|_{E}$ be a norm on $E$, $\|\cdot\|_{E,*}$ be the norm conjugate to $\|\cdot\|_{E}$, and let ${\cal L}$ be the norm of the linear mapping $y\to Ay:F\to E$ from the norm $\|\cdot\|$ on $F$ to the norm $\|\cdot\|_{E,*}$ on $E$, so that

Let $z_{y}(\zeta)=\operatorname{Prox}_{y}(\zeta):F\to Y.$ Function $f_{y}(\cdot)$ is continuously differentiable with the gradient

and this gradient is Lipschitz continuous:

Observe, first, that when summing up inequalities (36), we get

Second, for any $x_{t}\in X$, $1\leq t\leq N$, we have $\widehat{x}={1\over N}\sum_{t=1}^{N}x_{t}\in X$. Further, invoking (18) with $\lambda_{t}=N^{-1}$, $1\leq t\leq N$, and $z_{t}$ in the role of $y_{t}$ (which by (41) allows to set $\epsilon={\Omega^{2}[Y(X)]L\over 2N}$), we get

(we have used (39) and have taken into account that $z_{t}=z_{y_{t}}(\gamma[Gy_{t}-A^{*}x_{t}])$, see (39) and (35); recall that $\gamma=1/L$). Note that so far our conclusions were independent on how $x_{t}\in X$ are selected.

Relation (42) implies that when $x_{t}$ is a minimizer of $f_{y_{t}}(\cdot)$ on $X$, we have $\langle\nabla f_{y_{t}}(x_{t}),x_{t}-{x}\rangle\leq 0$ for all ${x}\in X$, and with this “ideal” for our purposes choice of $x_{t}$, (42) would imply

which is an $O(1/N)$ efficiency estimate, much better that the $O(1/\sqrt{N})$-efficiency estimate (32).

Of course, we cannot simply specify $x_{t}$ as a point from \mathop{\hbox{\rm Argmin\,}}_{X}f_{y_{t}}(x), since this would require solving precisely at every step of the MP recurrence (35) a large-scale convex optimization problem. What we indeed intend to do, is to solve this problem approximately. Specifically, given $y_{t}$ (so that $f_{y_{t}}(\cdot)$ is identified), we can apply the classical Conditional Gradient Algorithm (CGA) (which, as was explained in the introduction, is, basically, the only traditional algorithm capable to minimize a smooth convex function over an LMO-represented convex compact set) in order to generate an approximate solution $x_{t}$ to the problem $\min_{X}f_{y_{t}}(x)$ satisfying, for some prescribed $\epsilon>0$, the relation

By (42), this course of actions implies the efficiency estimate

2.3 Complexity analysis

Let us equip $E$ with a norm $\|\cdot\|_{E}$, the conjugate norm being $\|\cdot\|_{E,*}$, and let ${\cal L}$ be the operator norm of the mapping $y\mapsto Ay$ as defined in Lemma 3. Let, further, $R=R_{E}(X)$ be the radius of the smallest $\|\cdot\|_{E}$-ball, centered at the origin, which contains $X$. Taking into account (40) and applying the standard results on CGA (see section A.5), for every $\epsilon\in(0,{\cal L}R^{2})$ it takes at most $O(1){\cal L}R_{E}^{2}(X)/\epsilon$ CGA steps to generate a point $x_{t}$ with $\delta_{t}\leq\epsilon$; here and below $O(1)$’s are absolute constants. Specifying $\epsilon$ as ${\Omega^{2}[Y(X)]\over 2N}$, (44) becomes

while the computational effort to generate $\widehat{x}$ is dominated by the necessity to generate $x_{1},...,x_{N}$, which amounts to the total of

CGA steps. The effort per step is dominated by the necessity to compute the vector $g=\nabla f_{y}(x)$, given $y\in F$, $x\in E$, and to minimize the linear form $\langle g,u\rangle$ over $u\in X$. In particular, to ensure $\epsilon_{\hbox{\scriptsize\rm vi}}(\widehat{x}|\Phi,X)\leq\epsilon$, the total number of CGA steps should be proportional to $1/\epsilon^{2}$. We see that in terms of the theoretical upper bound on the number of calls to the LMO for $X$ needed to get an $\epsilon$-solution, our current scheme has no advantages as compared to the MD-based approach analyzed in Corollary 2. We, however, may hope that in practice the outlined MP-based scheme can be better than our basic MD-based one, provided that we apply CGA in a “smart” way, e.g., use CGA with memory, see .

Illustration

We apply our construction to the following problem (“matrix completion with spectral norm fit”):A more interesting for applications problem (cf. ) would be ${\mathop{\hbox{Opt}}}=\min_{v\in{\mathbf{R}}^{p_{v}\times q_{v}}}\left\{\|v\|_{\hbox{\scriptsize\rm nuc}}:\|{\cal A}v-b\|_{2,2}\leq\delta\right\};$ applying the approach from , this problem can be reduced to a “small series” of problems (45).

where ${\mathbf{R}}^{p\times q}$ is the space of $p\times q$ real matrices, $\|x\|_{\hbox{\scriptsize\rm nuc}}=\sum_{i}\sigma_{i}(x)$ is the nuclear norm on this space (sum of the singular values $\sigma_{i}(x)$ of $x$), $\|x\|_{2,2}=\max_{i}\sigma_{i}(x)$ is the spectral norm of $x$ (which is exactly the conjugate of the nuclear norm), and ${\cal A}$ is a linear mapping from ${\mathbf{R}}^{p_{v}\times q_{v}}$ to ${\mathbf{R}}^{p_{b}\times q_{b}}$. We are interested in the “large-scale” case, where the sizes of $p_{v},q_{v}$ of $v$ are large enough to make the full singular value decomposition of a $p_{v}\times q_{v}$ matrix prohibitively time consuming, what seemingly rules out the possibility to solve (45) by proximal type First Order algorithms. We assume, at the same time, that computing the leading singular vectors and the leading singular value of a $p_{v}\times q_{v}$ or a $p_{b}\times q_{b}$ matrix (which, computationally, is by far easier task than finding full singular value decomposition) still can be carried out in reasonable time.

We rewrite (45) as a bilinear saddle point problem

(from now on $\langle\cdot,\cdot\rangle_{\hbox{\scriptsize\rm Fro}}$ stands for Frobenius inner product, and $\|\cdot\|_{\hbox{\scriptsize\rm Fro}}$ – for the Frobenius norm on the space(s) of matrices). The domain $X$ of the problem is the direct product of two unit nuclear norm balls; minimizing a linear form over this domain reduces to minimizing, given $\xi$ and $\eta$, the linear forms ${\mathop{\hbox{\rm Tr}}}(v\xi^{T})$, ${\mathop{\hbox{\rm Tr}}}(w\eta^{T})$ over $\{v\in{\mathbf{R}}^{p_{v}\times q_{v}}:\|v\|\leq 1\}$, resp., $\{w\in{\mathbf{R}}^{p_{b}\times q_{b}}:\|w\|\leq 1\}$, which, in turn, reduces to computing the leading singular vectors and singular values of $\xi$ and $\eta$.

The monotone operator associated with (46) is affine and skew-symmetric:

From now on we assume that ${\cal A}$ is of spectral norm at most 1, i.e.,

(this always can be achieved by scaling).

We can represent the restriction of $\Phi$ on $X$ by the data

Indeed, in the notation from section 3.2, for $x=[v;w]\in X=\{[v;w]\in{\mathbf{R}}^{p_{v}\times q_{v}}\times{\mathbf{R}}^{p_{b}\times q_{b}}:\|v\|_{\hbox{\scriptsize\rm nuc}}\leq 1,\|w\|_{\hbox{\scriptsize\rm nuc}}\leq 1\}$, the solution $y(x)=[\xi(x);\eta(x)]$ to the linear system $A^{*}x=G(y)$ is given by $\eta(x)=-v$, $\xi(x)={\cal A}^{*}w$, so that both components of $y(x)$ are of Frobenius norm $\leq 1$ (recall that spectral norm of ${\cal A}$ is $\leq 1$), and therefore $y(x)\in Y$. Besides this,

We conclude that when $x=[v;w]\in X$, the just defined $y(x)$ meets all requirements from (10), and thus the data $F,A,a,y(\cdot),G(\cdot),Y$ given by (47) indeed represent the monotone operator $\Phi$ on $X$.

$\Psi$ given by the data $F,A,a,y(\cdot),G(\cdot),Y$ is

We use the Euclidean proximal setup for $Y$, i.e., we equip the space $F$ embedding $Y$ with the Frobenius norm and take, as the d.-g.f. for $Y$, the function

resulting in $\Omega[Y]=\sqrt{2}.$ Furthermore, from (48) and the fact that the spectral norm of ${\cal A}$ is bounded by 1 it follows that the monotone operator $\Theta(y)=-\Psi(y):Y\to F$ satisfies (22) with $M=2\sqrt{2}$ and (28) with $L=0$ and $M=4\sqrt{2}$.

Theorem 1 combines with Corollary 1 to imply that when converting an accuracy certificate ${\cal C}^{N}$ for the dual v.i. ${\hbox{\rm VI}}(-\Psi,Y)$ into a feasible solution $\widehat{x}^{N}$ to the primal v.i. ${\hbox{\rm VI}}(\Phi,X)$, we ensure that

with $f,V,W$ given by (46). In other words, in the representation $\widehat{x}^{N}=[\widehat{v}^{N};\widehat{w}^{N}]$, $\widehat{v}^{N}$ is a feasible solution to problem (45) (which is the primal problem associated with (46)), and $\widehat{w}^{N}$ is a feasible solution to the problem

(which is the dual problem associated with (46)) with the sum of non-optimalities, in terms of respective objectives, $\leq{\hbox{\rm Res}}({\cal C}^{N}|Y)$. Computing $\underline{f}(\widehat{w})$ (which, together with computing $\overline{f}(\widehat{v})$, takes a single call to LMO for $X$), we get a lower bound on ${\mathop{\hbox{Opt}}}(P)={\mathop{\hbox{Opt}}}(D)$ which certifies that $\overline{f}(\widehat{v})-{\mathop{\hbox{Opt}}}(P)\leq{\hbox{\rm Res}}({\cal C}^{N}|Y)$.

2 Numerical illustration

Here we report on some numerical experiments with problem (45). In these experiments, we used $p_{b}=q_{b}=:m$, $p_{v}=q_{v}=:n$, with $n=2m$, and the mapping ${\cal A}$ given by

In the first series of experiments, the dual v.i. ${\hbox{\rm VI}}(-\Psi,Y)$ is solved by the MD algorithm with $N=512$ steps for all but the largest instance, where $N=257$ is used. The MD is applied to the stationary sequence $H_{t}\equiv-\Psi$, $t=1,2,...$, of vector fields. The stepsizes $\gamma_{t}$ are proportional, with coefficient of order of 1, to those given by (23.$b$) with $\|\cdot\|\equiv\|\cdot\|_{*}=\|\cdot\|_{\hbox{\scriptsize\rm Fro}}$ and $\Omega[Y]=\sqrt{2}$ As we have already mentioned, with our proximal setup, the $\omega$-size of $Y$ is $\leq\sqrt{2}$, and (22) is satisfied with $M=2\sqrt{2}$.; the coefficient was tuned empirically in pilot runs on small instances and is never changed afterwards. We also use two straightforward “tricks”:

Instead of considering one accuracy certificate, ${\cal C}^{N}=\{y_{t},\;\lambda^{N}_{t}=1/N,\;-\Psi(y_{t})\}_{t=1}^{N}$, we build a “bunch” of certificates

where $\mu$ runs through a grid in $\{1,...,N\}$ (in this implementation, a 16-element equidistant grid), and $\nu\in\{\mu,\mu+1,...,N\}$ runs through another equidistant grid (e.g., for the largest problem instance, the grid $\{1,9,17,...,257\}$). We compute the resolutions of these certificates and identify the best (with the smallest resolution) certificate obtained so far. Every 8 steps, the best certificate is used to compute the current approximate solution to (46) along with the saddle point inaccuracy of this solution.

When applying MD to problem (46), the “dual iterates” $y_{t}=[\xi_{t};\eta_{t}]$ and the “primal iterates” $x_{t}:=x(y_{t})=[v_{t};w_{t}]$ are pairs of matrices, with $n\times n$ matrices $\xi_{t},\eta_{t},v_{t}$ and $m\times m$ matrices $w_{t}$ (recall that we are in the case of $p_{v}=q_{v}=n$, $p_{b}=q_{b}=m$). It is easily seen that with ${\cal A}$ given by (50), the matrices $\xi_{t},\eta_{t},v_{t}$ are linear combinations of rank 1 matrices $\alpha_{i}\beta_{i}^{T}$, $1\leq i\leq(k+1)t$, and $w_{t}$ are linear combinations of rank 1 matrices $\delta_{i}\epsilon_{i}^{T}$, $1\leq i\leq t$, with on-line computable vectors $\alpha_{i},\beta_{i},\delta_{i},\epsilon_{i}$. Every step of MD adds $k+1$ new $\alpha$- and $k+1$ new $\beta$-vectors, and a pair of new $\delta$- and $\epsilon$-vectors. Our matrix iterates were represented by the vectors of coefficients in the above rank 1 decompositions (let us call this representation incremental), so that the computations performed at a step of MD, including computing the leading singular vectors by straightforward power iterations, are as if the standard representations of matrices were used, but all these matrices were of the size (at most) $n\times[(k+1)N]$, and not $n\times n$ and $m\times m$, as they actually are. In our experiments, for $k=2$ and $N\leq 512$, this incremental representation of iterates yields meaningful computational savings (e.g., by factor of $6$ for $n=8192$) as compared to the plain representation of iterates by 2D arrays.

of our preliminary experiments are presented in Table 1. There ${\cal C}^{t}$ stands for the best certificate found in course of $t$ steps, and ${\hbox{\rm Gap}}({\cal C}^{t})$ denotes the saddle point inaccuracy of the solution to (46) induced by this certificate (so that ${\hbox{\rm Gap}}({\cal C}^{t})$ is a valid upper bound on the inaccuracy, in terms of the objective, to which the problem of interest (45) was solved in course of $t$ steps). The comments are as follows:

The results clearly demonstrate “nearly linear”, and not quadratic, growth of running time with $m,n$; this is due to the incremental representation of iterates.

When evaluating the “convergence patterns” presented in the table, one should keep in mind that we are dealing with a method with slow $O(1/\sqrt{N})$ convergence rate, and from this perspective, 50-fold reduction in resolution in 512 steps is not that bad.

A natural alternative to the proposed approach would be to solve the saddle point problem (46) “as it is,” by applying to the associated primal v.i. (where the domain is the product of two nuclear norm balls and the operator is Lipschitz continuous and even skew symmetric) a proximal type saddle point algorithm and computing the required prox-mappings via full singular value decompositions. The state-of-the-art MP algorithm when applied to this problem exhibits $O(1/N)$ convergence rate;For the primal v.i., (28) holds true for some $L>0$ and $M=0$. Moreover, with properly selected proximal setup for (45) the complexity bound (30) becomes ${\hbox{\rm Res}}({\cal C}^{N}|Y)\leq O(1)\sqrt{\ln(n)\ln(m)}/N$. yet, every step of this method would require 2 SVD’s of $n\times n$, and 2 SVD’s of $m\times m$ matrices. As applied to the primal v.i., MD exhibits $O(1/\sqrt{N})$ convergence rate, but the steps are cheaper – we need one SVD of $n\times n$, and one SVD of an $m\times m$ matrix, and we are unaware of a proximal type algorithm for the primal v.i. with cheaper iterations. For the sizes $m,n,k$ we are interested in, the computational effort required by the outlined SVD’s is, for all practical purposes, the same as the overall effort per step. Taking into account the actual SVD cpu times on the platform used in our experiments, the overall running times presented in Table 1, i.e., times required by 512 steps of MD as applied to the dual v.i., allow for the following iteration counts $N$ for MP as applied to the primal v.i.:

and for twice larger iteration counts for MD. From our experience, for $n=1024$ (and perhaps for $n=2048$ as well), MP algorithm as applied to the primal v.i. would yield solutions of better quality than those obtained with our approach. It, however, would hardly be the case, for both MP and MD, when $n=4096$, and definitely would not be the case for $n=8192$. Finally, with $n=16384$, CPU time used by the 257-step MD as applied to the dual v.i. is hardly enough to complete just one iteration of MD as applied to the primal v.i. We believe these data demonstrate that the approach developed in this paper has certain practical potential.

2.2 Experiments with the MP-based scheme

In this section we briefly report on the results obtained with the modified MP-based scheme presented in section 5. Same as above, we use the test problems and representation (47) of the monotone operator of interest (with the only difference that now $Y=F$), and the Euclidean proximal setup. Using the Euclidean setup on $Y=F$ makes prox-mappings and functions $f_{y}(\cdot)$, defined in (37), extremely simple:

When choosing $\|\cdot\|_{E}$ to be the Frobenius norm,

and taking into account that the spectral norm of ${\cal A}$ is $\leq 1$, it is immediately seen that the quantities $L$, $\gamma$, ${\cal L}$ introduced in section 5.2, can be set to 1, and what was called $R_{E}(X)$ in section 5.2.3, can be set to $\sqrt{2}$. As a result, by the complexity analysis of section 5.2.3, in order to find an $\epsilon$-solution to the problem of interest, we need $O(1)\epsilon^{-1}$ iterations of the recurrence (35), with $O(1)\epsilon^{-1}$ CGA steps of minimizing $f_{y_{t}}(\cdot)$ over $X$ per iteration, that is, the total of at most $O(1)\epsilon^{-2}$ calls to the LMO for $X=V\times W$. In fact, in our implementation $\epsilon$ is not fixed in advance; instead, we fix the total number $N=256$ of calls to LMO, and terminate CGA at iteration $t$ of the recurrence (35) when either a solution $x_{t}\in x$ with $\delta_{t}\leq 0.1/t$ is achieved, or the number of CGA steps reaches a prescribed limit (set to 32 in the experiment to be reported).