Algorithms and Applications to Weighted Rank-one Binary Matrix Factorization

4.1. Relation with Other Problems

The weighted rank-one BMF problem can be related to the maximum edge biclique problem [Peeters 2003], the maximum tile problem [Geerts et al. 2004], and the maximum edge weight biclique problem [Dawande et al. 1997].

A weighted rank-one BMF instance is defined as follows: Given a binary matrix A, weight parameters of w₀₁ and w₁₀, and an integer t, does there exist a pattern Y, such that ${\sum }_{i=1}^{m}I\left(A_{i.},Y\right)\ge t$?

Definition 4.1 (Maximum Edge Biclique). Given a bipartite graph G = (V₁ ∪ V₂, E) and a positive integer k, does G contain a biclique (a kind of bipartite graph where every vertex of the first set is connected to every vertex of the second set) with at least k edges [Peeters 2003]?

Statement 1. For any maximum edge biclique instance {(V₁ ∪ V₂, E), k}, we can construct an equivalent weighted rank-one BMF instance {A,w₀₁,w₁₀,t}, such that the answer to the maximum edge biclique instance is true, if and only if the answer to the weighted rank-one BMF instance is true.

The construction is as follows: (i) Let A be a m × n binary matrix, such that m = |V₁|, n = |V₂|, and A_ij = 1 if there is an edge in E, connecting vertices V₁(i) and V₂(j); (ii) let w₁₀ = 0 and w₀₁ ≥ max_i ‖A_i. ‖₁; (iii) let t be k. The mapping makes A be a binary matrix representation of the bipartite graph (V₁ ∪ V₂, E), as the vertices in V₁ correspond to rows in A, the vertices in V₂ correspond to columns in A, and each edge in E corresponds to a 1-entry in A. As w₀₁ ≥ max_i ‖A_i. ‖₁, 0-becoming-1 errors are penalized so severely that a pattern cannot cover 0-entries.So a pattern present in A naturally induces a biclique in (V₁ ∪ V₂, E). As w₁₀ = 0, 1-becoming-0 errors are not counted toward the pattern influence measure; in other words, the pattern influence measure is equal to the number of edges of the induced biclique. Therefore, the above statement is true.

Figure 1 illustrates the relation. The bipartite graph is the graph representation of the binary matrix on the left side in example 1.1. The nodes on the left side correspond to matrix rows, the nodes on the right side correspond to matrix columns, and edges (including both solid and dashed lines) correspond to 1’s cells in the matrix. The edges with solid lines correspond to the maximum edge biclique. It also corresponds to the optimal solution to the weighted rank-one BMF with w₁₀ being 0 and w₀₁ being a large enough number.

Definition 4.2 (Maximum Tile). Given a database D as a binary matrix A_m×n, find the tile with the largest area in D, where a tile corresponds to a row index subset R and a column index set C, such that A(i, j) = 1, ∀i ∈ R, j ∈ C [Geerts et al. 2004].

The maximum tile problem is identical to the maximum edge biclique problem, except that one is defined in terms of a matrix, whereas the other one is defined on an equivalent graph representation. Therefore the maximum tile problem can be reduced to the weighted rank-one BMF problem as well.

Definition 4.3 (Maximum Edge Weight Biclique Problem). Given a bipartite graph {V₁ ∪ V₂, E} and weights {w_ij} associated with edges {(V₁(i),V₂(j))}, respectively, find a biclique C, where the sum of the edge weights is maximum [Dawande et al. 1997].

Statement 2. The weighted rank-one BMF problem with w₁₀ of 0 is a special case of the maximum edge weight biclique problem.

For each instance of the weighted rank-one BMF problem with w₁₀ of 0, one can construct an equivalent instance of the maximum edge weight biclique problem. An instance of the weighted rank-one problem can be denoted by {A_m×n,w₀₁}. With w₁₀ being 0, the pattern influence measure changes to I (X,Y) = max{0, f₁₁(X,Y) − w₀₁f₀₁ (X,Y)}. So the weighted rank-one BMF problem becomes to find a pattern to cover the most 1-entries, while there is a penalty of w₀₁ for the coverage of each 0-entry. An instance of the maximum edge weight biclique problem can be denoted by {V₁ ∪ V₂, E,W}, where W is the weights of edges. An equivalent instance of the maximum edge weight biclique problem can be constructed as follows: (i) Let V₁ have m vertices and V₂ have n vertices; (ii) let {V₁ ∪ V₂, E} be a complete bipartite graph, i.e., ∀V₁(i) ∈ V₁,V₂(j) ∈ V₂, there is an edge connecting them; (iii) let W (i, j), the weight for the edge connecting nodes V₁(i) and V₂(j), be 1 if A(i, j) is 1; otherwise, −w₁₀. The construction maps each entry in A_m×n to an edge in the complete bipartite graph. As W (i, j) is 1 for A(i, j) of 1 and −w₁₀ for other edges, the maximum edge weight for the constructed biclique instance is equal to the maximum pattern influence value for the BMF instance. Therefore, the two instances are equivalent.

4.2. Optimization Models

In this section, we build optimization models, which allow us to take advantage of the latest advancement of optimization techniques.

4.3. Computational Hardness

The rank-one BMF problem is NP-hard, because the NP-hard maximum edge biclique problem can be reduced to the rank-one BMF problem [Peeters 2003].

4.4. Approximability

There exist 2-approximation algorithms for special cases of the weighted rank-one BMF problem.

Case 1: w₁₀ = 1 and w₀₁ ≥ max_i ‖A_i. ‖₁

The 2-approximation algorithm for such a case is based on Lemma 4.5.

Lemma 4.5. Any integer linear programming problem with n variables subject to m linear constraints with at most two variables per inequality, and with all variables bounded between 0 and U, has a 2-approximation algorithm that runs in polynomial time O(mnU ²log(Un²m)) [Hochbaum et al. 1992].

The weighted rank-one BMF problem with w₁₀ = 1 and w₀₁ ≥ max_i ‖A_i. ‖₁ can be formulated as the integer linear programming problem (20)–(24) with two variables per inequality as follows:

