Memory-Efficient Factorization Machines via Binarizing both Data and Model Coefficients

Even though SEFM improves the model expressiveness capacity of FM and therefore achieves higher accuracy than FM on highly nonlinear datasets, SEFM increases the memory cost of FM by $b$ times since the model parameters are $\mathbf{w}\in$ $\mathbb{R}^{p}$ and $\mathbf{V}\in\mathbb{R}^{p\times m}$ where $p$ is $b$ times larger than the original data dimensionality $d$. The high space complexity of SEFM will limit its applications in resource-constrained scenarios, such as model deployment on mobile phones or IoT devices.

In this section, we describe our proposed algorithm to solve (7) based on Straight Through Estimator (STE) (Bengio et al.,, 2013). STE is popularly used for training binary convolution neural networks (Alizadeh et al.,, 2018). The idea of STE is to learn a full-precision variable as the proxy of the binary variable during the training process. As in our proposed model, we introduce a full-precision $\mathbf{\widetilde{w}}$ as the proxy of binary $\mathbf{w}^{b}$ in (7) and a full-precision $\mathbf{\widetilde{V}}$ as the proxy of binary $\mathbf{V}^{b}$ in (7). Instead of directly learning binary variables $\mathbf{w}^{b}$ and $\mathbf{V}^{b}$, we learn their full-precision proxy variables $\mathbf{\widetilde{w}}$ and $\mathbf{\widetilde{V}}$. To estimate the binary variables $\mathbf{w}^{b}$ and $\mathbf{V}^{b}$ and scaling parameter $\alpha$, $\beta$, we propose to minimize $\|\mathbf{\widetilde{w}}-\alpha\mathbf{w}^{b}\|^{2}$ and $\|\mathbf{\widetilde{V}}-\beta\mathbf{V}^{b}\|^{2}$.

Here we present our detailed derivation of estimating $\mathbf{w}^{b}$, $\mathbf{V}^{b}$ and scaling parameter $\alpha$, $\beta$ based on the full-precision proxy variables $\mathbf{\widetilde{w}}$ and $\mathbf{\widetilde{V}}$. The derivation follows the similar procedures of optimizing binary convolution neural networks (Rastegari et al.,, 2016). In the STE framework, we use SGD to update full-precision proxy variables $\mathbf{\widetilde{w}}$ and $\mathbf{\widetilde{V}}$ instead of $\mathbf{w}^{b}$, $\mathbf{V}^{b}$. We propose to get an optimal approximation for $\mathbf{w}\approx\alpha\mathbf{w}^{b}$ and $\mathbf{V}\approx\beta\mathbf{V}^{b}$ by minimizing $\|\mathbf{\widetilde{w}}-\alpha\mathbf{w}^{b}\|^{2}$ and $\|\mathbf{\widetilde{V}}-\beta\mathbf{V}^{b}\|^{2}$.

Let us consider minimize $\|\mathbf{\widetilde{w}}-\alpha\mathbf{w}^{b}\|^{2}$ first. By expanding the equation$\|\mathbf{\widetilde{w}}-\alpha\mathbf{w}^{b}\|^{2}$, it can be rewritten as

Since $\mathbf{w}^{b}\in\{-1,1\}^{p}$, $(\mathbf{w}^{b})^{T}\mathbf{w}^{b}$ will equals to a constant $p$. The scaling parameter $\alpha>0$ and $\widetilde{w}$ is a known variable in (8). Therefore, to minimize (8) is equivalent to the following optimization problem,

It is straightforward to the optimal solution for (9), which is $\mathbf{w}^{b}=\text{sign}({\mathbf{\widetilde{w}}})$. To obtain the optimal solution for the scaling parameter $\alpha$, we take the derivative of (8) with respective to $\alpha$. It is $-2\mathbf{{\widetilde{w}}}^{T}\mathbf{w}^{b}+2\alpha p$. By setting it to 0, we will get $\alpha=\frac{\mathbf{{\widetilde{w}}}^{T}\mathbf{w}^{b}}{p}$. Since $\mathbf{w}^{b}=\text{sign}({\mathbf{\widetilde{w}}})$, the scaling parameter $\alpha$ will be

By following the same procedure, we can also get the optimal $\mathbf{V}^{b}=\text{sign}({\mathbf{\widetilde{V}}})$ and the optimal $\beta$ as

The optimal binary $\mathbf{w}^{b}$ is obtained by $\mathbf{w}^{b}=\text{sign}({\mathbf{\widetilde{w}}})$, where sign() is the element-wise sign function which return 1 if the element is larger or equal than zero and return $-1$ otherwise. Similarly, $\mathbf{V}^{b}$ can be obtained by $\mathbf{V}^{b}=\text{sign}({\mathbf{\widetilde{V}}})$. Since the sign() function are not differentiable, STE estimates the gradients with respect to $\mathbf{\widetilde{w}}$ and $\mathbf{\widetilde{V}}$ as if the non-differentiable function sign() is not present. In other word, STE will simply estimate $\frac{\partial\text{sign}(\tilde{w}_{j})}{\partial\tilde{w}_{j}}$ as $\frac{\partial\tilde{w}_{j}}{\partial\tilde{w}_{j}}=1$. In practice, we also employ the gradient clipping as in (Hubara et al.,, 2016). Then, the gradient for the non-differentiable sign function is

