Collaborative filtering is often applied in recommendation systems that primarily serve Internet services, such as e-commerce and video distribution platforms (Herlocker et al., 2000; Su and Khoshgoftaar, 2009). The essence of collaborative filtering lies in generating personalized recommendations based on the intrinsic similarities between users and items. Collaborative filtering relies on training data, in which users assign scores or ratings to various items. As it is common for users to omit ratings of specific items, leading to missing data, the central objective of collaborative filtering is to predict the scores for unrated items. Matrix factorization techniques, particularly nonnegative matrix factorization (NMF) (Lee and Seung, 1999), are frequently employed. When using NMF for collaborative filtering, the ranking matrix V, whose entries are nonnegative, is approximated as the product of two nonnegative matrices W and H, that is, V≈WH. The standard approach involves minimizing the difference between V and WH. In the optimization procedure, each element of W and H is constrained to be nonnegative. While the multiplicative update algorithm is the most prevalent approach for NMF (Lee and Seung, 2000), we focus on an alternative technique known as the nonnegative least-squares approach using the projected gradient method (PGM) (Lin, 2007). The convergence of the alternative update method for NMF was proved by Lin (2007). Such an alternative optimization method is essential for solving nonnegative/binary matrix factorization (NBMF), which is an extension of NMF. O'Malley et al. (2018) and Golden and O'Malley (2021) used D-Wave's quantum annealers to solve quadratic binary optimization problems involved in NBMF and demonstrated a speedup compared with two classical solvers.

In recent years, Ising machines, initially designed to solve combinatorial optimization problems efficiently, have found new applications in the field of machine learning, expanding their scope beyond their original purpose (Kitai et al., 2020; Willsch et al., 2020; Nath et al., 2021). Ising machines are special-purpose computers for solving combinatorial optimization problems. They are realized by several types of devices, such as quantum annealers (Johnson et al., 2011), digital processors based on simulated annealing (Yamaoka et al., 2016; Aramon et al., 2019; Yamamoto et al., 2020), digital processors based on simulated bifurcation (SB) (Goto et al., 2019; Hidaka et al., 2023), and coherent Ising machines (Inagaki et al., 2016; McMahon et al., 2016; Pierangeli et al., 2019). As Ising machines usually accept problems described by the Ising model or quadratic unconstrained binary optimization formulation, their application to machine learning requires hybrid methods that utilize both an Ising machine and a general-purpose computer (e.g., a CPU). In NBMF, the matrix elements of H are binary, whereas those of W are real and nonnegative. Therefore, an Ising machine is employed to accelerate the update of matrix H, whereas a general-purpose computer handles the update of matrix W. As the updates of matrices H and W are repeated alternately, NBMF inevitably involves a computation time overhead owing to the communication between the Ising machine and the CPU. The advantages and disadvantages of NMF and NBMF remain unclear in terms of solution quality, computation time, and applicability to sparse problems.

