Primal-Dual Block Frank-Wolfe

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.


Introduction
We consider optimization problems of the form: 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 x ∈ 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 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 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 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 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 2 regularizer, our algorithm yields an overall O(s(n + κ) log 1 ) time complexity to reach duality gap, where κ 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.

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 x max y {g(x) + y Ax − f (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 y Ax, we are able to exploit simple structure from the optimal solutions.

Notation
We briefly introduce the notation used throughout the paper. We use bold lower case letter to denote vectors, capital letter to represent matrices. · is 2 norm for vectors and Frobenius norm for matrices unless specified otherwise.
· * indicates the trace norm for a matrix.
We say a function f is α strongly . We use f * to denote the convex conjugate of f , i.e., f * (y) def = max x x, y − f (x). Some more parameters are problem-specific and are defined when needed.

A Theoretical Vignette
To elaborate the techniques we use to obtain the linear convergence for our Frank-Wolfe type algorithm, we consider the 1 norm constrained problem as an illustrating example: where f is L-smooth and µ-strongly convex. If we invoke the Frank Wolfe algorithm, we compute 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): We also call this new practice block Frank-Wolfe as in [1]. The proof of convergence can be completed within three lines.
(by convexity and µ-strong convexity of f ) Therefore, when η = µ 2L , h t+1 ≤ (1 − µ 4L ) t h 1 and the iteration complexity is O( L µ log(1/ )) to achieve 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 updatex 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 (x) = 1 2 x Ax, whose gradient is Ax for symmetric A. Asx is sparse, One can maintain the value of the gradient efficiently: Ax (t) ≡ (1 − η)Ax (t−1) + ηA I,:x , where I is the support set ofx. We therefore reduce the complexity of one iteration to O(sd), compared to O(d 2 ) with PGD. Similar benefits hold when we replace 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.

Primal-Dual Formulation
Note that the problem we are tackling is as follows We first focus on the setting where x ∈ R d is a vector and C is the 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: or its dual formulation: 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 ∇ x L and ∇ y L is dominated by computing A y and Ax respectively, sparse updates could reduce computational costs by a ratio of roughly O(d/s) for updating x and y while achieving good progress.

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 x is constrained over 1 norm ball, we conduct block Frank-Wolfe algorithm and achieve an s-sparse update. Meanwhile, for the dual variable 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 β-smooth, and is α strongly convex over some convex set (could be R), and linear otherwise.
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 L(x, y) is µ-strongly convex and L-smooth with respect to x, we set the primal learning rate η = µ 2L according to Section 3.2. Meanwhile, the dual learning rate δ 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 O(ns). Both primal and dual update could be viewed as roughly three steps: coordinate selection, variable update, and maintaining A T y or Ax. The coordinate selection as Eqn. (8)

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. Algorithm 1 Frank-Wolfe Block Primal-Dual Method for 1 Norm Ball 1: Input: Training data A ∈ R n×d , primal and dual step size η, δ > 0. 2: Initialize: Use Block Frank Wolfe to update the primal variable: x ← arg min

5:
Update w to maintain the value of Ax: 6: Consider the potential dual update:ỹ

7:
Choose greedily the dual coordinates to update: let I (t) be the top k coordinates that maximize Update the dual variable accordingly: 8: 9: end for 10: Output: Define the primal gap as ∆ is the primal optimal solution such that the dual D(y (t) ) = L(x (t) , y (t) ), and is sparse enforced by the 1 constraint. The dual gap is ∆ . We analyze the convergence rate of duality gap 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.
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. − Refer to Lemma B.5 for details. Primal Dual progress: The overall progress evolves as: primal hindrance In this way, we are able to connect the progress on duality gap with constant fraction of its value, and achieve linear convergence: that satisfies Assumption 4.1. Set s to upper bound the sparsity of the primal optimalx (t) , and learning rates η = µ 2L , δ = 1 k ( L µnβ + 5βR 2αµn 2 (1 + 4 L µ )) −1 . The duality gap Rβ nµ ) log 1 ) iterations to achieve error. The overall complexity is O(ns L µ (1 + β α Rβ nµ ) log 1 ). 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 β-smooth, and 2 regularizer g = µ 2 x 2 . Define the condition number κ = βR µ . Setting s upper bounds the sparsity of the primal optimalx (t) , and learning rates iterations to achieve error. The overall complexity is O(s(n + κ) log 1 ).
Our derivation of overall complexity implicitly requires ns ≥ d by setting k = sd/n ≥ 1. This is true for our considered applications like SVM. Otherwise we choose k = 1 and the complexity becomes O(max{d, ns}(1 + κ n ) log 1 ).
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.

