Higher-Order Factorization Machines

1 Introduction

Factorization machines (FMs) fm ; libfm are a supervised learning approach that can use second-order feature combinations efficiently even when the data is very high-dimensional. The key idea of FMs is to model the weights of feature combinations using a low-rank matrix. This has two main benefits. First, FMs can achieve empirical accuracy on a par with polynomial regression or kernel methods but with smaller and faster to evaluate models fm_icml . Second, FMs can infer the weights of feature combinations that were not observed in the training set. This second property is crucial for instance in recommender systems, a domain where FMs have become increasingly popular libfm ; fm_context . Without the low-rank property, FMs would fail to generalize to unseen user-item interactions.

Unfortunately, although higher-order FMs (HOFMs) were briefly mentioned in the original work of fm ; libfm , there exists to date no efficient algorithm for training arbitrary-order HOFMs. In fact, even just computing predictions given the model parameters naively takes polynomial time in the number of features. For this reason, HOFMs have, to our knowledge, never been applied to any problem. In addition, HOFMs, as originally defined in fm ; libfm

, model each degree in the polynomial expansion with a different matrix and therefore require the estimation of a large number of parameters.

In this paper, we propose the first efficient algorithms for training arbitrary-order HOFMs. To do so, we rely on a link between FMs and the so-called ANOVA kernel fm_icml . We propose linear-time dynamic programming algorithms for evaluating the ANOVA kernel and computing its gradient. Based on these, we propose stochastic gradient and coordinate descent algorithms for arbitrary-order HOFMs. To reduce the number of parameters, as well as prediction times, we also introduce two new kernels derived from the ANOVA kernel, allowing us to define new variants of HOFMs with shared parameters. We demonstrate the proposed approaches on four different link prediction tasks.

2 Factorization machines (FMs)

Second-order FMs. Factorization machines (FMs) fm ; libfm are an increasingly popular method for efficiently using second-order feature combinations in classification or regression tasks even when the data is very high-dimensional. Let $w \in R^{d}$ and $P \in R^{d \times k}$ , where $k \in N$ is a rank hyper-parameter. We denote the rows of $P$ by ${¯ p}_{j}$ and its columns by $p_{s}$ , for $j \in [d]$ and $s \in [k]$ , where $[d]\coloneqq{1,…,d}$ . Then, FMs predict an output $y \in R$

from a vector

x = [x_{1}, \dots, x_{d}]^{T}

^yFM(x)\coloneqq⟨w,x⟩+∑j′>j⟨¯pj,¯pj′⟩xjxj′.

(1)

An important characteristic of (1) is that it considers only combinations of distinct features (i.e., the squared features $x_{1}^{2}, \dots, x_{d}^{2}$ are ignored). The main advantage of FMs compared to naive polynomial regression is that the number of parameters to estimate is $O (d k)$ instead of $O (d^{2})$ . In addition, we can compute predictions in $O (2 d k)$ time¹¹1We include the constant factor for fair later comparison with arbitrary-order HOFMs. using

{^y}_{FM} (x) = w^{T} x + \frac{1}{2} (∥ P^{T} x ∥^{2} - k \sum s = 1 ∥ p_{s} \circ x ∥^{2}),

(2)

where $\circ$ indicates element-wise product convex_fm . Given a training set $X = [x_{1}, \dots, x_{n}] \in R^{d \times n}$ and $y = [y_{1}, \dots, y_{n}]^{T} \in R^{n}$ , $w$ and $P$ can be learned by minimizing the following non-convex objective

\frac{1}{n} n \sum i = 1 ℓ (y_{i}, {^y}_{FM} (x_{i})) + \frac{β_{1}}{2} ∥ w ∥^{2} + \frac{β_{2}}{2} ∥ P ∥^{2},

(3)

where $ℓ$

is a convex loss function and

β_{1} > 0, β_{2} > 0

are hyper-parameters. The popular libfm library libfm implements efficient stochastic gradient and coordinate descent algorithms for obtaining a stationary point of (3). Both algorithms have a runtime complexity of

O (2 d k n)

per epoch.

Higher-order FMs (HOFMs). Although no training algorithm was provided, FMs were extended to higher-order feature combinations in the original work of fm ; libfm . Let $P^{(t)} \in R^{d \times k_{t}}$ , where $t \in {2, \dots, m}$ is the order or degree of feature combinations considered, and $k_{t} \in N$ is a rank hyper-parameter. Let ${¯ p}_{j}^{(t)}$ be the $j^{th}$ row of $P^{(t)}$ . Then $m$ -order HOFMs can be defined as

^yHOFM(x)\coloneqq⟨w,x⟩+∑j′>j⟨¯p(2)j,¯p(2)j′⟩xjxj′+⋯+∑jm>⋯>j1⟨¯p(m)j1,…,¯p(m)jm⟩xj1xj2…xjm

(4)

where we defined $⟨¯p(t)j1,…,¯p(t)jt⟩\coloneqqsum(¯p(t)j1∘⋯∘¯p(t)jt)$ (sum of element-wise products). The objective function of HOFMs can be expressed in a similar way as for (3):

\frac{1}{n} n \sum i = 1 ℓ (y_{i}, {^y}_{HOFM} (x_{i})) + \frac{β_{1}}{2} ∥ w ∥^{2} + m \sum t = 2 \frac{β_{t}}{2} ∥ P^{(t)} ∥^{2},

(5)

where $β_{1}, \dots, β_{m} > 0$ are hyper-parameters. To avoid the combinatorial explosion of hyper-parameter combinations to search, in our experiments we will simply set $β_{1} = \dots = β_{m}$ and $k_{2} = \dots = k_{m}$ . While (4) looks quite daunting, fm_icml recently showed that FMs can be expressed from a simpler kernel perspective. Let us define the ANOVA²²2The name comes from the ANOVA decomposition of functions. wahba_book ; vapnik_book kernel vapnik_book of degree $2 \leq m \leq d$ by

Am(p,x)\coloneqq∑jm>⋯>j1m∏t=1pjtxjt.

(6)

For later convenience, we also define $A0(p,x)\coloneqq1$ and $A1(p,x)\coloneqq⟨p,x⟩$ . Then it is shown that