When w₁₀ = 1 and w₀₁ ≥ max_i ‖A_i. ‖₁, the weighted rank-one BMF problem is to find a pattern to cover the most 1-entries without covering one 0-entry. So the objective is to maximize Σ_{(i,j) |Ai j =1} Z_ij, where Z_ij = X_iY_j. Constraints (21) and (22) are used to enforce Z_ij = X_iY_j. As explained earlier, constraints −X_i − Y_j + 2Z_ij ≤ 0 and X_i + Y_j − Z_ij ≤ 1 enforce Z_ij = X_iY_j. But when the objective is to maximize Σ_{(i,j) |Ai j =1} Z_ij, the constraint X_i + Y_j − Z_ij ≤ 1 becomes redundant and can be dropped. Constraints Z_ij − X_i ≤ 0 and Z_ij − Y_j ≤ 0 are equivalent to the constraint −X_i −Y_j + 2Z_ij ≤ 0. Constraints (23) enforce no 0-entry being covered by the approximation.

Case 2: w₁₀ = 0 and w₀₁ = 1

The weighted rank-one BMF problem with w₁₀ = 0 and w₀₁ = 1 is equivalent to the conventional rank-one BMF problem. Conventional rank-one BMF considers observation X belongs to pattern Y, if ‖X‖₁ < ‖X − Y‖₁. As ‖X‖₁ = f₁₀(X,Y) + f₁₁(X,Y) and ‖X − Y‖₁ = f₁₀(X,Y) + f₀₁(X,Y), so ‖X‖₁ < ‖X − Y‖₁ induces f₁₁(X,Y) − f₀₁(X,Y) < 0, which is the special case of the weighted rank-one BMF problem with w₁₀ = 0 and w₀₁ = 1. A 2-approximation algorithm via linear programming relaxation is given in Shen et al. [2009].

4.5. Necessary Conditions for Observations of the Same Pattern

Weighted rank-one BMF essentially is to identify the dominant pattern with respect to given weights. Because of its inherent NP-hardness nature, it is difficult to find an efficient algorithm for weighted rank-one BMF on large datasets. However, if we can separate observations (records), which are determined not belonging to the same patterns, into disjoint groups, then we can perform weighted rank-one BMF on smaller datasets. To apply such a divide-and-conquer technique, we need an effective and efficient way to separate observations that cannot belong to the same pattern. The theorem below gives the necessary condition for two records to belong to the same pattern. In other words, if the condition is not satisfied, then two records can never belong to the same pattern. The theorem below also provides foundation for the practical iterative framework, described later in Algorithm 1 in Section 5.

Theorem 4.6. If A and B, the set representations of two binary vectors, belong to the same pattern, then the following are true:

where {z₁, z₂, z₃} denote {|A| − |A ∩ B|, |A ∩ B|, |B| − |A ∩ B|}, respectively.

However, we observe that given weights w₀₁ and w₁₀, some observations can never belong to the same pattern. So if we can, in other words, we utilize a divide-and-conquer approach. To this end, we need a property, which is easy to implement.

To prove the theorem, assume observations A and B belong to pattern C. Figure 2 illustrates the regular set intersection of the three sets A, B, and C.

According to the pattern presence definition, the following are true:

The above inequalities can be simplified as below by letting d₂ and d₇ be 0,

As such, the set intersection of A, B, and C is updated by Figure 3.

Since {z₁,z₂,z₃} denote {|A| − |A ∩ B|, |A ∩ B|, |B| − |A ∩ B|}, {z₁,z₂,z₃} equal {d₁ + d₄,d₅,d₃ + d₆}. So {d₄,d₅,d₆} equal {z₁ − d_1, z₂,z₃ − d₃}. Plug them into inequality (31) and (32) and obtain

Inequality (33) leads to

Plug it into inequality (34) and obtain

As d₁ is not less than 0, so if $w_{01}^2-{\left(1+w_{10}\right)}^2\le 0$, in order for the above inequality to satisfy, $\left(w_{01}^2-w_{10}-1\right)z_1-\left(w_{10}+w_{01}+1\right)z_2+w_{10}{w}_{01}{z}_3$ must be less than 0. Thus, the necessary condition (25) is proven.

Inequality (34) leads to

Plug it into inequality (33) and obtain

As d₃ is not less than 0, so if $w_{01}^2-{\left(1+w_{10}\right)}^2\le 0$, in order for the above inequality to satisfy, $w_{10}{w}_{01}{z}_1-\left(w_{10}+w_{01}+1\right)z_2+\left(w_{01}^2-w_{10}-1\right)z_3$ must be less than 0. The necessary condition (26) is proven.

Reconsider inequalities (31) and (32). {d₁,d₃,d₅} equal {z₁ − d₄,z₃ − d₆,z₂}. Plug them into inequalities (31) and (32) and obtain

Inequality (39) leads to

Plug it into inequality (40) and obtain

As d₆ is not less than 0, so if ${\left(1+w_{10}\right)}^2-w_{01}^2\le 0$, $w_{10}{w}_{01}{z}_1-\left(w_{10}+w_{01}+1\right)z_2+\left(w_{10}^2+w_{10}\right)z_3$ must be less than 0. The necessary condition (27) is proven.

Inequality (39) leads to

Plug it into inequality (40) and obtain

As d₄ is not less than 0, so if ${\left(1+w_{10}\right)}^2-w_{01}^2\le 0$, then $\left(w_{10}^2+w_{10}\right)z_1-\left(w_{10}+w_{01}+1\right)z_2+w_{10}{w}_{01}{z}_3$ must be less than 0. The necessary condition (28) is proven.

5.1. Practical Iterative Framework

The practical iterative framework, described in Algorithm 1, reduces a large weighted rank-one BMF problem to a number of small problems by taking advantage of the fact that some observations can never belong to the same pattern. At each iteration, the algorithm selects one observation from the remaining dataset and then applies necessary conditions (25)–(28) to identify observations that may belong to the same pattern with the selected observation. After that, employ either heuristics or an exact algorithm (if the identified group is very small) and execute weighted rank-one BMF on the identified observations group. The observations associated with the obtained pattern are deleted from the remaining dataset. The procedure continues until no observation remains.

However, we still need to solve a weighted rank-one BMF problem of reduced size at each iteration. Our experimental study finds that for a large dataset the original weighted rank-one BMF is indeed reduced to a number of problems of much smaller sizes. But due to notoriousness of binary optimization, it is still not practical to employ an exact algorithm to solve the reduced problems. Thus, efficient heuristics for weighted rank-one BMF are required.

