Primal-Dual Block Frank-Wolfe

TL;DR

This paper introduces a Primal-Dual Block Frank-Wolfe algorithm that efficiently solves sparse and low-rank optimization problems, including Elastic Net and regularized SVMs, with reduced per-iteration cost and proven linear convergence.

Contribution

It presents a novel Frank-Wolfe variant that lowers computational cost per iteration while maintaining convergence speed, applicable to various structured optimization problems.

Findings

Outperforms state-of-the-art methods on classification tasks

Reduces per-iteration computational cost based on solution structure

Maintains linear convergence rate

Abstract

We propose a variant of the Frank-Wolfe algorithm for solving a class of sparse/low-rank optimization problems. Our formulation includes Elastic Net, regularized SVMs and phase retrieval as special cases. The proposed Primal-Dual Block Frank-Wolfe algorithm reduces the per-iteration cost while maintaining linear convergence rate. The per iteration cost of our method depends on the structural complexity of the solution (i.e. sparsity/low-rank) instead of the ambient dimension. We empirically show that our algorithm outperforms the state-of-the-art methods on (multi-class) classification tasks.

Tables

Table 1. Table 1 : Time complexity comparisons on the setting of Corollary 4.1.1 . For clear comparison, we refer the per iteration cost as the time complexity of outer iterations.

