Polynomial Networks and Factorization Machines: New Insights and Efficient Training Algorithms

Introduction

Interactions between features play an important role in many classification and regression tasks. One of the simplest approach to leverage such interactions consists in explicitly augmenting feature vectors with products of features (monomials), as in polynomial regression. Although fast linear model solvers can be used (Chang et al., 2010; Sonnenburg & Franc, 2010), an obvious drawback of this kind of approach is that the number of parameters to estimate scales as $O(d^{m})$, where $d$ is the number of features and $m$ is the order of interactions considered. As a result, it is usually limited to second or third-order interactions.

Another popular approach consists in using a polynomial kernel so as to implicitly map the data via the kernel trick. The main advantage of this approach is that the number of parameters to estimate in the model is actually independent of $d$ and $m$. However, the cost of storing and evaluating the model is now proportional to the number of training instances. This is sometimes called the curse of kernelization (Wang et al., 2010). Common ways to address the issue include the Nyström method (Williams & Seeger, 2001), random features (Kar & Karnick, 2012) and sketching (Pham & Pagh, 2013; Avron et al., 2014).

Related work

2 Factorization machines

One of the simplest way to leverage feature interactions is polynomial regression (PR). For example, for second-order interactions, in this approach, we compute predictions by

Polynomial and ANOVA kernels

In this section, we show that the prediction functions used by polynomial networks and factorization machines can be written using (1) for a specific choice of kernel.

The polynomial kernel is a popular kernel for using combinations of features. The kernel is defined as

We thus see that $\mathcal{H}^{m}$ uses all monomials of degree $m$ (i.e., all combinations of features with replacement).

A much lesser known kernel is the ANOVA kernel (Stitson et al., 1997; Vapnik, 1998). Following (Shawe-Taylor & Cristianini, 2004, Section 9.2), the ANOVA kernel of degree $m$, where $2\leq m\leq d$, can be defined as

As a result, $\mathcal{A}^{m}$ uses only monomials composed of distinct features (i.e., feature combinations without replacement). For later convenience, we also define $\mathcal{A}^{0}(\boldsymbol{p},\boldsymbol{x})\coloneqq 1$ and $\mathcal{A}^{1}(\boldsymbol{p},\boldsymbol{x})\coloneqq\langle\boldsymbol{p},\boldsymbol{x}\rangle$.

With $\mathcal{H}^{m}$ and $\mathcal{A}^{m}$ defined, we are now in position to state the following lemma.

Let $\hat{y}_{\mathcal{K}}(\boldsymbol{x};\boldsymbol{\lambda},\boldsymbol{P})$ be defined as in (1). Then,

The relation easily extends to higher orders. This new view allows us to state results that will be very useful in the next sections. The first one is that $\mathcal{H}^{m}$ and $\mathcal{A}^{m}$ are homogeneous functions, i.e., they satisfy

Another key property of $\mathcal{A}^{m}(\boldsymbol{p},\boldsymbol{x})$ is multi-linearity.A function $f(\theta_{1},\dots,\theta_{k})$ is called multi-linear (resp. multi-convex) if it is linear (resp. convex) w.r.t. $\theta_{1},\dots,\theta_{k}$ separately.

Multi-linearity of $\mathcal{A}^{m}(\boldsymbol{p},\boldsymbol{x})$ w.r.t. $p_{1},\dots,p_{d}$

where $\boldsymbol{p}_{\neg j}$ denotes the $(d-1)$-dimensional vector with $p_{j}$ removed and similarly for $\boldsymbol{x}_{\neg j}$.

That is, everything else kept fixed, $\mathcal{A}^{m}(\boldsymbol{p},\boldsymbol{x})$ is an affine function of $p_{j}$, $\forall j\in[d]$. Proof is given in Appendix B.1.

Assuming $\boldsymbol{p}$ is dense and $\boldsymbol{x}$ sparse, the cost of naively computing $\mathcal{A}^{m}(\boldsymbol{p},\boldsymbol{x})$ by (8) is $O(n_{z}(\boldsymbol{x})^{m})$, where $n_{z}(\boldsymbol{x})$ is the number of non-zero features in $\boldsymbol{x}$. To address this issue, we will make use of the following lemma for computing $\mathcal{A}^{m}$ in nearly $O(mn_{z}(\boldsymbol{x}))$ time when $m\in\{2,3\}$.