Weighted rank-one BMF can be modeled by unconstrained binary quadratic programming. There have been a considerable number of studies on heuristics for unconstrained binary quadratic programming model, including greedy and local search [MerZ and Freisleben 2002], tabu search [Glover et al. 1998], simulated annealing [Beasley 1998], genetic algorithm [Merz and Freisleben 1999], and exact algorithm [Barahona et al. 1989]. As the purpose of weighted rank-one BMF is to facilitate large data analysis. We implemented and tested most mainstream heuristics. Experimental results show the lightweight local search algorithm in fact has the best overall performance in terms of running time and accuracy. Local search algorithms are heuristics that search in the neighborhood of the current solution for a better one until no further improvement can be made. Algorithm 2 provides the outline of a 1-opt local search algorithm, in which the neighborhood of a solution is defined by the solutions that can be obtained by flipping a single component in the binary vector. The algorithm iteratively examines all variables. If flipping the variable value improves the objective function, then make the move; otherwise, skip it.

The local search algorithm can be implemented efficiently. First, if the input matrix A_m×n is stored in sparse form such that each row of A_m×n is represented by a subset of {1, … ,n}, then the running time is reduced considerably for instances with low density. Second, as much of the running time of the heuristics is for evaluating the gain of the next move, the overall running time can be decreased significantly by evaluating a move’s gain on the basis of the current solution. Let the current solution be Y. Its influence on each observation A_i. is as follows:

Consider the move of flipping the value of the kth element of Y from 0 to 1. To evaluate the gain of the move, we only need to check whether A_i. contains the element k. If true, then the move increases f₁₁(A_i.,Y) by one and decreases f₁₀(A_i.,Y) by one. If false, then the move increases f₀₁(A_i.,Y) by one. Suppose the move is to flip the value of the kth element of Y from 1 to 0. If A_i. contains the element k, then the move decreases f₁₁(A_i.,Y) by one and increases f₁₀(A_i.,Y) by one. Otherwise, the move decreases f₀₁(A_i.,Y) by one. It takes O(log₂(n)) to check the existence of an element in an ordered list of n elements by binary search. So evaluation of a move’s gain takes O(mlog₂(n)).

5.2. Quick-and-Dirty

A quick-and-dirty solution is simply to randomly select one observation from the remaining dataset and take it as a tentative pattern at each iteration. Then, based on penalty weights identify observations belonging to the tentative pattern. Finally, take the rounded mean of the group of identified observations as a pattern. As such, rank-one BMF is avoided, while penalty weights are utilized. This random algorithm can easily be upgraded to be more robust. In this variant, at each iteration we randomly select log(n) observations as tentative patterns, where n is the number of total observations, and pick the one that has the most observations belonging to it. The variant would incur more computing time than the simple random algorithm but less time than the sophisticated algorithm employing BMF. The experimental results presented later show that the quick-and-dirty solutions have decent performance in terms of data compression, clustering, and pattern discovery and can be the preferred choice for extremely large datasets.

7.1. Effectiveness of Clustering and Pattern Discovery on Synthetic Data

This experiment is to compare weighted rank-one BMF with the state-of-the-art approaches for clustering and pattern discovery for binary data that include SVD [Andrews and Patterson 1976], k-means [Tan et al. 2005], conventional rank-one BMF [Koyuturk et al. 2005], denoted by conventional in this section, and Discrete Basis Problem (DBP) [Pauli et al. 2008], with respect to their effectiveness. We adopt the same synthetic data generator and experiment design employed in Koyuturk et al. [2005]. We generate 200 binary vectors of 80 dimensions by first sampling 5 binary vectors each with 50 components being 1 as underlying patterns, repeating each 40 times and then flipping the values of 5% entries randomly chosen in those records. The generated data can be represented by a 200 × 80 binary matrix and is plotted in white-black as shown in Figure 4(a), where a black point indicates an entry with the value of 1.

We run each approach mentioned above to extract the underlying patterns and find clusters by grouping records belonging to the same pattern. A good approach should discover all five underlying patterns and their associated clusters. We utilize the packages in Matlab to perform SVD and k-means, use the code in Pauli et al. [2008] to execute DBP, and run the algorithm described in Koyuturk et al. [2005] to implement conventional rank-one BMF. To carry out SVD, k-means, and DBP, the number of clusters (or matrix rank) is needed. Indeed, this experiment is designed to favor those approaches and assumes the number 5 of true clusters already available to those approaches, while in practice a prior step (e.g., some information criterion-based approach) is typically required to learn the number of clusters.

We first apply the five-way k-means clustering to discover five clusters, approximate each record by its associative cluster centroid, and plot the resultant approximation data in Figure 4(b). The result confirms the issue of k-means that was raised in Koyuturk et al. [2005] that k-means is not able to capture patterns with significant overlap. In Figure 4(a), the first two patterns are highly overlapped. However, k-means fails to identify them and falsely groups many records of the first two clusters into one big cluster. This example in fact manifests the necessity of applying weights on errors. If we increase penalty on discrepancy, then the records in the first two clusters would never be considered from the same cluster. Another issue with k-means that has been discussed in the literature is that it is highly susceptible to initial solutions and hence not very robust.

We run the algorithm described in Pauli et al. [2008] to execute the DBP to discover discrete patterns. To favor the DBP, the execution assumes that it is known there are five discrete patterns and each record is derived from one pattern. DBP is an algorithm that decomposes a binary matrix into the Boolean product of two binary matrices. One decomposed binary matrix can be interpreted as discrete basis or (patterns) while the other describes the data as a formation of the discrete basis. The DBP algorithm first constructs a set of candidate basis vectors, based on the confidence concept in association rule mining. A parameter r controls the confidence level. Then, a small set of candidate basis vectors are selected in a greedy way. Two weight parameters w⁺ and w⁻ control the selection of candidate basis vectors. The algorithm details can be found in Pauli et al. [2008]. We run the algorithm multiple times by varying r from 0.4 to 0.9 and w⁻ from 1 to 10 with w⁺ being 1 to select the best set of five basis vectors, and associate each record with the closest basis vector. We represent each record with its corresponding basis vector and plot the resultant approximation data in Figure 4(c). As it shows, the DBP fails to discover all five underlying patterns, and the resultant approximation significantly differs from the original data by 3,123 errors. One possible reason is that DBP is suitable for the pattern discovery problem where each record is generated from multiple patterns, while our synthetic data generator samples each record from only one pattern.

