Polynomial optimization methods formatrix factorization
Po-Wei Wang

Chun-Liang Li

J. Zico Kolter

Machine Learning Department
Carnegie Mellon University
Pittsburgh, PA 15213
poweiw@cs.cmu.edu

Machine Learning Department
Carnegie Mellon University
Pittsburgh, PA 15213
chunlial@cs.cmu.edu

School of Computer Science
Carnegie Mellon University
Pittsburgh, PA 15213
zkolter@cs.cmu.edu

Abstract
Matrix factorization is a core technique in many machinelearning problems, yet also presents a nonconvex and often
difficult-to-optimize problem. In this paper we present anapproach based upon polynomial optimization techniques thatboth improves the convergence time of matrix factorization
algorithms and helps them escape from local optima. Our
method is based on the realization that given a joint search
direction in a matrix factorization task, we can solve the “subspace search” problem (the task of jointly finding the steps
to take in each direction) by solving a bivariate quartic polynomial optimization problem. We derive two methods forsolving this problem based upon sum of squares moment relaxations and the Durand-Kerner method, then apply thesetechniques on matrix factorization to derive a direct coordinate descent approach and a method for speeding up existing
approaches. On three benchmark datasets we show the method
substantially improves convergence speed over state-of-the-art
approaches, while also attaining lower objective value.

1

Introduction

Matrix factorization has become a foundational operationin many machine learning tasks. Although there are several
variants of matrix factorization, one of the most common isthe collaborate filtering setting, where we wish to find a lowrank reconstruction of some matrix S ∈ RM ×N , where we
only observe the entries Sij at some relatively small numberof points (i, j) ∈ S. The goal is then to find some A ∈RM ×K and B ∈ RN ×K such thataTi bj ≈ Sij , ∀(i, j) ∈ S(1)where ai denotes the ith row of A and bj denotes the jth row
of B. Although many methods to perform this factorization
are possible, a natural approach is to solve the optimization
problemminimize L(A, B),
(2)
A,B

where
L(A, B) ≡X

(aTi bj − Sij )2 + λ kAk2F + kBk2F ,(i,j)∈S

in which A and B are the optimization variables and λ isa regularization parameter. This approach and extensions
Copyright c 2017, Association for the Advancement of Artificial
Intelligence (www.aaai.org). All rights reserved.

have been used in a number of collaborative filtering andrecommender system applications (Bell and Koren 2007;
Dror et al. 2012) and researchers have developed numerous
optimization approaches to solving these problems, discussed
in more detail below.
Despite it ubiquity, (2) is a non-convex optimization problem as it involves quartic terms in the A and B variables.
A number of approaches can still be used to optimize theobjective — two popular and straightforward approaches, forexample, are to use alternating minimization, since the objective is convex in either A or B separately, or simple gradientdescent in the non-convex objective — but the possibility of
local optima is real and a non-trivial reality of matrix factorization methods. A low-dimension example of (2) is given
in Appendix B, in which both methods stuck if not properly
initialized.
In this paper we propose an approach that both helps tomitigate the possibility of local optima and speeds up matrix
factorization approaches above the current state of the art.
The basic idea is straightforward: suppose we want to optimize a single pair of coefficient aik and bkj , or more generally suppose we have some joint search direction (∆A, ∆B)
over A and B (the case of a single set of coordinates simply
being when ∆A and ∆B have only one non-zero component). Then we can consider the problem of optimizing theobjective just over this subspace
minimize L(A + α∆A, B + β∆B)
α,β

(3)

where α, β ∈ R are our optimization variables. This problem
becomes a quartic polynomial in the two variables, α and β,which can be solved exactly using sum of squares methods (?,SOS,)]parrilo2003semidefinite,blekherman2013semidefinite;
alternatively, for this particular case we show the problem
can be reduced to finding the root of a quintic polynomial ina single variable, and solved exactly using the Durand-Kerner
method (also called the method of Weierstrass) (Durand 1959;
Petkovi et al. 1995) (we show, this approach is often several
orders of magnitude faster than the SOS-based methods).
Using this “exact subspace search” technique, we both derivea pure coordinate descent approach to matrix factorizationand modify existing optimization approaches to incorporate
this subspace search when adjusting the parameters. For
the case of the pure coordinate descent approach, we derive
an efficient caching mechanism that lets us efficiently