where we defined $\mathcal{D}^{{m}}(\boldsymbol{p},\boldsymbol{x})\coloneqq\sum_{j=1}^{d}(p_{j}x_{j})^{m}$ and $\mathcal{D}^{{m,n}}(\boldsymbol{p},\boldsymbol{x})\coloneqq\mathcal{D}^{{m}}(\boldsymbol{p},\boldsymbol{x})\mathcal{D}^{{n}}(\boldsymbol{p},\boldsymbol{x})$.

Direct approach

Multi-convexity of (14) when $\mathcal{K}=\mathcal{A}^{m}$

$D_{\mathcal{A}^{m}}$ is convex in $\boldsymbol{\lambda}$ and in each row of $\boldsymbol{P}$ separately.

Proof is given in Appendix B.3. As a corollary, the objective function of FMs of arbitrary order is thus multi-convex. Theorem 1 suggests that we can minimize (14) efficiently when $\mathcal{K}=\mathcal{A}^{m}$ by solving a succession of convex problems w.r.t. $\boldsymbol{\lambda}$ and the rows of $\boldsymbol{P}$. We next show that when $m$ is odd, we can just fix $\boldsymbol{\lambda}=\boldsymbol{1}$ without loss of generality.

When is it useful to fit $\boldsymbol{\lambda}$?

Let $\mathcal{K}=\mathcal{H}^{m}\text{ or }\mathcal{A}^{m}$. Then

The result stems from the fact that $\mathcal{H}^{m}$ and $\mathcal{A}^{m}$ are homogeneous functions. If we define $\boldsymbol{v}\coloneqq\operatorname*{sign}(\lambda)\sqrt[m]{|\lambda|}\boldsymbol{p}$, then we obtain $\lambda\mathcal{H}^{m}(\boldsymbol{p},\boldsymbol{x})=\mathcal{H}^{m}(\boldsymbol{v},\boldsymbol{x})~{}\forall\lambda$ if $m$ is odd, and similarly for $\mathcal{A}^{m}$. That is, $\lambda$ can be absorbed into $\boldsymbol{v}$ without loss of generality. When $m$ is even, $\lambda<0$ cannot be absorbed unless we allow complex numbers. Because FMs fix $\boldsymbol{\lambda}=\boldsymbol{1}$, Lemma 4 shows that the class of functions that FMs can represent is possibly smaller than our framework.

Lifted approach

We begin by rewriting the kernel definitions using rank-one tensors. For $\mathcal{H}^{m}(\boldsymbol{p},\boldsymbol{x})$, it is easy to see that

For $\mathcal{A}^{m}(\boldsymbol{p},\boldsymbol{x})$, we need to ignore irrelevant monomials. For convenience, we introduce the following notation:

We can now concisely rewrite the ANOVA kernel as

Our key insight is described in the following lemma.

Link between tensors and kernel expansions

If $\boldsymbol{\mathcal{W}}$ is decomposed as in (20), then from Lemma 5, we obtain $L_{\mathcal{K}}(\boldsymbol{\mathcal{W}})=D_{\mathcal{K}}(\boldsymbol{\lambda},\boldsymbol{P})$ for $\mathcal{K}=\mathcal{H}^{m}$ or $\mathcal{A}^{m}$. This suggests that we can convert the problem of learning $\boldsymbol{\lambda}$ and $\boldsymbol{P}$ to that of learning a symmetric tensor $\boldsymbol{\mathcal{W}}$ of (symmetric) rank $k$. Thus, the problem of finding a small number of bases $\boldsymbol{p}_{1},\dots,\boldsymbol{p}_{k}$ and their associated weights $\lambda_{1},\dots,\lambda_{k}$ is converted to that of learning a low-rank symmetric tensor. Following (Candès et al., 2013), we call this approach lifted. Intuitively, we can think of $\boldsymbol{\mathcal{W}}$ as a tensor that contains the weights for predicting $y$ of monomials of degree $m$. For instance, when $m=3$, $\boldsymbol{\mathcal{W}}_{i,j,k}$ is the weight corresponding to the monomial $x_{i}x_{j}x_{k}$.

