Recent Developments in Boolean Matrix Factorization

1 Introduction

Boolean matrix factorization (BMF) is a variant of the standard matrix factorization problem in the Boolean semiring: given a binary matrix, the task is to find two smaller binary matrices so that their product, taken over the Boolean semiring, is as close to the original matrix as possible. Because the matrix product is not done over a field but over a semiring, many standard matrix factorization techniques fail to work. Indeed, finding the best Boolean factorization is computationally hard.

The computational hardness of the problem has not prevented people from studying it. In psychometrics, some of the first algorithms appeared in the 1980's (see Bˇelohl´avek and Trnecka [2018]). Even before that, mathematicians studying combinatorics had studied the "Boolean linear algebra" [Kim, 1982; Monson et al., 1995]. More recently, Miettinen et al. [2006] introduced the problem to the data mining community, Vaidya et al. [2007] to the access control community, Bˇelohl´avek and Vychodil [2010] to formal concept analysts, and van den Broeck and Darwiche [2013] to lifted inference community, to name but a few examples. Even more recently, Ravanbakhsh et al. [2016] studied the problem in the framework of machine learning, increasing the interest of that community, and Chandran et al. [2016]; Ban et al. [2019]; Fomin et al. [2019] (re-)launched the interest in the problem in the theory community.

∗Contact Author

†The research leading to these results has received funding from the European Research Council under the European Community's Seventh Framework Programme (FP7/2007-2013) / ERC grant agreement No. 340506, and the Doctoral Programme "Vienna Graduate School on Computational Optimization" which is funded by the Austrian Science Fund (FWF, project no. W1260-N35).

The various fields in which BMF is used are connected by the desire to "keep the data binary". This can be because the factorization is used as a pre-processing step, and the subsequent methods require binary input, or because binary matrices are more interpretable in the application domain. In the latter case, Boolean algebra is often preferred over the field GF[2], as the XOR-operation can be harder to interpret.

The large number of areas where BMF or related problems are studied has given raise to many interesting approaches for solving this computationally hard problem. Unfortunately, it has also led to numerous re-inventions and repeated effort. Given the recent interest towards BMF in machine learning and AI communities, we believe that a timely survey will help to guide the research to novel directions and avoid reinventions of older ideas.

To that end, we present a concise survey of the recent theoretical results on the hardness of the problem in Section 4 and a survey on different types of algorithms proposed for the problem in Section 5. We will provide some lessons we have learned when developing and using some of these algorithms, as well as the main research directions and open problems we find most interesting. Before going further, however, we will provide the formal definitions of BMF in Section 2 and some applications thereof in Section 3.

2 Definitions and Related Formulations

In this section we will explain the basic definitions for Boolean matrix factorization (BMF) and related concepts. We will also cover the related formulations of the problem using bipartite graphs and set systems. These formulations are all equivalent, but they make connections to different problems more clear.

2.1 BooleanMatrix Factorization

The goal of matrix decompositions is to decompose (or factorize) a given input matrix into two (or more) smaller factor matrices whose product is close to the original matrix. Thus, to define BMF, we first define the corresponding product.

Definition 1. Given Boolean matrices B ∈ {0, 1}m×k and C ∈ {0, 1}k×n , the Boolean product of B and C is denoted B ◦ C ∈ {0, 1}m×n and it is defined element-wise as

The Boolean matrix product is like the standard matrix product of two binary-valued matrices, but the algebra is the Boolean semi-ring (the summation is defined as 1 + 1 = 1).

The Boolean rank of matrix A ∈ {0, 1}m×n is the least integer k such that there exists matrices B ∈ {0, 1}m×k and C ∈ {0, 1}k×n with A = B ◦ C. In practice, the Boolean rank is rarely that interesting; like normal matrix rank, it is not robust to noise, and many real-world matrices have almost-full rank (though exceptions exist, see Ene et al. [2008]).

Problem 1. Given a matrix A ∈ {0, 1}m×n and integer k, the goal of Boolean matrix factorization (BMF) is to find matrices B ∈ {0, 1}m×k and C ∈ {0, 1}k×n such that they minimize

In (2), subtraction, Frobenius norm ‖·‖F , and the summation are taken over the standard algebra while ⊕ stands for the element-wise XOR-operation. It might seem curious that the error is defined over different algebra than the matrix product - after all, the whole point of BMF is that it uses the Boolean matrix product - but it should be noted that the Boolean algebra does not allow for subtraction. On the other hand, as all elements of A and B ◦ C are binary, the (squared) Frobenius norm simplifies into the Hamming distance between A and B ◦ C, i.e. it is simply the number of entries with Aij ̸= (B ◦ C)ij .

