Frank-Wolfe Optimization for Symmetric-NMF under Simplicial Constraint

Han Zhao
Machine Learning Department
Carnegie Mellon University
Pittsburgh, PA 15213
han.zhao@cs.cmu.edu

Abstract
Symmetric nonnegative matrix factorization
has found abundant applications in various domains by providing a symmetric low-rank decomposition of nonnegative matrices. In thispaper we propose a Frank-Wolfe (FW) solver
to optimize the symmetric nonnegative matrix
factorization problem under a simplicial constraint, which has recently been proposed forprobabilistic clustering. Compared with existing solutions, this algorithm is simple to implement, and has no hyperparameters to be tuned.
Building on the recent advances of FW algorithms in nonconvex optimization, we provean O(1/ε2 ) convergence rate to ε-approximate
KKT points, via a tight bound Θ(n2 ) on the curvature constant, which matches the best known
result in unconstrained nonconvex setting using
gradient methods. Numerical results demonstrate the effectiveness of our algorithm. As aside contribution, we construct a simple nonsmooth convex problem where the FW algorithm
fails to converge to the optimum. This result
raises an interesting question about necessary
conditions of the success of the FW algorithm
on convex problems.

1

INTRODUCTION

Nonnegative matrix factorization (NMF) has found various applications in data mining, natural language processing, and computer vision [Xu et al., 2003, Liu et al., 2003,
Kuang et al., 2012], due to its ability to provide low rank
approximations and interpretable decompositions. The
decision version of NMF is known to be NP-hard [Vavasis,2009], which also implies that its optimization problem
is NP-hard. Recently, a variant of NMF where the input

Geoff Gordon
Machine Learning Department
Carnegie Mellon University
Pittsburgh, PA 15213
ggordon@cs.cmu.edu

matrix is constrained to be symmetric has become popularfor clustering [He et al., 2011, Kuang et al., 2012, 2015].
The problem is known as symmetric NMF (SymNMF),
and its goal is to minimize ||A − W W T ||2F under theconstraint that W ≥ 0 elementwise, where A is a symmetric matrix of cluster affinities. Compared with NMF,
SymNMF is applicable even when the algorithm does nothave direct access to the data instances, but only theirpairwise similarity scores. Note that in general the input
matrix A does not need to be nonnegative [Kuang et al.,2012]. SymNMF has been successfully applied in manydifferent settings and was shown to be competitive withstandard clustering algorithms; see [Kuang et al., 2012,2015] and the references therein for more details.
In this paper we investigate a constrained version of SymNMF where the input matrix is required to be both nonnegative and positive semidefinite. Furthermore, we requirethat W is normalized such that each row of W sums to
1. This problem has an interesting application in probabilistic clustering [Zhao et al., 2015] where the ith row of
W can be interpreted as the probability that the ith data
point lies in each clusters. Formally, we are interested in
the following optimization problem, which we name as
simplicial SymNMF (SSymNMF):
1||P − W W T ||2F4
subject to W ∈ Rn×k
W 1k = 1n+ ,minimizeW(1)where P ≥ 0 and P ∈ Sn+ is positive semidefinite. (1)was proposed as a formulation for probabilistic clustering [Zhao et al., 2015]. The input matrix P is interpreted
as the co-cluster affinity matrix, i.e., entry Pij corresponds
to the degree to which we encourage data instances xi
and xj to be in the same cluster. Each row of W thencorresponds to the probability distribution of instancexi being in different clusters. A similar simplicial constraint has been considered in NMF as well [Nguyen et al.,2013], where the goal is to seek a probabilistic part-based
decomposition for clustering.

Previous approaches to solve SymNMF or SSymNMF
use the penalty method to convert it into an unconstrained
optimization problem and then solve it iteratively, using either first-order or second-order methods [Kuang
et al., 2012, Zhao et al., 2015, Kuang et al., 2015]. Suchmethods usually include two loops where the outer loop
gradually increases the penalty coefficients and the inner
loop finds a stationary point of each fixed penalized objective, hence they are often computationally expensive andslow to converge. In this paper we first give an equivalent
geometric description of (1) and then propose a variantof the classic Frank-Wolfe (FW) algorithm [Frank andWolfe, 1956], a.k.a. the conditional gradient method [Levitin and Polyak, 1966], to solve it. We also provide a
non-asymptotic convergence guarantee of our algorithm
under an affine invariant stationarity measure (definedin Sec. 2). More specifically, for a given approximation
parameter ε > 0, we show that the algorithm converges toan ε-approximate KKT point of (1) in O(1/ε2 ) iterations.
This rate is analogous to the one derived by Nesterov[2013] for general unconstrained problems (potentiallynonconvex) using the gradient descent method, where
the measure of stationarity is given by the norm of thegradient. The O(1/ε2 ) rate has recently been shown to be
optimal [Cartis et al., 2010] in the unconstrained setting
for gradient methods, and it also matches the best known
rate to a stationary point with (accelerated) projected gradient descent [Ghadimi and Lan, 2016, Ghadimi et al.,2016] in the constrained smooth nonconvex setting.
Contributions. We first give a generalized definition of
the curvature constant [Jaggi, 2013, Lacoste-Julien et al.,2013, Lacoste-Julien, 2016] that works for both convex
and nonconvex functions, and we prove a tight bound of itin (1). We then propose a convergence measure in termsof the duality gap and show that the gap is 0 iff KKT conditions are satisfied. Using these two tools, we proposea FW algorithm to solve (1) and show that it has a nonasymptotic convergence rate O(1/ε2 ) to KKT points. On
the algorithmic side, we give a procedure that has the optimal linear time complexity and constant space complexity
to implement the linear minimization oracle (LMO) in
the FW algorithm. As a side contribution, we construct a
piecewise linear example where the FW algorithm fails
to converge to the optimum. Surprisingly, we can also
show that the FW algorithm works if we slightly change
the objective function, despite that the new function remains piecewise linear and has an unbounded curvature
constant. These two examples then raise an interestingquestion w.r.t. the necessary condition of the successof the FW algorithm. At the end, we conduct several
numerical experiments to demonstrate the efficiency ofthe proposed algorithm by comparing it with the penalty
method and projected gradient descent.

2

PRELIMINARY

2.1

SymNMF UNDER SIMPLICIALCONSTRAINT

One way to understand (1) is through its clustering based
explanation: the goal is to find a probabilistic clusteringof all the instances such that the given co-cluster affinity
Pij for a pair of instances (xi , xj ) is close to the trueprobability that xi and xj reside in the same cluster:Pr(xi ∼ xj ) ==kXPr(ci = cj = h)

h=1
k
X

Pr(ci = h) Pr(cj = h) = wiT wj

h=1