2 Multi-convex formulation

Due to multi-linearity of (25) w.r.t. $\boldsymbol{U}^{1},\dots,\boldsymbol{U}^{m}$, the objective function $L_{\mathcal{K}}$ is multi-convex in $\boldsymbol{U}^{1},\dots,\boldsymbol{U}^{m}$.

Computing predictions efficiently. When $\mathcal{K}=\mathcal{H}^{m}$, predictions are computed by $\langle\boldsymbol{\mathcal{W}},\boldsymbol{x}^{\otimes m}\rangle$. To compute them efficiently, we use the following lemma.

Symmetrization does not affect inner product

As a result, we never need to explicitly compute the symmetrized tensor. For the case $\mathcal{K}=\mathcal{A}^{2}$, cf. Appendix D.3.

Regularization

In some applications, the number of bases or the rank constraint are not enough for obtaining good generalization performance and it is necessary to consider additional form of regularization. For the lifted objective with $\mathcal{K}=\mathcal{H}^{2}$ or $\mathcal{A}^{2}$, we use the typical Frobenius-norm regularization

where $\beta>0$ is a regularization hyper-parameter. For the direct objective, we introduce the new regularization

This allows us to regularize $\boldsymbol{\lambda}$ and $\boldsymbol{P}$ with a single hyper-parameter. Let us define the following nuclear norm penalized objective:

We can show that (28), (29) and (30) are equivalent in the following sense.

Let $\mathcal{K}=\mathcal{H}^{2}$ or $\mathcal{A}^{2}$, then

Coordinate descent algorithms

Let us denote the elements of $\boldsymbol{P}$ by $p_{js}$. Then, our algorithm cyclically performs the following update for all $s\in[k]$ and $j\in[d]$:

The key challenge to use CD is computing $\frac{\partial\hat{y}_{i}}{\partial p_{js}}=\lambda_{s}\frac{\partial\mathcal{A}^{m}(\boldsymbol{p}_{s},\boldsymbol{x}_{i})}{\partial p_{js}}$ efficiently. Let us denote the elements of $\boldsymbol{X}$ by $x_{ji}$. Using Lemma 3, we obtain $\frac{\partial\mathcal{A}^{2}(\boldsymbol{p}_{s},\boldsymbol{x}_{i})}{\partial p_{js}}=\langle\boldsymbol{p}_{s},\boldsymbol{x}_{i}\rangle x_{ji}-p_{js}x_{ji}^{2}$ and $\frac{\partial\mathcal{A}^{3}(\boldsymbol{p}_{s},\boldsymbol{x}_{i})}{\partial p_{js}}=\mathcal{A}^{2}(\boldsymbol{p}_{s},\boldsymbol{x}_{i})x_{ji}-p_{js}x_{ji}^{2}\langle\boldsymbol{p}_{s},\boldsymbol{x}_{i}\rangle+p_{js}^{2}x_{ji}^{3}$. If for all $i\in[n]$ and for $s$ fixed, we maintain $\langle\boldsymbol{p}_{s},\boldsymbol{x}_{i}\rangle$ and $\mathcal{A}^{2}(\boldsymbol{p}_{s},\boldsymbol{x}_{i})$ (i.e., keep in sync after every update of $p_{js}$), then computing $\frac{\partial\hat{y}_{i}}{\partial p_{js}}$ takes $O(m)$ time. Hence the cost of one epoch, i.e. updating all elements of $\boldsymbol{P}$ once, is $O(mkn_{z}(\boldsymbol{X}))$. Complete details and pseudo code are given in Appendix D.1.