2.2 Covering by Bicliques

Any rectangular binary matrix A ∈ {0, 1}m×n can be viewed as the bi-adjacency matrix of a bipartite graph G = (U ∪V, E). The rows of A correspond to vertices in U, columns of A correspond to vertices in V , and Aij = 1 if {i, j} ∈ E. Vice versa, any (unweighted, undirected) bipartite graph G can be identified with its bi-adjacency matrix A(G).

For bipartite graphs, BMF is equivalent to biclique cover.

Problem 2. Given a bipartite graph G = (U ∪ V,E) and an integer k, the goal of biclique cover is to find k bicliques C1, . . . ,Ck, such that they minimize

A biclique is a complete bipartite subgraph Cℓ = (Uℓ ∪ Vℓ, E(Cℓ)), where Uℓ ⊆ U, Vℓ ⊆ V and E(Cℓ) = Uℓ × Vℓ, i.e. Cℓ contains edges between all pairs of vertices in Uℓ and Vℓ. Note that Problem 2 does not require that the edges of Cℓ are a subset of edges of G. Abusing the notation, in (3) the symbol ⊕ is denotes the symmetric difference between sets.

The connection between Problems 1 and 2 can observed as follows: Let bℓ be the ℓth column of B and cℓ be the ℓth row of C. Then we can write B ◦ C as the OR of outer products of bℓ and cℓ , i.e. B ◦C = ∨kℓ=1 bℓcℓ . Each of these rank-one matrices (called components) bℓcℓ defines a rectangle of 1s, i.e. (bℓcℓ )ij = 1 iff (bℓ)i = (cℓ )j = 1. Hence, the bi-adjacency matrix corresponding to bℓcℓ defines a biclique Cℓ = (Uℓ, Vℓ) with Uℓ = {i : (bℓ)i = 1} and Vℓ = {j : (cℓ )j = 1}. Then the bi-adjacency matrix of the bipartite graph with edges ⋃kℓ=1 E(Cℓ) is ∨kℓ=1 BiℓCℓj = B ◦ C.

2.3 Set Systems

A binary matrix A ∈ {0, 1}m×n can also be identified as a set system: the columns correspond to n items in universe U and the rows are incidence vectors of m sets in a collection of sets S = {Si ⊆ U : i = 1, . . . ,m}. Then the following problem is equivalent to BMF and will be called set basis.

Problem 3. Given a set system (U, S) and an integer k, the goal of set basis is to find a collection of k sets C = {Ci ⊆ U : ℓ = 1, . . . , k} such that for every set Si ∈ S there exists a subcollection of Bi ⊆ C that jointly minimize

In this formulation, the rows of A correspond to the sets Si with Si = {j : Aij = 1}. The matrix C encodes the incidence vectors of the sets Cℓ ∈ C, i.e. Cℓi = 1 iff i ∈ Cℓ. The choice of which sets to include in which subcollection Bi is given by the rows of B, i.e. Biℓ = 1 iff Cℓ ∈ Bi. Thus, Problem 3 is equivalent to Problem 1.

3 Applications

Let us now discuss applications of BMF. Wicker et al. [2012] use BMF to pre-process multi-label classification data. They decompose the binary label matrix and use the columns of the left factor matrix B as sort-of "super-labels." For a recent review of this and other approaches for dimensionality reduction for multi-label classification, see Siblini et al. [2019].

Somewhat similarly, van den Broeck and Darwiche [2013] use BMF to speed up lifted inference in Boolean data. They show that the inference is efficient if the data has low Boolean rank, and to transform the data, they compute low-rank Boolean decomposition before doing the inference.

Vaidya et al. [2007] propose using BMF for role mining. In role mining, we are given a binary matrix that represents users and permissions, and the goal is to assign users to (possibly overlapping) roles and to give these roles necessary permissions. Hence, if A is the user-permission matrix, B assigns users to k different roles and C assigns the permissions to these roles. Ene et al. [2008] showed that many real-world matrices in this domain actually have low Boolean rank.