update individual coordinates with minimal passes over the
training data. Our results show that our method improves
substantially over the existing state of the art in terms of
optimization speed, and also allows the methods to morerobustly escape from local optima.

2

Background and related work

Several solution methods have been proposed to solve the optimization problem (2) and similar problems. We discuss the
popular approaches briefly, and focus on two recent methods
that are considered as the state of the art by some metrics:
coordinate descent (Pilászy, Zibriczky, and Tikk 2010) andstochastic gradient descent (Koren, Bell, and Volinsky 2009).
As mentioned previously, problem (2) is nonconvex. However, when either A or B is fixed, the problem becomesan convex quadratic with an analytical solution. The alternating least square method (ALS)(Zhou et al. 2008) exploits this observation and iteratively solves the two leastsquares problems. The ALS is strictly decreasing but has an
O(|S|K 2 + (M + N )K 3 ) complexity per iteration, and thus
is not efficient on a large dataset with non-trivial K.
Another increasingly popular method is stochastic gradientdescent (SGD)(Koren, Bell, and Volinsky 2009). For each
update, the SGD randomly selects a pair (i, j) ∈ S andupdates the corresponding variables via
ai := ai −η(λai +Rij bj ),bj := bj −η(λbj +Rij ai ), (4)where Rij ≡ aTi bj − Sij . The SGD enjoys a O(|S|K +(M + N )K) complexity per data pass, and several parallel algorithms have been developed (Gemulla et al. 2011;
Recht and Ré 2013; Zhuang et al. 2013). However, it is notguaranteed to be strictly decreasing and requires fine-tuning
the additional step-size parameter η in practice.
Updating one variable at a time results in the coordinate descent algorithm (Pilászy, Zibriczky, and Tikk 2010).
With proper caching and update order, this method enjoys
O(|S|K + (M + N )K) complexity and is guaranteed to
be strictly decreasing. Surprisingly, the coordinate descentapproach for matrix factorization is fully parallelizable. To
see this, let us fix all the variables except for the variable aik ,the analytical solution for this variable isP
j|(i,j)∈S Rij/k bjk
P
aik := CCD(i, b:k , R:/k ) = −, (5)λ + j|(i,j)∈S b2jk
where Rij/k = Rij − aik bjk . Note that the update formula
is independent of aik for all i, because Rij/k has no dependence on aik terms. Thus, the update on aik for all i can be
performed concurrently. This extension is called CCD++ (Yuet al. 2012).
Background on polynomial optimization Although thenonconvexity of (2) can make it difficult to solve, not all nonconvex problems are intractable. Specifically, optimizing apolynomial function, including non-convex and multi-modal
polynomials, is a well-studied numerical problem, and there
exist multiple solution techniques that work well in practice. For example, one well-known approach to polynomial

optimization is via the sum of squares (SOS)(Parrilo 2003;
Blekherman, Parrilo, and Thomas 2013) and moment relaxation (Lasserre 2001; Laurent 2009) approach, which approximates the solution of a polynomial optimization problem
by a semidefinite program. For bivariate quartic polynomials (the precise type that will appear in our subspace searchmethod), the SOS relaxation is exact and we can recover
the exact solution to the polynomial optimization problem.
In addition, there are several approaches for root finding of
single variable polynomials, and we will exploit one of thesemethods, the Durand-Kerner approach (Durand 1959), to develop an even faster method for solving the bivariate quartic
optimization problem.

3

Polynomial optimization for matrix
factorization

First, we formally define the subspace search problem (3).
Consider an arbitrary pair of search directions U ≡ ∆A andV ≡ ∆B. In the subspace search, we want to solve theoptimization problem
minimize L(A + αU, B + βV ).
α,β

(6)

Expanding the definition of L in (2) and then shifting andscaling (see Appendix A for details), solving this problem is
equivalent tominimize F (α, β),(7)α,βwhere F (α, β) is defined as
1
1
1
C22 α2 β 2 + C11 αβ + C20 α2 + C10 α + C02 β 2 + C01 β,222in which C20 and C02 are strictly positive. In particular, as
we show in Appendix A, the coefficients on any possible α2 βand αβ 2 terms can be removed. When C22 = 0, the problem
can be solved by a linear equation by taking derivatives on α
and β. Otherwise, with proper scaling, we can assume that
C22 = 1 without loss of generality.