where $\eta\coloneqq\mu\sum_{i=1}^{n}\left(\frac{\partial\hat{y}_{i}}{\partial u_{js}^{t}}\right)^{2}+\beta$. The main difficulty is computing $\frac{\partial\hat{y}_{i}}{\partial u_{js}^{t}}=\prod_{t^{\prime}\neq t}\langle\boldsymbol{u}^{t^{\prime}}_{s},\boldsymbol{x}_{i}\rangle x_{ji}$ efficiently. If for all $i\in[n]$ and for $t$ and $s$ fixed, we maintain $\xi_{i}\coloneqq\prod_{t^{\prime}\neq t}\langle\boldsymbol{u}^{t^{\prime}}_{s},\boldsymbol{x}_{i}\rangle$, then the cost of computing $\frac{\partial\hat{y}_{i}}{\partial u_{js}^{t}}$ is $O(1)$. Hence the cost of one epoch is $O(mrn_{z}(\boldsymbol{X}))$, the same as SGD. Complete details are given in Appendix D.2.

Convergence. The above updates decrease the objective monotonically. Convergence to a stationary point is guaranteed following (Bertsekas, 1999, Proposition 2.7.1).

Inhomogeneous polynomial models

The algorithms presented so far are designed for homogeneous polynomial kernels $\mathcal{H}^{m}$ and $\mathcal{A}^{m}$. These kernels only use monomials of the same degree $m$. However, in many applications, we would like to use monomials of up to some degree. In this section, we propose a simple idea to do so using the algorithms presented so far, unmodified. Our key observation is that we can easily turn homogeneous polynomials into inhomogeneous ones by augmenting the dimensions of the training data with dummy features.

Next, we explain how to learn inhomogeneous polynomial models using $\mathcal{A}^{m}$. Using Lemma 2, we immediately obtain for $1\leq m\leq d$:

Experimental results

As explained in Section 4, there is no benefit to fitting $\boldsymbol{\lambda}$ when $m$ is odd, since $\mathcal{A}^{m}$ and $\mathcal{H}^{m}$ can absorb $\boldsymbol{\lambda}$ into $\boldsymbol{P}$. This is however not the case when $m$ is even: $\mathcal{A}^{m}$ and $\mathcal{H}^{m}$ can absorb absolute values but not negative signs (unless complex numbers are allowed for parameters). Therefore, when $m$ is even, the class of functions we can represent with models of the form (1) is possibly smaller if we fix $\boldsymbol{\lambda}=\boldsymbol{1}$ (as done in FMs).

To check that this is indeed the case, on the diabetes dataset, we minimized (29) with $m=2$ as follows:

minimize w.r.t. both $\boldsymbol{\lambda}$ and $\boldsymbol{P}$ alternatingly,

fix $\lambda_{s}=1$ for $s\in[k]$ and minimize w.r.t. $\boldsymbol{P}$,

fix $\lambda_{s}=\pm 1$ with proba. $0.5$ and minimize w.r.t. $\boldsymbol{P}$.

We initialized elements of $\boldsymbol{P}$ by $p_{js}\sim\mathcal{N}(0,0.01)$ for all $j\in[d]$, $s\in[k]$. Our results are shown in Figure 1. For $\mathcal{K}=\mathcal{A}^{2}$, we use CD and for $\mathcal{K}=\mathcal{H}^{2}$, we use L-BFGS. Note that since (29) is convex w.r.t. $\boldsymbol{\lambda}$, a) is insensitive to the initialization of $\boldsymbol{\lambda}$ as long as we fit $\boldsymbol{\lambda}$ before $\boldsymbol{P}$. Not surprisingly, fitting $\boldsymbol{\lambda}$ allows us to achieve a smaller objective value. This is especially apparent when $\mathcal{K}=\mathcal{H}^{2}$. However, the difference is much smaller when $\mathcal{K}=\mathcal{A}^{2}$. We give intuitions as to why this is the case in Section 10.

We emphasize that this experiment was designed to confirm that fitting $\boldsymbol{\lambda}$ does indeed improve representation power of the model when $m$ is even. In practice, it is possible that fixing $\boldsymbol{\lambda}=\boldsymbol{1}$ reduces overfitting and thus improves generalization error. However, this highly depends on the data.

2 Direct vs. lifted optimization

3 Recommender system experiment

To confirm the ability of the proposed framework to infer the weights of unobserved feature interactions, we conducted experiments on Last.fm and Movielens 1M, two standard recommender system datasets. Following (Rendle, 2012), matrix factorization can be reduced to FMs by creating a dataset of $(\boldsymbol{x}_{i},y_{i})$ pairs where $\boldsymbol{x}_{i}$ contains the one-hot encoding of the user and item and $y_{i}$ is the corresponding rating (i.e., number of training instances equals number of ratings). We compared four models:

