Conditional Gradient Algorithms for Norm-Regularized Smooth Convex Optimization

Introduction

We consider two norm-regularized convex optimization problems as follows:

In the large-scale case, first-order algorithms of proximal-gradient type are popular to tackle such problems, see for a recent overview. Among them, the celebrated Nesterov optimal gradient methods for smooth and composite minimization , and their stochastic approximation counterparts , are now state-of-the-art in compressive sensing and machine learning. These algorithms enjoy the best known so far theoretical estimates (and in some cases, these estimates are the best possible for the first-order algorithms). For instance, Nesterov’s algorithm for penalized minimization solves (2) to accuracy $\epsilon$ in $O(D_{0}\sqrt{L/\epsilon})$ iterations, where $L$ is the properly defined Lipschitz constant of the gradient of $f$, and $D_{0}$ is the initial distance to the optimal set, measured in the norm $\|\cdot\|$. However, applicability and efficiency of proximal-gradient algorithms in the large-scale case require from the problem to possess “favorable geometry” (for details, see [24, Section A.6]). To be more specific, consider proximal-gradient algorithm for convex minimization problems of the form

The comments to follow, with slight modifications, are applicable to problems such as (1) and (2) as well. In this case, a proximal-gradient algorithm operates with a “distance generating function” (d.g.f.) defined on the domain of the problem and $1$-strongly convex w.r.t. the norm $\|\cdot\|$. Each step of the algorithm requires minimizing the sum of the d.g.f. and a linear form. The efficiency estimate of the algorithm depends on the variation of the d.g.f. on the domain and on regularity of $f$ w.r.t. $\|\cdot\|$ i.e., the Lipschitz constant of $f$ w.r.t. $\|\cdot\|$ in the nonsmooth case, or the Lipschitz constant of the gradient mapping $x\mapsto f^{\prime}(x)$ w.r.t. the norm $\|\cdot\|$ on the argument and the conjugate of this norm on the image spaces in the smooth case.. As a result, in order for a proximal-gradient algorithm to be practical in the large scale case, two “favorable geometry” conditions should be met: (a) the outlined sub-problems should be easy to solve, and (b) the variation of the d.g.f. on the domain of the problem should grow slowly (if at all) with problem’s dimension. Both these conditions indeed are met in many applications; see, e.g., for examples. This explains the recent popularity of this family of algorithms.

However, sometimes conditions (a) and/or (b) are violated, and application of proximal algorithms becomes questionable. For example, for the case of $K=E$, (b) is violated for the usual $\|\cdot\|_{\infty}$-norm on ${\mathbf{R}}^{p}$ or, more generally, for $\|\cdot\|_{2,1}$ norm on the space of $p\times q$ matrices given by

where $\hbox{\rm Row}_{j}^{T}(x)$ denotes the $j$-th row of $x$. Here the variation of (any) d.g.f. on problem’s domain is at least $p$. As a result, in the case in question the theoretical iteration complexity of a proximal algorithm grows rapidly with the dimension $p$. Furthermore, for some high-dimensional problems which do satisfy (b), solving the sub-problem can be computationally challenging. Examples of such problems include nuclear-norm-based matrix completion, Total Variation-based image reconstruction, and multi-task learning with a large number of tasks and features. This corresponds to $\|\cdot\|$ in (1) or (2) being the nuclear norm or the TV-norm.

These limitations recently motivated alternative approaches, which do not rely upon favorable geometry of the problem domain and/or do not require to solve hard sub-problems at each iteration, and triggered a renewed interest in the Conditional Gradient (CndG) algorithm. This algorithm, also known as the Frank-Wolfe algorithm, which is historically the first method for smooth constrained convex optimization, originates from , and was extensively studied in the 70-s (see, e.g., and references therein). CndG algorithms work by minimizing a linear form on the problem domain at each iteration; this auxiliary problem clearly is easier, and in many cases – significantly easier than the auxiliary problem arising in proximal-gradient algorithms. Conditional gradient algorithms for collaborative filtering were studied recently , some variants and extensions were studied in . Those works consider constrained formulations of machine learning or signal processing problems, i.e., minimizing the discrepancy $f(x)$ under a constraint on the norm of the solution, as in (3). On the other hand, CndG algorithms for other learning formulations, such as norm minimization (1) or penalized minimization (2) remain open issues. An exception is the work of , where a Conditional Gradient algorithm for penalized minimization was studied, although the efficiency estimates obtained in that paper were suboptimal. In this paper, we present CndG-type algorithms aimed at solving norm minimization and penalized norm minimization problems and provide theoretical efficiency guarantees for these algorithms.