We run the SVD to find the rank-5 approximation to the original data and plot the approximation matrix as grayscale image in Figure 4(d). The SVD approximation is optimal in the sense of minimum least squares, with an approximation error of 29.78 in this example. However, the SVD clearly has the interpretability issue, as the retrieved patterns contain fractional and negative numbers, while the five true patterns are binary. For a fairer comparison, we quantize the result of SVD. We can quantize either the resultant approximation or the singular vectors. As analyzed and reported in Koyuturk et al. [2005], quantizing singular vectors obtains approximations that significantly deviate from the original data. So we only consider quantizing the resultant approximation. We round a value in the SVD approximation to 1 if it is greater than a threshold value r, else 0. We try different values of r and pick the one that minimizes the difference between the quantized result to the original data. The best result is depicted in Figure 4(e). The approximate matrix is very close to the original data without noise. However, it requires 10 patterns to describe the original data, of which five are false patterns.

We run the algorithm described in Koyuturk et al. [2005] to implement conventional rank-one BMF. The conventional rank-one BMF is to iteratively discover a binary presence vector X and a binary pattern vector Y that maximize ‖A − XY‖₂. It can be viewed as a special case of our weighted rank-one BMF with both error weights w₁₀ and w₀₁ of 1. We apply the conventional rank-one BMF to discover patterns, represent each record with their associative pattern and plot the approximation data in Figure 4(f). Two false patterns are retrieved, and the majority records are falsely clustered into one group. The reason that the conventional rank-one BMF is unable to discover the true patterns is the tolerance of errors, so it cannot distinguish patterns with large overlap as this example. If we increase penalty on differences, then more clusters would be revealed.

We finally apply our weighted rank-one BMF. Similarly to the conventional rank-one BMF, the weighted rank-one BMF is to iteratively discover a binary vector with the maximal pattern influence as defined in Equation (3), treat the discovered binary vector as an underlying pattern, and group the records with the presence of the pattern into a cluster. As it is a NP-hard problem to find the binary vector with the maximal pattern influence, we implement the weighted rank-one BMF approach with the three different heuristics proposed in Section 5. They are (1) local, which uses the local search heuristic to discover a pattern guaranteed to be a local optimum solution at each iteration; (2) random, where each iteration randomly selects a data record, groups its neighboring records into a cluster, and makes the rounded average of the cluster as pattern; and (3) upgraded random, which upgrades the random algorithm by repeating it log(n) times (n is the number of total records) random records as tentative patterns and then chooses the pattern with the most neighboring records. It is expected that local takes the most running time with the best performance, random takes the least running time at the cost of losing effectiveness, and upgraded random is in between the two. To execute the weighted rank-one BMF, we need to know the weights w₁₀ and w₀₁. Indeed, k-means, DBP, and SVD also need additional information, i.e., the number of clusters or rank, before they can be carried out. To favor those approaches, in the prior experiments, we simply assume the information is known to them. According to many experiments we conduct on the weighted rank-one BMF, we think it is reasonable to assume weight values in [1, 4]. For further convenience, we limit weight values to be increments of 0.5, i.e., {1.0, 1.5, … , 4.0}. That leaves 49 possible combinations of weight values. To determine the concrete weight values, we use the 10-fold cross-validation, a common approach for model selection, as discussed in Section 6. For each set of weight values, we partition the 200 records into 10 equal-size groups, use 9 groups as training data, apply the weights to learn patterns, associate each record in the remaining group with one of retrieved pattern that has the minimum Euclidean distance, and compute the approximation errors (i.e., the sum of Euclidean distances between records and their associative patterns). We repeat the above procedure 10 times by retaining a different group as testing data and calculate the average approximation errors. The set of weights that leads to the minimum average approximation errors is chosen. That is w₁₀ of 2 and w₀₁ of 2 in this case. Their induced patterns and clusters are the result of our weighted rank-one BMF approach. We represent each record with their associative pattern and plot the resultant approximations in Figure 4(g)–(i). The results show that all three implementations of the weighted rank-one BMF are able to discover the five true patterns and discover true underlying clusters. The experiment clearly demonstrates the effectiveness of weighted rank-one BMF in clustering and pattern discovery over other approaches.

7.2. Robustness of Clustering and Pattern Discovery on Synthetic Data

This experiment examines the robustness of the weighted rank-one BMF approach to noise, since robustness is a desired property for clustering and pattern discovery. The prior experiment has demonstrated the limitation of k-means in discerning patterns and clusters with large resemblance and the issue of SVD of discovering discrete patterns. So in this experiment, we compare our approach only with the conventional rank-one BMF and DBP. We use the same synthetic data generator to construct two synthetic datasets. One is added with 10% noise as plotted in Figure 5(a) and the other with 20% noise plotted in Figure 6(a).

The results of conventional are reported in Figures 5(b) and 6(b). For the data with 10% noise, conventional retrieves three false patterns with no resemblance to any true pattern and groups the majority records into one big cluster, as shown in Figure 5(b). The result on the data with 20% noise is not good either. Figure 6(b) shows that one true pattern and two false patterns are retrieved, and it fails to discover the five true clusters.

The results of DBP are reported in Figures 5(c) and 6(c). For DBP, we again assume the number of true clusters is known. By applying the rank-4 binary matrix decomposition, we retrieve four discrete basis vectors and represent each record with the closest discrete basis vector. For the data with 10% noise, Figure 5(c) shows that the DBP fails to discover the true patterns and the clustering result is incorrect. For the data with 20% noise, Figure 6(c) shows that the DBP performs even worse. The retrieved patterns hardly resemble any true pattern and the clustering result significantly deviates the underlying truth.