where we use the notation xi ∼ xj to mean “xi and xj
reside in the same cluster”; ci is the cluster assignmentof xi and wi denotes the ith row vector of W . Note thatthe second equation holds because of the assumption ofi.i.d. generation process of instances and their cluster
assignments.
As a first note, the optimal solution W ∗ to (1) is notunique: for any permutation πn over [n], an equivalent solution can be constructed by Wπ∗n = W ∗ Ππn , where Ππn
is a permutation matrix specified by πn . This corresponds
to an equivalence class of W by label switching. Hencefor any fixed k, there are at least k! optimal solutions to(1). The uniqueness of the solution to (1) up to permutation is still an open problem. Huang et al. [2014] studied
sufficient and necessary conditions for the uniqueness of
SymNMF, but they are NP-hard to check in general.
Zhao et al. [2015] proposed a penalty method to transform(1) into an unconstrained problem and solve it via sequential minimization. Roughly speaking, the penalty method
repeatedly solves an unconstrained problem, and enforces
the constraints in (1) by gradually increasing the coefficients of the penalty terms. This process iterates until asolution is both feasible and a stopping criterion w.r.t. theobjective function is met; see [Zhao et al., 2015, Algo. 1].
The penalty method contains 6 different hyperparameters
to be tuned, and it is not even clear whether it will converge to a KKT point of (1). To the best of our knowledge,
no other methods has been proposed to solve (1). To solve
SymNMF, Kuang et al. [2012] proposed a projected Newton method and Vandaele et al. [2016] developed a block
coordinate descent method. However, due to the couplingof columns of W introduced by the simplicial constraint,
it is not clear how to extend these two algorithms to solve(1). On the other hand, the simplicial constraint in (1)restricts the feasible set to be compact, which makes it
possible for us to apply the FW algorithm to solve it.

2.2

Frank-Wolfe ALGORITHM