The main body of the paper is organized as follows. In Section 2, we present detailed setting of problems (1), (2) along with basic assumptions on the “computational environment” required by the CndG-based algorithms we are developing. These algorithms and their efficiency bounds are presented in Sections 3 (problem (1)) and 5 (problem (2). In Section 6 we outline some applications, and in Section 7 present preliminary numerical results. All proofs are relegated to the appendix.

Problem statement

Throughout the paper, we shall assume that $K\subset E$ is a closed convex cone in Euclidean space $E$; we loose nothing by assuming that $K$ linearly spans $E$. We assume, further, that $\|\cdot\|$ is a norm on $E$, and $f:K\to{\mathbf{R}}$ is a convex function with Lipschitz continuous gradient, so that

where $\|\cdot\|_{*}$ denotes the norm dual to $\|\cdot\|$, whence

We consider two kinds of problems, detailed below.

To tackle (5), we consider the following parametric family of problems

Note that whenever (5) is feasible, which we assume from now on, we have

and both problems (5), (7) can be solved.

Given a tolerance $\epsilon>0$, we want to find an $\epsilon$-solution to the problem, that is, a pair $\rho_{\epsilon}$, $x_{\epsilon}\in K$ such that

where $X_{\rho}:=\{x\in E:\;x\in K,\;\|x\|\leq\rho\}$. Getting back to the problem of interest (5), $x_{\epsilon}$ is then “super-optimal” and $\epsilon$-feasible:

Penalized norm minimization.

We shall refer to (10) as the problem of composite optimization (CO). Given a tolerance $\epsilon>0$, our goal is to find an $\epsilon$-solution to (10), defined as a feasible solution $(x_{\epsilon},r_{\epsilon})$ to the problem satisfying $F([x_{\epsilon};r_{\epsilon}])-{\mathop{\hbox{Opt}}}\leq\epsilon$. Note that in this case $x_{\epsilon}$ is an $\epsilon$-solution, in the similar sense, to (9).

Special case.

In many applications where Problem (5) arise, (9) the function $f$ enjoys a special structure:

where $x\mapsto{\cal A}x-b$ is an affine mapping from $E$ to ${\mathbf{R}}^{m}$, and $\phi(\cdot):{\mathbf{R}}^{m}\to{\mathbf{R}}$ is a convex function with Lipschitz continuous gradient; we shall refer to this situation as to special case. In such case, the quantity $L_{f}$ can be bounded as follows. Let $\pi(\cdot)$ be some norm on ${\mathbf{R}}^{m}$, $\pi_{*}(\cdot)$ be the conjugate norm, and $\|{\cal A}\|_{\|\cdot\|,\pi}$ be the norm of the linear mapping $x\mapsto{\cal A}x$ induced by the norms $\|\cdot\|$, $\pi(\cdot)$ on the argument and the image spaces:

Let also $L_{\pi(\cdot)}[\phi]$ be the Lipschitz constant of the gradient of $\phi$ induced by the norm $\pi(\cdot)$, so that

Then, one can take as $L_{f}$ the quantity

Example 1: quadratic fit. In many applications, we are interested in $\|\cdot\|_{2}$-discrepancy between ${\cal A}x$ and $b$; the related choice of $\phi(\cdot)$ is $\phi(y)={1\over 2}y^{T}y$. Specifying $\pi(\cdot)$ as $\|\cdot\|_{2}$, we get $L_{\|\cdot\|_{2}}[\phi]=1$.

so that for $\beta=O(1)\ln(m)$ and $m$ large enough (specifically, such that $\beta\geq 2$), $\phi(y)$ is within absolute constant factor of ${1\over 2}\|y\|_{\infty}^{2}$. The latter situation can be interpreted as $\phi$ behaving as ${1\over 2}\|\cdot\|_{\infty}^{2}$). At the same time, with $\beta=O(1)\ln(m)$, $L_{\|\cdot\|_{\infty}}[\phi]\leq O(1)\ln(m)$ grows with $m$ logarithmically.

Another widely used choice of $\phi(\cdot)$ for this type of discrepancy is “logistic” function

For $\pi(\cdot)=\|\cdot\|_{\infty}$ we easily compute $L_{\|\cdot\|_{\infty}}[\phi]\leq\beta$ and $\|y\|_{\infty}\leq\phi(y)\leq\|y\|_{\infty}+\ln(2n)/\beta$.

Note that in some applications we are interested in “one-sided” discrepancies quantifying the magnitude of the vector $[{\cal A}x-b]_{+}=[\max[0,({\cal A}x-b)_{1}];...;\max[0,({\cal A}x-b)_{m}]]$ rather than the the magnitude of the vector ${\cal A}x-b$ itself. Here, instead of using \phi(y)=\mbox{\small\frac{1}{2}}\|y\|^{2}_{\beta} in the context of examples 1 and 2, one can use the functions $\phi_{+}(y)=\phi([y]_{+})$. In this case the bounds on $L_{\pi(\cdot)}[\phi_{+}]$ are exactly the same as the above bounds on $L_{\pi(\cdot)}[\phi]$. The obvious substitute for the two-sided logistic function is its “one-sided version:” $\phi_{+}(y)={1\over\beta}\ln\left(\sum_{i=1}^{m}\left[e^{\beta y_{i}}+1\right]\right)$ which obeys the same bound for $L_{\pi(\cdot)}[\phi_{+}]$ as its two-sided analogue.

First-order and Linear Optimization oracles.

We assume that $f$ is represented by a first-order oracle – a routine which, given on input a point $x\in K$, returns the value $f(x)$ and the gradient $f^{\prime}(x)$ of $f$ at $x$. As about $K$ and $\|\cdot\|$, we assume that they are given by a Linear Optimization (LO) oracle which, given on input a linear form $\langle\eta,\cdot\rangle$ on $E$, returns a minimizer $x[\eta]$ of this linear form on the set $\{x\in K:\|x\|\leq 1\}$. We assume w.l.o.g. that for every $\eta$, $x[\eta]$ is either zero, or is a vector of the $\|\cdot\|$-norm equal to 1. To ensure this property, it suffices to compute $\langle\eta,x[\eta]\rangle$ for $x[\eta]$ given by the oracle; if this inner product is 0, we can reset $x[\eta]=0$, otherwise $\|x[\eta]\|$ is automatically equal to 1.

Note that an LO oracle for $K$ and $\|\cdot\|$ allows to find a minimizer of a linear form of $z=[x;r]\in E^{+}:=E\times{\mathbf{R}}$ on a set of the form $K^{+}[\rho]=\{[x;r]\in E^{+}:x\in K,\|x\|\leq r\leq\rho\}$ due to the following observation:

Conditional Gradient algorithm

In this section, we present an overview of the properties of the standard Conditional Gradient algorithm, and highlight some memory-based extensions. These properties are not new. However, since they are key for the design of our proposed algorithms in the next sections, we present them for further reference.

Let $E$ be a Euclidean space and $X$ be a closed and bounded convex set in $E$ which linearly spans $E$. Assume that $X$ is given by a LO oracle – a routine which, given on input $\eta\in E$, returns an optimal solution $x_{X}[\eta]$ to the optimization problem

(cf. Section 2). Let $f$ be a convex differentiable function on $X$ with Lipschitz continuous gradient $f^{\prime}(x)$, so that

where $\|\cdot\|_{X}$ is the norm on $E$ with the unit ball $X-X$. We intend to solve the problem

A generic CndG algorithm is a recurrence which builds iterates $x_{t}\in X$, $t=1,2,...$, in such a way that

Basic implementations of a generic CndG algorithm are given by

in the sequel, we refer to them as CndGa and CndGb, respectively. As a byproduct of running generic CndG, after $t$ steps we have at our disposal the quantities

which, by convexity of $f$, are lower bounds on $f_{*}$. Consequently, at the end of step $t$ we have at our disposal a lower bound