Therefore, the gradient of (7) with respect to $\tilde{w}_{j}$ (i.e., the $j$-th element in the proxy variable $\mathbf{\tilde{w}}$) based on the $i$-th data sample can be estimated as

Any convex loss function can be used in (13) and the gradient $\frac{\partial\ell(y_{i},f(\mathbf{z}_{i}))}{\partial f(\mathbf{z}_{i})}$ can be easily obtained for different loss function.

Similarly, the gradient of (7) with respect to $\tilde{v}_{jf}$ (i.e., the $(j,f)$-th element in $\mathbf{\widetilde{V}}$) based on the $i$-th data sample can be estimated as

where $\frac{\partial f(\mathbf{z}_{i})}{v^{b}_{jf}}$ can be computed based a similar form as in (3). That is

In our objective optimization problem (7), the feature representation $\mathbf{Z}$ after binning is very sparse since it is generated by applying one-hot encoding to each original feature. Also, the sparsity rate for different features in $\mathbf{Z}$ will be very different, as shown in Figure 2. Therefore, parameters associated with features with high sparsity rate will be updated less frequently than parameters associated with features with low sparsity rate in SGD. If we directly applied the standard SGD where the same learning rate $\eta$ is set for all parameters, it could be either too small for the parameters associated with features with high sparsity rate or too large for the parameters associated with features with low sparsity rate.

To address this issue, we propose to use Adagrad (Duchi et al.,, 2011) to update proxy variables $\mathbf{\widetilde{w}}$ and $\mathbf{\widetilde{V}}$. Adagrad uses different learning rates for different model parameters and the learning rates can be adaptively adjusted during the training process. It has been shown that Adagrad can greatly improve the robustness of SGD (Ruder,, 2016). The basic intuition of Adagrad is to give low learning rates for parameters associated with features with low sparsity rate and to given high learning rates for parameters associated with features with high sparsity rate.

For a parameter (e.g., $\tilde{w}_{j}$ or $\tilde{v}_{jf}$) in our model, based on the idea of Adagrad, we set the learning rate in the $t$-th update $\eta_{t}$ to $\frac{\eta}{\sqrt{s_{t}}+\epsilon}$ where $s_{t}$ equals the sum of squares of previously observed gradients up to time $t$ and $\epsilon$ is a small constant that prevents division by 0.

Therefore, for updating the $j$-th element in $\mathbf{\tilde{w}}$, we can first obtain the gradient $g_{\tilde{w}_{j}}$ as in (13). Then, the updating rule for $\tilde{w}_{j}$ at the $t$-th time will be

where $s_{j}^{t}$ aggregates sum of squares of previously observed gradients up to time $t$. That is, $s_{j}^{t}=s_{j}^{t-1}+(g_{\tilde{w}_{j}})^{2}.$

Similarly, for updating the $(j,f)$-th parameter in $\mathbf{\tilde{V}}$, we first obtain the gradient $g_{\tilde{v}_{jf}}$ as in (14). Then, the updating rule for $\tilde{v}_{jf}$ at the $t$-th time will be

where $\tilde{s}_{jf}^{t}$ aggregates sum of squares of previously observed gradients with respect to $\tilde{v}_{jf}$ up to time $t$. That is, $\tilde{s}_{jf}^{t}=\tilde{s}_{jf}^{t-1}+(g_{\tilde{v}_{jf}})^{2}.$

As shown in (16) and (17), Adagrad adaptively adjusts the learning rate $\eta$ at each time step $t$ for each parameter in $\mathbf{\tilde{w}}$ and $\mathbf{\tilde{V}}$ based on their past gradients. In our experiment, we have shown that Adagrad can achieve lower training loss and is more robust than SGD.

We summarize our proposed method in Algorithm 1. The step 1 obtains the new feature representation $\mathbf{Z}$ by applying one-hot encoding on each feature in input data $\mathbf{X}$ which can be done in $O(nd)$ time (Lan and Geng,, 2019). The step 2 initializes the proxy variables $\mathbf{\tilde{w}}$, $\mathbf{\tilde{V}}$ and $\mathbf{s}\in\mathbb{R}^{p}=\mathbf{0}$, $\tilde{\mathbf{S}}\in\mathbb{R}^{p\times m}=\mathbf{0}$ which are used for adjusting learning rate in Adagrad. From step 3 to 21, we update the proxy parameters and the binary parameters using STE with Adagrad. The time complexity for learning the binarized parameter is $O(nnz(\mathbf{Z})ml)$, where $nnz(\mathbf{Z})$ denotes the number of non-zero values in $\mathbf{Z}$, $m$ is the low-rank parameter and $l$ is the number of passes over dataset $\mathbf{Z}$. Empirically, the algorithm converge very fast and $l$ is a small number in our experiments.

Similar to Lemma 3.1 in (Rendle,, 2010), the prediction score of our model for a given input $\mathbf{z}_{i}$ can be written as