3.1

Solving polynomial optimization

The subspace search problem (7) is in general not convex andadmits local minima. However, it can be solved efficiently
via polynomial optimization techniques. We here demonstrate two methods that solves the bivariate quartic problem
globally.
Moment relaxation SDP formulation The moment relaxation approach to solving the polynomial optimization problem, which is the dual of the SOS relaxation, is based upon
forming a positive semidefinite moment matrix that consistsof various polynomial degree terms of α and β. Specifically,>
define φ = (1 α β αβ) and let
1 y01 y10 y11
y02 y11 y12  y
Y = E[φφT ] =  01,
(8)
y10 y11 y20 y21 
y11 y12 y21 y22

where the expectation is over a properly constructed measure (Parrilo 2003). Then we can solve the semidefinite
program1
1minimize y22 + C11 y11 + C20 y20
Y 0
2
2
(9)
1
+ C10 y10 + C02 y02 + C01 y01 ,2with the guarantee that we are able to reconstruct α∗ andβ ∗ of (7) from the optimal solution Y ∗ of (9). This result is
important theoretically, because it provides a basis for the
global solvability of the such polynomial (bivariate quartic)
optimization problems.
Durand-Kerner method In practice, even for a relatively
small polynomial problem, solving the SDP as (9) is slow.
Therefore, instead of solving (7) directly through solving a
SDP, in this case we can transform it to a univariate polynomial and solve this directly using a root-finding method.
Specifically, considering the derivatives of the objective∇α F (α, β) = β(αβ + C11 ) + C20 α + C10(10)∇β F (α, β) = α(αβ + C11 ) + C02 β + C01(11)
We can perform a symmetric linear transform and write theseequations aspp
C02 · ∇α F (α, β) + C20 · ∇β F (α, β)(12)= D(t + D)x + D∆+ = 0,pp
C02 · ∇α F (α, β) − C20 · ∇β F (α, β)(13)= D(t − D)y + D∆− = 0,
where we introduce the new variables t, x, and yt = αβ + C11 ,
p
p
p
p
(14)
x = C20 α + C02 β, y = C20 α − C02 β,
and constants D, ∆+ , and ∆−p
D = C20 C02 ,
(15)
C10
C01
C10
C01
∆+ = √+√, ∆− = √
−√
.
C20
C02
C20
C02
Note that C02 and C20 are positive.
2−y 2
Starting now with the equality t − C11 = x 4D
, and
2
2
multiplying both sides by (t + D) (t − D) we get

1. Initialize ti , i = 1, . . . , 5 on a complex disk that willcontain all the roots.
2. While |f (ti )| ≥  for any i: repeatf (ti )
, i = 1, . . . , 5.
j6=i (ti − tj )ti ← ti − Q(18)

The Durand-Kerner method generates all 5 (complex) rootsof (17), including all the local minima. We choose the realroot that minimize the subspace search. To guarantee that our
initial disk contains all the roots requires potentially restartingthe Durand-Kerner method multiple times, but in practice
we are able to choose a disk size based upon the highest
magnitude coefficient C and have yet to observe divergence(and importantly, lack of divergence for the Durand-Kerner
method signifies that we have obtained the exact solution).

3.2

Polynomial coordinate descent

Using this exact bivariate subspace search procedure, wedevelop two approaches to solving the matrix factorization
problem. The first is the direct approach, where we extend
a coordinate descent approach to simultaneously optimizea pair of coefficients aik and bjk , where (i, j) will cycle
through the entries of S. In this setting, our optimization over
α and β becomes
F (α, β) =1 2 2
α β + Rij/k αβ
2
1 i/j
1 j/i
i/j
j/i
+ C20 α2 + C10 α + C02 β 2 + C01 β,
2
2
(19)

where
i/j

i
C10 = C10
− Rij/k bjk ,j/ii
C01 = C01
− Rij/k aik .

i
C20 = C20
− b2jk ,j
C02 = C02
− a2ik ,
i
C20
=λ+

X

b2qk ,