Finally, we define the approximate solution $\bar{x}_{t}$ found in course of $t=1,2,...$ steps as the best – with the smallest value of $f$ – of the points $x_{1},...,x_{t}$. Note that $\bar{x}_{t}\in X$.

The following statement summarizes the well known properties of CndG (to make the presentation self-contained, we provide in Appendix the proof).

For a generic CndG algorithm, in particular, for both CndGa, CndGb, we have

Some remarks regarding the conditional algorithm are in order.

Certifying quality of approximate solutions. An attractive property of CndG is the presence of online lower bound $f_{*}^{t}$ on $f_{*}$ which certifies the theoretical rate of convergence of the algorithm, see (20). This accuracy certificate, first established in , also provides a valuable stopping criterion when running the algorithm in practice.

CndG algorithm with memory. When computing the next search point $x_{t+1}$ the simplest CndG algorithm CndGa only uses the latest answer $x_{t}^{+}=x_{X}[f^{\prime}(x_{t})]$ of the LO oracle. Meanwhile, algorithm CndGb can be modified to make use of information supplied by previous oracle calls; we refer to this modification as CndG with memory (CndGM). Note that in the context of “classical” Frank-Wolfe algorithm – minimization of a smooth function over a polyhedral set – such modification is referred to as Restricted Simplicial Decomposition . Assume that we have already carried out $t-1$ steps of the algorithm and have at our disposal current iterate $x_{t}\in X$ (with $x_{1}$ selected as an arbitrary point of $X$) along with previous iterates $x_{\tau}$, $\tau<t$ and the vectors $f^{\prime}(x_{\tau})$, $x_{\tau}^{+}=x_{X}[f^{\prime}(x_{\tau})]$. At the step, we compute $f^{\prime}(x_{t})$ and $x^{+}_{t}=x_{X}[f^{\prime}(x_{t})]$. Thus, at this point in time we have at our disposal $2t$ points $x_{\tau},x^{+}_{\tau}$, $1\leq\tau\leq t$, which belong to $X$. Let $X_{t}$ be subset of these points, with the only restriction that the points $x_{t}$, $x^{+}_{t}$ are selected, and let us define the next iterate $x_{t+1}$ as

Clearly, it is again a generic CndG algorithm, so that conclusions in Theorem 1 are fully applicable to CndGM. Note that CndGb per se is nothing but CndGM with $X_{t}=\{x_{t},x_{t}^{+}\}$ and $M=2$ for all $t$.

CndGM: implementation issues. Assume that the cardinalities of the sets $X_{t}$ in CndGM are bounded by some $M\geq 2$. In this case, implementation of the method requires solving at every step an auxiliary problem (22) of minimizing over the standard simplex of dimension $\leq M-1$ a smooth convex function given by a first-order oracle induced by the first-oracle for $f$. When $M$ is a once for ever fixed small integer, the arithmetic cost of solving this problem within machine accuracy by, say, the Ellipsoid algorithm is dominated by the arithmetic cost of just $O(1)$ calls to the first-order oracle for $f$. Thus, CndGM with small $M$ can be considered as implementable Assuming possibility to solve (22) exactly, while being idealization, is basically as “tolerable” as the standard in continuous optimization assumption that one can use exact real arithmetic or compute exactly eigenvalues/eigenvectors of symmetric matrices. The outlined “real life” considerations can be replaced with rigorous error analysis which shows that in order to maintain the efficiency estimates from Theorem 1, it suffices to solve $t$-th auxiliary problem within properly selected positive inaccuracy, and this can be achieved in $O(\ln(t))$ computations of $f$ and $f^{\prime}$..

Note that in the special case (Section 2), where $f(x)=\phi(Ax-b)$, assuming $\phi(\cdot)$ and $\phi^{\prime}(\cdot)$ easy to compute, as is the case in most of the applications, the first-order oracle for the auxiliary problems arising in CndGM becomes cheap (cf. ). Indeed, in this case (22) reads

It follows that all we need to get a computationally cheap access to the first-order information on $g_{t}(\lambda^{t})$ for all values of $\lambda^{t}$ is to have at our disposal the matrix-vector products $Ax$, $x\in X_{t}$. With our construction of $X_{t}$, the only two “new” elements in $X_{t}$ (those which were not available at preceding iterations) are $x_{t}$ and $x^{+}_{t}$, so that the only two new matrix-vector products we need to compute at iteration $t$ are $Ax_{t}$ (which usually is a byproduct of computing $f^{\prime}(x_{t})$) and $Ax_{t}^{+}$. Thus, we can say that the “computational overhead,” as compared to computing $f^{\prime}(x_{t})$ and $x_{t}^{+}=x_{X}[f^{\prime}(x_{t})]$, needed to get easy access to the first-order information on $g_{t}(\cdot)$ reduces to computing the single matrix-vector product $Ax_{t}^{+}$.

Conditional gradient algorithm for parametric optimization

In this section, we describe a multi-stage algorithm to solve the parametric optimization problem (6), (7), using the conditional algorithm to solve inner sub-problems. (6), (7). The idea, originating from (see also ), is to use a Newton-type method for approximating from below the positive root $\rho_{*}$ of ${\mathop{\hbox{Opt}}}(\rho)$, with (inexact) first-order information on ${\mathop{\hbox{Opt}}}(\cdot)$ yielded by approximate solving the optimization problems defining ${\mathop{\hbox{Opt}}}(\cdot)$; the difference with the outlined references is that now we solve these problems with the CndG algorithm.

Our algorithm works stagewise. At the beginning of stage $s=1,2,...$, we have at hand a lower bound $\rho_{s}$ on $\rho_{*}$, with $\rho_{1}$ defined as follows:

We compute $f(0)$, $f^{\prime}(0)$ and $x[f^{\prime}(0)]$. If $f(0)\leq\epsilon$ or $x[f^{\prime}(0)]=0$, we are done — the pair ($\rho=0$, $x=0$) is an $\epsilon$-solution to (7) in the first case, and is an optimal solution to the problem in the second case (since in the latter case $0$ is a minimizer of $f$ on $K$, and (7) is feasible). Assume from now on that the above options do not take place (“nontrivial case”), and let