$\mathcal{K}=\mathcal{A}^{2}$ (linear combination): $\hat{y}=\langle\boldsymbol{w},\boldsymbol{x}\rangle+\hat{y}_{\mathcal{A}^{2}}(\boldsymbol{x})$,

$\mathcal{K}=\mathcal{H}^{2}$ (linear combination): $\hat{y}=\langle\boldsymbol{w},\boldsymbol{x}\rangle+\hat{y}_{\mathcal{H}^{2}}(\boldsymbol{x})$,

4 Low-budget non-linear regression experiment

In this experiment, we demonstrate the ability of the proposed framework to reach good regression performance with a small number of bases $k$. We compared:

Proposed with $\mathcal{K}=\mathcal{H}^{3}$ (with augmented features),

Proposed with $\mathcal{K}=\mathcal{A}^{3}$ (with augmented features),

Nyström method with $\mathcal{K}=\mathcal{P}^{3}_{1}$ and

Random Selection: choose $\boldsymbol{p}_{1},\dots,\boldsymbol{p}_{k}$ uniformly at random from training set and use $\mathcal{K}=\mathcal{P}^{3}_{1}$.

For a) and b) we used the lifted approach. For fair comparison in terms of model size (number of floats used), we set $r=k/3$. Results on the abalone, cadata and cpusmall datasets are shown in Figure 4. We see that i) the proposed framework reaches the same performance as kernel ridge regression with much fewer bases than other methods and ii) $\mathcal{H}^{3}$ tends to outperform $\mathcal{A}^{3}$ on these tasks. Similar trends were observed when using $\mathcal{K}=\mathcal{H}^{2}$ or $\mathcal{A}^{2}$.

Discussion

Ability to infer weights of unobserved interactions. In our view, one of the strengths of PNs and FMs is their ability to infer the weights of unobserved feature interactions, unlike traditional kernel methods. To see why, recall that in kernel methods, predictions are computed by $\hat{y}=\sum_{i=1}^{n}\alpha_{i}\mathcal{K}(\boldsymbol{x}_{i},\boldsymbol{x})$. When $\mathcal{K}=\mathcal{H}^{m}$ or $\mathcal{A}^{m}$, by Lemma 5, this is equivalent to $\hat{y}=\langle\widetilde{\boldsymbol{\mathcal{W}}},\boldsymbol{x}^{\otimes m}\rangle$ or $\langle\widetilde{\boldsymbol{\mathcal{W}}},\boldsymbol{x}^{\otimes m}\rangle_{>}$ if we set $\widetilde{\boldsymbol{\mathcal{W}}}\coloneqq\sum_{i=1}^{n}\alpha_{i}\boldsymbol{x}_{i}^{\otimes m}$. Thus, in kernel methods, the weight associated with $x_{j_{1}}\dots x_{j_{m}}$ can be written as a linear combination of the training data’s monomials:

Assuming binary features, the weights of monomials that were never observed in the training set are zero. In contrast, in PNs and FMs, we have $\boldsymbol{\mathcal{W}}=\sum_{s=1}^{k}\lambda_{i}\boldsymbol{p}_{s}^{\otimes m}$ and therefore the weight associated with $x_{j_{1}}\dots x_{j_{m}}$ becomes

Because parameters are shared across monomials, PNs and FMs are able to interpolate the weights of monomials that were never observed in the training set. This is the key property which makes it possible to use them on recommender system tasks. In future work, we plan to apply PNs and FMs to biological data, where this property should be very useful, e.g., for inferring higher-order interactions between genes.

and therefore the model is unable to predict negative values.

Empirically, we showed in Section 9.4 that $\mathcal{H}^{m}$ outperforms $\mathcal{A}^{m}$ for low-budget non-linear regression. In contrast, we showed in Section 9.3 that $\mathcal{A}^{m}$ outperforms $\mathcal{H}^{m}$ for recommender systems. The main difference between the two experiments is the nature of the features used: continuous for the former and binary for the latter. For binary features, squared features $x_{1}^{2},\dots,x_{d}^{2}$ are redundant with $x_{1},\dots,x_{d}$ and are therefore not expected to help improve accuracy. On the contrary, they might introduce bias towards first-order features. We hypothesize that the ANOVA kernel is in general a better choice for binary features, although this needs to be verified by more experiments, for instance on natural language processing (NLP) tasks.