BMF has many more applications and due to lack of space, we can only name a few. In bioinformatics, BMF is used for studying transcriptomic data [Liang et al., 2020] and for identifying functional interactions in brain networks [Haddad et al., 2018]. BMF was also used for adding missing attribute information to products in the Amazon catalog [Rukat et al., 2017b]. Further applications include approximate logic synthesis [Hashemi et al., 2018] and discovery of patterns in network data [Kocayusufoglu et al., 2018].

4 Theoretical Results

In this section we summarize computational lower bounds and upper bounds for computing BMF.

4.1 Computational Complexity

Let us review the computational complexity of Problem 1.

Orlin [1977] showed that, given a bipartite graph G and a parameter k, it is NP-complete to decide whether the edges of G can be covered with k bicliques. Via the discussion from Section 2.2 this implies that it is NP-complete to decide whether there exist Boolean matrices B and C of rank k such

that ‖A-B ◦ C‖2F = 0. Hence, Problem 1 is NP-complete. Chandran et al. [2016] showed that under a standard as-

sumption in computational complexity theory, one cannot distinguish between the case when the objective function in Problem 1 is 0 or at least 1 in time 22o(k) (mn)O(1) . This implies that either the running time of an algorithmmust be 22Ω(k) or it must be (mn)ω(1) , i.e. the algorithm's running time is either doubly exponential in k or superpolynomial in the size of the input matrix.

In fact, the above results hold for any α-approximate solution. Here, we say that matrices B′ ∈ {0, 1}m×k and C′ ∈ {0, 1}k×n provide an α-approximation of the optimal objective function value if

Observe that any α-approximation must be able to distinguish between whether minB,C ‖A-B ◦C‖2F = 0 or minB,C ‖A-B ◦C‖2F > 0. Hence, the above results imply that for any α ≥ 1, computing an α-approximation for Problem 1 is NP-complete and cannot be solved in time 22o(k) (mn)O(1) under a standard complexity assumption.

4.2 Theoretical Algorithms

Let us now look at algorithms with theoretical guarantees.

Gramm et al. [2008] and Fomin et al. [2020] showed that an exact solution for BMF can be computed in time 22O(k) (mn)O(1) , thus matching the above lower bound.

Furthermore, Fomin et al. [2019] provided an algorithm which computes a (1 + ε)-approximate BMF in time 22O(k) /ε2 ·lg2 (1/ε) · mn. Note that when k and ε are constants, then the running time of this algorithm becomes O(mn), i.e. the algorithm runs in time linear in the size of the m×n input matrix A. This is unlike the exact algorithms, where the exponent of the mn-term is strictly larger than 1, i.e. they require superlinear running time even when k is a fixed constant. Ban et al. [2019] obtained a similar algorithm with slightly worse running time, but their algorithm extends to any finite field.

Bhattacharya et al. [2019] extended ideas of Fomin et al. [2019] and Ban et al. [2019] to obtain a 4-pass streaming algorithm which computes a (1 + ε)-approximate BMF. Their algorithm never stores more than 2 ˜O(2k /ε2 ) · (lg n)2k rows of the matrix and it has running time 2 ˜O(2k /ε2 ) · (lgn)2k · mn. Note that when k is a fixed constant, then the algorithm only stores a polylogarithmic number of rows.

4.3 Discussion

From a theoretical perspective, BMF is relatively wellunderstood. Indeed, note that the lower bound results by Chandran et al. [2016] and the running times of the algorithms by Fomin et al. [2019] and Ban et al. [2019] match up to lower order terms. Hence, there is no hope in removing the 22O(k) -terms in the running times of the algorithms.

Unfortunately, all of the above algorithms are rather impractical. First of all, the 22O(k) -term in the running time grows extremely fast (even for k = 7, 227 > 3 · 1038 ). Secondly, the algorithms have in common that they use several exhaustive enumeration procedures which are extremely slow in practice.

There are two main open questions in this area. First, can one improve upon the results by Bhattacharya et al. [2019] by providing a streaming algorithm which performs less than 4 passes and returns a BMF. Second, it would be interesting whether there are practical algorithms for which one can obtain similar running time and approximation guarantees to the ones mentioned above.

5 Algorithms

In this section, we discuss several practical algorithms for computing BMF. From a high-level point of view, the algorithms can be placed into two categories: Those based on combinatorial optimization and those based on continuous optimization. We will now discuss each type of these algorithms and conclude the section with some open questions.

5.1 Combinatorial Optimization

Some of the earliest algorithms for BMF come from psychometrics community. One example is the 8M algorithm in the BMDP statistical software (see Bˇelohl´avek and Trnecka [2018] for a detailed description). The 8M algorithm follows a simple local update process: it starts by building initial matrices B and C, and then iteratively replaces the initial columns of B and rows of C with (locally) better ones.