Due to the origin of $x[\cdot]$, $d$ is positive, and $f(x)\geq f(0)+\langle f^{\prime}(0),x\rangle\geq f(0)-d\|x\|$ for all $x\in K$, which implies that $\rho_{*}\geq\rho_{1}:={f(0)\over d}>0.$

At stage $s$ we apply a generic CndG algorithm (e.g., CndGa,CndGb, or CndGM; in the sequel, we refer to the algorithm we use as to CndG) to the auxiliary problem

Note that the LO oracle for $K$, $\|\cdot\|$ induces an LO oracle for $K[\rho]$; specifically, for every $\eta\in E$, the point $x_{\rho}[\eta]:=\rho x[\eta]$ is a minimizer of the linear form $\langle\eta,x\rangle$ over $x\in K[\rho]$, see Lemma 1. $x_{\rho}[\cdot]$ is exactly the LO oracle utilized by CndG as applied to (23).

As explained above, after $t$ steps of CndG as applied to (23), the iterates being $x_{\tau}\in K[\rho_{s}]$, $1\leq\tau\leq t$ The iterates $x_{t}$, same as other indexed by $t$ quantities participating in the description of the algorithm, in fact depend on both $t$ and the stage number $s$. To avoid cumbersome notation when speaking about a particular stage, we suppress $s$ in the notation., we have at our disposal current approximate solution $\bar{x}_{t}\in\{x_{1},...,x_{t}\}$ such that $f(\bar{x}_{t})=\min_{1\leq\tau\leq t}f(x_{\tau})$ along with a lower bound $f_{*}^{t}$ on ${\mathop{\hbox{Opt}}}(\rho_{s})$. Our policy is as follows.

When $f(\bar{x}_{t})\leq\epsilon$, we terminate the solution process and output $\bar{\rho}=\rho_{s}$ and $\bar{x}=\bar{x}_{t}$;

When the above option is not met and $f_{*}^{t}<{3\over 4}f(\bar{x}_{t})$, we specify $x_{t+1}$ according to the description of CndG and pass to step $t+1$ of stage $s$;

Finally, when neither one of the above options takes place, we terminate stage $s$ and pass to stage $s+1$, specifying $\rho_{s+1}$ as follows: We are in the situation $f(\bar{x}_{t})>\epsilon$ and $f_{*}^{t}\geq{3\over 4}f(\bar{x}_{t})$. Now, for $k\leq t$ the quantities $f(x_{k})$, $f^{\prime}(x_{k})$ and $x[f^{\prime}(x_{k})]$ define affine function of $\rho\geq 0$

By Lemma 1 we have for every $\rho\geq 0$

is well defined and satisfies $\rho_{s}<r^{t}\leq\rho_{*}$. We compute $r^{t}$ (which is easy) and pass to stage $s+1$, setting $\rho_{s+1}=r^{t}$ and selecting, as the first iterate of the new stage, any point known to belong to $K[\rho]$ (e.g., the origin, or $\bar{x}_{t}$). The first iterate of the first stage is $0$.

The description of the algorithm is complete.

The complexity properties of the algorithm are given by the following proposition.

When solving a PO problem (6), (7) by the outlined algorithm,

(i) the algorithm terminates with an $\epsilon$-solution, as defined in Section 2 (cf. (8));

(ii) The number $N_{s}$ of steps at every stage $s$ of the method admits the bound

(iii) The number of stages before termination does not exceed the quantity

Conditional Gradient algorithm for Composite Optimization

In this section, we present a modification of the CndG algorithm capable to solve composite minimization problem (10). We assume in the sequel that $\|\cdot\|,K$ are represented by an LO oracle for the set $\{x\in K:\|x\|\leq 1\}$, and $f$ is given by a first order oracle. In order to apply CndG to the composite optimization problem (10), we make the assumption as follows:

Assumption A: There exists $D<\infty$ such that $\kappa r+f(x)\leq f(0)$ together with $\|x\|\leq r$, $x\in K$, imply that $r\leq D$.

We define $D_{*}$ as the minimal value of $D$ satisfying Assumption A, and assume that we have at our disposal a finite upper bound $D^{+}$ on $D_{*}$. An important property of the algorithm we are about to develop is that its efficiency estimate depends on the induced by problem’s data quantity $D_{*}$, and is independent of our a priori upper bound $D^{+}$ on this quantity, see Theorem 3 below.

We are about to present an algorithm for solving (10). Let $E^{+}=E\times{\mathbf{R}}$, and $K^{+}=\{[x;r]:\;x\in K,\,\|x\|\leq r\}$. From now on, for a point $z=[x;r]\in E^{+}$ we set $x(z)=x$ and $r(z)=r$. Given $z=[x;r]\in K^{+}$, let us consider the segment

Observe that by Lemma 1, for every $0\leq\rho\leq D^{+}$, the minimum of this form on $K^{+}[\rho]=\{[x;r]\in E^{+},x\in K,\|x\|\leq r\leq\rho\}$ is attained at a point of $\Delta(z)$ (either at $[\rho x[f^{\prime}(x)];\,\rho]$ or at the origin). A generic Conditional Gradient algorithm for composite optimization (COCndG) is a recurrence which builds the points $z_{t}=[x_{t};r_{t}]\in K^{+}$, $t=1,2,...$, in such a way that

Let $z_{*}=[x_{*};r_{*}]$ be an optimal solution to (10) (which under Assumption A clearly exists), and let $F_{*}=F(z_{*})$ (i.e., $F_{*}$ is nothing but ${\mathop{\hbox{Opt}}}$, see (9)).

A generic COCndG algorithm (24) maintains the inclusions $z_{t}\in K^{+}$ and is a descent algorithm: $F(z_{t+1})\leq F(z_{t})$ for all $t$. Besides this, we have

COCndG with memory.

The simplest implementation of a generic COCndG algorithm is given by the recurrence

Denoting $\widehat{z}_{\tau}:=D^{+}[x[f^{\prime}(x_{\tau})];1]$, the recurrence can be written