i/j

j/i

i
C10
=

q|i,q∈S
j
C02
=λ+

X
p|p,j∈S

X(20)

Riq/k bqk ,(21)

Rpj/k apk

(22)

q|i,q∈S

a2pk ,jC01=

X
p|p,j∈S(t − C11 )(t + D)2 (t − D)21
1=((t + D)2 x2 )(t − D)2 −((t − D)2 y 2 )(t + D)2 .
4D
4D
(16)
Applying (12), (13) to the above equality and rearrange theterms, we can finally write our gradient conditions in termsof a single quintic polynomial equation

Although the notation is somewhat cumbersome, these terms
just signify the total C coefficients with the coefficients involves the k terms subtracted out. The advantage to maintaining the coefficients in this form, is that after updating a single
aik , bjk pair, we can update all the coefficients in O(1) time,by simply adding back the new values. This bookkeeping
allows us to perform a pass over all (i, j) ∈ S (for a fixed k)in O(|S|) time. The full algorithm with the relevant updating
and downdating, is given in Algorithm 1.

f (t) ≡ (t − C11 )(t + D)2 (t − D)2
1 2
1 2
−(
∆ )(t − D)2 + (∆ )(t + D)2 = 0.
4D +
4D −
(17)
At this point, we can apply the Durand-Kerner method tosolve (17), which takes the form

Although the direct coordinate descent method is efficient
asymptotically, it still requires substantial computational
work, as it needs to perform polynomial root finding once
for each (i, j) ∈ S, and thus the constant term on the O(|S|)
complexity is quite large. An alternative which often works

3.3

Polynomial subspace search

for k = 1, . . . , K do
Rij/k = Rij − aik bjk , ∀(i, j) ∈ S;jj
Construct C01, C02
using (22) for all j;
for i = 1, . . . , M do
j/i
j/i
j
j
Update C01 , C02 from C01, C02
via (20) for
all j;
i
i
Construct C10, C20
by (21);for j | (i, j) ∈ S do
i/j
i/j
i
i
Update C10 , C20 from C10, C20
via (20);
Solve aik , bjk = arg minα,β F (α, β) with
i/j
i/j
j/i
j/i
coefficient C20 , C10 , C02 , C01 in (19);
i/j
i/j
i
i
Update C10, C20
from C10 , C20 via (20);
end
j/i
j/i
j
j
Update C01, C02
from C01 , C02 via (20) for
all j;
end
Rij = Rij/k + aik bjk , ∀(i, j) ∈ S;
end
Algorithm 1: PolyMF-CD

better in practice is to use an alternative matrix factorization algorithm to obtain search directions for an entire column/row of A and B respectively, and use the subspace
search to find the adjustment of A and B in this direction.
The subspace search could be applied to virtually any matrix factorization approach, including gradient descent, ALS,
SGD, and CCD++. In this paper, we adopt CCD++ for thesubspace it outperforms other algorithms in our benchmark
test. Concretely, for a given k we use CCD++ to compute thecoordinate descent direction of an entire row/column of A
and B respectively, and use our polynomial optimization tofind the exact minimum over this joint search direction. The
algorithm is described in Algorithm 2.
Finally, we note that a few technical issues come up when
applying this subspace search to CCD++. We keep trackingof the function decrease in inner iterations, and skip the inner
iteration early when the function decrease is too small, as
proposed in (Yu et al. 2014). Also, we skip the line search
procedure in the first iteration because the first few steps of
CCD++ is usually very large and avoiding subspace search
provides better initialization. The whole algorithm is also
fully parallelizable like CCD++, and we study the speedupof parallelism in Section 3.4.

3.4

Convergence

In this section we prove that the PolyMF-CD and PolyMFSS (with CCD++ search direction) algorithms converge toa stationary point. Although this a relatively weak form a
convergence, we know of no previous such guarantees for
coordinate descent methods for matrix factorization (similar
results exist for the cases of alternating least squares andmultiplicative updates).
Theorem 1. Suppose λ > 0. Then every limit point of the
PolyMF with coordinate descent direction and CCD direction
is a stationary point.