Since each element in $\mathbf{V}^{b}$ is either 1 or $-1$, ${v^{b}_{jf}}^{2}=1$ and $\sum_{j=1}^{p}{v^{b}_{jf}}^{2}z_{ij}^{2}=nnz(\mathbf{z}_{i})$ where $nnz(\mathbf{z}_{i})$ is the number of non-zero values in $\mathbf{z}_{i}$ which will be equal to $d$ because of one-hot encoding. Also note that $\mathbf{w}^{b}$ and $\mathbf{V}^{b}$ are binary weights and $\mathbf{z}_{i}$ is a sparse vector with $d$ nonzero values. The nonzero value in $z_{i}$ can only be 1. Therefore, the model inference as shown in (18) can be estimated by very efficiently only using addition and subtraction operations (without the need of multiplication operation). Also the arithmetic operations in (18) can be further replaced cheap XNOR and POPCOUNT bitwise operations (Rastegari et al.,, 2016).

During model inference, we do not need to keep the proxy variables $\mathbf{\tilde{w}}$ and $\mathbf{\tilde{V}}$. We only need the binary parameters $\mathbf{w}^{b}$, $\mathbf{V}^{b}$ and scaling parameter $\alpha$, $\beta$ to compute the prediction score for a test sample based on (18). We summarize the memory cost for different FM models in Table 1.

For each dataset, we randomly select 70% as training data and use the remaining 30% as test data. This procedure is repeated for 10 times and we report the average accuracy on test data. The low-rank parameter $m$ for FM, SEFM, DFM and Binarized FM is chosen from {16, 32, 64, 128}. The parameter $b$ (i.e., the number of bins) of SEFM and Binarized FM is chosen from {10, 20, 30, 40, 50}. The regularization parameter in Liblinear, FM, DFM and SEFM is chosen from $\{10^{-2},10^{-1},\dots,10^{2}\}$. The optimal parameter combination is selected by 5-fold cross-validation on training data. The suggested settings are used in NFM, but for webspam, the hidden factor and the number of layers are both decreased to make it successfully run in our computing environment.

The accuracy, prediction time and memory cost of all six algorithms are reported in Table 2. Note that the memory cost of model parameters is normalized by the memory cost of standard FM. For example, 10x means the memory cost is 10 times larger than standard FM. The best accuracy for each dataset is in bold. The second best accuracy is in italic.

Accuracy. As shown in Table 2, both SEFM and our proposed method Binarized FM achieve higher accuracy than Liblinear, FM and DFM on all eight datasets. This is due to the fact that the expressiveness capacities of Liblinear, FM and DFM are smaller than SEFM and our proposed binarized FM. Overall, the linear classifier (i.e., Liblinear) gets the lowest accuracy, while FM and DFM get higher accuracy than the linear classifier. DFM is less accurate than FM because of binarizing of feature interaction model parameters. Our proposed method Binarized FM achieves comparable results with SEFM with much lower memory cost.

Memory cost. With regard to memory cost, Liblinear and DFM have much lower memory cost than FM, as Liblinear is a linear classifier being linear in the parameters and DFM binarize the 2-dimensional parameter $\mathbf{V}$ of FM. As for our proposed method Binarized FM, which binarizes both the parameter $\mathbf{w}$ and $\mathbf{V}$ and transforms features by one-hot encoding, it has relatively lower memory cost on circles, moons, pendigits, ijcnn and webspam than FM and slightly higher memory costs on the other datasets. The memory costs of SEFM and NFM are much higher than FM and our proposed method on all of datasets.

We visually present the variation between memory cost and accuracy among four methods–FM, SEFM, DFM and Binarized FM on six different datasets as shown in Figure 3. It is clearly discovered that our proposed method Binarized FM can achieve higher accuracy than FM and DFM. With the same memory cost, Binarized FM has predominantly higher accuracy than FM. SEFM has the high accuracy and the highest memory cost, while the memory cost of our proposed method Binarized FM is obviously lower and its accuracy approaches that of SEFM.

The effect of scaling parameter $\alpha$ and $\beta$. The accuracy and prediction time of our algorithms with and without scaling parameter are reported in Table 3. After adding the scaling parameter $\alpha$ and $\beta$, our method achieves higher accuracy on most of datasets than without using them, especially on banana, segment, pendigits, ijcnn and webspam.

Illustration of Decision Boundaries on Synthetic Datasets. We first illustrate that our proposed method can produce highly nonlinear decision boundary as SEFM on the synthetic datasets. We show the decision boundaries of FM, SEFM, DFM and our proposed binarized FM on moons and banana datasets in Figure 4. As shown in this figure, both FM and DFM can not capture the highly nonlinear pattern on these two datasets since they only used polynomial expansion on the original data to model the nonlinear relationship. SEFM and our proposed binarized FM can produce complex piecewise axis perpendicular decision boundary.

Comparison of the Convergence of Adagrad and SGD in Binarized FM. In addition to demonstrating the performance of our proposed method, we illustrate our training process by presenting the variation of training loss during iteration to prove our optimization method–Adagrad has more efficient performance than SGD. As shown in Figure 5, Adagrad reaches the lower training loss fast and gradually tend to converge.