As for the CndG algorithm in section 3, the recurrence (26) admits a version with memory COCndGM still obeying (24) and thus sartisfying the conclusion of Theorem 3. Specifically, assume that we already have built $t$ iterates $z_{\tau}=[x_{\tau};r_{\tau}]\in K^{+}$, $1\leq\tau\leq t$, with $z_{1}=0$, along with the gradients $f^{\prime}(x_{\tau})$ and the points $x[f^{\prime}(x_{\tau})]$. Then we have at our disposal a number of points from $K^{+}$, namely, the iterates $z_{\tau}$, $\tau\leq t$, and the points $\widehat{z}_{\tau}=D^{+}[x[f^{\prime}(x_{\tau})];1]$. Let us select a subset $Z_{t}$ of the set $\{z_{\tau},\widehat{z}_{\tau},1\leq\tau\leq t\}$, with the only restriction that $Z_{t}$ contains the points $z_{t},\widehat{z}_{t}$, and set

Since $z_{t},\widehat{z}_{t}\in Z_{t}$, we have $\hbox{\rm Conv}\left(\Delta(z_{t})\cup\{z_{t}\}\right)\}\subset{\cal C}_{t}$, whence the procedure we have outlined is an implementation of generic COCndG algorithm. Note that the basic COCndG algorithm is the particular case of the COCndGM corresponding to the case where $Z_{t}=\{z_{t},\widehat{z}_{t}\}$ for all $t$. The discussion of implementability of CndGM in section 3 fully applies to COCndGM.

Let us outline several options which can be implemented in COCndGM; while preserving the theoretical efficiency estimates stated in Theorem 3 they can improve the practical performance of the algorithm. For the sake of definiteness, let us focus on the case of quadratic $f$: $f(x)=\|{\cal A}x-b\|_{2}^{2}$, with ${\rm Ker}{\cal A}=\{0\}$; extensions to a more general case are straightforward.

We lose nothing (and potentially gain) when extending ${\cal C}_{t}$ in (28) to the conic hull

of $Z_{t}$. When $K=E$, we can go further and replace (28) with

Note that the preceding “conic case” is obtained from (29) by adding to the constraints of the right hand side problem the inequalities $\lambda_{\zeta}\geq 0,\zeta\in Z_{t}$. Finally, when $\|\cdot\|$ is easy to compute, we can improve (29) to

(the definition of $\lambda^{*}$ assumes that $K=E$, otherwise the constraints of the problem specifying $\lambda^{*}$ should be augmented by the inequalities $\lambda_{\zeta}\geq 0,\zeta\in Z_{t}$).

In the case of quadratic $f$ and moderate cardinality of $Z_{t}$, optimization problems arising in (29) (with or without added constraints $\lambda_{\zeta}\geq 0$) are explicitly given low-dimensional “nearly quadratic” convex problems which can be solved to high accuracy “in no time” by interior point solvers. With this in mind, we could solve these problems for the given value of the penalty parameter $\kappa$ and also for several other values of the parameter. Thus, at every iteration we get feasible approximate solution to several instances of (9) for different values of the penalty parameter. Assume that we keep in memory, for every value of the penalty parameter in question, the best, in terms of the respective objective, of the related approximate solutions found so far. Then upon termination we will have at our disposal, along with the feasible approximate solution associated with the given value of the penalty parameter, provably obeying the efficiency estimates of Theorem 3, a set of feasible approximate solutions to the instances of (9) corresponding to other values of the penalty.

In the above description, $Z_{t}$ was assumed to be a subset of the set $Z^{t}=\{z_{\tau}=[x_{\tau};r_{\tau}],\widehat{z}_{\tau},\,1\leq\tau\leq t\}$ containing $z_{t}$ and $\widehat{z}_{t}$. Under the latter restriction, we lose nothing when allowing for $Z_{t}$ to contain points from $K^{+}\backslash Z^{t}$ as well. For instance, when $K=E$ and $\|\cdot\|$ is easy to compute, we can add to $Z_{t}$ the point $z_{t}^{\prime}=[f^{\prime}(x_{t});\|f^{\prime}(x_{t})\|]$. Assume, e.g., that we fix in advance the cardinality $M\geq 3$ of $Z_{t}$ and define $Z_{t}$ as follows: to get $Z_{t}$ from $Z_{t-1}$, we eliminate from the latter set several (the less, the better) points to get a set of cardinality $\leq M-3$, and then add to the resulting set the points $z_{t}$, $\widehat{z}_{t}$ and $z_{t}^{\prime}$. Eliminating the points according to the rule “first in – first out,” the projection of the feasible set of the optimization problem in (30) onto the space of $x$-variables will be a linear subspace of $E$ containing, starting with step $t=M$, at least $\lfloor M/3\rfloor$ (here $\lfloor a\rfloor$ stands for the largest integer not larger than $a$) of gradients of $f$ taken at the latest iterates, so that the method, modulo the influence of the penalty term, becomes a “truncated” version of the Conjugate Gradient algorithm for quadratic minimization. Due to nice convergence properties of Conjugate Gradient in the quadratic case, one can hope that a modification of this type will improve significantly the practical performance of COCndGM.

Application examples

In this section, we detail how the proposed conditional gradient algorithms apply to several examples. In particular, we detail the corresponding LO oracles, and how one could implement these oracles efficiently.

The first example where the proposed algorithms seem to be more attractive than the proximal methods are large-scale problems (5), (9) on the space of $p\times q$ matrices $E={\mathbf{R}}^{p\times q}$ associated with the nuclear norm $\|\sigma(x)\|_{1}$ of a matrix $x$, where $\sigma(x)=[\sigma_{1}(x);...;\sigma_{\min[p,q]}(x)]$ is the vector of singular values of a $p\times q$ matrix $x$. Problems of this type with $K=E$ arise in various versions of matrix completion, where the goal is to recover a matrix $x$ from its noisy linear image $y={\cal A}x+\xi$, so that $f=\phi({\cal A}x-y)$, with some smooth and convex discrepancy measure $\phi(\cdot)$, most notably, $\phi(z)={1\over 2}\|z\|_{2}^{2}$. In this case, $\|\cdot\|$ minimization/penalization is aimed at getting a recovery of low rank ( and references therein). Another series of applications relates to the case when $E={\mathbf{S}}^{p}$ is the space of symmetric $p\times p$ matrices, and $K={\mathbf{S}}^{p}_{+}$ is the cone of positive semidefinite matrices, with $f$ and $\phi$ as above; this setup corresponds to the situation when one wants to recover a covariance (and thus positive semidefinite symmetric) matrix from experimental data. Restricted from ${\mathbf{R}}^{p\times p}$ onto ${\mathbf{S}}^{p}$, the nuclear norm becomes the trace norm $\|\lambda(x)\|_{1}$, where $\lambda(x)\in{\mathbf{R}}^{p}$ is the vector of eigenvalues of a symmetric $p\times p$ matrix $x$, and regularization by this norm is, as above, aimed at building a low rank recovery.