To evaluate the weighted rank-one BMF, we again run the three implementations, namely local, random, and upgraded random. To determine the weights, we still use the 10-fold cross-validation as we did for the previous experiment. For the data with 10% noise, the learned weights are w₁₀ of 1 and w₀₁ of 2. For the data with 20% noise, the learned weights are w₁₀ of 1 and w₀₁ of 1. By applying the learned weights, we execute local, random, and upgraded random separately to retrieve patterns and represent each record with its cluster centroid. The approximation results on the data with 10% noise are plotted in Figure 5(d)–(f). As shown in Figure 5(d), local returns all five true underlying patterns and one false pattern. It successfully discovers the five underlying clusters with one mistake of separating one record from others and treating it as one unique cluster. As shown in Figure 5(e), random returns 15 patterns and the resultant approximation largely resembles the original data without 10% noise. The 15 clusters induced by the 15 patterns contain the five true clusters with 10 false clusters of few members that can be ruled out with some postprocessing step. Figure 5(f) shows that upgrade returns 10 patterns. The resultant approximation of upgrade is more distant from the original data than random. But upgrade returns less patterns, so in terms of precision on retrieving true patterns, upgrade performs slightly better than random. The approximation results on the data with 20% noise are plotted in Figure 6(d)–(f). As shown in Figure 6(d), local still performs well. It retrieves all five patterns and one false pattern. So in terms of effectiveness of pattern discovery, the local implementation of weighted rank-one BMF is robust to data noise. In terms of clustering, local performs slightly worse on the data with 20% noise than the data with 10% noise. Figure 6(d) shows 6 clusters with one false cluster of quite a few members, while the false cluster in Figure 5(d) has only one member. The results show that the amount of noise does affect the performance of weighted rank-one BMF. Figure 6(e) and (f) also verify the observation. Figure 6(e), random returns 32 patterns that include all five true patterns. In Figure 6(f), upgraded random returns 26 patterns that also include all true patterns. It is observed that the amount of noise does affect the performance of weighted rank-one BMF, as more data noise leads to more false patterns. It is also observed that weighted rank-one BMF performs significantly better than the other approaches. For the both noisy datasets, all three implementations of weighted rank-one BMF are able to retrieve all underlying patterns. The local implementation performs particularly well, as it only returns one false pattern. Another observation is that all three implementations of weighted rank-one BMF perform well in terms of clustering, as cluster centroids can well reconstruct the original data, while the approximations returned from conventional and DBP significantly deviate the original data. The experiment concludes that weighted rank-one BMF performs well in terms of robustness to data noise. It again verifies the advantage of weighted rank-one BMF in clustering and pattern discovery on binary data.

7.3. Running Time

This experiment studies the running time (complexity) of the weighted rank-one BMF algorithms. Specifically, we proposed two solutions. The first one, local, is the practical iterative framework combined with a local search heuristic. (Note that the iterative framework can be combined with other heuristics as well.) The second one, random, is a quick-and-dirty solution, which generates a random pattern at each iteration. The synthetic data generator works as follows. Denote the number of observations by m. First, generate 0.01 * m random binary vectors of size 200 with average density of 1-entries of 0.01 and use them as patterns. Then construct a m matrix by associating each row with a random pattern. Finally, randomly select 10% entries and flip their values. We generate six synthetic matrices with m having the following values: 10k, 20k, 40k, 60k, 80k, and 100k. Local and random with w₁₀ and w₀₁ of 1 are executed against the seven matrices. Their running times are reported in Figure 7. The running time of random grows almost linearly with data size. It takes less than 150 s for a dataset with size as large as 100k. The running time of local grows seemingly quadratically with data size. It takes about 350 s for a dataset with 100k. Both algorithms have good performance in terms of running time. However, for even larger datasets, random is preferred. It is possible to apply parallel algorithms to reduce running time by utilizing technologies like Hadoop MapReduce. Exploring that direction is a potential avenue for future work.

7.4. Application to Data Compression in Association Rule Mining

In this section, we illustrate the application of weighted rank-one BMF to data compression in the context of association rule mining. Data compression is an important topic in computer science and information theory. It involves encoding information using fewer bits than the original representation. Data compression is useful, because it helps reduce resource usage and has found numerous applications including audio transmission, video storage, wireless sensor network, and among others. A variety of data compression methods have been proposed. They are classified as either lossy or lossless, dependent on whether the original data can be reconstructed from the compressed data. By applying weighted rank-one BMF with appropriately chosen weight parameter values, one can extract patterns from the original binary data and associate patterns with data records. As one pattern typically represents many data records, so the extracted relatively few patterns and their associations to the original data records can be viewed as compressed data. Since original data records are not exactly the same as their associated patterns, the weighted rank-one BMF serves as a lossy compression method.

Data compression can benefit association rule learning. Association rule learning is a popular and well-researched method to identify strong rules in large databases. For example, the rule {milk,bread} → {butter} found in the sales data of a supermarket would indicate that if a customer buys milk and bread, he or she is likely to also buy butter. Association rules can be used as the basis for decisions about marketing activities such as promotional pricing or product placements. The core of association rule generation is to identify frequent itemsets. It is difficult to find all frequent itemsets in a database, since it involves searching all possible itemsets. Efficient algorithms, like apriori, exploit the anti-monotonicity property, which guarantees that for a frequent itemset, all its subsets are also frequent and thus for an infrequent itemset, all its supersets must also be infrequent. However, even for efficient algorithms, finding frequent itemsets from a large dataset is time-consuming.

To cope with the computational challenge in mining association rules from large datasets, data compression techniques have been used [Burdick et al. 2001; Koyuturk et al. 2005; Shenoy et al. 2000]. The general idea is to cut down data size so the mining algorithm can complete the mining task on large datasets. The conventional rank-one BMF [Koyuturk et al. 2005] is one of the recent applications of data compression in association rule mining. Its idea is to extract binary patterns from a large binary dataset and then performs the mining task on the extract patterns of much smaller size than the original data. Experimental studies [Koyuturk et al. 2005] show that frequent itemsets mined from compressed data (extracted patterns) are close to that mined from original data. However, the conventional rank-one BMF-based compression approach has one limitation. That is the size (number) of extracted patterns is unadjustable. In practice, one might want to control data compression ratio to balance efficiency and accuracy of the mining algorithm.

This experiment is designed to demonstrate that weighted rank-one BMF allows one to adjust data compression ratio and also achieves high accuracy mining results from compressed data. As what we demonstrated in the previous experiments, weighted rank-one BMF can be used to extract binary patterns from a large dataset. So a large dataset can be represented by a small number of extracted patterns and the sizes of their corresponding clusters. Then one can reduce the computational overhead of association rule mining on a large dataset by performing the mining task on the compressed data of much smaller size. As weighted rank-one BMF serves as a lossy data compression method, it is not necessary that the frequent itemsets from the compressed data are the same as those from the original data. So to validate effectiveness and flexibility of weighted rank-one BMF as a compression technique to support association rule mining, we conduct experiments on two FIMI workshop datasets, pumsb and connect.¹ The data descriptions are provided in Table 1. For each dataset, we first use the apriori algorithm, which proceeds by identifying the frequent individual items in the database and extending them to larger and larger item sets, to retrieve frequent itemsets. Then we apply weighted rank-one BMF to derive a small number of patterns and the sizes of their corresponding clusters (i.e., data compression). As both datasets are large, so we use the rand algorithm implementation for weighted rank-one BMF. Then we use the slightly modified apriori algorithm as described in Koyuturk et al. [2005] to mine frequent itemsets from the compressed dataset (i.e., a small number of patterns and the sizes of their corresponding clusters).