for k = 1, . . . , K do
Remove coordinate k by
Rij/k = Rij − aik bjk , ∀(i, j) ∈ S;
old_a:k = a:k , old_b:k = b:k ;max_func_decr = 0;
for iter=1, . . . , T do
foreach j = 1, . . . , N do
bjk := CCD(j, a:k , R:/k ) by (5) ;foreach i = 1, . . . , M do
aik := CCD(i, b:k , R:/k ) by (5) ;func_decr :=function decrease in this inner iteration;
max_func_decr :=max(max_func_decr, func_decr);
if func_decr < 10−8 max_func_decr then break;end
U = a:k − old_a:k , V = b:k − old_b:k ;Solve subspace search
minimizeα,β L(A + αU, B + βV );a:k := a:k + αU, b:k := b:k + βV ;Update and restore coordinate k by
Rij = Rij/k + aik bjk , ∀(i, j) ∈ S;
end
Algorithm 2: PolyMF-SS

Proof. Let Ai , B i be the solution at i-th outer iteration. Because the method is strictly decreasing, we haveλ(kAi k2 +kB i k2 ) ≤ L(Ai , B i ) ≤ L(A0 , B 0 ), ∀i = 0, 1, . . . .
(23)Thus, the sequence {(Ai , B i )} is bounded. By Weierstrass’s
sequential compactness theorem, every bounded sequencein a Euclidean space has an converging subsequence. Let(A∗ , B ∗ ) be any of the limit point of the converging subsequence.
Now we demonstrate how to use the optimality condition
for subproblem to construct a global condition for a stationary
point. For the PolyMF with coordinate descent direction,(A∗ , B ∗ ) is a fix point means that, taking derivative on (19),∇α F (a∗ik , b∗jk )
i/j

i/j

∗
∗
= a∗ik b∗2
jk + Rij/k bjk + C20 aik + C10 = 0, ∀(i, j) ∈ S.
(24)
i/j

i/j

Expanding the definition of C10 , C20 , and Rij/k , we obtainX
λa∗ik +
Riq b∗qk = 0, ∀i, k.
(25)
q|i,q∈S

By the same reasoning on ∇β F , there isX
λb∗jk +
Rpj a∗pk = 0, ∀j, k.

(26)

p|p,j∈S

Equations (25) and (26) means that the gradient of the fullproblem (2) vanishes on (A∗ , B ∗ ), because the definition of

Movielens10m
Netflix
Yahoo-MusicM
71,5672,649,429
1,000,990|S|
9,301,274
99,072,112
252,800,275

N
65,133
17,770624,961|Stest |698,780
1,408,395
4,003,960

k
5
10
10

λ
0.01
0.05
1

Table 1: The statistics and parameters for the datasets.
(b) Netflix(a) Movielens10m
10 10

10 5

10 0

10 -5

10 5

10 0

0

20

40

10 -5

60

10 15

LIBMF
CCD++
polyMF-SS

L − L∗LIBMF
CCD++
polyMF-SS

L − L∗

L − L∗10 10(c) Yahoo Music
LIBMF
CCD++
polyMF-SS

10 10

10 5

0

Time (s)

2000

4000

6000

10 0

0

Time (s)

2000

4000

6000

8000

Time (s)Figure 1: The comparison of L − L∗ in terms of time.
the gradient is∇ai L(A, B) = λai +X

Rij bj = 0, ∀i,(27)

Rij ai = 0, ∀j.

(28)

polyMF-CD This is the pair-wise coordinate descentmethod proposed in Algorithm 1, in which we perform an
optimal subspace search on each pair of coordinate (aik , bjk ).

j|(i,j)∈S∇bj L(A, B) = λbj +X
i|(i,j)∈S∗∗

Thus, (A , B ) is a stationary point. The same reasoning canbe also be applied to the proof for CCD++ steps.

4

Experiments