With the nuclear (or trace) norm in the role of $\|\cdot\|$, all known proximal algorithms require, at least in theory, computing at every iteration the complete singular value decomposition of $p\times q$ matrix $x$ (resp., complete eigenvalue decomposition of a symmetric $p\times p$ matrix $x$), which for large $p,q$ may become prohibitively time consuming. In contrast to this, with $K=E$ and $\|\cdot\|=\|\sigma(\cdot)\|_{1}$, LO oracle for $(K,\|\cdot\|=\|\sigma(\cdot)\|_{1})$ only requires computing the leading right singular vector $e$ of a $p\times q$ matrix $\eta$ (i.e., the leading eigenvector of $\eta^{T}\eta$): $x[\eta]=-\bar{f}\bar{e}^{T}$, where $\bar{e}=e/\|e\|_{2}$ and $\bar{f}=\eta e/\|\eta e\|_{2}$ for nonzero $\eta$ and $\bar{f}=0$, $\bar{e}=0$ when $\eta=0$. Computing the leading singular vector of a large matrix is, in most cases, much cheaper than computing the complete eigenvalue decomposition of the matrix. Similarly, in the case of $E={\mathbf{S}}^{p}$, $K={\mathbf{S}}^{p}_{+}$ and the trace norm in the role of $\|\cdot\|$, LO oracle requires computing the leading eigenvector $e$ of a matrix $\eta\in{\mathbf{S}}^{p}$: $x[-\eta]=\bar{e}\bar{e}^{T}$, where $\bar{e}=0$ when $e^{T}\eta e\geq 0$, and $\bar{e}=e/\|e\|_{2}$ otherwise. Here again, for a large symmetric $p\times p$ matrix, the required computation usually is much easier than computing the complete eigenvalue decomposition of such a matrix. As a result, in the situations under consideration, algorithms based on the LO oracle remain “practically implementable” in an essentially larger range of problem sizes than proximal methods.

An additional attractive property of the CndG algorithms we have described stems from the fact that since in the situations in question the matrices $x[\eta]$ are of rank 1, $t$-th approximate solution $x_{t}$ yielded by the CndG algorithms for composite minimization from Section 5 is of rank at most $t$. Similar statement holds true for $t$-th approximate solution $x_{t}$ built at a stage of a CndG algorithm for parametric optimization from Section 3, provided that the first iterate at every stage is the zero matrix. this property is an immediate corollary of the fact that in the situation in question, by description of the algorithms $x_{t}$ is a convex combination of $t$ points of the form $x[\cdot]$..

2 Regularization by Total Variation

Note that ${\hbox{\rm TV}}(\cdot)$ is a norm on the subspace $M^{n}_{0}$ of $M^{n}$ comprised of zero mean images $x$ (those with $\sum_{i,j}x(i,j)=0$) and vanishes on the orthogonal complement to $M_{0}^{n}$, comprised of constant images.

Consider the network (the oriented graph) $G$ with $n^{2}$ nodes $[i;j]\in\Gamma_{n,n}$ and $2n(n-1)$ arcs as follows: the first $n(n-1)$ arcs are of the form $([i+1;j],[i;j])$, $0\leq i<n-1$, $0\leq j<n$, the next $n(n-1)$ arcs are $([i;j+1],[i;j])$, $0\leq i<n$, $0\leq j<n-1$, and the remaining $2n(n-1)$ arcs (let us call them backward arcs) are the inverses of the just defined $2n(n-1)$ forward arcs. Let ${\cal E}$ be the set of arcs of our network, and let us equip all the arcs with unit capacities. Let us treat vectors from $E=M^{n}_{0}$ as vectors of external supplies for our network; note that the entries of these vectors sum to zero, as required from external supply. Now, given a nonzero vector $\eta\in M^{n}_{0}$, let us consider the network flow problem where we seek for the largest multiple $s\eta$ of $\eta$ which, considered as the vector of external supplies in our network, results in a feasible capacitated network flow problem. The problem in question reads

where $P$ is the incidence matrix of our network that is, the rows of $P$ are indexed by the nodes, the columns are indexed by the arcs, and in the column indexed by an arc $\gamma$ there are exactly two nonzero entries: entry 1 in the row indexed by the starting node of $\gamma$, and entry $-1$ in the row indexed by the terminal node of $\gamma$. and $\mathbf{e}$ is the all-ones vector. Now, problem (31) clearly is feasible, and its feasible set is bounded due to $\eta\neq 0$, so that the problem is solvable. Due to its network structure, this LP program can be solved reasonably fast even in the large scale case (say, when $n=512$ or $n=1024$, which already is of interest for actual imaging). Further, an intelligent network flow solver as applied to (31) will return not only the optimal $s=s_{*}$ and the corresponding flow, but also the dual information, in particular, the optimal vector $z$ of Lagrange multipliers for the linear equality constraints $Pr-s\eta=0$. Let $\bar{z}$ be obtained by subtracting from the entries of $z$ their mean; since the entries of $z$ are indexed by the nodes, $\bar{z}$ can be naturally interpreted as a zero mean image. It turns out that this image is nonzero, and the vector $x[\eta]=-\bar{z}/{\hbox{\rm TV}}(\bar{z})$ is nothing than a desired minimizer of $\langle\eta,\cdot\rangle$ on ${{\cal T}{\cal V}}$:

Let $\eta$ be a nonzero image with zero mean. Then (31) is solvable with positive optimal value, and the image $x[\eta]$, as defined above, is well defined and is a maximizer of $\langle\eta,\cdot\rangle$ on ${{\cal T}{\cal V}}$.