In this paper, we propose a novel approach for applying NBMF to collaborative filtering and demonstrate the advantages of utilizing a low-latency Ising machine to execute the proposed method. Previous studies have employed NBMF for image analysis that deals with dense data matrices, where the majority of matrix elements have nonzero values (O'Malley et al., 2018; Asaoka and Kudo, 2020, 2023). By contrast, collaborative filtering involves sparse data matrices, with most elements remaining undetermined. We propose a modified NBMF algorithm that masks undetermined elements within the data matrix to improve the prediction accuracy. In addition, we compare NBMF with NMF in terms of solution quality and computation time, and investigate the dependency of these characteristics on the sparsity and size of the problem. To accelerate the NBMF algorithm, we used an SB-based machine implemented with a field-programmable gate array (FPGA) (Goto et al., 2021; Hidaka et al., 2023) that supports up to 2,048 spins and has full spin-to-spin connectivity (no need for minor embedding techniques required for local-connectivity Ising machines). Incorporating an SB-based machine to update the binary matrix elements yields a substantial reduction in the overall computational time required for NBMF compared with NMF.

Furthermore, the low-latency characteristic of the SB-based machine is advantageous for executing the iterative method using a general-purpose computer (a CPU) and an Ising machine, alternatively, reducing the communication time between them. It is also possible to use a cloud-hosted Ising machine for executing the proposed method. While a high-performance cloud-hosted Ising machine can significantly reduce computation time, the communication time of accessing it may negate the benefits. Therefore, utilizing the low-latency system is crucial. This study presents the first empirical evidence that NBMF, when implemented with a low-latency Ising machine, surpasses NMF in terms of both solution quality and overall computational efficiency.

NBMF and NMF decompose a nonnegative n×m matrix V into an n×k matrix W and a k×m matrix H:

where W is a nonnegative real matrix. While H is a binary matrix in which each element is 0 or 1 for NBMF, it is a nonnegative real matrix for NMF. We assume that n>k and m>k, which implies that NBMF and NMF provide low-rank matrix approximations of V. The rank constraint is helpful to prevent overfitting. Moreover, NBMF can be more resilient to overfitting due to the binary nature of matrix H.

In the context of collaborative filtering, V is a rating matrix, where the (i, j) element v_ij represents user i's rating of item j. In matrix W, each row corresponds to a user, while each column represents a basis vector associated with user preferences. In other words, W consists of k basis vectors, with each dimension being n. Meanwhile, each column of H represents the coefficient vector related to the corresponding item. In NBMF, this coefficient vector indicates the combination of the selected basis vectors for the corresponding item. In general, the rating matrix contains numerous unrated entries. Matrix factorization techniques optimize W and H so that each rated entry in V is well approximated by the corresponding element of WH. Then, each unrated entry in V is estimated by the corresponding element of WH.

The comparison between NMF and other collaborative filtering techniques has already been extensively studied (Lee et al., 2014; Singh et al., 2024). Compared to user-based and item-based collaborative filtering techniques, matrix factorization techniques demonstrated better performance in recommendation systems on multi-criteria datasets (Singh et al., 2024). In particular, NMF is scalable to large datasets and can capture individual user preferences. However, there has been no direct comparison between NMF and NBMF. This paper focuses on comparing the two methods.

Figure 3 shows a comparison of RMSE and computation time of NBMF and NMF for 10 epochs. Each epoch involves updating matrix W followed by updating matrix H. The data points represent the averages of RMSE and computation time at each epoch, with some error bars too small to be observed. Figure 3 shows that the RMSE decays more rapidly in NBMF than in NMF for all datasets. Although the difference in the RMSE at each epoch between NBMF and NMF was negligible, the difference among the datasets was remarkable. The difference among the datasets originates from the frequency distribution of the ratings in each dataset. As elaborated later, when the distribution was sharp and the variance was small, the RMSE tended to be small. In contrast, when the variance in the frequency was large, the RMSE was relatively large.

The filling rate of a rating matrix, which is the proportion of rated elements, influences collaborative filtering. However, the filling rate dependence of the RMSE varied across the datasets, as shown in Figure 4. Here, the RMSE was calculated after 10 epochs and averaged over five trials.

NMF is expected to produce lower RMSE values than NBMF due to its higher resolution. However, in Figure 4, case (a) demonstrates that the RMSE values for both NBMF and NMF were similar when the filling rate was around 20%. Furthermore, in case (d), the RMSE for NBMF was smaller than that for NMF at the filling rate of ~20%. This inconsistent behavior suggests that the datasets have significant differences in their features.

Figure 5 shows the advantages of the masking procedure. The masking data for each dataset were the same as those used for NBMF in Figure 3. Notably, the RMSE without the masking procedure was about more than three times larger than the RMSE with the masking procedure for all datasets. Furthermore, the difference in computation time between the two procedures was minimal. These findings indicate that the masking procedure provides apparent benefits in collaborative filtering, despite a minor drawback in terms of computation time.

The results indicate that the RMSE reflects specific properties of the data. Here, we focus on the frequency distribution of ratings, as illustrated in Figure 6. The distribution represents the percentage of ratings (1, 2, 3, 4, and 5) among the rated elements in a rating matrix, with a filling rate of 20%. The distributions in (a) and (b) showed a broad peak, and the corresponding RMSE had a similar value at 10 epochs in Figures 3a, b. However, the distribution in (c) showed two peaks at 1 and 5, resulting in significant variance. The corresponding RMSE at 10 epochs in Figure 3c was larger than those of the other three data. By contrast, the distribution in (d) had a steep peak at 4, indicating that more than 80% of the rated elements had a value of 4. The corresponding RMSE at 10 epochs in Figure 3d was significantly smaller than those of the other three data. This observation indicates that a distribution with a sharp peak and small variance typically results in a smaller RMSE. However, a distribution with a broader peak and larger variance often results in a larger RMSE.

The computation time for NBMF was significantly shorter than that for NMF under the same problem setup, as shown in Figure 7, for all the datasets. The total time required to update matrix W over 10 epochs was the same for both NMF and NBMF. However, the time required to update matrix H in NMF was approximately six times longer than that in NBMF. This discrepancy suggests that the use of an SB-based machine accelerates the computation to update H. Additional time is required to minimize Equation 13 during the update of H. Executing minimization using the SB-based machine involves transforming the objective function into the Ising model form, as explained in Section 3.4. Using the SB-based machine causes the communication time between the CPU and the FPGA, although it is a small fraction of the total time. Nevertheless, the overall computation time for NBMF, including these additional factors, was shorter than that for NMF.

Throughout this study, the ratio of the number of features to the number of users was fixed at k/n = 0.5. In general, the RMSE tends to decrease as the ratio increases. However, the computation time increases with the ratio because the matrix sizes of H and W increase. Therefore, a moderate value needs to be selected for this ratio. As shown in Figure 8, the rate of improvement in the RMSE was slow for ratios of 0.5 or greater across all datasets. Considering this result, we chose k/n = 0.5 as an appropriate value for this ratio.

Our results support the computational advantages of NBMF. However, several limitations exist. Notably, the performance of NBMF is susceptible to the characteristics of datasets. NMF, which operates on continuous variables, shows comparable or superior accuracy in certain instances compared to NBMF. This higher accuracy is due to the greater resolution of continuous representations compared to binary ones. Furthermore, it is necessary to employ a low-latency system to realize the advantage of computation time. Even with a high-performance Ising machine, communication overhead between the CPU and the Ising machine can significantly impact overall performance. Therefore, utilizing a low-latency Ising machine is crucial for efficiently executing NBMF.

In summary, we proposed a novel approach that employs NBMF with masking for collaborative filtering, and our findings demonstrate a substantial improvement in estimation performance as a result of the masking procedure. Moreover, our results highlighted the computational advantage of employing an SB-based machine in NBMF. NBMF with masking can be applied in collaborative filtering across various datasets. This study reveals the potential of NBMF by utilizing an Ising machine for a wide range of applications.

The efficacy of hybrid concepts using both a general-purpose processor and an Ising machine can extend to other algorithms, indicating significant potential for further research in this area. Similar hybrid algorithms will become increasingly active in the future as low-latency Ising machines become more advanced and popular.