We compare the true frequent itemsets mined from the original dataset, denoted by A, and those retrieved from the compressed dataset, denoted by B. To evaluate a test’s accuracy, we use the F-measure. It considers both the precision and the recall of the test to compute the score. Precision is defined as $\frac{|A\cap B|}{|B|}$, where |A ∩ B| is the number of true frequent itemsets retrieved from the compressed dataset and |B| is the number of total frequent itemsets retrieved from the compressed dataset. Recall is defined as $\frac{|A\cap B|}{|B|}$, where |A| is the number of true frequent itemsets from the original dataset. F-measure is computed as: $\mathit{\text{F-measure }}=2\times \frac{\mathit{\text{precision}}\times \mathit{\text{recall}}}{\mathit{\text{precision}}+\mathit{\text{recall}}}$, which reaches its best value at 1 and worst score at 0. The compressed data consist of representative patterns and their associations to the original data records. As it takes negligible space to store the associations (one integer of pattern label for each record), so we measure the data compression ratio by $\frac{\mathit{\text{Number of Pat terns}}}{\mathit{\text{Number of Original Records}}}$.

As discussed in Section 6, greater values of weight parameters lead to lower compression and higher approximation, as more representative records are selected. But there is no deterministic relation among weight parameter values, compression, and approximation. So we explored different values and found that weight parameter values between 1 and 3 typically achieve a good balance between data compression and approximation. So we reported results of parameters varying from 1 to 3. Note that weighted rank-one BMF with w₁₀ and w₀₁ being 1 is in fact conventional rank-one BMF. The results for the pumsb transactions set are reported in Tables 2–4. We consider the frequent itemsets with support (the proportion of the transactions in a database that contain the set) over 0.7 and length (the number of items in the set) less than 6 in the original transactions set and compare them with the itemsets satisfying the same criterion, in the approximate transactions sets. From Table 2, we observe that the number of patterns required to represent the original transactions set increases as the values of w₁₀ and w₀₁ grow. The benefits of weighted rank-one BMF can be observed in Tables 2–4. The number of itemsets with the same support in the original dataset is 66,741. When w₁₀ and w₀₁ are 1, the data are compressed to about 0.003, and 184,372 frequent itemsets are retrieved from the approximate transactions set with a F-measure score of 0.5271, which could be interpreted as inaccurate. By increasing the values of w₁₀ and w₀₁, more data storage are needed to represent the original transactions set. However, the numbers are still much less than the original data size. In return, the accuracy of retrieved frequent itemsets is significantly improved. For instance, when w₁₀ and w₀₁ are 2.5, the data are compressed to one-tenth, while the retrieved frequent itemsets has high accuracy of the F-measure score of 0.9475. When w₁₀ and w₀₁ are greater than 2, the F-measure score is consistently greatly than 0.9, as opposed to 0.5271 when w₁₀ and w₀₁ are 1. We also observe that higher values of w₁₀ and w₀₁ do not necessary lead to better F-measure score, although the general trend exists. This is partially because the rand algorithm is heuristic oriented. To confirm the findings, we conduct another experiment on the connect transactions set. We compare the frequent itemsets of support over 0.45 and length less than 5 in the original transactions set with the frequent itemsets of the same characteristics derived from the approximate transactions sets. The number itemsets with the same support in the original dataset is 57,352. The results reported in Tables 5–7 are even better. As shown in Table 5, when both w₀₁ and w₁₀ increase to 3, the data are compressed to about 0.003, while the F-measure of the frequent itemsets derived from the compressed data is as high as 0.9459. The experimental results clearly demonstrate the benefits of weighted rank-one BMF as a data compression method.

Note again that the choice of weights w₁₀ and w₀₁ is application dependent. Greater values of weight parameters lead to lower compression and higher approximation, as more representative records are selected. Although our extensive experimental studies suggest weights between 1 and 3, users can always take a small data sample, try different values of weights, and choose the one that achieves their preferred compression ratio and accuracy.

We also compare our method with two recent algorithms, low-rank structured matrix factorization [Haeffele and Vidal 2019] and low-rank non-negative matrix factorization [Gligorijević et al. 2018]. Low-rank matrix factorization is commonly used for compressing matrix data. By decomposing A_m×n ≈ U_m×k × V_k×n, where k is significantly smaller than m and n, A is compressed to U and V. To reconstruct A, we compute the product of U and V and round values to 0/1. We apply the two matrix factorization methods with different ranks to the pumsb dataset. From their reconstructed binary matrices, we find frequent itemsets with support over 0.7 and compute the F-measure accuracy by comparing the returned frequent itemsets with that from the original dataset. We also apply the weighted rank-one BMF with weights of (2.5,2.5). The results are reported in Table 8. It shows that the traditional matrix factorization-based compression methods do not work well for preserving association rules for binary data, although the performance improves when the rank increases. Another limitation is that matrix decomposition methods are time-consuming for large matrices.

7.5. Application to Clustering in Information Retrieval

In this section, we illustrate the application of weighted rank-one BMF to clustering in the context of information retrieval. Information retrieval is finding material (usually documents) of an unstructured nature (usually text) that satisfies an information need from within large collections (usually stored on computers) [Manning et al. 2008]. In the world, hundreds of millions of people engage in information retrieval every day when they use a web search engine. Information retrieval is becoming the dominant form of information process. Due to its importance, information retrieval is an actively sought research area. One challenge in information retrieval is to process a large volume of data. There are a variety of methods to alleviate the computational overhead, such as indexing, data compression, and distributed computing. Clustering can benefit information retrieval in the sense that documents can be grouped into clusters such that documents in the same cluster behave similarly with respect to relevance to information needs. One application of clustering in information retrieval is search result clustering. Consider a search term jaguar of different word senses. There could be a large number of web sites relevant to the search term. A good search engine would provide clustering results, such as cars, club, and cat, to allow a user to choose what he/she is interested in.