When applying CndG algorithms to the TV-based problems (5), (9) with $E=M^{n}_{0}$ and $f(x)=\phi({\cal A}x-b)$, the efficiency estimates depend linearly on the associated quantity $L_{f}$, which, in turn, is readily given by the norm $\|{\cal A}\|_{{\hbox{\scriptsize TV}}(\cdot),\pi(\cdot)}$ of the mapping $x\mapsto{\cal A}x$, see the end of Section 2. Observe that in typical applications ${\cal A}$ is a simple operator (e.g., the discrete convolution), so that when restricting ourselves to the case when $\pi(\cdot)$ is $\|\cdot\|_{2}$ (quadratic fit), it is easy to find a tight upper bound on $\|{\cal A}\|_{\|\cdot\|_{2},\|\cdot\|_{2}}$. To convert this bound into an upper bound on $\|{\cal A}\|_{{\hbox{\scriptsize TV}}(\cdot),\|\cdot\|_{2}}$, we need to estimate the quantity

Bounding $Q_{n}$ is not a completely trivial question, and the answer is as follows:

$Q_{n}$ is nearly constant, specifically, $Q_{n}\leq O(1)\sqrt{\ln(n)}$ with a properly selected absolute constant $O(1)$.

Note that the result of Proposition 1 is in sharp contrast with one-dimensional case, where the natural analogy of $Q_{n}$ grows with $n$ as $\sqrt{n}$. We do not know whether it is possible to replace in Proposition 1 $O(1)\sqrt{\ln(n)}$ with $O(1)$, as suggested by Sobolev’s inequalities From the Sobolev embedding theorem it follows that for a smooth function $f(x,y)$ on the unit square one has $\|f\|_{L_{2}}\leq O(1)\|\nabla f\|_{1},$ $\|\nabla f\|_{1}:=\|f^{\prime}_{x}\|_{1}+\|f^{\prime}_{y}\|_{1}$, provided that $f$ has zero mean. Denoting by $f^{n}$ the restriction of the function onto a $n\times n$ regular grid in the square, we conclude that $\|f^{n}\|_{2}/{\hbox{\rm TV}}(f^{n})\to\|f\|_{L_{2}}/\|\nabla f\|_{1}\leq O(1)$ as $n\to\infty$. Note that the convergence in question takes place only in the 2-dimensional case.. Note that on inspection of the proof, Proposition extends to the case of $d$-dimensional, $d>2$, images with zero mean, in which case $Q_{n}\leq C(d)$ with appropriately chosen $C(d)$.

Numerical examples

We present here some very preliminary simulation results.

The goal of the first series of our experiments is to illustrate how the performance and requirements of CndG algorithm for parametric optimization, when applied to the matrix completion problem , scale with problem size. Specifically, we apply the algorithm of Section 4 to the problem of nuclear norm minimization

where $\sigma(x)$ is the singular spectrum of a $p\times q$ matrix $x$. In our experiments, the set $\Omega$ of observed entries $(i,j)\in\{1,...,p\}\times\{1,...,q\}$ of cardinality $m\ll pq$ was selected at random.

Note that the the implementation of the CndGM is especially simple for the problem (32) – at each method’s iteration it requires solving a simple quadratic problem with dimension of the decision variable which does not exceed the iteration count. This allows to implement efficiently the “full memory” version of CndGM (CndG algorithms with memory) (21), (22), in which the set $X_{t}$ contains $x_{t}$ and all the points $x^{+}_{\tau}$ for $1\leq\tau\leq t$.

We compare the performance of CndGM algorithms and of a “memoryless” version of the CndG. To this end we have conducted the following experiment:

For matrix sizes $p,q\in\times 10^{3}$ we generate $n=10$ sparse $p\times q$ matrices $y$ with density $d=0.1$ of non-vanishing entries as follows: we generate $p\times r$ matrix $U$ and $q\times r$ matrix $V$ with independent Gaussian entries $u_{ij}\sim{\cal N}(0,m^{-1}),\;v_{ij}\sim{\cal N}(0,n^{-1})$, and a $r\times r$ diagonal matrix $D={\rm diag}[d_{1},...,d_{r}]$ with $d_{i}$ drawn independently from a uniform distribution on $$. The non-vanishing entries of the sparse observation matrix$y$are obtained by sampling at random with probability$d$the entries of$x^{*}=UDV^{T}$, so that for every$i,j$,$y_{ij}$is, independently over$i,j$, set to$x^{*}_{ij}$with probability$d$and to$0$with probability$1-d$. Thus, the number of non-vanishing entries of$y$is approximately$m=dpq$. This procedure results in$m\sim 10^{5}$for the smallest matrices$y$($1000\times 1000$), and in$m\sim 10^{8}$for the largest matrices ($32000\times 32000$).

We apply to parametric optimization problem (32) MATLAB implementations of the CndGM with memory parameter $M=1$ (“memoryless” CndG), CndGM with $M=5$ and full memory CndGM. The parameter $\delta$ of (32) is chosen to be $\delta=0.001\|y\|^{2}_{\rm f}$ (here $\|y\|_{\rm f}=\left(\sum_{i,j}y^{2}_{ij}\right)^{1/2}$ stands for the Frobenius norm of $y$). The optimization algorithm is tuned to the relative accuracy $\varepsilon=1/4$, what means that it outputs an $\epsilon$-solution $\widehat{x}$ to (32), in the sense of (8), with absolute accuracy $\epsilon=\delta\varepsilon$.

For each algorithm (memoryless CndG, CndGM with memory $M=5$ and full memory CndGM) we present in table 1 the average, over algorithm’s runs on the (common for all algorithms) sample of $n=10$ matrices $y$ we have generated, 1) total number of iterations $N_{\rm it}$ necessary to produce an $\epsilon$-solution (it upper-bounds the rank of the resulting $\epsilon$-solutuion), 2) CPU time in seconds $T_{\rm cpu}$ and 3) MATLAB memory usage in megabytes $S_{\rm mem}$. This experiment was conducted on a Dell Latitude 6430 laptop equipped with Intel Core i7-3720QM CPU@2.60GHz and 16GB of RAM. Because of high memory requirements in our implementation of the full memory CndGM, this method was unable to complete the computation for the two largest matrix sizes.

