Collaborative filtering based on nonnegative/binary matrix factorization

Collaborative filtering generates recommendations by exploiting user-item similarities based on rating data, which often contains numerous unrated items. To predict scores for unrated items, matrix factorization techniques such as nonnegative matrix factorization (NMF) are often employed. Nonnegative/binary matrix factorization (NBMF), which is an extension of NMF, approximates a nonnegative matrix as the product of nonnegative and binary matrices. While previous studies have applied NBMF primarily to dense data such as images, this paper proposes a modified NBMF algorithm tailored for collaborative filtering with sparse data. In the modified method, unrated entries in the rating matrix are masked, enhancing prediction accuracy. Furthermore, utilizing a low-latency Ising machine in NBMF is advantageous in terms of the computation time, making the proposed method beneficial.

1 Introduction

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.

2 Problem formulation

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

$\begin{array}{ll}V\approx WH,& (1)\end{array}$

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.

3 Methods

3.1 Algorithm

The approach to conducting matrix factorization involves minimizing ||V−WH||_F, where the Frobenius norm is defined as $||A{||}_F=\sqrt{\sum _{i,j}{A}_{ij}^2}$, and A_ij is the (i, j) element of A. To achieve minimization, NBMF employs an iterative alternative update procedure as follows:

$\begin{array}{ll}W=\underset{X\in {ℝ}_{+}^{n\times k}}{\text{arg min}}\left(\Vert V-XH{\Vert }_F^2+{\lambda }_1\Vert X{\Vert }_F^2\right),& (2)\end{array}$

$\begin{array}{ll}H=\underset{X\in {\{0,1\}}^{k\times m}}{\text{arg min}}\left(\Vert V-WX{\Vert }_F^2+{\lambda }_2\Vert X{\Vert }_F^2\right),& (3)\end{array}$

where X is a matrix that corresponds to W and H in Equations 2, 3, respectively. Hyperparameters λ₁ and λ₂ are positive.

Matrix W is updated row-by-row. The objective function for the row vector x^⊤ of W is given by

$\begin{array}{ll}{f}_W(x)=\frac{1}{2}||v-H^{\top }x{||}^2+\frac{{\lambda }_1}{2}||x{||}^2,& (4)\end{array}$

where v^⊤ is the corresponding row vector in matrix V. We applied the PGM to minimize the objective function for each row vector, as detailed in Section 3.2. The PGM was executed using a general-purpose computer. In contrast, matrix H is updated column-by-column. The objective function for optimizing the column vector q (∈{0, 1}^k) of H is given by

$\begin{array}{ll}{f}_H(q)=\frac{1}{2}||u-Wq{||}^2+\frac{{\lambda }_2}{2}||q{||}^2,& (5)\end{array}$

where u is the corresponding column vector in matrix V. To minimize the objective function for each column, we employed an SB-based Ising machine, as Equation 5 can be reformulated in the Ising model form (see Section 3.4 for details).

The overall flow of NBMF is illustrated in Figure 1. The process of updating matrix W, followed by the update of matrix H, was repeated for 10 iterations in this paper.

Figure 1

Figure 1. Overall flow of NBMF. Matrix Ṽ is an approximation of matrix V.

In this study, we compared NBMF with NMF. In NMF, Equations 2, 3 are also used; however, X∈{0, 1}^k×m in Equation 3 is substituted by $X\in {ℝ}_{+}^{k\times m}$. Furthermore, each column vector q in Equation 5 is nonnegative. Equations 4, 5 were minimized using the PGM in NMF, and the computation was executed on a general-purpose processor (a CPU).

3.2 Projected gradient method

The PGM (Lin, 2007) for updating matrix W minimizes Equation 4, and the gradient is given by

$\begin{array}{ll}\nabla f_W=-H(v-H^{\top }x)+{\lambda }_{1}x.& (6)\end{array}$

The update rule for x is given by

$\begin{array}{ll}{x}^{t+1}=P\left[x^t-{\gamma }_t\nabla f_W(x^t)\right],& (7)\end{array}$

where the projection is defined as