Algorithm
Per

Extension to the Trace Norm Ball
We also extend our algorithm to matrix optimization over trace norm constraints: 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: Here 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: X ← arg min X * ≤λ,rank(X)≤s Algorithm 2 Frank-Wolfe Block Primal-Dual Method for Trace Norm Ball 1: Input: Training data A ∈ R n×d , primal and dual step size η, δ > 0. Target accuracy . 2: Initialize: for t = 1, 2, · · · , T do 4: Use Frank Wolfe to Update the primal variable: )-approximation of Eqn. (18).

5:
Update W to maintain the value of AX: 6: Consider the potential dual update:Ỹ 7: Choose greedily the rows of the dual variable to update: let I (t) be the top k coordinates that maximize Update the dual variable accordingly: 8: 9: end for 10: Output: However, an exact solution to (18) requires computing the top s left and right singular vectors of the matrix X (t−1) − 1 ηL (Z + ∇g(X (t−1) ) ∈ R d×c . Therefore we loosely compute an ( 1 2 , /2)-approximation, where is the target accuracy, based on the following definition: F be the objective function in (18), and let l * t = l t (X (t) ). Given parameters γ ≥ 0 and ≥ 0, a feasible solution V to (18) is 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) ≡ AX (t) requires complexity of nds + ncs, while updating Y (t) requires dck operations. Therefore, by setting k ≤ ns(1/c + 1/d), the iteration complexity's dependence on the data size becomes O(n(d + c)s) operations, instead of O(ndc) for conducting a full projected gradient step. Recall that s upper bounds the rank ofX (t) ≤ 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 β-smooth. Its convex conjugate f * i exists and satisfies 1 α -smooth on some convex set (could be R c ) and infinity otherwise.
• g is µ-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 L( With no much modification in the proof, we are able to derive similar convergence guarantees for the trace norm ball.
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 1 constraint to be 300 and the 2 regularize parameter to 10/n to achieve   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.

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 2 F and g to be 2 regularizer with µ = 10/n. The synthetic sensing matrix A ∈ R n×d is dense with n = 1000 and d = 800. Our observation B is of dimension 1000 × 600 and is generated by a ground truth matrix X 0 such that B = AX 0 . Here X 0 ∈ R 800×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.

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.

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 ∇f (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 (x) = f (Ax) + g(x), where A ∈ R n×d . We rephrase the connection for the Frank-Wolfe algorithm. Since the computation of gradient is essentially A ∇f |Ax + ∇g(x), to find its largest magnitude is to search maximum inner products among: whereã i ∈ R n is the i-th column of data matrix A, and ∇f |Ax is the gradient of f at Ax. In this way, we are able to select the greedy coordinates by conducting MIPS for a fixed R 2d×(n+1) matrix [A |I| − A | − I] and each newly generated vector [∇f |Ax |∇g i (x)]. Therefore when ∇g i is constant for linear function or ±λ for g(x) = λ 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 2 regularizer since the different coordinates in ∇ i g(x) creates d new vectors in each iteration and traditional MIPS could resolve it in time sublinear to d. Meanwhile, even with constant ∇ i g(x), it still requires at least 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 (Ax) + g(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.

B.1 Derivation of Primal-Dual Formulation
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 x (assuming function is convex to x) is bounded in a convex set, since we could change f (x, y) to f (x, y) + I C (x), where I C (x) = 0 if x ∈ C and ∞ otherwise. Then the property will be properly inherited.

B.2 Notation and simple facts
Recall primal, dual and Lagrangian forms: Similar to the definitions in [24], we introduce the primal gap defined as ∆ . Recall the assumptions: • f i is convex and β-smooth, and is α strongly convex over some convex set, and linear otherwise.
• g is µ-strongly convex and L-smooth.
To begin with, it is easy to verify that f * i is 1/β-strongly convex and is 1/α-smooth on a convex set and infinity otherwise (See Claim B.6). For simplicity we first assume α ≥ 1 2 β and then generalize the result. Claim B.1.
• Based on our update rule, ∃g ∈ ∂ y 1 n i f * i (y (t) ), such that y (t) And our update rule ensures that I (t) consists of indices i ∈ [n] that maximizes | 1 n a i x (t) − g i |.

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 y ∈ R n , denote y I = i∈I y i e i ∈ R k as the subarray of y indexed by I, with |I| = k. Recall Algorithm 1 selects the coordinates to update in the dual variable as I (t) .
Proof. Notice we have claimed that −D(y) is 1 β -strongly convex and for all g ∈ ∂ y 1 n n i f * i (y (t) ), Meanwhile since −L(x, y) is 1 α -smooth over its feasible set, Also, with the update rule of dual variables, we could make use of Eqn. (19) and re-write Eqn. (20) as: Together we get: Therefore we will connect the progress induced by − y (t) − y (t−1) and dual gap ∆ (t) d next.

B.5 Dual progress
Claim B.4. An α-strongly convex function f satisfies: This simply due to Since −D is 1 β -strongly convex, we get whereĪ is a set of size k that maximizes the values of A ix (t) − g i .
With the relation between 1 n A I (t) x (t) − g I (t) and y (t) − 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.
(from Eqn. (21) and (22)) Similarly to the previous setting, we get the whole primal-dual progress is bounded as follows: Therefore, when we set a proper k and δ such that β Rβ nµ ) log 1 ) iterations for the duality gap ∆ (t) to reach error.

B.7 Smooth Hinge Loss and Relevant Properties
Smooth hinge loss is defined as follows: Our loss function over a prediction p associated with a label i ∈ {±1} will be f i (p) = h(p i ). The derivative of smooth hinge loss h is: Its convex conjugate is: Notice since f i (p) = h( i p), f * i (p) = h * (p/ i ) = h * (p i ).
Claim B.6. For a convex and β-smooth scalar function f , if it is α strongly convex over some convex set, and linear otherwise, then its conjugate function f * is 1/β-strongly convex, and it is a 1/α-smooth function plus an indicator function over some interval [a, b].
Proof. To begin with, since f (x) ≤ β, ∀x, meaning f is β-smooth, then with duality we have f * is 1/β strongly convex [16]. Secondly, since f is α strongly convex over a convex set, meaning an interval for R, therefore f could only be linear on (−∞, a] or [b, ∞), and is α-strongly convex over the set [a, b] (Here for simplicity a < b could be ±∞).  y − a), ∀y ≤ a, and similarly

B.8 Convergence of Optimization over Trace Norm Ball
The convergence analysis for trace norm ball is mostly similar to the case of 1 ball. The most difference lies on the primal part, where our approximated update incur linear progress as well as some error.
Proof. Refer to the proof in [1] we have: L(X (t+1) , Y (t) ) − L(X (t) , Y (t) ) ≤ (1 − µ 8L ) L(X (t) , Y (t) )) − L(X (t) , Y (t) ) + µ 16L Now move the first term on the RHS to the left and rearrange we get: Therefore we get: Now back to the convergence guarantees on the trace norm ball. Figure 3: Convergence result comparison of different algorithms on smoothed hinge loss by varying the coefficient of the regularizer. The first row is the results ran on the rcv1.binary dataset, while the second row is the results ran on the news20.binary dataset. The first column is the result when the regularizer coeffcient µ is set to 1/n. The middle column is when µ = 10/n, and the right column is when µ = 100/n.

C More Results on Empirical Studies
C.1 More experiments with 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 µ that controls the strong convexity of the primal function. Experiments are shown in Figure 3.