We can make the following observation regarding the results summarized in table 1: CndG algorithm with memory consistently outperforms the standard – memoryless – version of CndG. The full memory CndGM requires the smallest number of iteration to produce an $\epsilon$-solution, which is of the smallest rank, as a result. On the other hand, the memory requirements of the full memory CndGM become prohibitive (at least, for the computer we used for this experiment and MATLAB implementation of the memory heap) for large matrices. On the other hand, a CndGM with memory $M=5$ appears to be a reasonable compromise in terms of numerical efficiency and memory demand.

2 CndG for composite optimization: multi-class classification with nuclear-norm regularization

We present here an empirical study of the CndG algorithm for composite optimization as applied to the machine learning problem of multi-class classification with nuclear-norm penalty. A brief description of the multi-class classification problem is as follows: we observe $N$ “feature vectors” $\xi_{i}\in{\mathbf{R}}^{q}$, each belonging to exactly one of $p$ classes $C_{1},...,C_{p}$. Each $\xi_{i}$ is augmented by its label $y_{i}\in\{1,...,p\}$ indicating to which class $\xi_{i}$ belongs. Our goal is to build a classifier capable to predict the class to which a new feature vector $\xi$ belongs. This classifier is given by a $p\times q$ matrix $x$ according to the following rule: given $\xi$, we compute the $p$-dimensional vector $x\xi$ and take, as the guessed class of $\xi$, the index of the largest entry in this vector.

In some cases (see ), when, for instance, one is dealing with a large number of classes, there are good reasons “to train the classifier” — to specify $x$ given the training sample $(\xi_{i},y_{i})$, $1\leq i\leq N$ — as the optimal solution to the nuclear norm penalized minimization problem

Below, we report on some experiments with this problem. Our goal was to compare two versions of CndG for composite minimization: the memoryless version defined in (24) and the version with memory defined in (28). To solve the corresponding sub-problems, we used the Center of Gravity method in the case of (24) and the Ellipsoid method in the case of (28) . In the version with memory we set $M=5$, as it appeared to be the best option from empirical evidence. We have considered the following datasets:

Simulated data: for matrix of sizes $p,q\in 10^{3}\times\{2^{s}\}_{s=1}^{4}$, we generate random matrices $x_{\star}=USV$, with $p\times p$ factor $U$, $q\times q$ factor $V$, and diagonal $p\times q$ factor $S$ with random entries sampled, independently of each other, from ${\cal N}(0,p^{-1})$ (for $U$), ${\cal N}(0,q^{-1})$ (for $V$), and the uniform distribution on $$(for diagonal entries in$S$). We use$N=20q$, with the feature vectors$\xi_{1},...,\xi_{N}$sampled, independently of each other, from the distribution${\cal N}(0,I_{q})$, and their labels$y_{i}$being the indexes of the largest entries in the vectors$x_{\star}\xi_{i}+\epsilon_{i}$, where$\epsilon_{i}\in{\mathbf{R}}^{p}$were sampled, independently of each other and of$\xi_{1},...,\xi_{N}$, from${\cal N}(0,{1\over 2}I_{p})$. The regularization parameter$\kappa$is set to$10^{-3}{\mathop{\hbox{\rm Tr}}}(x_{\star}x_{\star}^{T})$.

In both sets of experiments, the computations are terminated when the “$\epsilon$-optimality conditions”

were met, where $\|\sigma(\cdot)\|_{\infty}$ denotes the usual operator norm (the largest singular value). These conditions admit transparent interpretation as follows. For every $\bar{x}$, the function

underestimates $F_{\kappa}(x)$, see (33), whence ${\mathop{\hbox{Opt}}}(\kappa^{\prime})\geq f(\bar{x})-\langle f^{\prime}(\bar{x}),\bar{x}\rangle$ whenever $\kappa^{\prime}\geq\|\sigma(f^{\prime}(\bar{x}))\|_{\infty}$. Thus, whenever $\bar{x}=x_{t}$ satisfies the first relation in (34), we have ${\mathop{\hbox{Opt}}}(\kappa+\epsilon)\geq f(x_{t})-\langle f^{\prime}(x_{t}),x_{t}\rangle$, whence

We see that (34) ensures that $F_{\kappa}(x_{t})-{\mathop{\hbox{Opt}}}(\kappa+\epsilon)\leq\epsilon\;\|\sigma(x_{t})\|_{1}$, which, for small $\epsilon$, is a reasonable substitute for the actually desired termination when $F_{\kappa}(x_{t})-{\mathop{\hbox{Opt}}}(\kappa)$ becomes small. In our experiments, we use $\epsilon=0.001$.

In table 2 for each algorithm (memoryless CndG, CndGM with memory $M=5$) we present the average, over 20 collections of simulated data coming from 20 realizations of $x_{\star}$, of: 1) total number of iterations $N_{\rm it}$ necessary to produce an $\epsilon$-solution, 2) CPU time in seconds $T_{\rm cpu}$. The last row of the table corresponds to the real-world data. Experiments were conducted on a Dell R905 server equipped with four six-core AMD Opteron 2.80GHz CPUs and 64GB of RAM. A maximum of 32GB of RAM was used for the computations.

We draw the following conclusions from table 1: CndG algorithm with memory routinely outperforms the standard – memoryless – version of CndG. However, there is a trade-off between the algorithm progress at each iteration and the computational load of each iteration. Note that, for large $M$, solving the sub-problem (28) can be challenging.

3 CndG for composite optimization: TV-regularized image reconstruction

Here we report on experiments with COCndGM as applied to TV-regularized image reconstruction. Our problem of interest is of the form (9) with quadratic $f$, namely, the problem

In our experiments, the mapping $x\mapsto{\cal A}x$ is defined as follows: we zero-pad $x$ to extend it from $\Gamma_{n,n}$ to get a finitely supported function on ${\mathbf{Z}}^{2}$, then convolve this function with a finitely supported kernel $\alpha(\cdot)$, and restrict the result onto $\Gamma_{n,n}$. The observations $b\in M^{n}$ were generated at random according to

with mutually independent $\xi_{ij}$. The relative noise intensity $\sigma>0$, same as the convolution kernel $\alpha(\cdot)$, are parameters of the setup of an experiment.