{^y}_{HOFM} (x) = ⟨ w, x ⟩ + k_{2} \sum s = 1 A^{2} (p_{s}^{(2)}, x) + \dots + k_{m} \sum s = 1 A^{m} (p_{s}^{(m)}, x),

(7)

where $p_{s}^{(t)}$ is the $s^{th}$ column of $P^{(t)}$ . This perspective shows that we can view FMs and HOFMs as a type of kernel machine whose “support vectors” are learned directly from data. Intuitively, the ANOVA kernel can be thought as a kind of polynomial kernel that uses feature combinations without replacement (i.e., of distinct features). A key property of the ANOVA kernel is multi-linearity fm_icml :

A^{m} (p, x) = A^{m} (p_{\neg j}, x_{\neg j}) + p_{j} x_{j} A^{m - 1} (p_{\neg j}, x_{\neg j}),

(8)

where $p_{\neg j}$ denotes the $(d - 1)$ -dimensional vector with $p_{j}$ removed and similarly for $x_{\neg j}$ . That is, everything else kept fixed, $A^{m} (p, x)$ is an affine function of $p_{j} \forall j \in [d]$ . Although no training algorithm was provided, fm_icml showed based on (8) that, although non-convex, the objective function of arbitrary-order HOFMs is convex in $w$ and in each row of $P^{(2)}, \dots, P^{(m)}$ , separately.

Interpretability of HOFMs. An advantage of FMs and HOFMs is their interpretability. To see why this is the case, notice that we can rewrite (4) as