Conclusion

In this paper, we revisited polynomial networks (Livni et al., 2014) and factorization machines (Rendle, 2010, 2012) from a unified perspective. We proposed direct and lifted optimization approaches and showed their equivalence in the regularized case for $m=2$. With respect to PNs, we proposed the first CD solver with support for arbitrary integer $m\geq 2$. With respect to FMs, we made several novel contributions including making a connection with the ANOVA kernel, proving important properties of the objective function and deriving the first CD solver for third-order FMs. Empirically, we showed that the proposed algorithms achieve excellent performance on non-linear regression and recommender system tasks.

Acknowledgments

This work was partially conducted as part of “Research and Development on Fundamental and Applied Technologies for Social Big Data”, commissioned by the National Institute of Information and Communications Technology (NICT), Japan. We also thank Vlad Niculae, Olivier Grisel, Fabian Pedregosa and Joseph Salmon for their valuable comments.

References

Appendix A Symmetric tensors

In other words $\boldsymbol{\mathcal{M}}_{\boldsymbol{\sigma}}$ is a copy of $\boldsymbol{\mathcal{M}}$ with its axes permuted. This generalizes the concept of transpose to tensors.

where $k$ is called the symmetric rank of $\boldsymbol{\mathcal{W}}$. This generalizes the concept of eigendecomposition to tensors. These two concepts are illustrated in Figure 5.

A.2 Proof of Lemma 26

Appendix B Proofs related to ANOVA kernels

where we used $\mathcal{A}^{0}(\boldsymbol{p},\boldsymbol{x})=1$.

For $1<m\leq d$, first notice that we can rewrite (8) as

We can always permute the elements of $\boldsymbol{p}$ and $\boldsymbol{x}$ without changing $\mathcal{A}^{m}(\boldsymbol{p},\boldsymbol{x})$. It follows that

B.2 Efficient computation when m∈{2,3}𝑚23m\in\{2,3\}

Using the multinomial theorem, we can expand the homogeneous polynomial kernel as

To simplify notation, we define the shorthands

For $m=2$, the possible monomials are of the form $\rho_{j}^{2}$ for all $j$ and $\rho_{i}\rho_{j}$ for $j>i$. Applying (60), we obtain

This formula was already mentioned in (Stitson et al., 1997). It was also rediscovered in (Rendle, 2010, 2012), although the connection with the ANOVA kernel was not identified.

For $m=3$, the possible monomials are of the form $\rho_{j}^{3}$ for all $j$, $\rho_{i}\rho_{j}^{2}$ for $i\neq j$ and $\rho_{i}\rho_{j}\rho_{k}$ for $k>j>i$. Applying (60), we obtain

We can compute the second term efficiently by using

B.3 Proof of multi-convexity (Theorem 1)

Hence $\hat{y}_{\mathcal{A}^{m}}(\boldsymbol{x};\boldsymbol{\lambda},\boldsymbol{P})$ is an affine function of $\boldsymbol{\bar{p}}_{1},\dots,\boldsymbol{\bar{p}}_{d}$. The composition of a convex loss function and an affine function is convex. Therefore, (14) is convex in $\boldsymbol{\bar{p}}_{j}~{}\forall j\in[d]$. Convexity w.r.t. $\boldsymbol{\lambda}$ is obvious.

Appendix C Proof of equivalence between regularized problems (Theorem 2)

First, we are going to prove that the optimal solution of the nuclear norm penalized problem is a symmetric matrix. For that, we need the following lemma.

Upper-bound on nuclear norm of symmetrized matrix

Symmetry of optimal solution of nuclear norm penalized problem

Next, we recall the variational formulation of the nuclear norm based on the SVD.

Variational formulation of nuclear norm based on SVD

For a proof, see for instance (Mazumder et al., 2010, Section A.5).

Now, we give a specialization of the above for symmetric matrices, based on the eigendecomposition instead of SVD.