We evaluate the PolyMF algorithm (both the coordinate
descent and subspace search approaches) on several
benchmark problems in recommender systems. In particular,
we consider three datasets: Movielens10m, Netflix (Belland Koren 2007), and Yahoo Music (Dror et al. 2012),
from smaller to larger. See Table 1 for the dimensions
and parameters for each dataset. We first analyze theoptimization performance of our methods versus existingstate-of-the-art approaches in terms of minimizing the matrix
factorization objective, and we then discuss analyses ofdifferent additional aspects of the problem. We choose
regularization parameter λ as shown in Table 1, and focus
our analysis on the optimization performance compared to
the following algorithms:
LIBMF This package (Zhuang et al. 2013) implements a parallel SGD algorithm for matrix factorization withadaptive gradient steps (Duchi, Hazan, and Singer 2011). We
use its default parameters for learning rate and step-size because the package is already fine-tuned for the above datasets.
CCD++ This package (Yu et al. 2012) implements aparallel coordinate descent method for matrix factorization.
We use its default parameters because there is no need totune the step-size by the nature of coordinate descent.

polyMF-SS This is the subspace search method proposed in Algorithm 2, in which we perform an optimal
subspace-search on the CCD++ directions.
We conduct the experiments on a Intel Core i7-4790 machine with 32 GB memory. All experiments except the parallel speedup utilize 4 cores of the total 8 cores available.

4.1

Results on matrix factorization benchmarks

Figures 1 and 2 show the evolution of the objective function
L − L∗ (where L∗ is measured as the best objective foundby any of the methods) as a function of running time and iteration count. For all three data sets, the PolyMF-SS algorithm
outperforms LIBMF and CCD++ in terms of final objective
value and running time versus objective. For example, in theMovielens10m dataset, as shown in Figure 1(a), the PolyMFSS algorithm actually escapes the local optima that CCD++
converges to, and reaches a substantially lower final objective.
In all examples, the convergence to this objective is faster,with the largest gains present on the Netflix example. The
PolyMF-CD algorithm here is substantially slower by total
time (the time needed to perform a polynomial optimization
over each coordinate is prohibitively large), but we show
performance per iteration on the smaller Movielens10m data
set.

4.2

Comparative analysis

Finally, we analyze a few additional aspects of the algorithm,
such as the comparison between the SDP polynomial solver,
the data efficiency of PolyMF-CD and the parallel scaling
performance of the different algorithms.
Comparison of Kerner Method and SDP We test the twopolynomial solvers over 1 million iterations of the inner

(a) Movielens10m
10

10 5

10 0

10 -5

(c) Yahoo Music

10

10 5

10 0

0

50

100

10 -5

150

10

LIBMF
CCD++
polyMF-SS

L − L∗LIBMF
CCD++
polyMF-CD
polyMF-SS

L − L∗

L − L∗

10(b) Netflix

10

15

LIBMF
CCD++
polyMF-SS

10 10

10 5

0

100

Number of Iterations

200

300

Number of Iterations

400

10 0

0

50

100

150

200

250

Number of Iterations

Figure 2: The comparison of L − L∗ in terms of the number of iterations.
(a) Netflix(b) Yahoo Music

8

8
CCD++
polyMF-SS

7

SpeedupSpeedup

6

5
4
3
2
1

CCD++
polyMF-SS

7

6

5
4
3
2

2

4

6

Number of Cores

8

1

2

4

6

8

Number of Cores

Figure 3: Parallel Speedup

polynomial solve. The moment-based SDP approach takes
1.430 · 10−4 s ± 3.4611 · 10−6 s (average and standard deviation) to solve each problem, while the Durand-Kerner takes
2.0406 · 10−6 ± 5.1434 · 10−7 s, an improvement about 70
times.

we have shown that one can perform an exact subspace searchover joint update directions in both A and B terms, by solving (exactly) a quartic polynomial optimization problem. We
have discussed several approaches for applying these methods in either a direct coordinate descent fashion or in a subspace search fashion to speed up existing matrix factorization
approaches. We show that the resulting methods attain fasterobjective convergence, while also being able to escape local
optima that trap alternate approaches.
In future work, we aim to extend these approaches beyond the case of simple matrix factorization to richer settings
such as tensor factorization approaches (Kolda and Bader2009) or factorization machines (Rendle 2012). Both thesemethods have received significant attention in recent years,
yet they suffer from similar local optima problems as our
matrix factorization setting, and thus seem ripe for similar
approaches.