The iterative update approach is also used by Bˇelohl´avek and Trnecka [2018] in their GRECOND+ algorithm. GRECOND+ takes carefully into account the anti-monotonicity of Boolean algebra: if one component bℓcℓ introduces a 1 in location (i, j) where Aij = 0, this error is irrecoverable by the other components. Hence, GRECOND+ starts by generating a decomposition that does not over-cover (i.e. they generate a tiling) and then the components are iteratively updated to commit over-covering when that decreases the error. After one component is updated, other components are examined if they can be edited - or removed entirely - to improve the error.

The GRECOND+ algorithm belongs to a collection of BMF algorithms based on formal concept analysis [Ganter and Wille, 1999] (see also Bˇelohl´avek and Vychodil [2010]; Bˇelohl´avek and Trnecka [2015]). A formal concept, in the BMF framework, is a pair of row and column index sets (I, J) such that Aij = 1 for all i ∈ I and all j ∈ J and such that no new index can be added to neither I or J. Such a pair of indices corresponds to components bℓcℓ that do not over-cover (they also correspond to closed itemsets in frequent itemset mining). Most formal-concept based algorithms find non-overcovering factorizations as a consequence.

The tiling approach also motivated the PANDA+ algorithm [Lucchese et al., 2013]. However, instead of finding strictly non-overcovering initial components, as GRECOND+ does, PANDA+ finds dense initial components, that is, initial components that overcover only strictly limited amount of elements. In the second phase of the algorithm, these dense cores are extended, again provided that the extension does not increase the error. An important feature of PANDA+ is that the number of components, k, is not predefined. Rather, PANDA+ selects the number of components automatically, based on user-selected criteria for balancing the decrease in error and increase in the number of components.

Yet another approach is presented by the ASSO algorithm [Miettinen et al., 2008]. Unlike the previously-described algorithms, ASSO starts by generating a large collection of candidate columns for B. There will be a candidate column for every row ai of A. A candidate column is generated by comparing the association confidence from ai to aj , for any j = 1, . . . ,m, to a user-supplied threshold τ. The candidate column will have a 1 in the jth entry if the confidence, computed as ai (aj )T /(ai (ai )T ), is above the threshold τ .

ASSO then greedily selects the candidates one-by-one to be added to B by computing how much error each candidate would reduce. In order to compute the error, ASSO computes for each column a of A and for each candidate column d, how much the reconstruction error of a would change if d would be used to cover the column, that is, it computes ‖a -B ◦ c‖2F -

‖a - [B d] ◦ [ c1 ]‖2F , where B is the current left factor matrix and c is the column of the right factor matrix corresponding to a. The candidate is used in every column of A where it reduces the error, and the candidate that reduces the error the most over all columns of A is selected. The corresponding row of C is built based on the greedy computation: if the selected d reduces the error on column j, the corresponding column in C will be set to 1.

As mentioned, PANDA+ uses some measure to balance the reduction in error and the size of the factorization. There are different ways to do this, but one of the popular and principled ways for solving the model order selection problem is to use the minimum description length (MDL) principle [Rissanen, 1978]. Miettinen and Vreeken [2014] proposed to use two-part MDL for selecting the number of components in BMF. Twopart MDL computes the encoding length of the data L(A) as the length of the model (or hypothesis) L(M) plus the length of the data given the model L(A | M). In case of BMF, the model consists of the factor matrices, while the data given the model is the error the factorization causes. That is,

The length of encoding a binary matrix can be computed using standard methods for binary string encoding. Miettinen and Vreeken [2014] propose the use of so-called data-to-model encodings and differentiating over- and under-covering noise.

Miettinen and Vreeken [2014] use a variation of ASSO to find the decomposition that minimizes (5). Also PANDA+ can be used to find such decomposition. Karaev et al. [2015] also try to minimize (5) and their algorithm, called NASSAU, shares similarities with many of the above algorithms: it starts by finding a set of seed columns (akin to candidates in ASSO).

Unlike ASSO, though, NASSAU uses random walks in the bipartite graph corresponding to A to find the seeds. In short, they consider the probability of a random walk starting from a vertex u ∈ U to be in another vertex u′ ∈ U. The seed column corresponding to vertex u ∈ U has 1-entries in the rows corresponding to the vertices u′ where the random walk would be with a sufficiently high probability. The seed columns are then turned into rank-1 components using a similar greedy strategy as used by ASSO. After a new rank-1 component is added, the existing ones are iteratively updated, similarly to 8M, potentially dropping some if they become redundant.