Algorithm	Per Iteration Cost	Iteration Complexity
Frank Wolfe	$𝒪\left(nd\right)$	$𝒪\left(\frac{1}{ϵ}\right)$
Accelerated PGD [29]	$𝒪\left(nd\right)$	$𝒪\left(\sqrt{\kappa }\log \frac{1}{ϵ}\right)$
SVRG [15]	$𝒪\left(nd\right)$	$𝒪\left(\left(1+\kappa /n\right)\log \frac{1}{ϵ}\right)$
SCGS [23]	$𝒪\left({\kappa }^2\frac{\#{\text{iter}}^3}{{ϵ}^2}d\right)$	$𝒪\left(\frac{1}{ϵ}\right)$
STORC [13]	$𝒪\left({\kappa }^{2}d+nd\right)$	$𝒪\left(\log \frac{1}{ϵ}\right)$
Primal Dual FW (ours)	$𝒪\left(ns\right)$	$𝒪\left(\left(1+\kappa /n\right)\log \frac{1}{ϵ}\right)$

Table 2. Table 2 : Summary of the properties of the datasets.

Dataset Name	# Features	# Samples	# Non-Zero	Solution Sparsity (Ratio)
duke breast-cancer [4]	7,129	44	313,676	423 (5.9%)
rcv1 [4]	47,236	20,242	1,498,952	1,169 (2.5%)
news20.binary [4]	1,355,191	19,996	9,097,916	1,365 (0.1%)
MNIST.RB 0 VS 9 [4, 32]	894,499	11,872	1,187,200	8,450 (0.9%)
ijcnn.RB [4, 32]	58,699	49,990	14,997,000	715 (1.2%)
cob-rna.RB [4, 32]	81,398	59,535	5,953,500	958 (1.2%)

Equations

$$\underset{\mathbf{x}\in C}{\min }:\sum _i{f}_i({\mathbf{a}}_i^{\top }\mathbf{x})+g(\mathbf{x}),$$

$$\underset{\mathbf{x}\in {\mathbb{R}}^d,\Vert \mathbf{x}{\Vert }_1\le \tau }{\mathrm{argmin}}f(\mathbf{x}),$$

$${\mathbf{x}}^{(t)}\leftarrow (1-\eta ){\mathbf{x}}^{(t-1)}+\eta \tilde{\mathbf{x}},\text{ where }\tilde{\mathbf{x}}\leftarrow \underset{\Vert \mathbf{x}{\Vert }_1\le \tau }{\mathrm{argmin}}⟨\nabla f({\mathbf{x}}^{(t-1)}),\mathbf{x}⟩.$$

$${\mathbf{x}}^{(t)}\leftarrow (1-\eta ){\mathbf{x}}^{(t-1)}+\eta \tilde{\mathbf{x}},\text{ where }\tilde{\mathbf{x}}\leftarrow \underset{\Vert \mathbf{x}{\Vert }_1\le \tau ,\Vert \mathbf{x}{\Vert }_0\le s}{\mathrm{argmin}}\{⟨\nabla f({\mathbf{x}}^{(t-1)}),\mathbf{x}⟩+\frac{L}{2}\eta \Vert {\mathbf{x}}^{(t-1)}-\mathbf{x}{\Vert }_2^2\}.$$

$$h_t$$

$$\le h_{t-1}+\eta ⟨\nabla f({\mathbf{x}}^{(t-1)}),\tilde{\mathbf{x}}-{\mathbf{x}}^{(t-1)}⟩+\frac{L}{2}{\eta }^2\Vert \tilde{\mathbf{x}}-{\mathbf{x}}^{(t-1)}{\Vert }^2$$

$$\le h_{t-1}+\eta ⟨\nabla f({\mathbf{x}}^{(t-1)}),{\mathbf{x}}^{\ast }-{\mathbf{x}}^{(t-1)}⟩+\frac{L}{2}{\eta }^2\Vert {\mathbf{x}}^{\ast }-{\mathbf{x}}^{(t-1)}{\Vert }^2$$

$$\le (1-\eta +\frac{L}{\mu }{\eta }^2)h_{t-1}$$

$$\underset{\mathbf{x}\in C}{\min }\left\{P(\mathbf{x})\equiv \frac{1}{n}\sum _{i=1}^n{f}_i({\mathbf{a}}_i^{\top }\mathbf{x})+g(\mathbf{x})\right\},$$

$$\text{eqref}eqn:prima{l}_{c}omposite$$

$$\text{eqref}eqn:prima{l}_{c}omposite\iff \underset{\mathbf{y}}{\max }\left\{D(\mathbf{y}):=\underset{\mathbf{x}\in C}{\min }\left\{g(\mathbf{x})+\frac{1}{n}⟨\mathbf{y},A\mathbf{x}⟩\right\}-\frac{1}{n}\sum _{i=1}^n{f}_i^{\ast }(y_i)\right\}.$$

$$\tilde{\mathbf{x}}\leftarrow \underset{\Vert \mathbf{x}{\Vert }_1\le \lambda ,\Vert \mathbf{x}{\Vert }_0\le s}{\mathrm{argmin}}\{⟨\frac{1}{n}{\mathbf{z}}^{(t-1)}+\nabla g({\mathbf{x}}^{(t-1)}),\mathbf{x}⟩+\frac{L}{2}\eta \Vert \mathbf{x}-{\mathbf{x}}^{(t-1)}{\Vert }^2\}$$

$${\mathbf{x}}^{(t)}\leftarrow (1-\eta ){\mathbf{x}}^{(t-1)}+\eta \tilde{\mathbf{x}}$$

$${\mathbf{w}}^{(t)}\leftarrow (1-\eta ){\mathbf{w}}^{(t-1)}+\eta A\Delta \mathbf{x}$$

$$\tilde{\mathbf{y}}=\underset{{\mathbf{y}}^{\prime }}{\mathrm{argmax}}\left\{\frac{1}{n}⟨{\mathbf{w}}^{(t)},{\mathbf{y}}^{\prime }⟩-f^{\ast }({\mathbf{y}}^{\prime })-\frac{1}{2\delta }\Vert {\mathbf{y}}^{\prime }-{\mathbf{y}}^{(t-1)}{\Vert }^2\right\}.$$

$$|{\tilde{y}}_i-y_i^{(t-1)}|,i\in [n].$$

$$y_{i}^{(t)}\leftarrow\left\{\begin{array}[]{lc}\tilde{y}_{i}&\text{ if }i\in I^{(t)}\\ y_{i}^{(t-1)}&\text{ otherwise.}\end{array}\right.$$

$${\mathbf{z}}^{(t)}\leftarrow {\mathbf{z}}^{(t-1)}+A_{:,I^{(t)}}^{\top }({\mathbf{y}}^{(t)}-{\mathbf{y}}^{(t-1)})$$

$$\mathcal{L}({\mathbf{x}}^{(t+1)},{\mathbf{y}}^{(t)})-\mathcal{L}({\mathbf{x}}^{(t)},{\mathbf{y}}^{(t)})\le -\frac{\eta }{2}{\Delta }_p^{(t)}.$$

$$-\Vert {\mathbf{y}}^{(t)}-{\mathbf{y}}^{(t-1)}{\Vert }^2\le -\frac{k\delta }{n\beta }{\Delta }_d^{(t)}+\frac{k\delta }{n^2}R\Vert {\bar{\mathbf{x}}}^{(t)}-{\mathbf{x}}^{(t)}{\Vert }_2^2$$

$${\Delta }^{(t)}-{\Delta }^{(t-1)}\le \stackrel{\text{primal progress}}{\overbrace{\mathcal{L}({\mathbf{x}}^{(t+1)},{\mathbf{y}}^{(t)})-\mathcal{L}({\mathbf{x}}^{(t)},{\mathbf{y}}^{(t)})}}-\frac{1}{4\delta }\stackrel{\text{dual progress}}{\overbrace{\Vert {\mathbf{y}}^{(t)}-{\mathbf{y}}^{(t-1)}{\Vert }^2}}+\frac{3\delta Rk}{2n^2}\stackrel{\text{primal hindrance}}{\overbrace{\Vert {\bar{\mathbf{x}}}^{(t)}-{\mathbf{x}}^{(t)}{\Vert }^2}}.$$

$$X^{(t)}\leftarrow (1-\eta )X^{(t-1)}+\eta \tilde{X},\text{ where }\tilde{X}\leftarrow (\frac{1}{2},\frac{ϵ}{8})\text{-approximation of Eqn. eqrefeqn:deltaX}.$$

$$W^{(t)}\leftarrow (1-\eta )W^{(t-1)}+\eta A\tilde{X}$$

$${\tilde{Y}}^{(t)}\leftarrow \underset{Y}{\mathrm{argmax}}\left\{⟨W,Y⟩-f^{\ast }(Y)-\frac{1}{2\delta }\Vert Y-Y^{(t-1)}{\Vert }^2\right\}$$

$${\Vert {\tilde{Y}}_{i,:}-Y_{i,:}^{(t-1)}\Vert }_2,i\in [n].$$

$$Y_{i,:}^{(t)}\leftarrow\left\{\begin{array}[]{lc}\tilde{Y}_{i,:}&\text{ if }i\in I^{(t)}\\ Y_{i,:}^{(t-1)}&\text{ otherwise.}\end{array}\right.$$

$$Z^{(t)}\leftarrow Z^{(t-1)}+A^{\top }(Y^{(t)}-Y^{(t-1)})$$

$$\underset{\Vert X{\Vert }_{\ast }\le \lambda ,X\in {\mathbb{R}}^{d\times c}}{\min }\left\{\frac{1}{n}\sum _{i=1}^n{f}_i({\mathbf{a}}_i^{\top }X)+g(X)\right\}.$$

$$\underset{\Vert X{\Vert }_{\ast }\le \lambda ,X\in {\mathbb{R}}^{d\times c}}{\min }\underset{Y\in {\mathbb{R}}^{n\times c}}{\max }\left\{\mathcal{L}(X,Y)\equiv g(X)+\frac{1}{n}⟨AX,Y⟩-\frac{1}{n}\sum _{i=1}^n{f}_i^{\ast }({\mathbf{y}}_i)\right\}.$$

$$\tilde{X}\leftarrow \underset{\Vert X{\Vert }_{\ast }\le \lambda ,\text{rank}(X)\le s}{\mathrm{argmin}}\left\{⟨\frac{1}{n}Z+\nabla g(X^{(t-1)}),X⟩+\frac{L}{2}\eta \Vert X-X^{(t-1)}{\Vert }^2\right\},Z\equiv A^{\top }{Y}^{(t-1)}.$$

$$\pm ⟨[{\tilde{\mathbf{a}}}_i^{\top }|1],[\nabla f_{|A\mathbf{x}}^{\top }|{\nabla }_{i}g(\mathbf{x})]⟩,\text{i.e.}\pm \left({\tilde{\mathbf{a}}}_i^{\top }\nabla f_{|A\mathbf{x}}+{\nabla }_{i}g(\mathbf{x})\right),$$

$$\underset{\mathbf{x}\in C}{\min }P(\mathbf{x})=\frac{1}{n}\sum _{i=1}^n{f}_i({\mathbf{a}}_i^{\top }\mathbf{x})+g(\mathbf{x})$$

Full text

1 Introduction

We consider optimization problems of the form:

[TABLE]

directly motivated by regularized and constrained Empirical Risk Minimization (ERM). Particularly, we are interested in problems whose solution has special “simple” structure like low-rank or sparsity. The sparsity constraint applies to large-scale multiclass/multi-label classification, low-degree polynomial data mapping [5], random feature kernel machines [32], and Elastic Net [39]. Motivated by recent applications in low-rank multi-class SVM, phase retrieval, matrix completion, affine rank minimization and other problems (e.g., [9, 31, 2, 3]), we also consider settings where the constraint $\bm{x}\in C$ (e.g., trace norm ball) while convex, may be difficult to project onto. A wish-list for this class of problems would include an algorithm that (1) exploits the function finite-sum form and the simple structure of the solution, (2) achieves linear convergence for smooth and strongly convex problems, (3) does not pay a heavy price for the projection step.

We propose a Frank-Wolfe (FW) type method that attains these three goals. This does not come without challenges: Although it is currently well-appreciated that FW type algorithms avoid the cost of projection [14, 1], the benefits are limited to constraints that are hard to project onto, like the trace norm ball. For problems like phase retrieval and ERM for multi-label multi-class classification, the gradient computation requires large matrix multiplications. This dominates the per-iteration cost, and the existing FW type methods do not asymptotically reduce time complexity per iteration, even without paying the expensive projection step. Meanwhile, for simpler constraints like the $\ell_{1}$ norm ball or the simplex, it is unclear if FW can offer any benefits compared to other methods. Moreover, as is generally known, FW suffers from sub-linear convergence rate even for well-conditioned problems that enjoy strong convexity and smoothness.

In this paper we tackle the challenges by exploiting the special structure induced by the constraints and FW steps. We propose a variant of FW that we call Primal-Dual Block Frank Wolfe. The main advantage is that the computational complexity depends only on the sparsity of the solution, rather than the ambient dimension, i.e. it is dimension free. This is achieved by conducting partial updates in each iteration, i.e., sparse updates for $\ell_{1}$ and low-rank updates for the trace norm ball. While the benefits of partial updates seem unclear for the original problem, they significantly benefit a primal-dual reformulation. This reduces the per iteration cost to roughly a ratio of $\frac{s}{d}$ compared to naive Frank-Wolfe, where $s$ is the sparsity (or rank) of the optimal solution, and $d$ is the feature dimension. Meanwhile, the per iteration progress of our proposal is comparable to a full gradient descent step, thus retaining linear convergence rate.

For strongly convex and smooth $f$ and $g$ we show that our algorithm achieves linear convergence with per-iteration cost $sn$ over $\ell_{1}$-norm ball, where $s$ upper bounds the sparsity of the primal optimal. Specifically, for sparse ERM with smooth hinge loss or quadratic loss with $\ell_{2}$ regularizer, our algorithm yields an overall ${\mathcal{O}}(s(n+\kappa)\log\frac{1}{\epsilon})$ time complexity to reach $\epsilon$ duality gap, where $\kappa$ is the condition number (smoothness divided by strong convexity). Our theory has minimal requirements on the data matrix $A$.

Experimentally we observe our method yields significantly better performance compared to prior work, especially when the data dimension is large and the solution is sparse.

2 Related Work

We review relevant algorithms that improve the overall performance of Frank-Wolfe type methods. Such improvements are roughly obtained for two reasons: the enhancement on convergence speed and the reduction on iteration cost. Very few prior works benefit in both.

Nesterov’s acceleration has proven effective as in Stochastic Condition Gradient Sliding (SCGS) [23] and other variants [36, 26, 10]. Restarting techniques dynamically adapt to the function geometric properties and fills in the gap between sublinear and linear convergence for FW method [18]. Some variance reduced algorithms obtain linear convergence as in [13], however, the number of inner loops grows significantly and hence the method is not computationally efficient.

Linear convergence has been obtained specifically for polytope constraints like [27, 20], as well as the work proposed in [21, 11] that use the Away-step Frank Wolfe and Pair-wise Frank Wolfe, and their stochastic variants. One recent work [1] focuses on trace norm constraints and proposes a FW-type algorithm that yields similar progress as projected gradient descent per iteration but is almost projection free. However, in many applications where gradient computation dominates the iteration complexity, the reduction on projection step doesn’t necessarily produce asymptotically better iteration costs.

The sparse update introduced by FW steps was also appreciated by [22], where they conducted dual updates with a focus on SVM with polytope constraint. Their algorithm yields low iteration costs but still suffer from sub-linear convergence.

On the other hand, the recently popularized primal-dual formulation $\min_{{\bm{x}}}\max_{{\bm{y}}}\{g({\bm{x}})+{\bm{y}}^{\top}A{\bm{x}}-f({\bm{y}})\}$ has proven useful for different machine learning tasks like reinforcement learning, ERM, and robust optimization [8]. Especially for the ERM related problems, the primal-dual formulation still inherits the finite-sum structure from the primal form, and could be used to reduce variance [38, 35] or reduces communication complexity in the distributed setting [37]. However, there has been a lack of studies for constrained problems for this form. Meanwhile, as we emphasize in this paper, since the most time-consuming part of this formulation comes from the bilinear term ${\bm{y}}^{\top}A{\bm{x}}$, we are able to exploit simple structure from the optimal solutions.

3 Setup

3.1 Notation

We briefly introduce the notation used throughout the paper. We use bold lower case letter to denote vectors, capital letter to represent matrices. $\|\cdot\|$ is $\ell_{2}$ norm for vectors and Frobenius norm for matrices unless specified otherwise. $\|\cdot\|_{*}$ indicates the trace norm for a matrix.

We say a function $f$ is $\alpha$ strongly convex if $f({\bm{y}})\geq f({\bm{x}})+\langle{\bm{g}},{\bm{y}}-{\bm{x}}\rangle+\frac{\alpha}{2}\|{\bm{y}}-{\bm{x}}\|^{2},$ where ${\bm{g}}\in\partial f({\bm{x}})$ is any sub-gradient of $f$. Similarly, $f$ is $\beta$-smooth when $f({\bm{y}})\leq f({\bm{x}})+\langle{\bm{g}},{\bm{y}}-{\bm{x}}\rangle+\frac{\beta}{2}\|{\bm{y}}-{\bm{x}}\|^{2}.$ We use $f^{*}$ to denote the convex conjugate of $f$, i.e., $f^{*}({\bm{y}})\overset{\text{def}}{=}\max_{{\bm{x}}}\langle{\bm{x}},{\bm{y}}\rangle-f({\bm{x}})$. Some more parameters are problem-specific and are defined when needed.

3.2 A Theoretical Vignette

To elaborate the techniques we use to obtain the linear convergence for our Frank-Wolfe type algorithm, we consider the $\ell_{1}$ norm constrained problem as an illustrating example:

[TABLE]

where $f$ is $L$-smooth and $\mu$-strongly convex. If we invoke the Frank Wolfe algorithm, we compute

[TABLE]

Even when the function $f$ is smooth and strongly convex, (2) converges sublinearly. As inspired by [1], if we assume the optimal solution is $s$-sparse, we can enforce a sparse update while maintaining linear convergence by a mild modification on (2):

[TABLE]

We also call this new practice block Frank-Wolfe as in [1]. The proof of convergence can be completed within three lines. Let $h_{t}=f({\bm{x}}^{(t)})-f^{*}$.

[TABLE]

Therefore, when $\eta=\frac{\mu}{2L}$, $h_{t+1}\leq(1-\frac{\mu}{4L})^{t}h_{1}$ and the iteration complexity is ${\mathcal{O}}(\frac{L}{\mu}\log(1/\epsilon))$ to achieve $\epsilon$ error.

Although we achieve linear convergence, the advantage of overall complexity against classical methods (e.g. Projected Gradient Descend (PGD)) is not shown yet. Luckily, with the update $\tilde{{\bm{x}}}$ being sparse, it is possible to improve the iteration complexity, while maintaining the linear convergence rate. In order to differentiate, we name the sparse update nature of (3) as partial update.

Consider quadratic functions as an example:$f({\bm{x}})=\frac{1}{2}{\bm{x}}^{\top}A{\bm{x}}$, whose gradient is $A{\bm{x}}$ for symmetric $A$. As $\tilde{{\bm{x}}}$ is sparse, One can maintain the value of the gradient efficiently: $A{\bm{x}}^{(t)}\equiv(1-\eta)A{\bm{x}}^{(t-1)}+\eta A_{I,:}\tilde{{\bm{x}}}$, where $I$ is the support set of $\tilde{{\bm{x}}}$. We therefore reduce the complexity of one iteration to ${\mathcal{O}}(sd)$, compared to ${\mathcal{O}}(d^{2})$ with PGD. Similar benefits hold when we replace ${\bm{x}}$ by a matrix $X$ and conduct a low-rank update on $X$.

The benefit of partial update is not limited to quadratic functions. Next we show that for a class of composite function, we are able to take the full advantage of the partial update, by taking a primal-dual re-formulation.

4 Methodology

4.1 Primal-Dual Formulation

Note that the problem we are tackling is as follows

[TABLE]

We first focus on the setting where ${\bm{x}}\in\mathbb{R}^{d}$ is a vector and $C$ is the $\ell_{1}$-norm ball. This form covers general classification or regression tasks with $f_{i}$ being some loss function and $g$ being a regularizer. Extension to matrix optimization over a trace norm ball is introduced in Section 4.4.

Even with the constraint, we could reform (5) as a primal-dual convex-concave saddle point problem:

[TABLE]

or its dual formulation:

[TABLE]

Notice (7) is not guaranteed to have an explicit form. Therefore some existing FW variants like [22] that optimizes over (7) may not apply. Instead, we directly solve the convex concave problem (6) and could therefore solve more general problems, including complicated constraint like trace norm.

Since the computational cost of the gradient $\nabla_{{\bm{x}}}{\mathcal{L}}$ and $\nabla_{{\bm{y}}}{\mathcal{L}}$ is dominated by computing $A^{\top}{\bm{y}}$ and $A{\bm{x}}$ respectively, sparse updates could reduce computational costs by a ratio of roughly ${\mathcal{O}}(d/s)$ for updating ${\bm{x}}$ and ${\bm{y}}$ while achieving good progress.

4.2 Primal-Dual Block Frank-Wolfe

With the primal-dual formulation, we are ready to introduce our algorithm. The idea is simple: since the primal variable ${\bm{x}}$ is constrained over $\ell_{1}$ norm ball, we conduct block Frank-Wolfe algorithm and achieve an $s$-sparse update. Meanwhile, for the dual variable ${\bm{y}}$ we conduct greedy coordinate ascent method to select and update $k$ coordinates ($k$ determined later). We selected coordinates that allow the largest step, which is usually referred as a Gauss-Southwell rule denoted by GS-r [30]. Our algorithm is formally presented in Algorithm 1. The parameters to set in the algorithm request setting assumptions on $f$ and $g$:

Assumption 4.1**.**

We assume the functions satisfy the following properties:

• •

Each loss function $f_{i}$ is convex and $\beta$-smooth, and is $\alpha$ strongly convex over some convex set (could be ${\mathbb{R}}$), and linear otherwise.

• •

$\max_{i}\|{\bm{a}}_{i}\|_{2}^{2}\leq R$. Therefore $\frac{1}{n}\sum_{i=1}^{n}f_{i}({\bm{a}}_{i}^{\top}{\bm{x}})$ is $\beta R$-smooth.

• •

$g$* is $\mu$-strongly convex and $L$-smooth.*

Suitable loss functions $f_{i}$ include smooth hinge loss [33] and quadratic loss function. Relevant applications covered are Support Vector Machine (SVM) with smooth hinge loss, elastic net [39], and linear regression problem with quadratic loss.

As ${\mathcal{L}}({\bm{x}},{\bm{y}})$ is $\mu$-strongly convex and $L$-smooth with respect to ${\bm{x}}$, we set the primal learning rate $\eta=\frac{\mu}{2L}$ according to Section 3.2. Meanwhile, the dual learning rate $\delta$ is set to balance its effect on the dual progress as well as the primal progress. We specify it in the theoretical analysis part.

The computational complexity for each iteration in Algorithm 1 is ${\mathcal{O}}(ns)$. Both primal and dual update could be viewed as roughly three steps: coordinate selection, variable update, and maintaining $A^{T}{\bm{y}}$ or $A{\bm{x}}$. The coordinate selection as Eqn. (8) for primal and the choice of $I^{(t)}$ for dual variable respectively take ${\mathcal{O}}(d)$ and ${\mathcal{O}}(n)$ on average if implemented with the quick selection algorithm. Meanwhile, the variable update costs ${\mathcal{O}}(s)$ and ${\mathcal{O}}(k)$. The dominating cost is maintaining $A{\bm{x}}$ as in Eqn. (9) that takes ${\mathcal{O}}(ns)$, and maintaining $\mathcal{A}^{\top}{\bm{y}}$ as in Eqn. (12) that takes ${\mathcal{O}}(dk)$. To balance the time budget for primal and dual step, we set $k=ns/d$ and achieve an overall complexity of ${\mathcal{O}}(ns)$ per iteration.

4.3 Theoretical Analysis

We derive convergence analysis under Assumption 4.1. The derivation consists of the analysis on the primal progress, the balance of the dual progress, and their overall effect.

Define the primal gap as $\Delta^{(t)}_{p}\overset{\text{def}}{=}{\mathcal{L}}({\bm{x}}^{(t+1)},{\bm{y}}^{(t)})-{\mathcal{L}}(\bar{{\bm{x}}}^{(t)},{\bm{y}}^{(t)})$, where $\bar{{\bm{x}}}^{(t)}$ is the primal optimal solution such that the dual $D({\bm{y}}^{(t)})={\mathcal{L}}(\bar{{\bm{x}}}^{(t)},{\bm{y}}^{(t)})$, and is sparse enforced by the $\ell_{1}$ constraint. The dual gap is $\Delta^{(t)}_{d}\overset{\text{def}}{=}D^{*}-D({\bm{y}}^{(t)})$. We analyze the convergence rate of duality gap $\Delta^{(t)}\equiv\max\{1,(\beta/\alpha-1)\}\Delta^{(t)}_{p}+\Delta^{(t)}_{d}$.

**Primal progress: ** Firstly, similar to the analysis in Section 3.2, we could derive that primal update introduces a sufficient descent as in Lemma B.2.

[TABLE]

**Dual progress: ** With the GS-r rule to carefully select and update the most important $k$ coordinates in the dual variable in (10), we are able to derive the following result on dual progress that diminishes dual gap as well as inducing error.

[TABLE]

Refer to Lemma B.5 for details.

Primal Dual progress: The overall progress evolves as:

[TABLE]

In this way, we are able to connect the progress on duality gap with constant fraction of its value, and achieve linear convergence:

Theorem 4.1**.**

Given a function $P({\bm{x}})=\sum_{i=1}^{n}f_{i}({\bm{a}}_{i}^{\top}{\bm{x}})+g({\bm{x}})$ that satisfies Assumption 4.1. Set $s$ to upper bound the sparsity of the primal optimal $\bar{{\bm{x}}}^{(t)}$, and learning rates $\eta=\frac{\mu}{2L},\delta=\frac{1}{k}(\frac{L}{\mu n\beta}+\frac{5\beta R}{2\alpha\mu n^{2}}(1+4\frac{L}{\mu}))^{-1}$. The duality gap $\Delta^{(t)}=\max\{1,\frac{\beta}{\alpha}-1\}\Delta^{(t)}_{p}+\Delta^{(t)}_{d}$ generated by Algorithm 1 takes ${\mathcal{O}}(\frac{L}{\mu}(1+\frac{\beta}{\alpha}\frac{R\beta}{n\mu})\log\frac{1}{\epsilon})$ iterations to achieve $\epsilon$ error. The overall complexity is ${\mathcal{O}}(ns\frac{L}{\mu}(1+\frac{\beta}{\alpha}\frac{R\beta}{n\mu})\log\frac{1}{\epsilon})$.

For our target applications like elastic net, or ERM with smooth hinge loss, we are able to connect the time complexity to the condition number of the primal form.

Corollary 4.1.1**.**

Given a smooth hinge loss or quadratic loss $f_{i}$ that is $\beta$-smooth, and $\ell_{2}$ regularizer $g=\frac{\mu}{2}\|{\bm{x}}\|^{2}$. Define the condition number $\kappa=\frac{\beta R}{\mu}$. Setting $s$ upper bounds the sparsity of the primal optimal $\bar{{\bm{x}}}^{(t)}$, and learning rates $\eta=\frac{1}{2},\delta=\frac{1}{k}(\frac{1}{n\beta}+\frac{25R}{2\mu n^{2}})^{-1}$, the duality gap $\Delta^{(t)}$ takes ${\mathcal{O}}((1+\frac{\kappa}{n})\log\frac{1}{\epsilon})$ iterations to achieve $\epsilon$ error. The overall complexity is ${\mathcal{O}}(s(n+\kappa)\log\frac{1}{\epsilon})$.

Our derivation of overall complexity implicitly requires $ns\geq d$ by setting $k=sd/n\geq 1$. This is true for our considered applications like SVM. Otherwise we choose $k=1$ and the complexity becomes ${\mathcal{O}}(\max\{d,ns\}(1+\frac{\kappa}{n})\log\frac{1}{\epsilon})$.

In Table 1, we briefly compare the time complexity of our algorithm with some benchmark algorithms: (1) Accelerated Projected Gradient Descent (PGD) (2) Frank-Wolfe algorithm (FW) (3) Stochastic Variance Reduced Gradient (SVRG) [15] (4) Stochastic Conditional Gradient Sliding (SCGS) [23] and (5) Stochastic Variance-Reduced Conditional Gradient Sliding (STORC) [13]. The comparison is not thorough but intends to select constrained optimization that improves the overall complexity from different perspective. Among them, accelerated PGD improves conditioning of the problem, while SCGS and STORC reduces the dependence on number of samples. In the experimental session we show that our proposal outperforms the listed algorithms under various conditions.

4.4 Extension to the Trace Norm Ball

We also extend our algorithm to matrix optimization over trace norm constraints:

[TABLE]

This formulation covers multi-label multi-class problems, matrix completion, affine rank minimization, and phase retrieval problems (see reference therein [3, 1]). Equivalently, we solve the following primal-dual problem:

[TABLE]

Here ${\bm{y}}_{i}$ is the $i$-th row of the dual matrix $Y$. For this problem, the partial update we enforced on the primal matrix is to keep the update matrix low rank:

[TABLE]

However, an exact solution to (18) requires computing the top $s$ left and right singular vectors of the matrix $X^{(t-1)}-\frac{1}{\eta L}(Z+\nabla g(X^{(t-1)})\in{\mathbb{R}}^{d\times c}$. Therefore we loosely compute an $(\frac{1}{2},\epsilon/2)$-approximation, where $\epsilon$ is the target accuracy, based on the following definition:

Definition 4.2** (Restated Definition 3.2 in [1]).**

Let $l_{t}(V)=\langle\nabla_{X}{\mathcal{L}}(X^{(t)},Y^{(t)}),V-X^{(t)}\rangle+\frac{L}{2}\eta\|V-X^{(t)}\|_{F}^{2}$ be the objective function in (18), and let $l^{*}_{t}=l_{t}(\bar{X}^{(t)})$. Given parameters $\gamma\geq 0$ and $\epsilon\geq 0$, a feasible solution $V$ to (18) is called $(\gamma,\epsilon)$-approximate if it satisfies $l(V)\leq(1-\gamma)l^{*}_{t}+\epsilon$.

The time dependence on the data size $n,c,d,s$ is $ncs+s^{2}(n+c)$ [1], and is again independent of $d$. Meanwhile, the procedures to keep track of $W^{(t)}\equiv AX^{(t)}$ requires complexity of $nds+ncs$, while updating $Y^{(t)}$ requires $dck$ operations. Therefore, by setting $k\leq ns(1/c+1/d)$, the iteration complexity’s dependence on the data size becomes ${\mathcal{O}}(n(d+c)s)$ operations, instead of ${\mathcal{O}}(ndc)$ for conducting a full projected gradient step. Recall that $s$ upper bounds the rank of $\bar{X}^{(t)}\leq\min\{d,c\}$.

The trace norm version mostly inherits the convergence guarantees for vector optimization. Refer to the Appendix for details.

Assumption 4.2**.**

We assume the following property on the primal form (17):

• •

$f_{i}$* is convex, and $\beta$-smooth. Its convex conjugate $f_{i}^{*}$ exists and satisfies $\frac{1}{\alpha}$-smooth on some convex set (could be ${\mathbb{R}}^{c}$) and infinity otherwise.*

• •

Data matrix $A$ satisfies $R=\max_{|I|\leq k,I\subset[n]}\sigma_{\max}^{2}(A_{I,:})$ ($\leq\|A\|_{2}^{2}$). Here $\sigma_{\max}(X)$ denotes the largest singular value of $X$.

• •

$g$* is $\mu$-strongly convex and $L$-smooth.*

The assumptions also cover smooth hinge loss as well as quadratic loss. With the similar assumptions, the convergence analysis for Algorithm 2 is almost the same as Algorithm 1. The only difference comes from the primal step where approximated update produces some error:

**Primal progress: ** With the primal update rule in Algorithm 2, it satisfies ${\mathcal{L}}(X^{(t+1)},Y^{(t)})-{\mathcal{L}}(X^{(t)},Y^{(t)})\leq-\frac{\mu}{8L}\Delta_{p}^{(t)}+\frac{\epsilon}{16}$. (See Lemma B.7.) With no much modification in the proof, we are able to derive similar convergence guarantees for the trace norm ball.

Theorem 4.3**.**

Given a function $\frac{1}{n}\sum_{i=1}^{n}f_{i}({\bm{a}}_{i}^{\top}X)+g(X)$ that satisfies Assumption 4.2. Setting $s\geq$rank$(\bar{X}^{(t)})$, and learning rate $\eta=\frac{\mu}{2L},\delta\leq\frac{1}{k}(\frac{L}{\mu n\beta}+\frac{5\beta R}{2\alpha\mu n^{2}}(1+8\frac{L}{\mu}))^{-1}$, the duality gap $\Delta^{(t)}$ generated by Algorithm 2 satisfies $\Delta^{(t)}\leq\frac{k\delta}{k\delta+8\beta n}\Delta^{(t-1)}+\frac{\epsilon}{16}.$ Therefore it takes ${\mathcal{O}}(\frac{L}{\alpha}(1+\frac{\beta}{\alpha}\frac{R\beta}{n\mu})\log\frac{1}{\epsilon})$ iterations to achieve $\epsilon$ error.

We also provide a brief analysis on the difficulty to extend our algorithm to polytope-type constraints in the Appendix B.9.

5 Experiments

We evaluate the primal-dual Frank-Wolfe algorithm by its performance on binary classification with smoothed hinge loss. We refer the readers to Appendix B.7 for details about smoothed hinge loss.

We compare the proposed algorithm against five benchmark algorithms: (1) Accelerated Projected Gradient Descent (Acc PG) (2) Frank-Wolfe algorithm (FW) (3) Stochastic Variance Reduced Gradient (SVRG) [15] (4) Stochastic Conditional Gradient Sliding (SCGS) [23] and (5) Stochastic Variance-Reduced Conditional Gradient Sliding (STORC) [13]. We presented the time complexity for each algorithm in Table 1. Three of them (FW, SCGS, STORC) are projection-free algorithms, and the other two (Acc PG, SVRG) are projection-based algorithms. Algorithms are implemented in C++, with the Eigen linear algebra library [12].

The six datasets used here are summarized in Table 2. All of them can be found in LIBSVM datasets [4]. We augment the features of MNIST, ijcnn, and cob-rna by random binning [32], which is a standard technique for kernel approximation. Data is normalized. We set the $\ell_{1}$ constraint to be $300$ and the $\ell_{2}$ regularize parameter to $10/n$ to achieve reasonable prediction accuracy. We refer the readers to the Appendix C.1 for results of other choice of parameters. These datasets have various scale of features, samples, and solution sparsity ratio.

The results are shown in Fig 1. To focus on the convergence property, we show the decrease of loss function instead of prediction accuracy. From Fig 1, our proposed algorithm consistently outperforms the benchmark algorithms. The winning margin is roughly proportional to the solution sparsity ratio, which is consistent with our theory.

5.1 Experiments with trace norm ball on synthetic data

For trace norm constraints, we also implemented our proposal Primal Dual Block Frank Wolfe to compare with some prior work, especially Block FW [1]. Since prior work were mostly implemented in Matlab to tackle trace norm projections, we therefore also use Matlab to show fair comparisons. We choose quadratic loss $f(AX)=\|AX-B\|_{F}^{2}$ and $g$ to be $\ell_{2}$ regularizer with $\mu=10/n$. The synthetic sensing matrix $A\in{\mathbb{R}}^{n\times d}$ is dense with $n=1000$ and $d=800$. Our observation $B$ is of dimension $1000\times 600$ and is generated by a ground truth matrix $X_{0}$ such that $B=AX_{0}$. Here $X_{0}\in{\mathbb{R}}^{800\times 600}$ is constructed with low rank structure. We vary its rank $s$ to be $10,20,$ and $100$. The comparisons with stochastic FW, blockFW [1], STORC [13], SCGS [23], and projected SVRG [15] are presented in Figure 2, which verifies that our proposal PDBFW consistently outperforms the baseline algorithms.

6 Conclusion

In this paper we consider a class of problems whose solutions enjoy some simple structure induced by the constraints. We propose a FW type algorithm to exploit the simple structure and conduct partial updates, reducing the time cost for each update remarkably while attaining linear convergence. For a class of ERM problems, our running time depends on the sparsity/rank of the optimal solutions rather than the ambient feature dimension. Our empirical studies verify the improved performance compared to various state-of-the-art algorithms.

Appendix A Efficient Coordinate Selections

The modified Block Frank-Wolfe step in Eqn. (3) achieves an $s$-sparse update of the iterates and could be computed efficiently when one knows which $s$ coordinates to update. However, in order to find the $s$ coordinates, one needs to compute the full gradient $\nabla f({\bm{x}})$ with naive implementation. This phenomenon reminds us of greedy coordinate descent.

Even with the known fact that coordinate descent converges faster with greedy selection than with random order[30], there have been hardness to propogate this idea because of expensive greedy selections since the arguments that GCD converges similarly with RCD in [28], except for special cases [25, 24, 6, 17]. This is also probability why the partial updates nature of FW steps is less exploited before.

We investigate some possible tricks to boost GCD method that could be possibly applied to FW methods. A recent paper [17], Karimireddy et al. make connections between the efficient choice of the greedy coordinates with the problem of Maximum Inner Product Search (MIPS) for a composite function $P({\bm{x}})=f(A{\bm{x}})+g({\bm{x}})$, where $A\in{\mathbb{R}}^{n\times d}$. We rephrase the connection for the Frank-Wolfe algorithm. Since the computation of gradient is essentially $A^{\top}\nabla f_{|A{\bm{x}}}+\nabla g({\bm{x}})$, to find its largest magnitude is to search maximum inner products among:

[TABLE]

where $\tilde{{\bm{a}}}_{i}\in{\mathbb{R}}^{n}$ is the $i$-th column of data matrix $A$, and $\nabla f_{|A{\bm{x}}}$ is the gradient of $f$ at $A{\bm{x}}$. In this way, we are able to select the greedy coordinates by conducting MIPS for a fixed ${\mathbb{R}}^{2d\times(n+1)}$ matrix $[A^{\top}|I|-A^{\top}|-I]^{\top}$ and each newly generated vector $[\nabla f_{|A{\bm{x}}}^{\top}|\nabla g_{i}({\bm{x}})]$. Therefore when $\nabla g_{i}$ is constant for linear function or $\pm\lambda$ for $g({\bm{x}})=\lambda\|{\bm{x}}\|_{1}$, we could find the largest magnitude of the gradient in sublinear time. Still, the problems it could conquer is very limited. It doesn’t even work for $\ell_{2}$ regularizer since the different coordinates in $\nabla_{i}g({\bm{x}})$ creates $d$ new vectors in each iteration and traditional MIPS could resolve it in time sublinear to $d$. Meanwhile, even with constant $\nabla_{i}g({\bm{x}})$, it still requires at least ${\mathcal{O}}((2d)^{c}\log(d))$ times of inner products of dimension $n+1$ for some constant $c$ [34].

However, we have shown that for general composite form $f(A{\bm{x}})+g({\bm{x}})$ with much more relaxed requirements on the regularizer $g$, we are able to select and update each coordinate with constant times of inner products on average while achieving linear convergence. Therefore the usage of these tricks applied on FW method (MIPS as well as the nearest neighbor search [6]) is completely dominated by our contribution and we omit them in the main text of this paper.

Appendix B Proofs

B.1 Derivation of Primal-Dual Formulation

[TABLE]

We use Von Neumann-Fan minimax theorem for the whole derivation when swapping each min-max formula [7]. For the last equality, there is a convex constraint in the minimization part. Although the original Von Neumann-Fan doesn’t have constraints, it naturally applies to the case when ${\bm{x}}$ (assuming function is convex to ${\bm{x}}$) is bounded in a convex set, since we could change $f({\bm{x}},{\bm{y}})$ to $f({\bm{x}},{\bm{y}})+I_{C}({\bm{x}})$, where $I_{C}({\bm{x}})=0$ if ${\bm{x}}\in C$ and $\infty$ otherwise. Then the property will be properly inherited.

B.2 Notation and simple facts

Recall primal, dual and Lagrangian forms:

[TABLE]

Similar to the definitions in [24], we introduce the primal gap defined as $\Delta^{(t)}_{p}\overset{\text{def}}{=}{\mathcal{L}}({\bm{x}}^{(t+1)},{\bm{y}}^{(t)})-D({\bm{y}}^{(t)})$, and dual gap $\Delta^{(t)}_{d}\overset{\text{def}}{=}D^{*}-D({\bm{y}}^{(t)})$. Recall the assumptions:

• •

$f_{i}$ is convex and $\beta$-smooth, and is $\alpha$ strongly convex over some convex set, and linear otherwise.

• •

$R=\max_{i}\|{\bm{a}}_{i}\|_{2}^{2},\forall i\in[n]$.

• •

$g$ is $\mu$-strongly convex and $L$-smooth.

To begin with, it is easy to verify that $f_{i}^{*}$ is $1/\beta$-strongly convex and is $1/\alpha$-smooth on a convex set and infinity otherwise (See Claim B.6). For simplicity we first assume $\alpha\geq\frac{1}{2}\beta$ and then generalize the result.

Claim B.1**.**

• •

Since $D({\bm{y}})=\min_{{\bm{x}}\in C}\{g({\bm{x}})+\frac{1}{n}{\bm{y}}^{\top}A{\bm{x}}\}-\frac{1}{n}\sum_{i=1}^{n}f^{*}_{i}({\bm{y}})$, $-D({\bm{y}})$ is at least $\frac{1}{\beta}$-strongly convex.

• •

Based on our update rule, $\exists{\bm{g}}\in\partial_{{\bm{y}}}\frac{1}{n}\sum_{i}f_{i}^{*}({\bm{y}}^{(t)})$, such that

[TABLE]

And our update rule ensures that $I^{(t)}$ consists of indices $i\in[n]$ that maximizes $|\frac{1}{n}{\bm{a}}_{i}^{\top}{\bm{x}}^{(t)}-g_{i}|$.

B.3 Primal Progress

Lemma B.2**.**

(Primal Progress)

[TABLE]

Or equivalently,

[TABLE]

Proof.

Simply replace $h_{t}$ as ${\mathcal{L}}({\bm{x}}^{(t)},{\bm{y}}^{(t)})-D({\bm{y}}^{(t)})$ and $h_{t+1}$ as ${\mathcal{L}}({\bm{x}}^{(t+1)},{\bm{y}}^{(t)})-D({\bm{y}}^{(t)})$ in Inequality (4). We could conclude that $h_{t+1}\leq(1-\eta+\eta^{2}\frac{L}{\mu})h_{t}$. Therefore when $\eta\leq\frac{\mu}{2L}$, $h_{t+1}\leq(1-\frac{\eta}{2})h_{t}$ and the first part of Lemma B.2 is true. Some simple rearrangement suffices the second part of the lemma. ∎

B.4 Primal Dual Progress

In order to get a clue on how to analyze the dual progress, we first look at how the primal and dual evolve through iterations.

For an index set $I$ and a vector ${\bm{y}}\in{\mathbb{R}}^{n}$, denote ${\bm{y}}_{I}=\sum_{i\in I}y_{i}{\bm{e}}_{i}\in{\mathbb{R}}^{k}$ as the subarray of ${\bm{y}}$ indexed by $I$, with $|I|=k$. Recall Algorithm 1 selects the coordinates to update in the dual variable as $I^{(t)}$.

Lemma B.3**.**

(Primal-Dual Progress).

[TABLE]

Proof.

Notice we have claimed that $-D({\bm{y}})$ is $\frac{1}{\beta}$-strongly convex and for all ${\bm{g}}\in\partial_{{\bm{y}}}\frac{1}{n}\sum_{i}^{n}f^{*}_{i}({\bm{y}}^{(t)})$,

[TABLE]

Meanwhile since $-{\mathcal{L}}({\bm{x}},{\bm{y}})$ is $\frac{1}{\alpha}$-smooth over its feasible set,

[TABLE]

Also, with the update rule of dual variables, we could make use of Eqn. (19) and re-write Eqn. (20) as:

[TABLE]

Together we get:

[TABLE]

∎

Therefore we will connect the progress induced by $-\|{\bm{y}}^{(t)}-{\bm{y}}^{(t-1)}\|$ and dual gap $\Delta_{d}^{(t)}$ next.

B.5 Dual progress

Claim B.4**.**

An $\alpha$-strongly convex function $f$ satisfies:

[TABLE]

This simply due to $f({\bm{x}})-f^{*}\leq\langle\nabla f({\bm{x}}),{\bm{x}}-\bar{{\bm{x}}}\rangle-\frac{\alpha}{2}\|{\bm{x}}-\bar{{\bm{x}}}\|_{2}^{2}\leq\frac{1}{2\alpha}\|\nabla f({\bm{x}})\|^{2}+\frac{\alpha}{2}\|{\bm{x}}-\bar{{\bm{x}}}\|^{2}-\frac{\alpha}{2}\|{\bm{x}}-\bar{{\bm{x}}}\|^{2}=\frac{1}{2\alpha}\|\nabla f({\bm{x}})\|^{2}$.

Since $-D$ is $\frac{1}{\beta}$-strongly convex, we get

[TABLE]

where $\bar{I}$ is a set of size $k$ that maximizes the values of $A_{i}^{\top}\bar{{\bm{x}}}^{(t)}-g_{i}$.

Lemma B.5** (Dual Progress).**

[TABLE]

Proof of Lemma B.5.

Define $\Delta=\frac{1}{n}A(\bar{{\bm{x}}}^{(t)}-{\bm{x}}^{(t)})$. Since

[TABLE]

With the relation between $\frac{1}{n}A_{I^{(t)}}^{\top}{\bm{x}}^{(t)}-{\bm{g}}_{I^{(t)}}$ and ${\bm{y}}^{(t)}-{\bm{y}}^{(t-1)}$ we finish the proof. ∎

B.6 Convergence on Duality Gap

Now we are able to merge the primal/dual progress to get the overall progress on the duality gap.

Proof of Theorem 4.1.

We simply blend Lemma B.2 and Lemma B.5 with the primal-dual progress (Lemma B.3):

[TABLE]

When setting $\frac{k\delta}{2n\beta}=(1-\frac{5R\delta k}{\mu n^{2}})\frac{\mu}{4L}-\frac{5R\delta k}{\mu n^{2}}$, we get that $\Delta^{(t)}\leq\frac{1}{1+a}\Delta^{(t-1)}$, where $1/a={\mathcal{O}}(\frac{L}{\mu}(1+\frac{R\beta}{n\mu}))$. Therefore it takes ${\mathcal{O}}(\frac{L}{\mu}(1+\frac{R\beta}{n\mu})\log\frac{1}{\epsilon})$ for $\Delta^{(t)}$ to reach $\epsilon$.

When $\beta>2\alpha$, we could redefine the primal-dual process as $\Delta^{(t)}:=(\frac{\beta}{\alpha}-1)\Delta^{(t)}_{d}+\Delta^{(t)}_{p}$ and rewrite some of the key steps, especially for the overall primal-dual progress.

[TABLE]

Similarly to the previous setting, we get the whole primal-dual progress is bounded as follows:

[TABLE]

Therefore, when we set a proper $k$ and $\delta$ such that $\frac{\beta}{4\alpha}\frac{k\delta}{n\beta}=(\frac{\beta}{\alpha}-1)\left((1-\frac{5\beta R\delta k}{2\alpha\mu n^{2}})\frac{\mu}{4L}-\frac{5\beta R\delta k}{2\alpha\mu n^{2}}\right)$, and since $\frac{\beta}{\alpha}-1\geq\frac{\beta}{2\alpha}$, we get $\delta=\frac{1}{k}(\frac{L}{\mu n\beta}+\frac{5\beta R}{2\alpha\mu n^{2}}(1+4\frac{L}{\mu}))^{-1}$. And we have $\Delta^{(t)}-\Delta^{(t-1)}\leq-1/a\Delta^{(t)}$, where $a={\mathcal{O}}(\frac{L}{\mu}(1+\frac{\beta}{\alpha}\frac{R\beta}{n\mu}))$. Therefore it takes $t={\mathcal{O}}(\frac{L}{\mu}(1+\frac{\beta}{\alpha}\frac{R\beta}{n\mu})\log\frac{1}{\epsilon})$ iterations for the duality gap $\Delta^{(t)}$ to reach $\epsilon$ error. ∎

B.7 Smooth Hinge Loss and Relevant Properties

Smooth hinge loss is defined as follows:

[TABLE]

Our loss function over a prediction $p$ associated with a label $\ell_{i}\in\{\pm 1\}$ will be $f_{i}(p)=h(p\ell_{i})$. The derivative of smooth hinge loss $h$ is:

[TABLE]

Its convex conjugate is:

[TABLE]

Notice since $f_{i}(p)=h(\ell_{i}p)$, $f_{i}^{*}(p)=h^{*}(p/\ell_{i})=h^{*}(p\ell_{i})$.

Claim B.6**.**

For a convex and $\beta$-smooth scalar function $f$, if it is $\alpha$ strongly convex over some convex set, and linear otherwise, then its conjugate function $f^{*}$ is $1/\beta$-strongly convex, and it is a $1/\alpha$-smooth function plus an indicator function over some interval $[a,b]$.

Proof.

To begin with, since $f^{\prime\prime}(x)\leq\beta,\forall x$, meaning $f$ is $\beta$-smooth, then with duality we have $f^{*}$ is $1/\beta$ strongly convex [16]. Secondly, since $f$ is $\alpha$ strongly convex over a convex set, meaning an interval for ${\mathbb{R}}$, therefore $f$ could only be linear on $(-\infty,a]$ or $[b,\infty)$, and is $\alpha$-strongly convex over the set $[a,b]$ (Here for simplicity $a<b$ could be $\pm\infty$). We denote $f^{\prime}(-\infty):=\lim_{x\rightarrow-\infty}f^{\prime}(x)$ and $f^{\prime}(-\infty)$ likewise. It’s easy to notice that $f^{\prime}(-\infty)\leq f^{\prime}(a)<f^{\prime}(b)\leq f^{\prime}(\infty)$ since $f$ is convex overall and strongly convex over $[a,b]$. Therefore $f(y)>f(a)+f^{\prime}(a)(y-a)$ when $y>a$ and $f(y)=f(a)+f^{\prime}(a)(y-a)$ when $y\leq a$.

Now since $f^{*}(x^{*})\equiv\max_{x}\{x^{*}x-f(x)\}$, it’s easy to verify that when $x^{*}<f^{\prime}(a)$, $x^{*}x-f(x)=x^{*}x-f(a)-f^{\prime}(a)(x-a)=-(f^{\prime}(a)-x^{*})x-f(a)+f^{\prime}(a)a\rightarrow\infty$ when $x\rightarrow-\infty$. Similarly, when $x^{*}>f^{\prime}(b)$, $f^{*}(x^{*})=\infty$. On the other hand, when $x^{*}\in[f^{\prime}(a),f^{\prime}(b)]$, $f^{*}(x^{*})=\max_{x}\{x^{*}x-f(x)\}=\max_{x\in[a,b]}\{x^{*}x-f(x)\}$. This is because $x^{*}a-f(a)\geq x^{*}y-f(y)=x^{*}y-f(y)-f^{\prime}(a)(y-a),\forall y\leq a$, and similarly $x^{*}b-f(b)\geq x^{*}y-f(y)\forall y>b$. Therefore $f^{*}$ is $1/\alpha$ smooth over the interval $[f^{\prime}(a),f^{\prime}(b)]$, where $-\infty\leq f^{\prime}(a)<f^{\prime}(b)\leq\infty$.

∎

B.8 Convergence of Optimization over Trace Norm Ball

The convergence analysis for trace norm ball is mostly similar to the case of $\ell_{1}$ ball. The most difference lies on the primal part, where our approximated update incur linear progress as well as some error.

Lemma B.7** (Primal Progress for Algorithm 2).**

Suppose rank $\bar{X}^{(t)}\leq s$ and $\epsilon>0$. If each $\tilde{X}$ computed in our algorithm is a $(\frac{1}{2},\frac{\epsilon}{8})$-approximate solution to (18), then for every $t$, it satisfies ${\mathcal{L}}(X^{(t+1)},Y^{(t)})-{\mathcal{L}}(X^{(t)},Y^{(t)})\leq-\frac{\mu}{8L}\Delta_{p}^{(t)}+\frac{\epsilon}{16}$.

Proof.

Refer to the proof in [1] we have:

[TABLE]

Now move the first term on the RHS to the left and rearrange we get:

[TABLE]

Therefore we get:

[TABLE]

∎

Now back to the convergence guarantees on the trace norm ball.

Proof of Theorem 4.3.

We again define $\Delta=\frac{1}{n}A(\bar{X}^{(t)}-X^{(t)})$. $G=\nabla_{Y}{\mathcal{L}}(X^{(t)},Y^{(t)})$ such that $Y^{(t)}_{I^{(t)},:}-Y^{(t-1)}_{I^{(t)},:}=\delta(\frac{1}{n}\langle A_{I^{(t)},:}X^{(t)}\rangle-G_{I^{(t)},:})$. Again we get $\|\Delta\|_{F}^{2}\leq\frac{R}{n^{2}}\|\bar{X}^{(t)}-X^{(t)}\|_{F}^{2}.$

[TABLE]

Other parts are exactly the same and we get:

[TABLE]

Therefore when $\delta\leq\frac{1}{k}(\frac{L}{\mu n\beta}+\frac{5\beta R}{2\alpha\mu n^{2}}(1+8\frac{L}{\mu}))^{-1}$, it satisfies $\Delta^{(t)}-\Delta^{(t-1)}\leq-\frac{k\delta}{2\beta n}\Delta^{(t)}+\frac{\epsilon}{16}.$ Therefore denote $a=\frac{2\beta n}{k\delta}$, we get $\Delta^{(t)}\leq\frac{a}{a+1}(\Delta^{(t-1)}+\frac{\epsilon}{16})$. Therefore we get $\Delta^{(t)}\leq(\frac{a}{a+1})^{t}\Delta^{(0)}+\frac{\epsilon}{16}\sum_{i=1}^{t}(\frac{a}{a+1})^{i}\leq(\frac{c}{c+1})^{t}\Delta^{(0)}+\epsilon/16$. Since $(\frac{a}{a+1})^{t}\leq e^{-t/a}$, it takes around $a={\mathcal{O}}(\frac{L}{\mu}(1+\frac{\beta}{\alpha}\frac{R\beta}{n\mu})\log\frac{1}{\epsilon})$ iterations for the duality gap to get $\epsilon$-error. ∎

B.9 Difficulty on Extension to Polytope Constraints

Another important type of constraint we have not explored in this paper is the polytope constraint. Specifically,

[TABLE]

where ${\mathcal{A}}\subset{\mathbb{R}}^{d},|{\mathcal{A}}|=m$ is a finite set of vectors that is usually referred as atoms. It is worth noticing that this linear minimization oracle (LMO) for FW step naturally chooses a single vector in ${\mathcal{A}}$ that minimizes the inner product with ${\bm{x}}$. Again, this FW step creates some "partial update" that could be appreciated in many machine learning applications. Specifically, if our computation of gradient is again dominated by a matrix-vector (data matrix versus variable ${\bm{x}}$) inner product, we could possibly pre-compute each value of ${\bm{v}}_{i}:=A{\bm{x}}_{i},{\bm{x}}_{i}\in{\mathcal{A}}$, and simply use ${\bm{v}}_{i}$ to update the gradient information when ${\bm{x}}_{i}$ is the greedy direction provided by LMO.

When connecting to our sparse update case, we are now looking for a $k$-sparse update, $k\ll m=|{\mathcal{A}}|$, with the basis of ${\mathcal{A}}$, i.e., $\tilde{{\bm{x}}}=\sum_{i=1}^{k}\lambda_{i}{\bm{x}}_{n_{i}},{\bm{x}}_{n_{i}}\in{\mathcal{A}}$. In this way, when we update ${\bm{x}}^{+}\leftarrow(1-\eta){\bm{x}}+\eta\tilde{{\bm{x}}}$, we will only need to compute $\sum_{i=1}^{k}{\bm{v}}_{n_{i}}$ which is ${\mathcal{O}}(kd)$ time complexity.

However, to enforce such update that is "sparse" on ${\mathcal{A}}$ is much harder. To migrate our algorithms with $\ell_{1}$ ball or trace norm ball, we will essentially be solving the following problem:

[TABLE]

where $\Delta^{m}$ is the $m$ dimensional simplex, and ${\bm{g}}$ is the current gradient vector.

Unlike the original sparse recovery problem that could be relaxed with an $\ell_{1}$ constraint to softly encourage sparsity, it’s generally much harder to find the $k$ sparse $\Lambda$ in this case. Actually, it is as hard as the lattice problem [19] and is NP hard in general.

Therefore we are not able to achieve linear convergence with cheap update with polytope-type constraints. Nonetheless, the naive FW with primal dual formulation should still be computational efficient in terms of per iteration cost, where a concentration on SVM on its dual form has been explored by [22].

Appendix C More Results on Empirical Studies

C.1 More experiments with $\ell_{1}$ norm

To investigate more on how our algorithms perform with different choices of parameters, we conducted more empirical studies with different settings of condition numbers. Specifically, we vary the parameter $\mu$ that controls the strong convexity of the primal function. Experiments are shown in Figure 3.