References
Data efficiency of PolyMF-CD The implementation of
CCD++ requires 4 data pass per iteration. Our subspace
search method, polyMF-SS, only requires 1 more data passthan CCD++ but saves more iterations. However, both methods require maintaining the transpose of the score matrix, or
more precisely, iterating over rows and columns of S. Thismight be cubersome when the data is too large to be store
in the main memory. On the other hand, the PolyMF-CD
algorithm requires similar number of data pass and do notneed to maintain the transpose matrix.
Parallel performance Finally, we emphasize that because
PolyMF-SS spends the bulk of its running time in computing
the search directions U and V (using any desired method,
though here we use CCD++), the PolyMF-SS inherits thecomputational advantages of the underlying search direction
algorithm. Thus, when applying our PolyMF-SS algorithm
to the Netflix and Yahoo Music domains, the speedup of themethods with additional cores is roughly identical to that of
CCD++, as shown in Figure 3 (results on Movielens10m are
not shown, as the dataset is too small to obtain substantial
speedup).

5

Conclusion and future work

In this paper, we have derived and presented a polynomial
optimization approach to matrix factorization. Specifically,

Bell, R. M., and Koren, Y. 2007. Lessons from the netflix
prize challenge. SIGKDD Explorations.
Blekherman, G.; Parrilo, P. A.; and Thomas, R. R. 2013.
Semidefinite optimization and convex algebraic geometry.
SIAM.
Dror, G.; Labs, Y.; Koenigstein, N.; Koren, Y.; and Weimer,
M. 2012. The yahoo! music dataset and kddcup’11. In
Proceedings of KDD Cup 2011 Competition.
Duchi, J.; Hazan, E.; and Singer, Y. 2011. Adaptive subgradient methods for online learning and stochastic optimization.
Journal of Machine Learning Research 12(Jul):2121–2159.
Durand, É. 1959. Solutions numeriques des equations algebriques. tome i: Equations du type f (x).
Gemulla, R.; Nijkamp, E.; Haas, P. J.; and Sismanis, Y. 2011.
Large-scale matrix factorization with distributed stochastic
gradient descent. In Proceedings of the 17th ACM SIGKDDInternational Conference on Knowledge Discovery and Data
Mining.
Kolda, T. G., and Bader, B. W. 2009. Tensor decompositions
and applications. SIAM review 51(3):455–500.
Koren, Y.; Bell, R.; and Volinsky, C. 2009. Matrix factorization techniques for recommender systems. Computer.
Lasserre, J. B. 2001. Global optimization with polynomials
and the problem of moments. SIAM Journal on Optimization.

Laurent, M. 2009. Sums of squares, moment matrices andoptimization over polynomials. In Emerging applications ofalgebraic geometry.
Parrilo, P. A. 2003. Semidefinite programming relaxations
for semialgebraic problems. Mathematical programming.
Petkovi, M. S.; Carstensen, C.; Trajkov, M.; et al. 1995.
Weierstrass formula and zero-finding methods. Numerische
Mathematik.
Pilászy, I.; Zibriczky, D.; and Tikk, D. 2010. Fast alsbased matrix factorization for explicit and implicit feedback
datasets. In Proceedings of the fourth ACM conference on
Recommender systems.
Recht, B., and Ré, C. 2013. Parallel stochastic gradientalgorithms for large-scale matrix completion. Mathematical
Programming Computation.
Rendle, S. 2012. Factorization machines with libfm. ACM
Transactions on Intelligent Systems and Technology (TIST)
3(3):57.
Yu, H.-F.; Hsieh, C.-J.; Si, S.; and Dhillon, I. S. 2012. Scalable coordinate descent approaches to parallel matrix factorization for recommender systems. In IEEE International
Conference of Data Mining.
Yu, H.-F.; Hsieh, C.-J.; Si, S.; and Dhillon, I. S. 2014. Parallel
matrix factorization for recommender systems. Knowledge
and Information Systems.
Zhou, Y.; Wilkinson, D.; Schreiber, R.; and Pan, R. 2008.
Large-scale parallel collaborative filtering for the netflix prize.
In Algorithmic Aspects in Information and Management.
Zhuang, Y.; Chin, W.-S.; Juan, Y.-C.; and Lin, C.-J. 2013. A
fast parallel sgd for matrix factorization in shared memory
systems. In Proceedings of the 7th ACM conference on
Recommender systems.

A