The FW algorithm [Frank and Wolfe, 1956, Levitin andPolyak, 1966] is a popular first-order method to solve constrained convex optimization problems of the following
form:
minimizex∈D f (x)(2)
where f : Rd → R is a convex and continuously differentiable function over the convex and compact domain D. The FW method has recently attracted a surge
of interest in machine learning due to the fact that itnever requires us to project onto the constraint set D
and its ability to cheaply exploit structure in D [Clarkson, 2010, Jaggi, 2011, 2013, Lacoste-Julien et al., 2013,
Lan, 2013]. Compared with projected gradient descent,it provides arguably wider applicability since projection
can often be computationally expensive (e.g., for the general `p ball), and in some case even computationally intractable [Collins et al., 2008]. At each iteration, the FW
algorithm finds a feasible search corner s by minimizing
a linear approximation at the current iterate x over D:
w(x) := min f (x) + ∇f (x)T (ŝ − x)
ŝ∈D

(3)

The linear minimization oracle (LMO) at x is defined as:T

s = LMO(x) := arg min f (x) + ∇f (x) (ŝ − x) (4)ŝ∈DGiven s = LMO(x), the next iterate is then updated as aconvex combination of s and the current iterate x. If theFW algorithm starts with a corner as the initial point, thenthis property implies that at the t-th iteration, the current
iterate is a convex combination of at most t + 1 corners.
The difference between f (x) and w(x) is known as theFrank-Wolfe gap (FW-gap):
g(x) := f (x) − w(x) = max ∇f (x)T (x − ŝ)
ŝ∈D

(5)

which turns out to be a special case of the general Fenchel
duality gap when we transform (2) into an unconstrained
problem by adding an indicator function over D into theobjective function [Lacoste-Julien et al., 2013, Appendix
D]. Due to the convexity of f , we have the following
inequality: ∀x, y ∈ D, w(y) ≤ f ∗ ≤ f (x), where f ∗is the globally optimal value of f . Specifically, ∀x ∈ D,
w(x) ≤ f (x) and as a result g(x) can be used as an upperbound for the optimality gap: ∀x ∈ D, f (x) − f ∗ ≤g(x), so that g(x) is a certificate for the approximation
quality of the current solution. Furthermore, the duality
gap g(x) can be computed essentially for free in each
iteration: g(x) = ∇f (x)T (x − s).
It is well known that for smooth convex optimization problems the FW algorithm converges to an ε-approximate
solution in O(1/ε) iterations [Frank and Wolfe, 1956,

Dunn and Harshbarger, 1978]. This result has recently
been generalized to the setting where the LMO is solved
only approximately [Clarkson, 2010, Jaggi, 2013]. The
analysis of convergence depends on a crucial concept
known as the curvature constant Cf defined for a convex
function f :
Cf :=

sup
x,s∈D,
γ∈(0,1],y=x+γ(s−x)

2f (y) − f (x) − ∇f (x)T (y − x)2γ(6)The curvature constant measures the relative deviation ofa convex function f from its linear approximation. It isclear that Cf ≥ 0. Furthermore, the definition of Cf only
depends on the inner product of the underlying space,which makes it affine invariant. In fact, the curvature
constant, the FW algorithm, and its convergence analysis
are all affine invariant [Lacoste-Julien et al., 2013]. Later
we shall generalize the curvature constant to a function fthat is not necessarily convex and state our result in termsof the generalized definition. The above definition of thecurvature constant still works for nonconvex functions,but for a concave function f , Cf = 0, which loses its
geometric interpretation as measuring the curvature of f .
To proceed our discussion, we first establish some standard terminologies used in this paper. A continuously
differentiable function f is called L-smooth w.r.t. thenorm || · || if ∇f is L-Lipschitz continuous w.r.t. the norm|| · ||: ∀x, y ∈ D, ||∇f (x) − ∇f (y)|| ≤ L||x − y||.
Throughout this paper, we will use || · || to mean the
Euclidean norm, i.e., the `2 norm for vectors and the
Frobenius norm for matrices if not explicitly specified.
It is standard to assume f to be L-smooth in the convergence analysis of the FW algorithm [Clarkson, 2010,
Jaggi, 2013]. As we will see below, the smoothness of f
also implies the boundedness of Cf on a compact set.

3Frank-Wolfe FOR SymNMF UNDERSIMPLICIAL CONSTRAINT

3.1

A GEOMETRIC PERSPECTIVE

In this section we complement our discussion of (1) in
Sec. 2 with a geometric interpretation, which allows us
to make a connection between the decision version of (1)to the well known problem of completely positive matrix
factorization [Berman and Shaked-Monderer, 2003, Dickinson, 2013, Dickinson and Gijben, 2014]. The decisionversion of the optimization problem in (1) is formulated
as follows:
Definition 1 (SSymNMF). Given a matrix P ∈ Rn×n∩+
Sn+ , can it be factorized as P = W W T for some integer
k such that W ∈ Rn×k
and W 1k = 1n ?
+

Clearly, for P ∈ Sn+ , we can decompose it as P = U U T ,where U ∈ Rn×r , r = rank(P ). Let ui be the ith row
vector of U ; then P is the Gram matrix of a set of rdimensional vectors U = {u1 , . . . , un } ⊆ Rr . Similarly,
W W T can be understood as the Gram matrix of a set ofk dimensional vectors W = {w1 , . . . , wn } ⊆ Rk , where
wi is the ith row vector of W . The simplex constraintin (1) further restricts each wi ∈ ∆k−1 , i.e., wi resides
in the k − 1 dimensional probability simplex. Henceequivalently, SSymNMF asks the following question:
Definition 2. Given a set of n instances U = {ui }ni=1 ⊆Rr , does there exist an integer k and an embedding T :Rr → ∆k−1 , such that inner product is preserved under
T , i.e., ∀i, j ∈ [n], hui , uj i = hT (ui ), T (uj )i ?
An affirmative answer to SSymNMF will give a certificate of the existence of such embedding T : T (ui ) =wi , ∀i ∈ [n]. The goal of (1) can thus be understood as
follows: find an embedding into the probability simplex
such that the discrepancy of inner products between theimage space and the original space is minimized. The
fact that inner product is preserved immediately implies
that distances between every pair of instances are alsopreserved. If k = r, such an embedding is also known as
an isometry. In this case T is unitary. Note in the above
definition of SSymNMF we do not restrict that k = r,
and T does not have to be linear.
SSymNMF is closely connected to the strong membership problem for completely positive matrices [Berman
and Plemmons, 1994, Berman and Shaked-Monderer,
2003]. A completely positive matrix is a matrix P thatcan be factorized as P = W W T where W ∈ Rn×k
for+some k. The set of completely positive matrices forms
a convex cone, and is known as the completely positive
cone (CP cone). From this definition we can see that thedecision version of SymNMF corresponds to the strong
membership problem of the CP cone, which has recently
been shown by Dickinson and Gijben [2014] to be NPhard. Geometrically, the strong membership problem for
CP matrices asks the following question:
Definition 3 (SMEMCP). Given a set of n instances
U = {ui }ni=1 ⊆ Rr , does there exist an integer k andan embedding T : Rr → Rk+ , such that inner product is preserved under T , i.e., ∀i, j ∈ [n], hui , uj i =hT (ui ), T (uj )i ?
Note that the only difference between SSymNMF and
SMEMCP lies in the range of the image space: the former
asks for an embedding in ∆k−1 while the latter only asks
for embedding to reside in Rk+ . SSymNMF is therefore
conjectured to be NP-hard as well, but this assertion has
not been formally proved yet. A simple reduction from
SMEMCP to SSymNMF does not work: a “yes” answer
to the former does not imply a “yes” answer to the latter.

3.2

ALGORITHM

We list the pseudocode of the FW method in Alg. 1 anddiscuss how Alg. 1 can be efficiently applied and implemented to solve SSymNMF. We start by deriving thegradient of f (W ) = 41 ||P − W W T ||2F :∇f (W ) = (W W T − P )W ∈ Rn×k

(7)

The normalization constraint keeps the feasible set decomposable (even though the objective function f is notdecomposable): it can be equivalently represented as the
product of n probability simplices ∆k−1 × · · · × ∆k−1 =:Πn ∆k−1 . Hence at the t-th iteration of the algorithm
we solve the following linear optimization problem overΠn ∆k−1 :minimizeStr S T ∇f (W (t) )(8)subject to S ∈ Πn ∆k−1
Because of the special structure of the constraint set in(8), given ∇f (W ), we can efficiently compute the two
key quantities x(t+1) and gt in Alg. 1 in O(nk) time andO(1) space. The pseudocode is listed in Alg. 2. The key
observation that allows us to achieve this efficient implementation is that at each iteration t ∈ [T ], LMO(x(t) )is guaranteed to be a sparse matrix that contains exactly
one 1 in each row. The time complexity of Alg. 2 is 3nk,
which can be further reduced to 2nk by additional O(n)space to store the index of the nonzero element at eachrow of LMO(x(t) ). Alg. 1, together with Alg. 2 as a subprocedure, is very efficient while at the same time being
simple to implement. Furthermore, it does not have any
hyperparameter to be tuned: this is in sharp contrast withthe penalty method.
Algorithm 1 Frank-Wolfe algorithm (non-convex variant)
Input: Initial point x(0) ∈ D, approximation parameter
ε>0
1: for t = 0 to T do:=
2:Compute s(t)=LMO(x(t) )T(t)arg mins∈D s ∇f (x )3:Compute update direction dt := s(t) − x(t)4:Compute FW-gap gt := −∇f (x(t) )T dt
5:
if gt ≤ ε then return x(t)6:Compute γt := min{gt /C, 1} for any C ≥ C̄f(defined in (10))
7:
Update x(t+1) := x(t) + γt dt
8: end for
9: return x(T )Algorithm 2 Compute next iterate and the gap function
(t)

(t)

Input: ∇f (x ), x
1: gt := h∇f (x(t) ), x(t) i
2: for i = 1 to n do
3:
gt ← gt − minj∈[k] ∇f (x(t) )ij
4: end for
5: Compute γt := min{gt /C, 1}
6: x(t+1) := (1 − γt )x(t)7: for i = 1 to n do
8:
ji := arg minj∈[k] ∇f (x(t) )ij
(t+1)

(t+1)

9:
xiji ← xiji
10: end for
11: return x(t+1) , gt

3.3

+ γt

CONVERGENCE ANALYSIS

In this section we provide a non-asymptotic convergence
rate for the FW algorithm for solving (1) and derive a tightbound for its curvature constant. Our analysis is based on
the recent work [Lacoste-Julien, 2016], where we redefine
the curvature constant so that the new definition works inboth convex and nonconvex settings. Due to space limit,
we only provide partial proofs for theorems and lemmas
derived in this paper, and refer readers to supplementarymaterial for all detailed proofs.
When applied to smooth convex constrained optimization
problems, the FW algorithm is known to converge to an
ε-approximate solution in O(1/ε) iterations. However
the convergence of global optimality is usually unrealistic to hope for in nonconvex optimization, where evenchecking local optimality itself can be computationally intractable [Murty and Kabadi, 1987]. Clearly, to talk about
convergence, we need to first establish a convergence criterion. For unconstrained nonconvex problems, the normof the gradient ||∇f || has been used to measure convergence [Ghadimi et al., 2016, Ghadimi and Lan, 2016],
since limt→∞ ||∇f (x(t) )|| = 0 means every limit pointof the sequence {x(t) } is a stationary point. But such a
convergence criterion is not appropriate for constrained
problems because a stationary point can lie on the boundary of the feasible region while not having a zero gradient.
To address this issue, we will use the FW-gap g(x) as ameasure of convergence. Note that the gap g(x) works as
an optimality gap only if the original problem is convex.
To see why this is also a good measure of convergencefor constrained nonconvex problems, we first prove thefollowing theorem:
Theorem 1. Let f be a differentiable function and
D be a convex compact domain. Define g(x) :=maxŝ∈D ∇f (x)T (x − ŝ). Then ∀x ∈ D, g(x) ≥ 0 andg(x) = 0 iff x is a Karush–Kuhn–Tucker (KKT) point.

Proof. To see g(x) ≥ 0, we have:g(x) := max ∇f (x)T (x − ŝ) ≥ ∇f (x)T (x − x) = 0ŝ∈DReformulate the original constrained problem in the following way:
minimize f (x)xsubject to ID (x) ≤ 0(9)where ID (x) is the indicator function of D which takes
value 0 iff x ∈ D otherwise ∞. Define ND (x) as thenormal cone at x in a convex set D:
ND (x) := {z | zT x ≥ zT y, ∀y ∈ D}
and realize the fact that the subdifferential of the indicatorfunction when x ∈ D is precisely the normal cone at x,i.e., ∂ ID (x) = ND (x), we have:g(x) = 0 ⇔ max ∇f (x)T (x − y) = 0y∈D⇔ ∀y ∈ D, ∇f (x)T (y − x) ≥ 0 ⇔ − ∇f (x) ∈ ND (x)Note that the normal cone for a convex set is a convex
cone, which implies that ∃λ ≥ 0, s.t.,∇f (x) + λ∂ ID (x) = 0 ⇔ ∇x L(x, λ) = 0
where L(x, λ) := f (x) + λ ID (x) is the Lagrangian of(9). By construction, λ ≥ 0 satisfies the dual feasibility
condition and x ∈ D satisfies the primal feasibility condition, which also means the complementary slackness is
satisfied, i.e., λ ID (x) = 0. It follows that g(x) = 0 iff x
is a KKT point of (9).
Remark. Note that in Thm. 1 we do not assume f to beconvex. In fact we can also relax the differentiability of f ,as long as the subgradient exists at every point of f . The
proof relies on the fact that g(x) = 0 implies −∇f (x) isin the normal cone at x. Thm. 1 justifies the use of g(x)as a convergence measure: if limt→∞ g(x(t) ) = 0, thenby the continuity of g, every limit point of {x(k) } is a
KKT point of f . Since being a KKT point is also a sufficient condition for optimality when f is convex, Thm. 1
also recovers the case where g(x) is used as convergencemeasure for convex problems.
For a continuously differentiable function f , we now extend the definition of curvature constant as follows:C̄f :=

supx,s∈D,
γ∈(0,1],y=x+γ(s−x)

2f (y) − f (x) − ∇f (x)T (y − x)γ2(10)Clearly C̄f ≥ 0 and the new definition reduces to Cf
when f is a convex function. In general we have C̄f ≥Cf for any function f . Again, C̄f measures the relative
deviation of f from its linear approximation, and is stillaffine invariant. The difference of C̄f from the original
Cf becomes clear when f is concave: in this case Cf = 0,
but C̄f > 0 and C̄f = C−f . On the downside, the factthat C̄f ≥ Cf means our asymptotic bound is of constanttimes larger than the original one proved by LacosteJulien [2016], but they share the same asymptotic rate interms of the given precision parameter ε. Finally, C̄f = 0
iff f is affine. As in [Jaggi, 2013, Lacoste-Julien et al.,2013], for a smooth function f with Lipschitz constant L
over compact set D, we can bound C̄f in terms of L:
Lemma 1. Let f be a L-smooth function over a
convex compact domain D, and define diam(D) :=supx,y∈D ||x − y||. Then C̄f ≤ diam2 (D)L.

The proof of Lemma 1 does not require f to be convex.
Furthermore, f does not need to be second-order differentiable — being smooth is sufficient. We proceed to derivea convergence bound for Alg. 1 using our new C̄f , which
better reflects the geometric nonlinearity of nonconvex
functions. The main idea of the proof is to bound thedecrease of the gap function by minimizing a quadratic
function iteratively.
Theorem 2. Consider the problem (2) where f is a continuously differentiable function that is potentially nonconvex, but has a finite curvature constant C̄f as definedby (10) over the compact convex domain D. Consider
running Frank-Wolfe (Algo. 1), then the minimal FW gap
g̃T := min0≤t≤T gt encountered by the iterates during
the algorithm after T iterations satisfies:pmax{2h0 C̄f , 2h0 C̄f }
√
g̃T ≤
,
T +1

∀T ≥ 0(11)

where h0 := f (x(0) ) − minx∈D f (x) is the initial global
suboptimality. It thus takes at most O(1/ε2 ) iterations tofind an approximate KKT point with gap smaller than ε.

Lemma 2 (Nesterov [2013]). Let f be twice differentiable. If ||∇2 f (x)||2 ≤ L for all x in the domain, then f
is L-smooth.
Note that (7) is a matrix function of a matrix variable,
whose Hessian is a matrix of order nk×nk. Although one
can compute all the elements of the Hessian by computing
all the partial derivatives ∂∇f (W )ij /∂Wst separately,
this approach is tedious and may hide the structure of
the Hessian matrix. Instead we apply matrix differential
calculus [Magnus and Neudecker, 1985, Magnus, 2010]
to derive the Hessian:Lemma 3. Let f (W ) = 41 ||P − W W T ||2F and define∇2 f (W ) := ∂ vec ∇f (W )/∂ vec W . Then:∇2 f (W ) = W T W ⊗ In + Ik ⊗ (W W T − P )+ (W T ⊗ W )Knk

(12)

where Knk is a commutation matrix such thatKnk vec W = vec W T .
The derivation of the above lemma uses the matrix differential calculus with basic properties of tensors and commutation matrices. Before we proceed, we first present alemma that will be useful to bound the spectral norm of
the above Hessian matrix:
Lemma 4. sup W ≥0, ||W T W ||2 = n.
W 1k =1n

We are now ready to bound the spectral norm of the Hessian ∇2 f (W ) and use it to bound the smoothness constantof f .
Lemma 5. Let c := ||P ||2 . f = 14 ||P − W W T ||2F is(3n + c)-smooth on D = {W ∈ Rn×k
| W 1k = 1n }.
+The above lemma follows from Lemma 2 where we use
Lemma 4 to help bound the spectral norm of the Hessian
matrix in (12).
Bound on diameter of D. The following lemma can beeasily shown:

We comment that the O(1/ε2 ) convergence rate is thesame as the one provided by Lacoste-Julien [2016] using
the original curvature constant Cf for smooth nonconvex
functions, and it is also analogous to the ones derived forgradient descent for unconstrained smooth problems and(accelerated) projected gradient descent for constrainedsmooth problems. The convergence rate of Alg. 1 depends
on the curvature constant of f over D. We now bound thesmoothness constant L and the diameter of the feasible
set in (1).

Lemma 6. Let D = {W ∈ Rn×k
| W 1k = 1n }. Then+
diam2 (D) = 2n with respect to the Frobenius norm.

Bound on smoothness constant. The objective functionin (1) is second-order differentiable, hence we can boundthe smoothness constant by bounding the spectral normof the Hessian instead:

Lemma 7. inf W ≥0, ||∇2 f (W )||2 ≥ n/k 2 − c.

Combining Lemma 6 with Lemma 5 and assuming c isa constant that does not depend on n, we immediately
have C̄f ≤ 2n(3n + c) = O(n2 ) by Lemma 1. A natural
question to ask is: can we get better dependency on n in
the upper bound for C̄f given the special structure thatD = Πn ∆k−1 ? The answer is negative, as we can prove
the following lower bound on the Hessian:
W 1k =1nUsing Lemma 7, we can prove a tight bound on our curvature constant C̄f :Theorem 3. The curvature constant C̄f for f = 14 ||P −W W T ||2F on D = {W ∈ Rn×k
| W 1k = 1n } satisfies:+
2n(n/k 2 − c) ≤ C̄f ≤ 2n(3n + c)
where c := ||P ||2 . Specifically, we have C̄f = Θ(n2 ).
Proof. The upper bound part is clear by combining
Lemma 1 and Lemma 5. We only need to show the lower
bound. Since f is twice-differentiable, we have:C̄f :=

supx,s∈D,
γ∈(0,1],y=x+γ(s−x)=

supx,s∈D,
γ∈(0,1],y=x+γ(s−x)

≥

≥

sup
x,s∈D,
γ∈(0,1],y=x+γ(s−x)

2f (y) − f (x) − ∇f (x)T (y − x)γ22 1|(y − x)T ∇2 f (ξ)(y − x)|
γ2 2
1||y − x||22 · inf ||∇2 f (ξ)||2
ξ∈D
γ2

n

n
1 2
2
γ
diam
(D)
−
c
= 2n( 2 − c)γ2k2kThe second equality is due to the mean-value theorem,
and the third inequality holds by the definition of inf. The
last inequality follows from Lemma 6 and 7.
Combining all the analysis above and using Alg. 2 toimplement the linear minimization oracle in Alg. 1, we
can bound the time complexity of the FW algorithm tosolve (1):Corollary 1. The FW algorithm (Alg. 1) achieves anε-approximate KKT point of (1) in O(n3 k/ε2 ) time.
Proof. In each iteration Alg. 1 takes O(nk) time to compute the gap function as well as the next iterate. Basedon Thm. 2, the iteration complexity to achieve an εapproximate KKT point is O(C̄f /ε2 ). The result follows
from Thm. 3 showing that C̄f = Θ(n2 ).
4

A FAILURE CASE OF Frank-Wolfe
ON A NONSMOOTH CONVEXPROBLEMFor convex problems, the theoretical convergence guarantee of Frank-Wolfe algorithm and its variants depends
crucially on the smoothness of the objective function:[Freund and Grigas, 2016, Garber and Hazan, 2015, Harchaoui et al., 2015, Nesterov et al., 2015] require thegradient of the objective function to be Lipschitz continuous or Hölder continuous; the analysis in [Jaggi, 2013]
requires a finite curvature constant Cf . The Lipschitz

continuous gradient condition is sufficient for the convergence analysis, but not necessary: Odor et al. [2016]shows that Frank-Wolfe works for a Poisson phase retrieval problem, where the gradient of the objective is not
Lipschitz continuous. As we discuss in the last section,
Lipschitz continuous gradient implies a finite curvature
constant. Hence an interesting question to ask is: does
there exist a constrained convex problem with unboundedcurvature constant such that the Frank-Wolfe algorithm
can find its optimal solution? Surprisingly, the answer
to the above question is affirmative. But before we givethe example, we first construct an example where the FW
algorithm fails when the objective function is convex but
nonsmooth.
We first construct a very simple example where the objective function is a piecewise linear and the constraint set isa polytope. We will show that when applied to this simple
problem, the FW algorithm, along with its line search variant and its fully corrective variant, do not even converge
to the optimal solution. In fact, as we will see shortly,
the limit point of the sequence can be arbitrarily far awayfrom the global optimum. Consider the following convex,
constrained optimization problem:
minimize
x1 ,x2

max{5x1 + x2 , −5x1 + x2 }(13)

subject to x2 ≤ 3, 3x1 + x2 ≥ 0, −3x1 + x2 ≥ 0
We plot the objective function of this example in the left
figure of Fig. 1. The unique global optimum for thisproblem is given by x∗ = (0, 0) with f (x∗ ) = 0. The
feasible set D contains three vertices: (−1, 3), (1, 3) and(0, 0). If we apply the FW algorithm to this problem, it isstraightforward to verify that the LMO(x) is given by:
x1 > 0(−1, 3)LMO(x) = (1, 3)
x1 < 0{(−1, 3), (1, 3), (0, 0)} x1 = 0
Note that when x1 = 0, the function is not differentiable,
and the subdifferential is given by conv{(5, 1), (−5, 1)}.
In this case the FW algorithm chooses arbitrary subgradient from the subdifferential and computes the corresponding LMO. From the fundamental theorem of linear
programming, it is easy to see that when x1 = 0, LMO(x)can be any of the three corners depending on the choiceof subgradient at x. Now suppose the FW algorithm stops(0)(0)
in T iterations. For any initial point x(0) = (x1 , x2 )(0)
where x1 6= 0, the final point output by the FW algorithm will be a convex combination of (−1, 3), (1, 3),(0, 0) and x(0) . Let {γt }Tt=1 be the sequence of step sizes
chosen by the FW algorithm. Then we can easily check
QT
(T )
(0)
that x2 ≥ 3 − t=1 (1 − γt )(1 − x2 /3) → 3 asmax{5x1 + x2 , −5x1 + x2 }

3
2.5

2.5
2
x2

2
x2

max{0.5x1 + x2 , −0.5x1 + x2 }

3

1.5

1.5

1

1

0.5

0.5

0

−1.5

0

−1

−0.5

0
x1

2

0.5

4

1

6

0

1.5

8

−1.5

0

0.5

−1

−0.5

1

1.5

0
x1

0.5

2

1

2.5

1.5

3

3.5

Figure 1: Two convex examples of piecewise linear objective functions on a convex polytope. The FW algorithm
fails on the left example but works on the right one. The
objective functions of both examples have unboundedcurvature constants.

It is not hard to see that FW converges to the global optimum, and the curvature constant Cf = ∞ as well for thisnew problem. Combining with example (14), we can seethat Cf < ∞ is not a necessary condition for the successof FW algorithm on convex problems, either. Piecewise
linear functions form a rich class of objectives that are
frequently encountered in practice, while depending on
the structure of the problem, the FW algorithm may or
may not work for them. This thus raises an interesting
problem: can we develop a necessary condition for the
success of the FW algorithm on convex problems? Another interesting question is, can we develop sufficientconditions for piecewise linear functions under which the
FW algorithm converges to global optimum?

5
T → ∞. Note that we can readily change this example
by extending D so that the distance between x(T ) andthe optimum x∗ becomes arbitrarily large. Furthermore,both the line search and the fully-corrective variants failsince the vertices picked by the algorithm remains thesame: {(−1, 3), (1, 3)}. Finally, for any regular probability distribution that is absolutely continuous w.r.t. theLebesgue measure, with probability 1 the initial points
sampled from the distribution will converge to suboptimal
solutions. As a comparison, it can be shown that the subgradient method works for this problem since the function
f itself is Lipschitz continuous [Nesterov, 2013].
Pick x = (−ε, δ), s−x = (1, 0) and y = x+γ(s−x) =(γ − ε, δ), where γ > ε, and plug them in the definitionof the curvature constant Cf . We have:
Cf ≥

lim

2

ε→0, γ>ε γ 2

= lim

γ→0+

(f (y) − f (x) − ∇f (x)T (y − x))

20
=∞
γ

i.e., the curvature constant of this piecewise linear function is unbounded. The problem for this failure case of
FW lies in the fact that the curvature constant is infinity.
On the other hand, we can also show that FW works evenwhen Cf = ∞ by slightly changing the objective function
while keeping the constraint set:1
1
max{ x1 + x2 , − x1 + x2 }(14)2
2
subject to x2 ≤ 3, 3x1 + x2 ≥ 0, −3x1 + x2 ≥ 0
minimize
x1 ,x2The objective function of the second example is shown in
the right figure of Fig. 1. Still, the unique global optimum
for the new problem is given by x∗ = (0, 0) with f (x∗ ) =0, but now the LMO(x) is:LMO(x) = (0, 0),∀x ∈ D

NUMERICAL RESULTS

We evaluate the effectiveness of the FW algorithm(Alg. 1) in solving (1) by comparing it with the penalty
method [Zhao et al., 2015] and the projected gradientdescent method (PGD) on 4 datasets (Table 1). These
datasets are standard for clustering analysis: two of them
are used in [Zhao et al., 2015], and we add two moredatasets with various sizes to make a more comprehensive evaluation. The instances in each dataset are associated with true class labels. For each data set, we use thenumber of classes as the true number of clusters.
Table 1: Statistics about datasets.
Datasetblood
yeast
satimagependigits# inst. (n)
748
1,484
4,435
10,992# feats. (p)
4
8
36
16

# clusters (k)
10
10
6
100

Given a data matrix X ∈ Rn×p , we use a Gaussian kernel
with fixed bandwidth 1.0 to construct the co-cluster matrix
P as Pij = exp(−||xi − xj ||22 ). For each dataset, we use
its number of clusters as the rank of the decomposition,i.e., W ∈ Rn×k . We implement the penalty method
based on [Zhao et al., 2015], where we set the maximum
number of inner loops to be 50, and choose the step factor
for the coefficients of the two penalty terms to be 2. Ineach iteration of PGD, we use backtracking line search
to choose the step size. For all three algorithms, the stop
conditions are specified as follows: if the difference offunction values in two consecutive iterations is smaller
than the fixed gap ε = 10−3 , or the number of (outer)
iterations exceed 50, then the algorithms will stop. Also,for each dataset, all three algorithms share the same initialpoint.
We plot the convergence speed to local optimum of these

8

12

10

8

pendigits

Penalty Method
Projected GD
Frank-Wolfe

14

log function value

9

Penalty Method
Projected GD
Frank-Wolfe

log function value

Penalty Method
Projected GD
Frank-Wolfe

log function value

log function value

satimage

yeast

blood
10

13
12
11

−4

−3

−2

−1

0

−4

1

log time (seconds)

12

11

10

7
−5

Penalty Method
Projected GD
Frank-Wolfe

13

−3

−2

−1

0

1

2

−1

3

log time (seconds)

0

1

2

3

4

5

2

log time (seconds)

4

6

log time (seconds)

Figure 2: Convergence speed of three algorithms. The x-axis measures log-seconds and the y-axis measures log of
objective function value. Note that the intermediate iterates of the penalty method are not feasible solutions, so weshould only compare the convergent point of the penalty method with the other two.
yeast

blood
Projected GD
Frank-Wolfe

Projected GD
Frank-Wolfe

pendigits

Projected GD
Frank-Wolfe

14

14

Projected GD
Frank-Wolfe

12
12

4

log(f − f ∗ )10

6log(f − f ∗ )8log(f − f ∗ )log(f − f ∗ )10satimage

10

10

5

8

8
6

6
2
4

0
−5

−4

−3

−2

−1

0

log time (seconds)

−4

−3

−2

−1

0

1

log time (seconds)

2

4
−1

0

1

2

log time (seconds)

3

1

2

3

4

5

log time (seconds)

Figure 3: Convergence speed comparison between the projected gradient descent method and the FW algorithm. The
x-axis measures log-seconds and the y-axis measures log(f − f ∗ ), where f is the objective function value and f ∗ is the
local optimum achieved by the algorithm. Note that the local optimum f ∗ achieved by PGD and FW can be different.
three algorithms in Fig. 2. Clearly, the penalty method isorders of magnitude slower than both PGD and the FW
algorithm, and it usually converges to a worse solution.
On the other hand, although the FW algorithm tends tohave more iterations before it converges, due to its cheap
computation in each iteration, it consistently takes lesstime than PGD. Another distinction between PGD and
FW is that PGD, when implemented with backtracking
line search, is a monotone descent method, while FW
is not. For a better visualization to compare between
the PGD and the FW algorithms, we omit the penalty
method in Fig. 3, and draw the log-gap plot of the four
datasets. We can confirm from Fig. 3 that both PGD
and the FW algorithm have roughly the same order ofconvergence speed for solving (1), which is consistent
with the theoretical result proved in the previous section(Thm. 2). However, the FW algorithm often converges
faster than the PGD method, in terms of the gap between
the objective function value and local optimum.

6

CONCLUSION

We propose a FW algorithm to solve the SymNMF problem under a simplicial constraint. Compared with existing solutions, the proposed algorithm enjoys a non-asymptotic convergence guarantee to KKT points, is simple to implement, contains no hyperparameter to be tuned,
and is also demonstrated to be much more efficient inpractice. Theoretically, we establish a close connectionof this problem to the famous completely positive matrix factorization by providing an equivalent geometricdescription. We also derive a tight bound on the curvature
constant of this problem. As a side contribution, we givea pair of nonsmooth convex examples where the FW algorithm converges or fails to converge to its optimum. This
result raises an interesting question w.r.t. the necessary
condition of the success of the FW algorithm.
Acknowledgements
HZ would like to thank Simon Lacoste-Julien for his
insightful comments and helpful discussions. HZ and
GG gratefully acknowledge support from ONR, award
number N000141512365.

References
A. Berman and R. J. Plemmons. Nonnegative matrices in
the mathematical sciences. SIAM, 1994.
A. Berman and N. Shaked-Monderer. Completely positive
matrices. World Scientific, 2003.

C. Cartis, N. I. Gould, and P. L. Toint. On the complexity of steepest descent, newton’s and regularized
newton’s methods for nonconvex unconstrained optimization problems. SIAM Journal on Optimization, 20(6):2833–2852, 2010.
K. L. Clarkson. Coresets, sparse greedy approximation,and the Frank-Wolfe algorithm. ACM Transactions on
Algorithms (TALG), 6(4):63, 2010.
M. Collins, A. Globerson, T. Koo, X. Carreras, and P. L.
Bartlett. Exponentiated gradient algorithms for conditional random fields and max-margin markov networks.
Journal of Machine Learning Research, 9(Aug):1775–
1822, 2008.
P. J. Dickinson and L. Gijben. On the computational
complexity of membership problems for the completely
positive cone and its dual. Computational optimization
and applications, 57(2):403–415, 2014.
P. J. C. Dickinson. The copositive cone, the completely
positive cone and their generalisations. Citeseer, 2013.
J. C. Dunn and S. Harshbarger. Conditional gradientalgorithms with open loop step size rules. Journal of
Mathematical Analysis and Applications, 62(2):432–
444, 1978.

algorithm for symmetric decomposition. IEEE Transactions on Signal Processing, 62(1):211–224, 2014.
M. Jaggi. Sparse convex optimization methods for machine learning. PhD thesis, ETH Zürich, 2011.
M. Jaggi. Revisiting Frank-Wolfe: Projection-Free Sparse
Convex Optimization. In ICML, pages 427–435, 2013.
D. Kuang, C. Ding, and H. Park. Symmetric nonnegative matrix factorization for graph clustering. In SDM,
pages 106–117. SIAM, 2012.
D. Kuang, S. Yun, and H. Park. SymNMF: nonnegative
low-rank approximation of a similarity matrix for graphclustering. Journal of Global Optimization, 62(3):545–
574, 2015.
S. Lacoste-Julien.
Convergence rate of FrankWolfe for non-convex objectives. arXiv preprint
arXiv:1607.00345, 2016.
S. Lacoste-Julien, M. Jaggi, M. Schmidt, and P. Pletscher.
Block-Coordinate Frank-Wolfe Optimization for Structural SVMs. In ICML, pages 53–61, 2013.
G. Lan. The complexity of large-scale convex programming under a linear optimization oracle. arXiv preprint
arXiv:1309.5550, 2013.

M. Frank and P. Wolfe. An algorithm for quadratic programming. Naval research logistics quarterly, 3(1-2):
95–110, 1956.

E. S. Levitin and B. T. Polyak. Constrained minimization methods. USSR Computational mathematics andmathematical physics, 6(5):1–50, 1966.

R. M. Freund and P. Grigas. New analysis and results forthe frank–wolfe method. Mathematical Programming,
155(1-2):199–230, 2016.

W. Liu, N. Zheng, and X. Lu. Non-negative matrix factorization for visual coding. In Acoustics, Speech, and Signal Processing, 2003. Proceedings.(ICASSP’03). 2003
IEEE International Conference on, volume 3, pages
III–293. IEEE, 2003.

D. Garber and E. Hazan. Faster Rates for the Frank-Wolfe
Method over Strongly-Convex Sets. In ICML, pages
541–549, 2015.
S. Ghadimi and G. Lan. Accelerated gradient methods
for nonconvex nonlinear and stochastic programming.
Mathematical Programming, 156(1-2):59–99, 2016.
S. Ghadimi, G. Lan, and H. Zhang. Mini-batch stochastic
approximation methods for nonconvex stochastic composite optimization. Mathematical Programming, 155(1-2):267–305, 2016.
Z. Harchaoui, A. Juditsky, and A. Nemirovski. Conditional gradient algorithms for norm-regularized smoothconvex optimization. Mathematical Programming, 152(1-2):75–112, 2015.
Z. He, S. Xie, R. Zdunek, G. Zhou, and A. Cichocki. Symmetric nonnegative matrix factorization: Algorithms
and applications to probabilistic clustering. IEEE
Transactions on Neural Networks, 22(12):2117–2131,
2011.
K. Huang, N. D. Sidiropoulos, and A. Swami. Nonnegative matrix factorization revisited: Uniqueness and

J. R. Magnus. On the concept of matrix derivative. Journal of Multivariate Analysis, 101(9):2200–2206, 2010.
J. R. Magnus and H. Neudecker. Matrix differential calculus with applications to simple, hadamard, and kronecker products. Journal of Mathematical Psychology,
29(4):474–492, 1985.
K. G. Murty and S. N. Kabadi. Some np-complete problems in quadratic and nonlinear programming. Mathematical programming, 39(2):117–129, 1987.
Y. Nesterov. Introductory lectures on convex optimization: A basic course, volume 87. Springer Science &Business Media, 2013.
Y. Nesterov et al. Complexity bounds for primal-dual
methods minimizing the model of objective function.
Center for Operations Research and Econometrics,
CORE Discussion Paper, 2015.
D. K. Nguyen, K. Than, and T. B. Ho. Simplicial nonnegative matrix factorization. In Computing and Communication Technologies, Research, Innovation, and Visionfor the Future (RIVF), 2013 IEEE RIVF International
Conference on, pages 47–52. IEEE, 2013.
G. Odor, Y.-H. Li, A. Yurtsever, Y.-P. Hsieh, Q. TranDinh, M. El Halabi, and V. Cevher. Frank-Wolfe works
for non-Lipschitz continuous gradient objectives: Scalable poisson phase retrieval. In Acoustics, Speech andSignal Processing (ICASSP), 2016 IEEE International
Conference on, pages 6230–6234. IEEE, 2016.
A. Vandaele, N. Gillis, Q. Lei, K. Zhong, and I. Dhillon.
Efficient and non-convex coordinate descent for symmetric nonnegative matrix factorization. IEEE Transactions on Signal Processing, 64(21):5571–5584, 2016.
S. A. Vavasis. On the complexity of nonnegative matrix
factorization. SIAM Journal on Optimization, 20(3):1364–1377, 2009.
W. Xu, X. Liu, and Y. Gong. Document clustering basedon non-negative matrix factorization. In Proceedings ofthe 26th annual international ACM SIGIR conference
on Research and development in informaion retrieval,
pages 267–273. ACM, 2003.
H. Zhao, P. Poupart, Y. Zhang, and M. Lysy. Sof: Softcluster matrix factorization for probabilistic clustering.
In AAAI, pages 3188–3195, 2015.

A

MISSING PROOFS

Lemma 1. Let f be a L-smooth function over a convex compact domain D, and define diam(D) := supx,y∈D ||x − y||.
Then C̄f ≤ diam2 (D)L.
Proof. Let ∀x, s ∈ D, γ ∈ (0, 1], and y = x + γ(s − x). The smoothness of f implies that f is continuously
differentiable, hence we have:f (y) − f (x) − ∇f (x)T (y − x)
Z 1
T=(∇f (x + t(y − x)) − ∇f (x)) (y − x) dt
≤
≤
≤

0
Z 1
0

Z 1
0

Z 1
0

(Mean-value theorem)T(∇f (x + t(y − x)) − ∇f (x)) (y − x) dt(Triangle inequality)||∇f (x + t(y − x)) − ∇f (x)|| · ||y − x|| dt(Cauchy-Schwarz inequality)tLγ 2 ||s − x||2 dt ≤Lγ 2
diam2 (D)2

(Smoothness assumption of f )It immediately follows thatC̄f ≤2 Lγ 2
diam2 (D) = diam2 (D)Lγ2 2Theorem 2. Consider the problem (2) where f is a continuously differentiable function that is potentially nonconvex,
but has a finite curvature constant C̄f as defined by (10) over the compact convex domain D. Consider runningFrank-Wolfe (Algo. 1), then the minimal FW gap g̃T := min0≤t≤T gt encountered by the iterates during the algorithm
after T iterations satisfies:pmax{2h0 C̄f , 2h0 C̄f }√g̃T ≤, ∀T ≥ 0(11)
T +1
where h0 := f (x(0) ) − minx∈D f (x) is the initial global suboptimality. It thus takes at most O(1/ε2 ) iterations to find
an approximate KKT point with gap smaller than ε.
Proof. Let y := x + γd, where d := s − x is the update direction found by the LMO in Alg. 1. Using the definition ofC̄f , we have:f (y) = f (y) − f (x) − γ∇f (x)T d + f (x) + γ∇f (x)T d≤ f (x) + γ∇f (x)T d + f (y) − f (x) − γ∇f (x)T d≤ f (x) + γ∇f (x)T d +

γ2
C̄f
2

Now using the definition of the FW gap g(x) and for ∀C ≥ C̄f , we get:f (y) ≤ f (x) − γg(x) +

γ2
C̄f ,
2

∀γ ∈ (0, 1](15)

Depending on whether C > 0 or C = 0, the R.H.S. of (15) is a either a quadratic function with positive second ordercoefficient or an affine function. In the first case, the optimal γ ∗ that minimizes the R.H.S. is γ ∗ = g(x)/C. In thesecond case, γ ∗ = 1. Combining the constraint that γ ∗ ≤ 1, we have γ ∗ = min{1, g(x)/C}. Thus we obtain:f (y) ≤ f (x) − min

 2



g (x)
C
, g(x) −Ig(x)>C2C
2(16)

(16) holds for each iteration in Alg. 1. A cascading sum of (16) through iteration step 1 to T + 1 shows that:f (x(T +1)) ≤ f (x

(0)

T
X


 2 (t) Cg (x )(t))−, g(x ) −Ig(x(t) )>C
min
2C
2
t=0

(17)

Define g̃T := min0≤t≤T g(x(t) ) be the minimal FW gap in T + 1 iterations. Then we can further bound inequality (17)
as:
 2 


g̃T
C
f (x(T +1) ) ≤ f (x(0) ) − (T + 1) min, g̃T −
Ig̃T >C(18)
2C
2
We discuss two subcases depending on whether g̃T > C or not. The main idea is to get an upper bound on g̃T byshowing that g̃T cannot be too large, otherwise the R.H.S. of (18) can be smaller than the global minimum of f , which
is a contradiction. For the ease of notation, define h0 := f (x(0) ) − minx∈D f (x), i.e., the initial gap to the global
minimum of f .
g̃ 2

TCase I. If g̃T > C and g̃T − C2 ≤ 2C
, from (18), then:0 ≤ f (x(T +1) ) − min f (x) ≤ f (x(0) ) − min f (x) − (T + 1)(g̃T −
x∈D

x∈D

C
C
) = h0 − (T + 1)(g̃T − )22which implies
C < g̃T ≤

h0
C
2h0 C
+
⇒ g̃T ≤= O(1/T )
T +1
2
T +1

On the other hand, solving the following inequality:
C−

C
C
g̃ 2
4h20 C 2 1≤ g̃T −
≤ T ≤
2
2
2C
(T + 1)2 2C

we getT + 1 ≤ 2h0
This means that g̃T decreases in rate O(1/T ) only for at most the first 2h0 iterations.
g̃ 2

TCase II. If g̃T ≤ C or g̃T − C2 > 2C
. Similarly, from (18), we have:

g̃ 2
g̃ 2
0 ≤ f (x(T +1) ) − min f (x) ≤ f (x(0) ) − min f (x) − (T + 1) T = h0 − (T + 1) T
x∈D
x∈D
2C
2C
which yields
r
g̃T ≤

2h0 C
T +1

0C
Combining the two cases together, we get g̃T ≤ 2h
T +1 if T + 1 ≤ 2h0 ; otherwise g̃T ≤√T + 1 ≤ T + 1, thus we can further simplify the upper bound of g̃T as:√
max{2h0 C, 2h0 C}
√
g̃T ≤
T +1

q

2h0 C
T +1 . Note that for T ≥ 0,Lemma 3. Let f (W ) = 14 ||P − W W T ||2F and define ∇2 f (W ) := ∂ vec ∇f (W )/∂ vec W . Then:∇2 f (W ) = W T W ⊗ In + Ik ⊗ (W W T − P )+ (W T ⊗ W )Knk

where Knk is a commutation matrix such that Knk vec W = vec W T .

(12)

Proof. Using the theory of matrix differential calculus, the Hessian of a matrix-valued matrix function is defined as:∂ vec ∇f (W )∂ vec W
Using the differential notation, we can compute the differential of ∇f (W ) as:∇2 f (W ) :=d∇f (W ) = d(W W T − P )W = (dW )W T W + W (dW )T W + W W T dW − P dW

Vectorize both sides of the above equation and make use of the identity that vec(ABC) = (C T ⊗ A) vec B for A, B, C
with appropriate shapes, we get:vec d ∇f (W ) = (W T W ⊗ In ) vec dW + (W T ⊗ W ) vec dW T + (Ik ⊗ (W W T − P )) vec dW

Let Knk be a commutation matrix such that Knk vec W = vec W T . We can further simplify the above equation as:
vec d∇f (W ) = W T W ⊗ In + (W T ⊗ W )Knk + Ik ⊗ (W W T − P ) vec dW(19)
It then follows from the first identification theorem [Magnus and Neudecker, 1985, Thm. 6] that the Hessian is given by∇2 f (W ) = W T W ⊗ In + Ik ⊗ (W W T − P ) + (W T ⊗ W )Knk ∈ Rnk×nk

As a sanity check, the first two terms in ∇2 f (W ) are clearly symmetric. The third term can be verified as symmetric as−1
T
well by realizing that Knk
= Knk
, and
W ⊗ W T = Knk (W T ⊗ W )Knk

Lemma 4. sup W ≥0, ||W T W ||2 = n.
W 1k =1n

Proof. ∀W ≥ 0, if W 1k = 1n , then by the Courant-Fischer theorem:||W T W ||2 := max ||W T W v||2

(Courant-Fischer theorem)= max ||W T W v||2

(Perron-Frobenius theorem)≤ max ||W T W v||2

(B2 (0, 1) ⊆ B∞ (0, 1))= ||W T 1n ||2(W ≥ 0, W 1k = 1n )v∈Rk ,||v||2 =1v∈Rk+,||v||2 =1v∈Rk+,||v||∞ ≤1

T≤ ||W 1n ||1 = n

To achieve this upper bound, consider W = 1n eT1 , where e1 is the first column vector of the identity matrix Ik .
In this case W T W = e1 1Tn 1n eT1 = ne1 eT1 , which is a rank one matrix with a positive eigenvalue n. Hencesup ||W T W ||2 = n.
Lemma 5. Let c := ||P ||2 . f = 14 ||P − W W T ||2F is (3n + c)-smooth on D = {W ∈ Rn×k
| W 1k = 1n }.
+Proof. Recall that the spectral norm || · ||2 is sub-multiplicative and the spectrum of A ⊗ B is the product of thespectrums of A and B. Using (12), we have:||∇2 f (W )||2 = ||W T W ⊗ In + Ik ⊗ (W W T − P ) + (W T ⊗ W )Knk ||2≤ ||W T W ⊗ In ||2 + ||Ik ⊗ (W W T − P )||2 + ||(W T ⊗ W )Knk ||2
T

T

T= ||W W ||2 ||In ||2 + ||Ik ||2 ||W W − P ||2 + ||W ⊗ W ||2 ||Knk ||2
T

T

T= ||W W ||2 + ||W W − P ||2 + ||W ⊗ W ||2
T≤ 3||W W ||2 + ||P ||2≤ 3n + c

The result then follows from Lemma 2.

(Triangle inequality)(submultiplicativity of || · ||2 )(||In ||2 = ||Ik ||2 = ||Knk ||2 = 1)
(Triangle inequality)(Lemma 4)Lemma 6. Let D = {W ∈ Rn×k
| W 1k = 1n }. Then diam2 (D) = 2n with respect to the Frobenius norm.
+
Proof.
diam2 (D) = sup ||W − Z||2F
W,Z∈D

= sup

X
X
2
(Wij − Zij )2 = supWij2 + Zij− 2Wij Zij

W,Z∈D ij

≤ sup

W,Z∈D

X
W,Z∈D

W,Z∈D ij2
Wij2 + Zij
≤ sup

W,Z∈D

X

Wij + Zij

W,Z∈D

= 2n
Note that choosing W = 1eT1 and Z = 1n eT2 make all the equalities hold in the above inequalities. Hence diam2 (D) =2n.
Lemma 7. inf W ≥0, ||∇2 f (W )||2 ≥ n/k 2 − c.
W 1k =1n

Proof. For a matrix A, we will use σi (A) to mean the ith largest singular value of A and λmax (A), λmin (A) to mean thelargest and smallest eigenvalues of A, respectively. Recall ∇2 f (W ) = W T W ⊗In +Ik ⊗(W W T −P )+(W T ⊗W )Knk .
For W ≥ 0, W 1k = 1n , let r = rank(W ). Clearly r ≥ 1. We have the following inequalities hold:||∇2 f (W )||2 = ||W T W ⊗ In + Ik ⊗ (W W T − P ) + (W T ⊗ W )Knk ||2≥ λmax W W T ⊗ In + (W T ⊗ W )Knk + λmin Ik ⊗ (W W T − P )≥ λmax (W W T ⊗ In ) + λmin (W T ⊗ W )Knk + λmin Ik ⊗ (W W T − P )(Weyl’s inequality)= λmax (W W T ) + λmin (W T ⊗ W ) + λmin (W W T − P )≥ λmax (W W T ) + λmin (W T ⊗ W ) + λmin (W W T ) − λmax (P )= σ12 (W ) + 2σr2 (W ) − λmax (P )≥ σ12 (W ) − c
1≥ ||W ||2F − cr1≥ ||W ||2F − cknk
1 XX 2=W −ck i=1 j=1 ij
n

1X
≥
k
k i=1
n
= 2 −ck

Pk

j=1 Wijk(||P ||2 ≤ ||P ||F )(r · σ12 (W ) ≥ ||W ||2F )(rank(W ) ≤ k)!2−cwhere the first three inequalities all follow from Weyl’s inequality.

(Cauchy ineq.)