Weighted rank-one BMF can be used as a document clustering method in information retrieval. A document can be represented as a document-term binary matrix, where 1’s cell denotes the corresponding document contains the corresponding term. By iteratively applying weighted rank-one BMF with appropriately chosen weight parameters, one can derive dominant patterns and associate patterns with documents. As each data record is linked to one pattern, documents are naturally partitioned into clusters of relevant documents. Each pattern consists of a cluster of terms, from which semantic information can be derived and fine grain classification of terms is enabled. To validate our idea, we conduct experiments on the 20 newsgroup dataset, which is originally collected by Lang [1995], and compare our method with k-means, the classic text clustering method. The 20 newsgroups collection is a popular dataset for experiments in text classification and text clustering. The data contain 20 different newsgroups, each corresponding to a different topic, e.g., comp.graphics, rec.autos, and sci.crypt. Some of the newsgroups are very closely related to each other, while others are not. The data have 18,774 documents and 61,188 words in total.² Since stop words do not contain important significance to be used in search queries, we filter out a list of 119 English stop words³ from the dataset.

The experimental design is as follows. We vary weight parameters and for each weight parameter setting run the rand algorithm to discover clusters. We then use the same number of returned clusters and apply the k-means algorithm. As the newsgroup dataset is quite large, so we use a fast implementation of k-means algorithm [Sculley 2010] that is based on mini-batch optimization. To compare the clustering results of k-means and our methods, we consider two measures, cohesion and normalized mutual information (NMI). Cohesion measures how closely related objects in a cluster are and is defined as ${\sum }_{i=1}^k{\sum }_{x\in C_i}\mathit{\text{dist}}\left(x,c_i\right)$, where c_i denotes the centroid of cluster i, C_i denotes the points set of cluster i, k denotes the number of clusters, and dist (x,c_i) denotes the distance between point x and the cluster centroid c_i. There are several ways of measuring distance of two binary vectors. We use L₁ distance, L₂ distance, and cosine distance. We report the results of $\frac{\mathit{\text{cohesion of BMF clustering result}}}{\mathit{\text{cohes ion of }}k-\mathit{\text{means clustering result}}}$ in Table 9. Normalized mutual information determines the similarity of two different clusterings of a dataset. There are various forms of normalized mutual information. We adopt the form used in Kvalseth [1987]. Let U = {U₁,U₂, … ,U_R} and V = {V₁,V₂, … ,V_C} be two clustering (partitional) results of a set of N data items. The NMI value of two clustering results is computed as $\frac{2I(\mathbf{U},\mathbf{V})}{H(\mathbf{U})+H( \mathbf{V})}$, where $I(\mathbf{U},\mathbf{V})=\sum _{i=1}^R\sum _{j=1}^C\frac{|U_i\cap U_j|}{N}log\frac{|U_i\cap U_j|N}{|U_j||V_j|}$, $H(\mathbf{U})=-\sum _{i=1}^R\frac{|U_i|}{N}log\frac{|U_i|}{N}$, and $H(\mathbf{V})=-\sum _{i=1}^C\frac{|V_i|}{N}log\frac{|V_i|}{N}$. The newsgroup dataset is extremely sparse. The overlapping terms of two documents are much less than their non-overlapping terms. We deliberately set weight parameters to be small values, so that the number of resulting clusters is reasonable. Note that there is no deterministic relation between the number of clusters and weight parameters. The selection of parameters is done by exploratory search. We report the results in Table 9 with weight parameters ranging from 0.01 to 0.07.

We observe that with respect to cohesion, weighted rank-one BMF is comparable to k-means in any distance measure, as all values in columns 4–6 of Table 9 are around 1. With respect to NMI, values in column 7 suggest no significant similarity between the BMF and k-means clustering results. We suspect the randomness nature of both algorithms causes the fact. So we run the k-means algorithm one more time and study the similarity between the two k-means clustering results. The results are reported in the last column of Table 9. We observe that the values in the last column are close to their counterparts in the last second column. It confirms our suspect on the limitation of BMF and k-means that they both have randomness in nature and their results are susceptible to random choices made in the course of their heuristics.

We also compare our methods with the two recent k-means advancements, i.e., online k-means [Liberty et al. 2016] and local k-means [Gupta et al. 2017]. Online k-means assumes that we receive one point at each time. In this case a point is either assigned to an existing cluster or a new cluster. Local k-means incorporates an additional “cluster” to hold outliers. We compare our method with online k-means and local k-means with respect to cohesion ratio. The results are reported in Table 10. It shows that BMF performs better than both online k-means and local k-means. The main reason is that these two advanced k-means algorithms are not suitable for this application setting. Online k-means performs well when data points arrive in an online fashion. Local k-means is good for the dataset with many extreme outliers. However, neither works well with binary data.

As numbers do not tell the whole story, so we further look into the patterns (centroids) returned from BMF and k-means to see if they do provide any insight on the dataset. We find that BMF does return semantically meaningful patterns that were not discovered by k-means. For instance, BMF discovers a pattern of terms “canada,” “russia,” “austria,” “sweden,” “switzerland,” and “italy,” which exists in a few documents in the rec.sport.hockey category. The pattern can be interpreted as “World Hockey.” Another returned pattern includes terms of “playoff,” “hawks,” “farenebt,” “farenell,” “binghamton,” “moncton,” “breton,” “candiens,” “adirondack,” and others. The pattern exists in many documents in the rec.sport.hockey category and can be interpreted as “North American Hockey.” BMF also discovers many interesting sub-topics in other categories. The examples illustrate that weighted rank-one BMF is able to classify documents and terms at an adjustable resolution that is much finer than the available categorical information. As an exploratory research task, clustering is to group similar objects in attempt to draw knowledge to make conceptual distinction or posit an exploratory relationship. In that sense, the weighted rank-one BMF provided with fast and effective heuristics serves as a useful cluster analysis tool.

7.6. Application to Pattern Discovery in Role-based Access Control