Variational formulation of nuclear norm based on eigendecomposition

Now, computing $\frac{1}{2}(\sum_{s}\|\boldsymbol{u}_{s}\|^{2}+\|\boldsymbol{v}_{s}\|^{2})$ gives $\sum_{s=1}^{k}|\lambda_{s}|~{}\|\boldsymbol{p}_{s}\|^{2}$. The minimum value $\|\boldsymbol{M}\|_{*}=\|\boldsymbol{\lambda}\|_{1}$ follows from the fact that $\boldsymbol{P}$ is orthonormal and hence $\|\boldsymbol{p}_{s}\|^{2}=1~{}\forall s\in[k]$. $\square$

We now have all the tools to prove our result. The equivalence between (28) and (30) when $r=\operatorname*{rank}(\boldsymbol{M}^{*})$ is a special case of (Mazumder et al., 2010, Theorem 3). From Lemma 74, we know that the optimal solution of (30) is symmetric. This allows us to substitute (75) with (76), and therefore, (29) is equivalent to (30) with $k=\operatorname*{rank}(\boldsymbol{M}^{*})$. As discussed in (Mazumder et al., 2010), the result also holds when $r=k$ is larger than $\operatorname*{rank}(\boldsymbol{M}^{*})$.

Appendix D Efficient coordinate descent algorithms

The fact that the second derivative is null is a consequence of the multi-linearity of $\mathcal{A}^{m}$.

For an efficient implementation, we need to maintain $\hat{y}_{i}~{}\forall i\in[n]$ and statistics that depend on $\boldsymbol{p}_{s}$. For the former, we need $O(n)$ memory. For the latter, we need $O(kmn)$ memory for an implementation with full cache. However, this requirement is not realistic for a large training set. In practice, the memory requirement can be reduced to $O(mn)$ if we recompute the quantities then sweep through $p_{1s},\dots,p_{ds}$ for $s$ fixed. Overall the cost of one epoch is $O(kn_{z}(\boldsymbol{X}))$. A similar implementation technique is described for factorization machines with $m=2$ in (Rendle, 2012).

The first and second coordinate-wise derivatives are given by

We consider the following regularized objective function

For an efficient implementation, the two quantities we need to maintain are $\hat{y}_{i}~{}\forall i\in[n]$ and $\prod_{t^{\prime}\neq t}\langle\boldsymbol{u}^{t^{\prime}}_{s},\boldsymbol{x}_{i}\rangle~{}\forall i\in[n],s\in[r],t\in[m]$. For the former, we need $O(n)$ memory. For the latter, we need $O(rmn)$ memory for an implementation with full cache. However, this requirement is not realistic for a large training set. In practice, the memory requirement can be reduced to $O(mn)$ if we recompute the quantity then sweep through $u^{t}_{1s},\dots,u^{t}_{ds}$ for $t$ and $s$ fixed. Overall the cost of one epoch is $O(mrn_{z}(\boldsymbol{X}))$.

For $\langle\cdot,\cdot\rangle_{>}$, efficient computations are more involved since we need to ignore irrelevant monomials. Nevertheless, we can also compute the predictions directly without explicitly symmetrizing the model. For $m=2$, it suffices to subtract the effect of squared features. It is easy to verify that we then obtain

where $\circ$ indicates element-wise product. The coordinate-wise derivatives are given by

Generalizing this to arbitrary $m$ is a future work.

Appendix E Datasets

For regression experiments, we used the following public datasets.

The diabetes dataset is available in scikit-learn (Pedregosa et al., 2011). Other datasets are available from http://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/.

For recommender system experiments, we used the following two public datasets.

The design matrix $\boldsymbol{X}$ was constructed following (Rendle, 2010, 2012). Namely, for each rating $y_{i}$, the corresponding $\boldsymbol{x}_{i}$ is set to the concatenation of the one-hot encodings of the user and item indices. Hence the number of samples $n$ is the number of ratings and the number of features is equal to the sum of the number of users and items. Each sample contains exactly two non-zero features. It is known that factorization machines are equivalent to matrix factorization when using this representation (Rendle, 2010, 2012).

We split samples uniformly at random between 75% for training and 25% for testing.