{^y}_{HOFM} (x) = ⟨ w, x ⟩ + \sum j^{'} > j W_{j, j^{'}}^{(2)} x_{j} x_{j^{'}} + \dots + \sum j_{m} > \dots > j_{1} W_{j_{1}, \dots, j_{m}}^{(m)} x_{j_{1}} x_{j_{2}} \dots x_{j_{m}},

(9)

where we defined $W(t)\coloneqq∑kts=1p(t)s⊗⋯⊗p(t)st times$ . Intuitively, $W^{(t)} \in R^{d^{t}}$ is a low-rank $t$

-way tensor which contains the weights of feature combinations of degree

t

. For instance, when

t = 3

W_{i, j, k}^{(3)}

is the weight of

x_{i} x_{j} x_{k}

. Similarly to the ANOVA decomposition of functions, HOFMs consider only combinations of distinct features (i.e.,

x_{j_{1}} x_{j_{2}} \dots x_{j_{m}}

for

j_{m} > \dots > j_{2} > j_{1}

This paper. Unfortunately, there exists to date no efficient algorithm for training arbitrary-order HOFMs. Indeed, computing (6) naively takes $O (d^{m})$ , i.e., polynomial time. In the following, we present linear-time algorithms. Moreover, HOFMs, as originally defined in fm ; libfm require the estimation of $m - 1$ matrices $P^{(2)}, \dots, P^{(m)}$ . Thus, HOFMs can produce large models when $m$ is large. To address this issue, we propose new variants of HOFMs with shared parameters.

3 Linear-time stochastic gradient algorithms for HOFMs

The kernel view presented in Section 2 allows us to focus on the ANOVA kernel as the main “computational unit” for training HOFMs. In this section, we develop dynamic programming (DP) algorithms for evaluating the ANOVA kernel and computing its gradient in only $O (d m)$ time.

Evaluation. The main observation (see also (kernel_book, , Section 9.2)) is that we can use (8) to recursively remove features until computing the kernel becomes trivial. Let us denote a subvector of $p$ by $p_{1 : j} \in R^{j}$ and similarly for $x$ . Let us introduce the shorthand $aj,t\coloneqqAt(p1:j,x1:j)$ . Then, from (8),

a_{j, t} = a_{j - 1, t} + p_{j} x_{j} a_{j - 1, t - 1} \forall d \geq j \geq t \geq 1.

(10)

For convenience, we also define $a_{j, 0} = 1 \forall j \geq 0$ since $A^{0} (p, x) = 1$ and $a_{j, t} = 0 \forall j < t$ since there does not exist any $t$ -combination of features in a $j < t$ dimensional vector.

	$j = 0$	$j = 1$	$j = 2$	$\dots$	$j = d$
$t = 0$	1	$1$	$1$	$1$	$1$
$t = 1$	0	$a_{1, 1}$	$a_{2, 1}$	$\dots$	$a_{d, 1}$
$t = 2$	0	$0$	$a_{2, 2}$	$\dots$	$a_{d, 2}$
$⋮$	$⋮$	$⋮$	$⋮$	$⋱$	$⋮$
$t = m$	0	$0$	$0$	$\dots$	$a_{d, m}$

Table 1: Example of DP table

The quantity we want to compute is $A^{m} (p, x) = a_{d, m}$ . Instead of naively using recursion (10), which would lead to many redundant computations, we use a bottom-up approach and organize computations in a DP table. We start from the top-left corner to initialize the recursion and go through the table to arrive at the solution in the bottom-right corner. The procedure, summarized in Algorithm 1, takes $O (d m)$ time and memory.

Gradients. For computing the gradient of $A^{m} (p, x)$ w.r.t. $p$ , we use reverse-mode differentiation autodiff

(a.k.a. backpropagation in a neural network context), since it allows us to compute the entire gradient in a single pass. We supplement each variable

a_{j, t}

in the DP table by a so-called adjoint

~aj,t\coloneqq∂ad,m∂aj,t

, which represents the sensitivity of

a_{d, m} = A^{m} (p, x)

w.r.t.

a_{j, t}

. From recursion (10), except for edge cases,

a_{j, t}

influences

a_{j + 1, t + 1}

and

a_{j + 1, t}

. Using the chain rule, we then obtain

{~ a}_{j, t} = \frac{\partial a_{d, m}}{\partial a_{j + 1, t}} \frac{\partial a_{j + 1, t}}{\partial a_{j, t}} + \frac{\partial a_{d, m}}{\partial a_{j + 1, t + 1}} \frac{\partial a_{j + 1, t + 1}}{\partial a_{j, t}} = {~ a}_{j + 1, t} + p_{j + 1} x_{j + 1} {~ a}_{j + 1, t + 1} \forall d - 1 \geq j \geq t \geq 1.

(11)

Similarly, we introduce the adjoint $~pj\coloneqq∂ad,m∂pj ∀j∈[d]$ . Since $p_{j}$ influences $a_{j, t} \forall t \in [m]$ , we have

{~ p}_{j} = m \sum t = 1 \frac{\partial a_{d, m}}{\partial a_{j, t}} \frac{\partial a_{j, t}}{\partial p_{j}} = m \sum t = 1 {~ a}_{j, t} a_{j - 1, t - 1} x_{j} .

(12)

We can run recursion (11) in reverse order of the DP table starting from ${~ a}_{d, m} = \frac{\partial a_{d, m}}{\partial a_{d, m}} = 1$ . Using this approach, we can compute the entire gradient $\nabla A^{m} (p, x) = [{~ p}_{1}, \dots, {~ p}_{d}]^{T}$ w.r.t. $p$ in $O (d m)$ time and memory. The procedure is summarized in Algorithm 2.

Input:

p \in R^{d}

x \in R^{d}

a_{j, t} \leftarrow 0 \forall t \in [m], j \in [d] \cup {0}

a_{j, 0} \leftarrow 1 \forall j \in [d] \cup {0}

for

t\coloneqq1,…,m

for

j\coloneqqt,…,d

a_{j, t} \leftarrow a_{j - 1, t} + p_{j} x_{j} a_{j - 1, t - 1}

end for

Output:

A^{m} (p, x) = a_{d, m}

Algorithm 1 Evaluating

A^{m} (p, x)

O (d m)

Input:

p \in R^{d}

x \in R^{d}

{a_{j, t}}_{j, t = 0}^{d, m}

{~ a}_{j, t} \leftarrow 0 \forall t \in [m + 1], j \in [d]

{~ a}_{d, m} \leftarrow 1

for

t\coloneqqm,…,1

for

j\coloneqqd−1,…,t

{~ a}_{j, t} \leftarrow {~ a}_{j + 1, t} + {~ a}_{j + 1, t + 1} p_{j + 1} x_{j + 1}

end for

~pj\coloneqq∑mt=1~aj,taj−1,t−1xj ∀j∈[d]

Output:

\nabla A^{m} (p, x) = [{~ p}_{1}, \dots, {~ p}_{d}]^{T}

Algorithm 2 Computing

\nabla A^{m} (p, x)

O (d m)

Stochastic gradient (SG) algorithms. Based on Algorithm 1 and 2, we can easily learn arbitrary-order HOFMs using any gradient-based optimization algorithm. Here we focus our discussion on SG algorithms. If we alternatingly minimize (5) w.r.t $P^{(2)}, \dots, P^{(m)}$ , then the sub-problem associated with degree $m$ is of the form

F(P)\coloneqq1nn∑i=1ℓ(yi,k∑s=1Am(ps,xi)+oi)+β2∥P∥2,

(13)

where $o_{1}, \dots, o_{n} \in R$ are fixed offsets which account for the contribution of degrees other than $m$ to the predictions. The sub-problem is convex in each row of $P$ fm_icml . A SG update for (13) w.r.t. $p_{s}$ for some instance $x_{i}$ can be computed by $p_{s} \leftarrow p_{s} - η ℓ^{'} (y_{i}, {^y}_{i}) \nabla A^{m} (p_{s}, x_{i}) - η β p_{s}$ , where $η$ is a learning rate and where we defined $^yi\coloneqq∑ks=1Am(ps,xi)+oi$ . Because evaluating $A^{m} (p, x)$ and computing its gradient both take $O (d m)$ , the cost per epoch, i.e., of visiting all instances, is $O (m d k n)$ . When $m = 2$ , this is the same cost as the SG algorithm implemented in libfm.

Sparse data. We conclude this section with a few useful remarks on sparse data. Let us denote the support of a vector $x = [x_{1}, \dots, x_{d}]^{T}$ by $supp(x)\coloneqq{j∈[d]:xj≠0}$ and let us define $xS\coloneqq[xj:j∈S]T$ . It is easy to see from (8) that the gradient and $x$ have the same support, i.e., $supp (\nabla A^{m} (p, x)) = supp (x)$ . Another useful remark is that $A^{m} (p, x) = A^{m} (p_{supp (x)}, x_{supp (x)})$ , provided that $m \leq n_{z} (x)$ , where $n_{z} (x)$ is the number of non-zero elements in $x$ . Hence, when the data is sparse, we only need to iterate over non-zero features in Algorithm 1 and 2. Consequently, their time and memory cost is only $O (n_{z} (x) m)$ and thus the cost per epoch of SG algorithms is $O (m k n_{z} (X))$ .

4 Coordinate descent algorithm for arbitrary-order HOFMs

We now describe a coordinate descent (CD) solver for arbitrary-order HOFMs. CD is a good choice for learning HOFMs because their objective function is coordinate-wise convex, thanks to the multi-linearity of the ANOVA kernel fm_icml . Our algorithm can be seen as a generalization to higher orders of the CD algorithms proposed in libfm ; fm_icml .

An alternative recursion. Efficient CD implementations typically require maintaining statistics for each training instance, such as the predictions at the current iteration. When a coordinate is updated, the statistics then need to be synchronized. Unfortunately, the recursion we used in the previous section is not suitable for a CD algorithm because it would require to store and synchronize the DP table for each training instance upon coordinate-wise updates. We therefore turn to an alternative recursion:

A^{m} (p, x) = \frac{1}{m} m \sum t = 1 (- 1)^{t + 1} A^{m - t} (p, x) D^{t} (p, x),

(14)

where we defined $Dt(p,x)\coloneqq∑dj=1(pjxj)t$ . Note that the recursion was already known in the context of traditional kernel methods (c.f., (vapnik_book, , Section 11.8)) but its application to HOFMs is novel. Since we know that $A^{0} (p, x) = 1$ and $A^{1} (p, x) = ⟨ p, x ⟩$ , we can use (14) to compute $A^{2} (p, x)$ , then $A^{3} (p, x)$ , and so on. The overall evaluation cost for arbitrary $m \in N$ is $O (m d + m^{2})$ .

Coordinate-wise derivatives. We can apply reverse-mode differentiation to recursion (14) in order to compute the entire gradient (c.f., Appendix C). However, in CD, since we only need the derivative of one variable at a time, we can simply use forward-mode differentiation:

\frac{\partial A^{m} (p, x)}{\partial p_{j}} = \frac{1}{m} m \sum t = 1 (- 1)^{t + 1} [\frac{\partial A^{m - t} (p, x)}{\partial p_{j}} D^{t} (p, x) + A^{m - t} (p, x) \frac{\partial D^{t} (p, x)}{\partial p_{j}}],

(15)

where $\frac{\partial D^{t} (p, x)}{\partial p_{j}} = t p_{j}^{t - 1} x_{j}^{t}$ . The advantage of (15) is that we only need to cache $D^{t} (p, x)$ for $t \in [m]$ . Hence the memory complexity per sample is only $O (m)$ instead of $O (d m)$ for (10).

Use in a CD algorithm. Similarly to fm_icml , we assume that the loss function $ℓ$ is $μ$ -smooth and update the elements $p_{j, s}$ of $P$ in cyclic order by $p_{j, s} \leftarrow p_{j, s} - η_{j, s}^{- 1} \frac{\partial F (P)}{\partial p_{j, s}}$ , where we defined

ηj,s\coloneqqμnn∑i=1(∂Am(ps,xi)∂pj,s)2+βand∂F(P)∂pj,s=1nn∑i=1ℓ′(yi,^yi)∂Am(ps,xi)∂pj,s+βpj,s.

(16)

The update guarantees that the objective value is monotonically non-increasing and is the exact coordinate-wise minimizer when $ℓ$ is the squared loss. Overall, the total cost per epoch, i.e., updating all coordinates once, is $O (τ (m) k n_{z} (X))$ , where $τ (m)$ is the time it takes to compute (15). Assuming $D^{t} (p_{s}, x_{i})$ have been previously cached, for $t \in [m]$ , computing (15) takes $τ (m) = m (m + 1) / 2 - 1$ operations. For fixed $m$ , if we unroll the two loops needed to compute (15), modern compilers can often further reduce the number of operations needed. Nevertheless, this quadratic dependency on $m$ means that our CD algorithm is best for small $m$ , typically $m \leq 4$ .

5 HOFMs with shared parameters

HOFMs, as originally defined in fm ; libfm , model each degree with separate matrices $P^{(2)}, \dots, P^{(m)}$ . Assuming that we use the same rank $k$ for all matrices, the total model size of $m$ -order HOFMs is therefore $O (k d m)$ . Moreover, even when using our $O (d m)$ DP algorithm, the cost of computing predictions is $O (k (2 d + \dots + m d)) = O (k d m^{2})$ . Hence, HOFMs tend to produce large, expensive-to-evaluate models. To reduce model size and prediction times, we introduce two new kernels which allow us to share parameters between each degree: the inhomogeneous ANOVA kernel and the all-subsets kernel. Because both kernels are derived from the ANOVA kernel, they share the same appealing properties: multi-linearity, sparse gradients and sparse-data friendliness.

5.1 Inhomogeneous ANOVA kernel

It is well-known that a sum of kernels is equivalent to concatenating their associated feature maps (kernel_book, , Section 3.4). Let $θ = [θ_{1}, \dots, θ_{m}]^{T}$ . To combine different degrees, a natural kernel is therefore

A1→m(p,x;θ)\coloneqqm∑t=1θtAt(p,x).

(17)

The kernel uses all feature combinations of degrees $1$ up to $m$ . We call it inhomogeneous ANOVA kernel, since it is an inhomogeneous polynomial of $x$ . In contrast, $A^{m} (p, x)$ is homogeneous. The main difference between (17) and (7) is that all ANOVA kernels in the sum share the same parameters. However, to increase modeling power, we allow each kernel to have different weights $θ_{1}, \dots, θ_{m}$ .

Evaluation. Due to the recursive nature of Algorithm 1, when computing $A^{m} (p, x)$ , we also get $A^{1} (p, x), \dots, A^{m - 1} (p, x)$ for free. Indeed, lower-degree kernels are available in the last column of the DP table, i.e., $A^{t} (p, x) = a_{d, t} \forall t \in [m]$ . Hence, the cost of evaluating (17) is $O (d m)$ time. The total cost for computing $^y = \sum_{s = 1}^{k} A^{1 \to m} (p_{s}, x; θ)$ is $O (k d m)$ instead of $O (k d m^{2})$ for ${^y}_{HOFM} (x)$ .

Learning. While it is certainly possible to learn $P$ and $θ$ by directly minimizing some objective function, here we propose an easier solution, which works well in practice. Our key observation is that we can easily turn $A^{m}$ into $A^{1 \to m}$ by adding dummy values to feature vectors. Let us denote the concatenation of $p$ with a scalar $γ$ by $[γ, p]$ and similarly for $x$ . From (8), we easily obtain

A^{m} ([γ_{1}, p], [1, x]) = A^{m} (p, x) + γ_{1} A^{m - 1} (p, x) .

(18)

Similarly, if we apply (8) twice, we obtain:

(19)

Applying the above to $m = 2$ and $m = 3$ , we obtain

A^{2} ([γ_{1}, p], [1, x]) = A^{1 \to 2} (p, x; [γ_{1}, 1]) and A^{3} ([γ_{1}, γ_{2}, p], [1, 1, x]) = A^{1 \to 3} (p, x; [γ_{1} γ_{2}, γ_{1} + γ_{2}, 1]) .

(20)

More generally, by adding $m - 1$ dummy features to $p$ and $x$ , we can convert $A^{m}$ to $A^{1 \to m}$ . Because $p$ is learned, this means that we can automatically learn $γ_{1}, \dots, γ_{m - 1}$ . These weights can then be converted to $θ_{1}, \dots, θ_{m}$ by “unrolling” recursion (8). Although simple, we show in our experiments that this approach works favorably compared to directly learning $P$ and $θ$ . The main advantage of this approach is that we can use the same software unmodified (we simply need to minimize (13) with the augmented data). Moreover, the cost of computing the entire gradient by Algorithm 2 using the augmented data is just $O (d m + m^{2})$ compared to $O (d m^{2})$ for HOFMs with separate parameters.

5.2 All-subsets kernel

We now consider a closely related kernel called all-subsets kernel (kernel_book, , Definition 9.5):

S(p,x)\coloneqqd∏j=1(1+pjxj).

(21)

The main difference with the traditional use of this kernel is that we learn $p$ . Interestingly, it can be shown that $S (p, x) = 1 + A^{1 \to d} (p, x; 1) = 1 + A^{1 \to n_{z} (x)} (p, x; 1)$ , where $n_{z} (x)$ is the number of non-zero features in $x$ . Hence, the kernel uses all combinations of distinct features up to order $n_{z} (x)$ with uniform weights. Even if $d$ is very large, the kernel can be a good choice if each training instance contains only a few non-zero elements. To learn the parameters, we simply substitute $A^{m}$ with $S$ in (13). In SG or CD algorithms, all it entails is to substitute $\nabla A^{m} (p, x)$ with $\nabla S (p, x)$ . For computing $\nabla S (p, x)$ , it is easy to verify that $S (p, x) = S (p_{\neg j}, x_{\neg j}) (1 + p_{j} x_{j}) \forall j \in [d]$ and therefore we have

(22)

Therefore, the main advantage of the all-subsets kernel is that we can evaluate it and compute its gradient in just $O (d)$ time. The total cost for computing $^y = \sum_{s = 1}^{k} S (p_{s}, x)$ is only $O (k d)$ .

6 Experimental results

6.1 Application to link prediction

Problem setting. We now demonstrate a novel application of HOFMs to predict the presence or absence of links between nodes in a graph. Formally, we assume two sets of possibly disjoint nodes of size $n_{A}$ and $n_{B}$ , respectively. We assume features for the two sets of nodes, represented by matrices $A \in R^{d_{A} \times n_{A}}$ and $B \in R^{d_{B} \times n_{B}}$ . For instance, $A$ can represent user features and $B$ movie features. We denote the columns of $A$ and $B$ by $a_{i}$ and $b_{j}$ , respectively. We are given a matrix $Y \in {0, 1}^{n_{A} \times n_{B}}$ , whose elements indicate presence (positive sample) or absence (negative sample) of link between two nodes $a_{i}$ and $b_{j}$ . We denote the number of positive samples by $n_{+}$ . Using this data, our goal is to predict new associations. Datasets used in our experiments are summarized in Table 2. Note that for the NIPS and Enzyme datasets, $A = B$ .

Dataset	$n_{+}$	Columns of $A$	$n_{A}$	$d_{A}$	Columns of $B$	$n_{B}$	$d_{B}$
NIPS nips_data	4,140	Authors	2,037	13,649
Enzyme yamanishi_2005	2,994	Enzymes	668	325
GD imc	3,954	Diseases	3,209	3,209	Genes	12,331	25,275
Movielens 100K movielens	21,201	Users	943	49	Movies	1,682	29

Table 2: Datasets used in our experiments. For a detailed description, c.f. Appendix A.

Conversion to a supervised problem. We need to convert the above information to a format FMs and HOFMs can handle. To predict an element $y_{i, j}$ of $Y$ , we simply form $x_{i, j}$ to be the concatenation of $a_{i}$ and $b_{j}$ and feed this to a HOFM in order to compute a prediction ${^y}_{i, j}$ . Because HOFMs use feature combinations in $x_{i, j}$ , they can learn the weights of feature combinations between $a_{i}$ and $b_{j}$ . At training time, we need both positive and negative samples. Let us denote the set of positive and negative samples by $Ω$ . Then our training set is composed of $(x_{i, j}, y_{i, j})$ pairs, where $(i, j) \in Ω$ .

Models compared.

[topsep=0pt,itemsep=1ex,partopsep=0ex,parsep=0.3ex,leftmargin=2ex]
HOFM: ${^y}_{i, j} = {^y}_{HOFM} (x_{i, j})$ as defined in (4) and as originally proposed in fm ; libfm . We minimize (5) by alternating minimization of (13) for each degree.
HOFM-shared: ${^y}_{i, j} = \sum_{s = 1}^{k} A^{1 \to m} (p_{s}, x_{i, j}; θ)$ . We learn $P$ and $θ$ using the simple augmented data approach described in Section 5.1 (HOFM-shared-augmented). Inspired by SimpleMKL simple_mkl , we also report results when learning $P$ and $θ$ directly by minimizing $\frac{1}{| Ω |} \sum_{(i, j) \in Ω} ℓ (y_{i, j}, {^y}_{i, j}) + \frac{β}{2} ∥ P ∥^{2}$ subject to $θ \geq 0$ and $⟨ θ, 1 ⟩ = 1$ (HOFM-shared-simplex).
All-subsets: ${^y}_{i, j} = \sum_{s = 1}^{k} S (p_{s}, x_{i, j})$ . As explained in Section 5.2, this model is equivalent to the HOFM-shared model with $m = n_{z} (x_{i, j})$ and $θ = 1$ .
Polynomial network: ${^y}_{i, j} = \sum_{s = 1}^{k} (γ_{s} + ⟨ p_{s}, x_{i, j} ⟩)^{m}$ . This model can be thought as factorization machine variant that uses a polynomial kernel instead of the ANOVA kernel (c.f., livni ; fm_icml ; tensor_machines ).
Low-rank bilinear regression: ${^y}_{i, j} = a_{i} U V^{T} b_{j}$ , where $U \in R^{d_{A} \times k}$ and $V \in R^{d_{B} \times k}$ . Such model was shown to work well for link prediction in link_prediction and imc . We learn $U$ and $V$ by minimizing $\frac{1}{| Ω |} \sum_{(i, j) \in Ω} ℓ (y_{i, j}, {^y}_{i, j}) + \frac{β}{2} (∥ U ∥^{2} + ∥ V ∥^{2})$ .

Experimental setup and evaluation. In this experiment, for all models above, we use CD rather than SG to avoid the tuning of a learning rate hyper-parameter. We set $ℓ$

to be the squared loss. Although we omitted it from our notation for clarity, we also fit a bias term for all models. We evaluated the compared models using the area under the ROC curve (AUC), which is the probability that the model correctly ranks a positive sample higher than a negative sample. We split the

n_{+}

positive samples into

50 %

for training and

50 %

for testing. We sample the same number of negative samples as positive samples for training and use the rest for testing. We chose

β

from

10^{- 6}, 10^{- 5}, \dots, 10^{6}

by cross-validation and following link_prediction we empirically set

k = 30

. Throughout our experiments, we initialized the elements of

P

randomly by

N (0, 0.01)

	NIPS	Enzyme	GD	Movielens 100K
HOFM ( $m = 2$ )	0.856	0.880	0.717	0.778
HOFM ( $m = 3$ )	0.875	0.888	0.717	0.786
HOFM ( $m = 4$ )	0.874	0.887	0.717	0.786
HOFM ( $m = 5$ )	0.874	0.887	0.717	0.786
HOFM-shared-augmented ( $m = 2$ )	0.858	0.876	0.704	0.778
HOFM-shared-augmented ( $m = 3$ )	0.874	0.887	0.704	0.787
HOFM-shared-augmented ( $m = 4$ )	0.836	0.824	0.663	0.779
HOFM-shared-augmented ( $m = 5$ )	0.824	0.795	0.600	0.621
HOFM-shared-simplex ( $m = 2$ )	0.716	0.865	0.721	0.701
HOFM-shared-simplex ( $m = 3$ )	0.777	0.870	0.721	0.709
HOFM-shared-simplex ( $m = 4$ )	0.758	0.870	0.721	0.709
HOFM-shared-simplex ( $m = 5$ )	0.722	0.869	0.721	0.709
All-subsets	0.730	0.840	0.721	0.714
Polynomial network ( $m = 2$ )	0.725	0.879	0.721	0.761
Polynomial network ( $m = 3$ )	0.789	0.853	0.719	0.696
Polynomial network ( $m = 4$ )	0.782	0.873	0.717	0.708
Polynomial network ( $m = 5$ )	0.543	0.524	0.648	0.501
Low-rank bilinear regression	0.855	0.694	0.611	0.718

Table 3: Comparison of area under the ROC curve (AUC) as measured on the test sets.

Results are indicated in Table 3. Overall the two best models were HOFM and HOFM-shared-augmented, which achieved the best scores on 3 out of 4 datasets. The two models outperformed low-rank bilinear regression on 3 out 4 datasets, showing the benefit of using higher-order feature combinations. HOFM-shared-augmented achieved similar accuracy to HOFM, despite using a smaller model. Surprisingly, HOFM-shared-simplex did not improve over HOFM-shared-augmented except on the GD dataset. We conclude that our augmented data approach is convenient yet works well in practice. All-subsets and polynomial networks performed worse than HOFM and HOFM-shared-augmented, except on the GD dataset where they were the best. Finally, we observe that HOFM were quite robust to increasing $m$ , which is likely a benefit of modeling each degree with a separate matrix.

6.2 Solver comparison

[scale=0.19]fig/nips2

(a) Convergence when

m = 2

[scale=0.19]fig/nips3

(b) Convergence when

m = 3

[scale=0.19]fig/nips4

m = 4

[scale=0.19]fig/scalability_m

(d) Scalability w.r.t. degree

m

Figure 1: Solver comparison for minimizing (13) when varying the degree

m

on the NIPS dataset with

β = 0.1

and

k = 30

. Results on other datasets are in Appendix B.

We compared AdaGrad adagrad , L-BFGS and coordinate descent (CD) for minimizing (13) when varying the degree $m$ on the NIPS dataset with $β = 0.1$ and $k = 30$ . We constructed the data in the same way as explained in the previous section and added $m - 1$ dummy features, resulting in $n = 8, 280$ sparse samples of dimension $d = 27, 298 + m - 1$ . For AdaGrad and L-BFGS, we computed the (stochastic) gradients using Algorithm 2. All solvers used the same initialization.

Results are indicated in Figure 1. We see that our CD algorithm performs very well when $m \leq 3$ but starts to deteriorate when $m \geq 4$ , in which case L-BFGS becomes advantageous. As shown in Figure 1 d), the cost per epoch of AdaGrad and L-BFGS scales linearly with $m$ , a benefit of our DP algorithm for computing the gradient. However, to our surprise, we found that AdaGrad is quite sensitive to the learning rate $η$ . AdaGrad diverged for $η \in {1, 0.1, 0.01}$ and the largest value to work well was $η = 0.001$ . This explains why AdaGrad did not outperform CD despite the lower cost per epoch. In the future, it would be useful to create a CD algorithm with a better dependency on $m$ .

7 Conclusion and future directions

In this paper, we presented the first training algorithms for HOFMs and introduced new HOFM variants with shared parameters. A popular way to deal with a large number of negative samples is to use an objective function that directly maximize AUC link_prediction ; rendle_bpr . This is especially easy to do with SG algorithms because we can sample pairs of positive and negative samples from the dataset upon each SG update. We therefore expect the algorithms developed in Section 3 to be especially useful in this setting. Recently, difacto proposed a distributed SG algorithm for training second-order FMs. It should be straightforward to extend this algorithm to HOFMs based on our contributions in Section 3. Finally, it should be possible to integrate Algorithm 1 and 2

into a deep learning framework such as TensorFlow

tensorflow , in order to easily compose ANOVA kernels with other layers (e.g., convolutional).

References

(1) M. Abadi et al.
TensorFlow: Large-scale machine learning on heterogeneous systems, 2015.
(2) A. G. Baydin, B. A. Pearlmutter, and A. A. Radul. Automatic differentiation in machine learning: a survey. arXiv preprint arXiv:1502.05767, 2015.
(3) M. Blondel, A. Fujino, and N. Ueada. Convex factorization machines. In Proceedings of European Conference on Machine Learning and Principles and Practice of Knowledge Discovery in Databases (ECML PKDD), 2015.
(4) M. Blondel, M. Ishihata, A. Fujino, and N. Ueada. Polynomial networks and factorization machines: New insights and efficient training algorithms. In Proceedings of International Conference on Machine Learning (ICML), 2016.
(5) J. Duchi, E. Hazan, and Y. Singer. Adaptive subgradient methods for online learning and stochastic optimization. Journal of Machine Learning Research, 12:2121–2159, 2011.
(6) GroupLens. http://grouplens.org/datasets/movielens/, 1998.
(7) M. Li, Z. Liu, A. Smola, and Y.-X. Wang. Difacto – distributed factorization machines. In Proceedings of International Conference on Web Search and Data Mining (WSDM), 2016.
(8) R. Livni, S. Shalev-Shwartz, and O. Shamir. On the computational efficiency of training neural networks. In Advances in Neural Information Processing Systems, pages 855–863, 2014.
(9) A. K. Menon and C. Elkan. Link prediction via matrix factorization. In Machine Learning and Knowledge Discovery in Databases, pages 437–452. 2011.
(10) N. Natarajan and I. S. Dhillon. Inductive matrix completion for predicting gene–disease associations. Bioinformatics, 30(12):i60–i68, 2014.
(11) V. Y. Pan. Structured Matrices and Polynomials: Unified Superfast Algorithms. Springer-Verlag New York, Inc., 2001.
(12) A. Rakotomamonjy, F. Bach, S. Canu, and Y. Grandvalet. Simplemkl. Journal of Machine Learning Research, 9:2491–2521, 2008.
(13) S. Rendle. Factorization machines. In Proceedings of International Conference on Data Mining, pages 995–1000. IEEE, 2010.
(14) S. Rendle. Factorization machines with libfm. ACM Transactions on Intelligent Systems and Technology (TIST), 3(3):57–78, 2012.
(15) S. Rendle, C. Freudenthaler, Z. Gantner, and L. Schmidt-Thieme. Bpr: Bayesian personalized ranking from implicit feedback. In
Proceedings of the twenty-fifth conference on uncertainty in artificial intelligence
, pages 452–461, 2009.
(16) S. Rendle, Z. Gantner, C. Freudenthaler, and L. Schmidt-Thieme. Fast context-aware recommendations with factorization machines. In SIGIR, pages 635–644, 2011.
(17) S. Roweis. http://www.cs.nyu.edu/~roweis/data.html, 2002.
(18) J. Shawe-Taylor and N. Cristianini. Kernel Methods for Pattern Analysis. Cambridge University Press, 2004.
(19) V. Vapnik. Statistical learning theory. Wiley, 1998.
(20) G. Wahba. Spline models for observational data, volume 59. Siam, 1990.
(21) Y. Yamanishi, J.-P. Vert, and M. Kanehisa. Supervised enzyme network inference from the integration of genomic data and chemical information. Bioinformatics, 21:i468–i477, 2005.
(22) J. Yang and A. Gittens. Tensor machines for learning target-specific polynomial features. arXiv preprint arXiv:1504.01697, 2015.

Appendix A Dataset descriptions

[topsep=0pt,itemsep=1ex,partopsep=0ex,parsep=1ex,leftmargin=2ex]
NIPS: co-author graph of authors at the first twelve editions of NIPS, obtained from [17]. For this dataset, as well as the Enzyme dataset below, we have $A = B$ . The co-author graph comprises $n_{A} = n_{B} = 2, 037$ authors represented by bag-of-words vectors of dimension $d_{A} = d_{B} = 13, 649$ (words used by authors in their publications). The number of positive samples is $n_{+} = 4, 140$ .
Enzyme: metabolic network obtained from [21]. The network comprises $n_{A} = n_{B} = 668$ enzymes represented by three sets of features: a $157$ -dimensional vector of phylogenetic information, a $145$ -dimensional vector of gene expression information and a $23$ -dimensional vector of gene location information. We concatenate the three sets of information to form feature vectors of dimension $d_{A} = d_{B} = 325$ . Original enzyme similarity scores are between $0$ and $1$
. We binarize the scores using
$0.95$ as threshold. The resulting number of positive samples is $n_{+} = 2, 994$ .
GD: human gene-disease association data obtained from [10]. The bipartite graph is comprised of $n_{A} = 3, 209$ diseases and $n_{B} = 12, 331$ genes. We represent each disease using a vector of $d_{A} = 3, 209$ dimensions, whose elements are similarity scores obtained from MimMiner. The study [10] also used bag-of-words vectors describing each disease but we found these to not help improve performance both for baselines and proposed methods. We represent each gene using a vector of $d_{B} = 25, 275$ features, which are the concatenation of $12, 331$ similarity scores obtained from HumanNet and $12, 944$ gene-phenotype associations from 8 other species. See [10] for a detailed description of the features. The number of positive samples is $n_{+} = 3, 954$ .
Movielens 100K: recommender system data obtained from [6]. The bipartite graph is comprised of $n_{A} = 943$ users and $n_{B} = 1, 682$
movies. For users, we convert age, gender, occupation and living area (first digit of zipcode) to a binary vector using a one-hot encoding. For movies, we use the release year and genres. The resulting vectors are of dimension
$d_{A} = 49$ and $d_{B} = 29$ , respectively. Original ratings are between $1$ and $5$ . We binarize the ratings using $5$ as threshold, resulting in $n_{+} = 21, 201$ positive samples.

Appendix B Additional experiments

b.1 Solver comparison

We also compared AdaGrad, L-BFGS and coordinate descent (CD) on the Enzyme, Gene-Disease (GD) and Movielens 100K datasets. Results are indicated in Figure 2, 3 and 4, respectively.

[scale=0.14]fig/metabolic2

(a)

m = 2

[scale=0.14]fig/metabolic3

(b)

m = 3

[scale=0.14]fig/metabolic4

(c)

m = 4

Figure 2: Solver comparison for minimizing (13) on the Enzyme dataset. We set

β

to the values with best test-set performance, which were

β = 0.1

β = 0.1

and

β = 0.01

, respectively. We set

k = 30

[scale=0.14]fig/gd2

(a)

m = 2

[scale=0.14]fig/gd3

(b)

m = 3

[scale=0.14]fig/gd4

(c)

m = 4

Figure 3: Solver comparison for minimizing (13) on the GD dataset. We set

β

to the values with best test-set performance, which were

β = 0.01

β = 0.01

and

β = 0.0001

, respectively. We set

k = 30

[scale=0.14]fig/movielens2

(a)

m = 2

[scale=0.14]fig/movielens3

(b)

m = 3

[scale=0.14]fig/movielens4

(c)

m = 4

Figure 4: Solver comparison for minimizing (13) on the Movielens 100K dataset. We set

β

to the values with best test-set performance, which were

β = 10^{- 3}

β = 10^{- 4}

and

β = 10^{- 6}

, respectively. We set

k = 30

b.2 Recommender system experiments

As we explained in Section 5.2, the all-subsets kernel can be a good choice if the number of non-zero elements per sample is small. To verify this assumption, we ran experiments on two recommender system tasks: Movielens 1M and Last.fm. We used the exact same setting as in [4, Section 9.3]. For each rating $y_{i}$ , the corresponding $x_{i}$ was set to the concatenation of the one-hot encodings of the user and item indices. We compared the following models:

FM: ${^y}_{i} = ⟨ w, x_{i} ⟩ + \sum_{s = 1}^{k} A^{2} (p_{s}, x_{i})$
FM-augmented: ${^y}_{i} = \sum_{s = 1}^{k} A^{2} (p_{s}, {~ x}_{i})$ where ${~ x}_{i}^{T} = [1, x_{i}^{T}]$
All-subsets: ${^y}_{i} = \sum_{s = 1}^{k} S (p_{s}, x_{i})$
Polynomial networks: ${^y}_{i} = {~ x}_{i} U V^{T} {~ x}_{i}$ (c.f. [4] for more details)

Results are indicated in Figure 5. We see that All-subsets performs relatively well on these tasks.

[scale=0.22]fig/recsys_m1m

(a) Movielens 1M

[scale=0.22]fig/recsys_lastfm

(b) Last.FM

Figure 5: Model comparison on two recommender system datasets.

Appendix C Reverse-mode differentiation on the alternative recursion

We now describe how to apply reverse-mode differentiation to the alternative recursion (14) in order to compute the entire gradient efficiently. Let us introduce the shorthands $at\coloneqqAt(p,x)$ and $dt\coloneqqDt(p,x)$ . We can then write the recursion as

a_{m} = \frac{1}{m} m \sum t = 1 (- 1)^{t + 1} a_{m - t} d_{t} .

(23)

For concreteness, let us illustrate the recursion for $m = 3$ . We have

a_{1} = a_{0} d_{1}, a_{2} = \frac{1}{2} (a_{1} d_{1} - a_{0} d_{2}) and a_{3} = \frac{1}{3} (a_{2} d_{1} - a_{1} d_{2} + a_{0} d_{3}) .

(24)

We see that $a_{2}$ influences $a_{3}$ , and $a_{1}$ influences $a_{2}$ and $a_{3}$ . Likewise, $d_{3}$ influences $a_{3}$ , $d_{2}$ influences $a_{2}$ and $a_{3}$ , and $d_{1}$ influences $a_{1}$ , $a_{2}$ and $a_{3}$ . Let us denote the adjoints $~at\coloneqq∂am∂at$ and $~dt\coloneqq∂am∂dt$ . For general $m$ , summing over quantities that influences $a_{t}$ and $d_{t}$ , we obtain

{~ a}_{t} = m \sum s = t + 1 \frac{(- 1)^{s - t + 1}}{s} {~ a}_{s} d_{s - t} and {~ d}_{t} = (- 1)^{t + 1} m \sum s = t \frac{1}{s} {~ a}_{s} a_{s - t} .

(25)

Let us denote the adjoint of $p_{j}$ by $~pj\coloneqq∂am∂pj$ . We know that $p_{j}$ directly influences only $d_{1}, \dots, d_{m}$ and therefore

{~ p}_{j} = m \sum t = 1 \frac{\partial a_{m}}{\partial d_{m}} \frac{\partial d_{m}}{\partial p_{j}} = m \sum t = 1 {~ d}_{t} t p_{j}^{t - 1} x_{j}^{t} .

(26)

Assuming that $d_{1}, \dots, d_{m}$ and $a_{1}, \dots, a_{m}$ have been previously computed, which takes $O (d m + m^{2})$ , the procedure for computing the gradient can be summarized as follows:

Initialize ${~ a <}$

Higher-Order Factorization Machines

Authors

Higher-Order Ranking and Link Prediction: From Closing Triangles to Closing Higher-Order Motifs

Adaptive Factorization Network: Learning Adaptive-Order Feature Interactions

Factorization Machines with Regularization for Sparse Feature Interactions

DS-FACTO: Doubly Separable Factorization Machines

Strongly Hierarchical Factorization Machines and ANOVA Kernel Regression

Stabilizing Linear Prediction Models using Autoencoder

Efficiently Finding Higher-Order Mutants

Code Repositories

tffm

GSOCFactorizationMachinesSystemML2

1 Introduction

2 Factorization machines (FMs)

3 Linear-time stochastic gradient algorithms for HOFMs

4 Coordinate descent algorithm for arbitrary-order HOFMs

5 HOFMs with shared parameters

5.1 Inhomogeneous ANOVA kernel

5.2 All-subsets kernel

6 Experimental results

6.1 Application to link prediction

6.2 Solver comparison

7 Conclusion and future directions

References

Appendix A Dataset descriptions

Appendix B Additional experiments

b.1 Solver comparison

b.2 Recommender system experiments

Appendix C Reverse-mode differentiation on the alternative recursion

Higher-Order Factorization Machines

Authors

Higher-Order Ranking and Link Prediction: From Closing Triangles to Closing Higher-Order Motifs

Adaptive Factorization Network: Learning Adaptive-Order Feature Interactions

Factorization Machines with Regularization for Sparse Feature Interactions

DS-FACTO: Doubly Separable Factorization Machines

Strongly Hierarchical Factorization Machines and ANOVA Kernel Regression

Stabilizing Linear Prediction Models using Autoencoder

Efficiently Finding Higher-Order Mutants

Code Repositories

tffm

GSOCFactorizationMachinesSystemML2

This week in AI

1 Introduction

2 Factorization machines (FMs)

3 Linear-time stochastic gradient algorithms for HOFMs

4 Coordinate descent algorithm for arbitrary-order HOFMs

5 HOFMs with shared parameters

5.1 Inhomogeneous ANOVA kernel

5.2 All-subsets kernel

6 Experimental results

6.1 Application to link prediction

6.2 Solver comparison

7 Conclusion and future directions

References

Appendix A Dataset descriptions

Appendix B Additional experiments

b.1 Solver comparison

b.2 Recommender system experiments

Appendix C Reverse-mode differentiation on the alternative recursion