In this section, we illustrate the application of weighted rank-one BMF to pattern discovery in the context of role-based access control (RBAC) [Ferraiolo et al. 2001]. RBAC has been regarded as the de facto access control model and deployed by many organizations and companies. As opposed to the conventional permission-based access control model (i.e., an access control list specifies which users or system processes are granted access to objects, as well as what operations are allowed on given objects), RBAC groups a large number of permissions into a few roles and associates roles with users. This substantially simplifies the complexity of users and permissions management. Another advantage of RBAC is that the use of roles as authorization subjects, instead of users, avoids having to revoke and re-grant authorizations whenever users change their positions and/or duties within the organization. Given these benefits, companies and organizations want to migrate from their old permission-based access control systems to the RBAC model. Roles can be defined by carefully analyzing and decomposing business processes into smaller units in a functionally independent manner. These functional units are then associated with permissions on information systems. But such a top-down approach is not viable when an organization involves more than thousands users and operations. In contrast, the bottom-up approach utilizes the existing permission assignments and discovers roles by extracting data patterns, which is commonly referred to as role mining [Kuhlmann et al. 2003]. To offer justifications for the identified roles and provide insightful analysis, the number of roles has been widely considered as a measure for the goodness/interestingness of a role set. Discovering a minimum role set is referred to as the basic RMP [Vaidya et al. 2007].

From the pattern discovery perspective, the basic RMP can be viewed as an overlapping clustering problem. This is because a user can be assigned to multiple roles, e.g., a person can have both STUDENT and TEACHING ASSISTANT roles. Weighted rank-one BMF can be applied to role mining due to its lineage to Boolean matrix decomposition, which has been used to model the basic RMP. Let binary matrices A, X, and Y denote user-permission assignments, user-role assignments, and role-permission assignments, respectively. For example, A(i, j) = 1 denotes that user i has permission j, while being 0 indicates not; X (i, j) = 1 means user i is assigned role j, while being 0 indicates not. The basic RMP can be formulated as a Boolean matrix decomposition problem, represented by A_m×n = X_m×k ⊗ Y_k×n, where the ⊗ sign denotes the Boolean matrix product, where 1 plus 1 is 1, as opposed to being 2 in standard matrix product. The use of Boolean matrix product is to reflect the fact that when a user is assigned to multiple roles containing a common permission, he/she gets that permission assignment only once. Given an input binary matrix A, there exist many possible decompositions.

The basic idea of applying weighted rank-one BMF to role mining is to iteratively perform weighted rank-one BMF to discover dominant patterns (roles) and assign them to records (users). Different from the previous studied applications, where each record belongs to one pattern, the basic RMP allows a user to be assigned to multiple roles (patterns). To reflect that, after each iteration, we keep those records that have been assigned to the previously identified roles. We mark the permission assignments (1’s elements) that have been covered by the identified roles. For the subsequent weighted rank-one BMF operations, we do not count those marked assignments. 0-becoming-1 (over-assignment) errors are strictly forbidden in role mining, as an over-assignment error may lead to abuse of rights and cause severe security problems. Therefore, for weighted rank-one BMF, we set w₁₀, the penalty weight on 0-becoming-1 errors, to be a sufficiently large number. Since we do not expect one role to fulfil the complete assignments of a user, therefore we set w₀₁ as 0 and w₁₁ as 1. In other words, at each iteration, an identified role is assigned to a user if the assignment does not introduce 0-becoming-1 errors. However, when we count the influence of a pattern, we count the number of newly covered assignments. In other words, at each iteration we try to discover the pattern, which covers the most remaining 1’s elements. To implement this, we use the local search algorithm as described in Algorithm 2. We generate a random pattern (role) and then search its neighborhood by iteratively flipping its component till no better solution can be found. Note that in this scenario (i) the influence of a pattern is computed as the number of newly covered assignments and (ii) a pattern is considered to exist in a record if the record contains the pattern.

We test the algorithm on the real access control datasets collected by Ene et al. [2008]. The description of datasets are provided in Table 11. We compare our algorithm with the benchmark algorithm DBP, which discovers discrete basis (patterns) from a binary matrix. The DBP implementation allows to set preferences for 1-becoming-0 and 0-becoming-1 errors. As 0-becoming-1 errors are prohibited in the role mining scenario, so for the DBP implementation, we set the weight of 1-becoming-0 errors as 1 and the weight of 0-becoming-1 errors as a sufficiently large number. The DBP implementation requires the number of roles (basis) as an input. We set the number of roles as the number of records, as the minimum number of roles cannot be greater than the number of records. The DBP implementation returns a basis matrix (Y) and a usage matrix (X) that describes how to reconstruct records (A) with the discovered basis. If a column of the usage matrix (X) consists of all zeros, then the corresponding basis vector (role) is not used. Then we delete that basis vector and its usage from X and Y accordingly and only keep those basis vectors (roles) that are in use. Results are reported in Table 11. Note that the reported errors are essentially 1-becoming-0 errors, because we configured both DBP and BMF to avoid 0-becoming-1 errors, because the 0-becoming-1 error type is strictly prohibited in the role mining stetting. The result shows that DBP is not able to fully cover permission assignments for all datasets. For instance, for the emea dataset, 32 roles are discovered and after assigning these roles to users, there are still 52 original user-permission assignments left uncovered. The issue cannot be solved by increasing the number of roles, as the additional roles returned by the DBP cannot cover those user-permission assignments either. The underlying reason is that the DBP algorithm has a fixed pool of candidate roles, which are derived from the association of attributes (permissions) of records. In contrast, BMF is able to successfully cover all user-permission assignments, as a role is discovered iteratively and there is no fixed candidate role pool. We can also observe that in addition to full coverage, the number of roles returned by BMF is very close to that of DBP. For the customer dataset, BMF finds a set of 276 roles with full coverage, while DBP returns 275 roles and still have 241 permission assignments uncovered. In some cases, the result is even better. For example, BMF finds a set of 456 roles to fully cover the apj dataset, while the DBP returns a set of 465 roles, which still leave 40 permission assignments uncovered. Thus, the BMF approach outperforms the DBP approach in the role mining setting, too.

As mentioned, one main advantage of weighted BMF over traditional matrix decomposition methods is that some real applications require the decomposed values to be binary. If we forcefully apply traditional methods, then it would return odd results. For instance, we could apply nonnegative matrix factorization to the binary setting by forcing the decomposed values to be within [0, 1] and then rounding values to 0/1. To compare with weighted BMF, we decompose a user-permission binary matrix into the two small non-negative matrices with the same rank as weighted BMF. We round decomposed values to 0/1 and compute their Boolean product to reconstruct the original binary matrix and count errors. The results are reported at the last two columns of Table 11. It shows that the reconstructed data are far different from the original one, which suggests the necessity of binary-valued matrix methods.