Many of the above algorithms first build some initial solution or candidate set, and then refine it. The BMAD framework [Tyukin et al., 2014] formalizes this process and allows the user to define the approach for candidate generation, the approach for selecting the candidates, and the approach for building the factorization from the candidates.

5.2 Continuous Optimization

Next, let us look at continuous optimization algorithms.

Hess et al. [2017] propose to use a proximal alternating linearized minimization technique for finding the MDLoptimizing BMF. The idea is to replace BMF with a surrogate continuous optimization problem with a two-part goal: minimizing the error and making the factor matrices binary. This approach allows alternating between gradient descent optimization of the factor matrices and computing the proximal function. Hess et al. [2018] use a similar technique for finding Boolean decompositions that limit the false discovery rate, i.e. the probability that the found patterns arise from noise.

Next, we discuss a sequence of works which try to find factor matrices B and C with maximum likelihood under the given input matrix A. This includes works of Ravanbakhsh et al. [2016]; Rukat et al. [2017a]; Rukat et al. [2018]; Liang and Lu [2019]; Frolov et al. [2014]. In order to compute the likelihoods of the factor matrices, they make some prior assumptions on the data-generating process. We will only discuss the model of Ravanbakhsh et al. [2016] in detail.

From a high-level point of view, Ravanbakhsh et al. [2016] assume that the input matrix A for Problem 1 is generated as follows. First, two ground-truth factor matrices B ∈ {0, 1}m×k and C ∈ {0, 1}k×n are generated randomly according to a certain distribution. This gives a matrix Z = B◦C of Boolean rank k. Now the matrix A is generated from Z by setting A = Z and then flipping the bit in each entry Aij according to some distribution.

Let us now look consider this model in more detail. For all entries of the factor matrices B ∈ {0, 1}m×k and C ∈ {0, 1}k×n there exist priors pBir (·) and pCrj (·) for all i = 1, . . . ,m, r = 1, . . . , k and j = 1, . . . , n. In particular, it is assumed that Bir = 1 (Crj = 1) with probability pBir (Bir) (pCrj (Crj)) and that for both factor matrices the prior takes the following separable product form:

Next, given B and C as above, consider the product matrix Z = B ◦ C ∈ {0, 1}m×n and observe that Z only depends on the randomness of the entries in B and C. It is assumed that for all i = 1, . . . , m and j = 1, . . . , n there exists another prior pAij (Aij | Zij), i.e. the probability of flipping the bit in entry Aij can depend on i, j and Zij. It is further assumed that for the matrix A it holds that

Recall that in Problem 1 the input is the matrix A and we are actually interested in computing the (unknown) factor matrices B and C. The approach of Ravanbakhsh et al. [2016] is to solve the maximum a posteriori inference problem

To solve the above inference problem, they introduce a graphical model for the posterior and use belief propagation to find the matrices B and C which maximize the posterior.

Frolov et al. [2014] and Liang and Lu [2019] used a similar generative model and expectation maximization instead of belief propagation. Rukat et al. [2017a] used a Metropolised Gibbs sampler for the posterior inference. Later, Rukat et al. [2018] extended their algorithm for Boolean tensors.

5.3 Open Problems and Lessons Learned

The algorithms based on combinatorial optimization use various heuristics and iterative update approaches to avoid the problem that their optimization landscape is not at all smooth. Indeed, most iterative update algorithms will converge in just few rounds (see, e.g. the experiments by Miettinen [2009]). No currently-published algorithm combines all of the approaches (many candidates, preference for under-covering, iterative local updates, etc.). Whether a combination of these approaches would improve the results will remain as a topic of further research, although our experience indicates that combining multiple approaches has diminishing returns in most cases.

At the same time it should also be noticed that while it is easy to build examples where greedy heuristics will fail spectacularly, they often tend to work very well in real-world applications. In case of BMF, the experiments by Miettinen [2009] seem to confirm this at least when it comes to the problem of selecting the candidates to build the factorization.

When the data is very sparse, some algorithms based on continuous optimization outperform all combinatorial algorithms we are aware of. Specifically, the experiments by Ravanbakhsh et al. [2016] show that on synthetic data, their algorithm recovers planted structures close to what is possible (they show that the reconstruction error of their algorithm is close to an information theoretic lower bound by Davenport et al. [2014]). This is also supported by the experiments by Neumann [2018].