$\begin{array}{ll}P[x_i]=\{\begin{array}{ll}0\hfill & (x_i\le 0),\hfill \\ x_i\hfill & (0<x_i<x_{\text{max}}),\hfill \\ x_{\text{max}}\hfill & (x_{\text{max}}\le x_i).\hfill \end{array}& (8)\end{array}$

In this study, we set x_max = 1 as the upper bound of x_i. The learning rate γ_t was adjusted at each step t to satisfy the following inequality:

$\begin{array}{ll}{f}_W(x^{t+1})-f_W(x^t)\le \sigma \nabla f_W{(x^t)}^{\top }(x^{t+1}-x^t),& (9)\end{array}$

where σ = 0.01 in our experiments. Initially, we assigned γ_t−1 to γ_t (γ₀ = 1). If γ_t satisfies Equation 9, it is repeatedly divided by β, where we set β = 0.1 in our experiments, while the inequality holds. If γ_t does not satisfy Equation 9, it is repeatedly multiplied by β until the inequality is satisfied. Following this adjustment, we calculated x^t+1 using Equation 7. This procedure is repeated until ||x^t+1−x^t||≪ϵ, where ϵ = 10⁻⁷ in our experiments.

3.3 Masking procedure

Given that the rating matrix is typically sparse, the handling of unrated entries has a significant impact on the performance of collaborative filtering. A straightforward approach is to assign a rating of zero to unrated entries, which is a simple and practical choice. Another method for handling unrated entries is to introduce a mask matrix of the same size as matrix V after assigning them a zero rating. The elements of the mask matrix M are defined as follows:

$\begin{array}{ll}{M}_{ij}=\{\begin{array}{ll}1\hfill & (V_{ij}\ne 0),\hfill \\ 0\hfill & (V_{ij}=0).\hfill \end{array}& (10)\end{array}$

For collaborative filtering, we propose a modified NBMF method in which the masked matrix is decomposed as

$\begin{array}{ll}MᵒV\approx Mᵒ(WH),& (11)\end{array}$

where ° denotes the Hadamard product (MᵒV)_ij = M_ijV_ij.

In the modified NBMF algorithm, the objective function for updating matrix W, as defined by Equation 4, is replaced with

$\begin{array}{ll}{f}_W(x)=\frac{1}{2}||\tilde{v}-{\tilde{H}}^{\top }x{||}^2+\frac{{\lambda }_1}{2}||x{||}^2.& (12)\end{array}$

When updating the ith row, the jth element consists of ṽ_j = M_ijV_ij and ${({\tilde{H}}^{\top }x)}_j=\sum _l{M}_{ij}{H}_{lj}{x}_l$. Similarly, the objective function for updating matrix H, as expressed in Equation 5, is replaced by

$\begin{array}{ll}{f}_H(q)=\frac{1}{2}||\tilde{u}-\tilde{W}q{||}^2+\frac{{\lambda }_2}{2}||q{||}^2.& (13)\end{array}$

When updating the jth column, the ith element consists of ũ_i = M_ijV_ij and ${(\tilde{W}q)}_i=\sum _l{M}_{ij}{W}_{il}{q}_l$.

3.4 Ising formulation

The Ising machine (the SB-based machine in this study) seeks spin configurations that minimize the energy of the Ising model defined by

$\begin{array}{ll}E=-\frac{1}{2}{\displaystyle \sum _{i,j}}{J}_{ij}{s}_i{s}_j+{\displaystyle \sum _i}{h}_i{s}_j.& (14)\end{array}$

Here, s_i = ±1 represents the ith spin, J_ij is the coupling coefficient between the ith and jth spins, and h_i is the local field on the ith spin. For minimizing Equation 5, J_ij and h_i are given as follows:

$\begin{array}{ll}{J}_{ij}=\{\begin{array}{ll}-\frac{1}{2}{\displaystyle {\sum }_r{W}_{ri}}{W}_{rj}\hfill & (i\ne j),\hfill \\ 0\hfill & (i=j),\hfill \end{array}& (15)\end{array}$

$\begin{array}{ll}{h}_i=\frac{1}{2}\left({\displaystyle \sum _r{W}_{ri}}\left({\displaystyle \sum _j{W}_{rj}}-2u_r\right)+{\lambda }_2\right).& (16)\end{array}$

For minimizing Equation 13 to update the jth row of Mᵒ(WH), W_rl in Equations 15, 16 are replaced with ${(\tilde{W})}_{rl}=M_{rj}{W}_{rl}$.

3.5 Simulated-bifurcation-based Ising machine

The SB algorithm, which is based on the adiabatic evolution in classical nonlinear systems that exhibit bifurcation, was introduced to accelerate combinatorial optimization (Goto et al., 2019; Tatsumura et al., 2020; Goto et al., 2021). The SB algorithm has several variants, including adiabatic, ballistic, and discrete SB. In this study, we employed the ballistic SB method, whose update rule is described below (Goto et al., 2021):

$\begin{array}{ll}{y}_i(t_{k+1})=y_i(t_k)& (17)\end{array}$

$\begin{array}{ll}+\left\{-\left[a_0-a(t_k)\right]x_i(t_k)-\eta h_i+c_0{\displaystyle \sum _j}{J}_{ij}{x}_j(t_k)\right\}{\Delta }_t,& (18)\end{array}$

$\begin{array}{ll}{x}_i(t_{k+1})=x_i(t_k)+a_0{y}_i(t_{k+1}){\Delta }_t,& (19)\end{array}$

where x_i and y_i are real numbers corresponding to the ith spin; a₀, c₀, and η are positive constants; and a(t) is a control parameter that increases from zero to a₀. The time increment is Δ_t; thus, t_k+1 = t_k+Δ_t. After updating x_i at each time step, if |x_i| > 1, we replace x_i with sgn(x_i) = ±1 and set y_i = 0.

In our experiments, we employed a device with the FPGA implementation of the SB algorithm (SB-based Ising machine) to minimize Equations 5, 13. The SB-based Ising machine (Figure 2) can solve fully-connected 2,048-spin Ising problems (the machine size M is 2,048), featuring a computational precision of 32-bit floating points and a system clock frequency of 259 MHz. As shown in Figure 2a, the FPGA (Intel Stratix 10 SX 2800 FPGA) on the board (Intel FPGA PAC D5005 accelerator card) is connected to a CPU (Intel Core i9-9900K, 3.60 GHz) via a PCIe bus (Gen 3 × 16, peak bandwidth of 15.75 GB/s). The NBMF process is executed using a CPU; however, the Ising problems described in Equations 15, 16 are transferred/solved (offloaded) to/using the SB-based Ising machine. The computation times shown in Section 4 include the processing times of the CPU and FPGA and the data transfer times (overhead times) between them. The NMF process was executed only on the CPU (no data transfer time). The column update problems involved in updating matrix H (Equation 5), each formulated as an Ising problem of size k (Equations 15, 16), are independent and thus can be processed simultaneously. By packing the multiple-column update problems as a large Ising problem, as shown in Figure 2b (placing the small J matrices on the diagonal line with zero padding to the remaining off-diagonal components), we solve ⌊M/k⌋ column update problems simultaneously using the SB-based Ising machine with size M, where ⌊A⌋ is the floor function of A∈ℝ.

Figure 2

Figure 2. SB-based Ising machine. (a) System configuration. (b) Packing of multiple small Ising problems as a large Ising problem.

3.6 Data preparation

In this study, we used the MovieLens 1M dataset (Harper and Konstan, 2015); Netflix Prize data¹; Yahoo! Music user ratings of songs with artist, album, and genre meta information²; and the CiaoDVD dataset (Guo et al., 2014). These datasets were sparse, as shown in Table 1. The numbers of users and items presented in Table 1 are the dataset sizes imported for the calculation in this study. The filling rate, which is the proportion of rated entries, differs among the datasets. To compare the results of these datasets, we extracted data from them to create a rating matrix with a specified filling rate.

Table 1

Table 1. Dataset sizes (the numbers of users and items) and filling rates used in this study.

The method for extracting data at a specified filling rate is as follows. First, we sorted the columns in descending order by the percentage of filled elements in each column and then sorted the rows similarly. Next, we selected an n×m matrix whose (1, 1) element coincides with the first-row and first-column element of the sorted table, and calculated the filling rate of the matrix. By shifting the (1, 1)-element location by one row and one column in the sorted table, we repeated the calculation of the filling rate. The n×m matrix with the closest filling rate to the desired filling rate was selected as the rating matrix.

3.7 Parameter settings and evaluation

By extracting data from each dataset, we constructed a rating matrix in which 20% of the elements were rated unless otherwise specified. The numbers of users (rows) and items (columns) in the rating matrix are n = 250 and m = 500, respectively, and the number of features is set to half the number of users, that is, k = n/2, unless otherwise specified. For the learning process, which involved the execution of NBMF/NMF, we concealed 20% of the rated elements together with the unrated ones. To evaluate the performance, we used the root mean squared error (RMSE) of the rated elements:

$\begin{array}{ll}\sqrt{\frac{1}{\left|\mathcal{D}\right|}{\displaystyle \sum _{(i,j)\in \mathcal{D}}}{(v_{ij}-r_{ij})}^2},& (20)\end{array}$

where $\mathcal{D}$ is the set of rated elements, and $\left|\mathcal{D}\right|$ is the number of rated elements. v_ij is user i's rating for item j, and r_ij denotes the corresponding predicted rating.

We set the hyperparameters in Equations 2, 3 as ${\lambda }_1=10^{-2}$ and ${\lambda }_2=10^{-5}$, which were tuned for the MovieLens dataset using a grid search. Parameters in Equations 9, 17 were also tuned for the case of the MovieLens dataset. Although the optimal values may depend on the dataset and matrix size, we used the fixed values for simplicity.

In our experiments, we applied NBMF and NMF to the same rating matrix. To ensure a comprehensive evaluation, we divided the rated elements into five distinct sets and performed five trials, masking one set at a time. The average was calculated for five trials unless otherwise specified.

4 Results and discussion

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.

Figure 3

Figure 3. RMSE and computation time at each epoch, averaged over five trials, for (a) MovieLens, (b) Netflix, (c) Yahoo, and (d) CiaoDVD datasets. The error bars denote the standard deviation.

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.

Figure 4

Figure 4. Filling rate dependence of RMSE for (a) MovieLens, (b) Netflix, (c) Yahoo, and (d) CiaoDVD datasets. Training data with filling rates of more than 45%, 25%, and 20% could not be extracted for (b–d), respectively. The data were averaged over five trials, with error bars denoting the standard deviation.

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.

Figure 5

Figure 5. RMSE and computation time at each epoch with and without the masking procedure for NBMF for (a) MovieLens, (b) Netflix, (c) Yahoo, and (d) CiaoDVD datasets. The data were averaged over five trials, with error bars denoting the standard deviation.

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.

Figure 6

Figure 6. Frequency distributions of ratings expressed as percentages for (a) MovieLens, (b) Netflix, (c) Yahoo, and (d) CiaoDVD datasets.

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.

Figure 7

Figure 7. Computation time for 10 epochs in NBMF and NMF methods for the (a) MovieLens, (b) Netflix, (c) Yahoo, and (d) CiaoDVD datasets. The blue and green bars represent the total times spent on updating matrices W and H, respectively.

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.

Figure 8

Figure 8. RMSE as a function of the ratio of the number of features k to that of users n for (a) MovieLens, (b) Netflix, (c) Yahoo, and (d) CiaoDVD datasets. The data were averaged over five trials, with error bars denoting the standard deviation.

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.

5 Conclusions

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.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions

YT: Formal analysis, Investigation, Writing – review & editing. YI: Formal analysis, Writing – review & editing, Investigation. YH: Investigation, Writing – review & editing, Software. KT: Software, Investigation, Writing – review & editing. KK: Supervision, Writing – original draft, Funding acquisition, Writing – review & editing.

Funding

The author(s) declare that financial support was received for the research and/or publication of this article. This study was partially supported by JSPS KAKENHI (grant numbers JP23H04499 and JP25H01522) and Murata Science and Education Foundation.

Conflict of interest

KT was an inventor on a Japanese patent application related to this work filed by the Toshiba Corporation (no. P2019-164742, filed 10 September 2019). This study was conducted as collaborative research between Ochanomizu University and Toshiba Corporation.

YH and KT were employed at the Toshiba Corporation.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that no Gen AI was used in the creation of this manuscript.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

1. ^Netflix Prize data, https://www.kaggle.com/datasets/netflix-inc/netflix-prize-data.

2. ^Verizon Media Webscope, Yahoo! Music User Ratings of Songs with Artist, Album, and Genre Meta Information, v. 1.0, https://library.mcmaster.ca/data/yahoo-webscope-program-datasets.

References

Aramon, M., Rosenberg, G., Valiante, E., Miyazawa, T., Tamura, H., Katzgraber, H. G., et al. (2019). Physics-inspired optimization for quadratic unconstrained problems using a digital annealer. Front. Phys. 7:48. doi: 10.3389/fphy.2019.00048

Crossref Full Text | Google Scholar

Asaoka, H., and Kudo, K. (2020). Image analysis based on nonnegative/binary matrix factorization. J. Phy. Soc. Jpn. 89:085001. doi: 10.7566/JPSJ.89.085001

PubMed Abstract | Crossref Full Text | Google Scholar

Asaoka, H., and Kudo, K. (2023). Nonnegative/binary matrix factorization for image classification using quantum annealing. Sci. Rep. 13:16527. doi: 10.1038/s41598-023-43729-z

PubMed Abstract | Crossref Full Text | Google Scholar

Golden, J., and O'Malley, D. (2021). Reverse annealing for nonnegative/binary matrix factorization. PLoS ONE 16:e0244026. doi: 10.1371/journal.pone.0244026

PubMed Abstract | Crossref Full Text | Google Scholar

Goto, H., Endo, K., Suzuki, M., Sakai, Y., Kanao, T., Hamakawa, Y., et al. (2021). High-performance combinatorial optimization based on classical mechanics. Sci. Adv. 7:eabe7953. doi: 10.1126/sciadv.abe7953

PubMed Abstract | Crossref Full Text | Google Scholar

Goto, H., Tatsumura, K., and Dixon, A. R. (2019). Combinatorial optimization by simulating adiabatic bifurcations in nonlinear Hamiltonian systems. Sci. Adv. 5:eaav2372. doi: 10.1126/sciadv.aav2372

PubMed Abstract | Crossref Full Text | Google Scholar

Guo, G., Zhang, J., Thalmann, D., and Yorke-Smith, N. (2014). “ETAF: an extended trust antecedents framework for trust prediction," in Proceedings of the 2014 International Conference on Advances in Social Networks Analysis and Mining (ASONAM) (Beijing IEEE). doi: 10.1109/ASONAM.2014.6921639

Crossref Full Text | Google Scholar

Harper, F. M., and Konstan, J. A. (2015). The movielens datasets: history and context. ACM Trans. Interact. Intell. Syst. 5:19. doi: 10.1145/2827872

Crossref Full Text | Google Scholar

Herlocker, J. L., Konstan, J. A., and Riedl, J. (2000). “Explaining collaborative filtering recommendations," in Proceedings of the 2000 ACM Conference on Computer Supported Cooperative Work, CSCW '00 (New York, NY: Association for Computing Machinery), 241–250. doi: 10.1145/358916.358995

Crossref Full Text | Google Scholar

Hidaka, R., Hamakawa, Y., Nakayama, J., and Tatsumura, K. (2023). Correlation-diversified portfolio construction by finding maximum independent set in large-scale market graph. IEEE Access 11, 142979–142991. doi: 10.1109/ACCESS.2023.3341422

Crossref Full Text | Google Scholar

Inagaki, T., Haribara, Y., Igarashi, K., Sonobe, T., Tamate, S., Honjo, T., et al. (2016). A coherent Ising machine for 2000-node optimization problems. Science 354, 603–606. doi: 10.1126/science.aah4243

PubMed Abstract | Crossref Full Text | Google Scholar

Johnson, M. W., Amin, M. H. S., Gildert, S., Lanting, T., Hamze, F., Dickson, N., et al. (2011). Quantum annealing with manufactured spins. Nature 473, 194–198. doi: 10.1038/nature10012

PubMed Abstract | Crossref Full Text | Google Scholar

Kitai, K., Guo, J., Ju, S., Tanaka, S., Tsuda, K., Shiomi, J., et al. (2020). Designing metamaterials with quantum annealing and factorization machines. Phys. Rev. Res. 2:013319. doi: 10.1103/PhysRevResearch.2.013319

Crossref Full Text | Google Scholar

Lee, D. D., and Seung, H. S. (1999). Learning the parts of objects by non-negative matrix factorization. Nature 401, 788–791. doi: 10.1038/44565

PubMed Abstract | Crossref Full Text | Google Scholar

Lee, D. D., and Seung, H. S. (2000). “Algorithms for non-negative matrix factorization," in Proceedings of the 13th International Conference on Neural Information Processing Systems, NIPS'00 (Cambridge, MA: MIT Press), 535–541.

Google Scholar

Lee, J., Bengio, S., Kim, S., Lebanon, G., and Singer, Y. (2014). “Local collaborative ranking," in Proceedings of the 23rd International Conference on World Wide Web, WWW '14 (New York, NY: Association for Computing Machinery), 85–96. doi: 10.1145/2566486.2567970

Crossref Full Text | Google Scholar

Lin, C. J. (2007). Projected gradient methods for nonnegative matrix factorization. Neural Comput. 19, 2756–2779. doi: 10.1162/neco.2007.19.10.2756

PubMed Abstract | Crossref Full Text | Google Scholar

McMahon, P. L., Marandi, A., Haribara, Y., Hamerly, R., Langrock, C., Tamate, S., et al. (2016). A fully programmable 100-spin coherent Ising machine with all-to-all connections. Science 354, 614–617. doi: 10.1126/science.aah5178

PubMed Abstract | Crossref Full Text | Google Scholar

Nath, R. K., Thapliyal, H., and Humble, T. S. (2021). A review of machine learning classification using quantum annealing for real-world applications. SN Comput. Sci. 2:365. doi: 10.1007/s42979-021-00751-0

Crossref Full Text | Google Scholar

O'Malley, D., Vesselinov, V. V., Alexandrov, B. S., and Alexandrov, L. B. (2018). Nonnegative/binary matrix factorization with a D-Wave quantum annealer. PLoS ONE 13:e0206653. doi: 10.1371/journal.pone.0206653

PubMed Abstract | Crossref Full Text | Google Scholar

Pierangeli, D., Marcucci, G., and Conti, C. (2019). Large-scale photonic Ising machine by spatial light modulation. Phys. Rev. Lett. 122:213902. doi: 10.1103/PhysRevLett.122.213902

PubMed Abstract | Crossref Full Text | Google Scholar

Singh, R., Dwivedi, P., and Kant, V. (2024). Comparative analysis of collaborative filtering techniques for the multi-criteria recommender systems. Multimed. Tools Appl. 83, 64551–64571. doi: 10.1007/s11042-024-18164-5

PubMed Abstract | Crossref Full Text | Google Scholar

Su, X., and Khoshgoftaar, T. M. (2009). A survey of collaborative filtering techniques. Adv. Artif. Intell. 2009:e421425. doi: 10.1155/2009/421425

Crossref Full Text | Google Scholar

Tatsumura, K., Hidaka, R., Yamasaki, M., Sakai, Y., and Goto, H. (2020). “A currency arbitrage machine based on the simulated bifurcation algorithm for ultrafast detection of optimal opportunity," in 2020 IEEE International Symposium on Circuits and Systems (ISCAS) (Seville: IEEE), 1–5. doi: 10.1109/ISCAS45731.2020.9181114

Crossref Full Text | Google Scholar

Willsch, D., Willsch, M., De Raedt, H., and Michielsen, K. (2020). Support vector machines on the D-wave quantum annealer. Comput. Phys. Commun. 248:107006. doi: 10.1016/j.cpc.2019.107006

Crossref Full Text | Google Scholar

Yamamoto, K., Ando, K., Mertig, N., Takemoto, T., Yamaoka, M., Teramoto, H., et al. (2020). “Statica: a 512-spin 0.25M-weight full-digital annealing processor with a near-memory all-spin-updates-at-once architecture for combinatorial optimization with complete spin-spin interactions," in 2020 IEEE International Solid-State Circuits Conference- (ISSCC) (San Francisco, CA: IEEE), 138–140. doi: 10.1109/ISSCC19947.2020.9062965

Crossref Full Text | Google Scholar

Yamaoka, M., Yoshimura, C., Hayashi, M., Okuyama, T., Aoki, H., Mizuno, H., et al. (2016). A 20k-spin Ising chip to solve combinatorial optimization problems with CMOS annealing. IEEE J. Solid-State Circ. 51, 303–309. doi: 10.1109/JSSC.2015.2498601

PubMed Abstract | Crossref Full Text | Google Scholar