Proof on transformation of (7)In this section, we will show that any subspace search problem can be reduced to (7), with the coefficients of the α2and β 2 terms being positive. Consider subspace search on
arbitrary pair of direction U and V ,minimize L(A + αU, B + βV ).
(29)
α,β

Expand the definition of L in (2), we have1
L(A + αU, B + βV )
2
M
N
X
λ X
= (
kai + αui k2 +kbj + βvj k2 )
2 i=1
j=1
+

1 X
((ai + αui )> (bj + βvj ) − Sij )22(30)This way, the problem has been reduced to (7). Now we
examine the positivity of the coefficient for α̂2 . Observe that
P
PM
>
2
2
(i,j)∈S (ui bj ) + λ
i=1 kui k
2
PĈ20 − Ĉ21 =>2(i,j)∈S (ui vj )!2P>>(i,j)∈S (ui vj )(ui bj )
P
−
.
(34)
>
2
(i,j)∈S (ui vj )By Cauchy inequality, we have



X
X
2 
2

(u>
(u>
i bj )
i vj )(i,j)∈S(31)(i,j)∈S

1= C22 α2 β 2 + C21 α2 β + C12 αβ 2 + C11 αβ211+ C20 α2 + C10 α + C02 β 2 + C01 β,
2
2
where
X
>
2
C22 =
S (ui vj ) ,
X
>
>
C21 =
S (ui vj )(ui bj ),
X
>
>
>
C11 =
S (ui vj )Rij + (ui bj )(vj ai ),
X
>
>
C12 =
S (ui vj )(vj ai ),
C20 =

X

>
2
S (ui bj ) + λ

M
X

≥ 

(32)

X

C02 =

X

C01 =

X

>
S Rij (ui bj ) + λ

>
 .
(u>
i vj )(ui bj )

Thus,
PM

2
i=1 kui k> 0.
>2(i,j)∈S (ui vj )λ2Ĉ20 − Ĉ12
≥P

(36)

2
Similarly, the coefficient Ĉ02 − Ĉ12
for β̂ 2 is also positive.

B

Example of local minima

Consider a toy problem of 2 scores, {10, −10}, with A ∈R1×1 and B ∈ R2×1 .
L(A, B) = (a2 + b21 + b22 )/2 + (ab1 − 10)2 + (ab2 + 10)2 .
Fixing b2 = 1, we can see that the function plot for (a, b1 )contains exactly two local minima (Figure 4), and only oneof which is the global minimum.

i=1

C10 =

(35)

2
X

(i,j)∈S

kui k2 ,MX(i,j)∈Su>i ai ,i=1>
2
S (vj ai ) + λNX

kvj k2 ,j=1
NX

vj> bj .

400

(33)

j=1

P
P
Note that we denote
S for simplicity. If(i,j)∈S as
C22 = 0, we have u>
v
=
0
for
all
(i,
j) ∈ S. This case,i j
C21 = C12 = 0 and the minimization of L can be solved by alinear equation. Otherwise, we show that we can remove C21
and C12 by some transformation, and the resulting equation
will satisfy the constraint of (7).
Let α̂ = α + C12 , β̂ = β + C21 . Normalize every constantC·,· by C22 and denote them as Ĉ·,· , we have11
L = α̂2 β̂ 2 + (Ĉ11 − 2Ĉ21 C12 )α̂β̂
2C22
2
1
1
2
2
+ (Ĉ20 − Ĉ21
)α̂2 + (Ĉ02 − Ĉ12
)β̂ 2
2
2
2
+ (Ĉ10 + 2Ĉ12 Ĉ21
− Ĉ11 Ĉ21 − Ĉ12 Ĉ20 )α̂2+ (Ĉ01 + 2Ĉ21 Ĉ12
− Ĉ11 Ĉ12 − Ĉ21 Ĉ02 )β̂+ constant.

L>
S Rij (vj ai ) + λ

200
2
0

−2

0

2

−2

a

b1

Figure 4: plot of L(a, (b1 , 1))
Further, the difference between the two minima can bearbitrary large by scaling the scores. Thus, if not properly
initialized, both CCD++ and SGD may stuck at the local
minima, but our polyMF method will always converge toglobal optima because we did an exact subspace search. Thisis one scenario that our method performs better than the otherstate-of-the-art methods.