A major open question is to make existing algorithms more scalable. Indeed, existing algorithms do not easily scale to matrices with tens or hundreds of thousands rows and columns. For the combinatorial optimization algorithms, a major bottleneck is that they need to generate a lot of candidate factor matrices and, hence, need to cover the input matrix many times. For the continuous optimization algorithms, a major bottleneck is that usually they maintain dense real-valued versions of the factor matrices; this makes them slow and memory-intensive.

However, these methods may benefit from the recent GPU development and methods such as stochastic gradient descent; indeed, such techniques are used by Hess et al. [2017]. In practice, the scalability problem is often bypassed by pruning all rows and columns from the input matrix which have less than a user-specified threshold of non-zero entries.

None of the above algorithms comes with provable guarantees to return a solution within a certain approximation ratio. While some of themwill converge to an optimal solution, it is not clear how many iterations will be necessary.

6 RelatedWork

BMF has been extended (and constrained) in various forms to work better in different applications.

Tatti and Miettinen [2019] aim at finding BMF where the rows and columns of the data and factor matrices can be permuted such that the 1s in the columns of B and in the rows of C all appear consecutively. This is used for sorting vertices when graphs are drawn as edge bundles.

Many datasets are multi-view: they present different data about the same entities. Sometimes such data can be expressed as two (or more) matrices, and they can be jointly decomposed. Miettinen [2012] solves joint BMF: given A and B whose rows represent the same entities, find W, U, V , H, and L such that A ≈ [W U] ◦ H and B ≈ [W V ] ◦ L. Similarly, Hess and Morik [2017] find decompositions of type A ≈ (W∨U)◦H and B ≈ (W∨V )◦L (note that the dimensions of the matrices are not the same in these two definitions).

One of the most actively studied extension of BMF is Boolean tensor decompositions (see, e.g. Miettinen [2011]; Erd˝os and Miettinen [2013]; Rukat et al. [2018]). Briefly, tensors extend the matrix (2-way) data to multi-way data. For example, a Boolean rank-r CP decomposition of a three-way tensor A would decompose it into three factor matrices, B, C, and D, minimizing

While BMF can be seen as covering a bipartite graph with (potentially overlapping) bicliques, many algorithms have been proposed for the problem of covering a bipartite graph with disjoint bicliques. These problems are often studied under names such as biclustering [Lim et al., 2015], coclustering [Dhillon, 2001] or bipartite graph partitioning [Zha et al., 2001]. Kumar et al. [2019] provide a constant factor approximation algorithm for covering a bipartite graph with disjoint bicliques with running time 2O(k2 ) (mn)O(1) .

Covering random bipartite graphs with disjoint bicliques has recently received some attention [Lim et al., 2015; Xu et al., 2014; Zhou and Amini, 2018; Zhou and Amini, 2019; Ndaoud et al., 2019]. In these bipartite stochastic block models one assumes that the input bipartite graph was randomly generated based on some hidden bicliques. The goal is to state conditions under which the hidden bicliques can be recovered when the input only consists of the random bipartite graph. These random graph models are similar to the datagenerating models mentioned in Section 5.2 (assumptions on the hidden bicliques correspond to priors on factor matrices). Neumann [2018] showed that in such random bipartite graphs a simple two-step algorithm can recover even very tiny clusters, which is not possible in general graphs.

7 Future Directions and Conclusion

In this survey, we have summarized the current state of the art for BMF algorithms including several recent developments based on results from the machine learning and theoretical computer science community.

From a practical point of view, perhaps the most important open problem is to make the existing algorithms more scalable.

From a theoretical point of view, the lower bounds from Section 4.1 essentially rule out any practical algorithms which provably approximate BMF. Nonetheless, we have discussed in Section 5.3 that there are algorithms which appear to recover very high quality results in practice. Hence, we believe that one should try to analyze BMF algorithms not using the current worst-case analysis approach, but using other complexity measures. In our opinion, one promising approach would be to generalize the results for bipartite stochastic block models we mentioned in Section 6 from disjoint biclique covers to overlapping biclique covers. This would imply practical algorithms for BMF with provable quality guarantees. Since some of the bipartite stochastic block model algorithms work well in practice and their running times do not carry large overheads in some parameters as the ones from Section 4.2, we believe that there is hope to obtain such results.