Boolean Matrix Factorization (BMF)

• Boolean Matrix Factorization is the process of decomposing a binary matrix into the Boolean product of lower-rank binary matrices, capturing latent structures with high interpretability.
• It employs a range of methods—from formal concept analysis and combinatorial heuristics to SAT and probabilistic models—to minimize Hamming error under NP-hard constraints.
• BMF finds applications in bioinformatics, recommendation systems, and logic synthesis, effectively balancing precise reconstruction with scalable algorithmic solutions.

Boolean matrix factorization (BMF) is the problem of approximating a binary matrix $X \in \{0,1\}^{m \times n}$ by the Boolean product of two lower-rank binary matrices, typically $W \in \{0,1\}^{m \times r}$ and $H \in \{0,1\}^{r \times n}$, with product defined entrywise by $$(W \odot H)_{ij} = \bigvee_{k=1}^r (W_{ik} \wedge H_{kj})$$. In the standard formulation, one fixes a target rank $r$ and minimizes Hamming reconstruction error, or equivalently the squared Frobenius norm for binary data. The field is motivated by the high interpretability of binary factors and constrained by substantial computational hardness: exact factorization is NP-complete, the fixed-rank optimization problem is NP-hard, and practical work therefore spans combinatorial heuristics, formal-concept methods, SAT and integer-programming formulations, probabilistic models, and continuous relaxations (Miettinen et al., 2020).

1. Algebraic formulation and representational semantics

The Boolean product differs fundamentally from ordinary arithmetic multiplication. For binary factors $B \in \{0,1\}^{m \times r}$ and $C \in \{0,1\}^{r \times n}$, the product $B \odot C$ or $B \circ C$ is defined by a disjunction of conjunctions, so each reconstructed 1-entry indicates that at least one latent factor simultaneously activates the corresponding row and column coordinates. In this setting, the standard optimization problem is

$$\min_{W \in \{0,1\}^{m \times r},\, H \in \{0,1\}^{r \times n}} \|X-(W \odot H)\|_F^2,$$

with exact factorization yielding the Boolean rank of $W \in \{0,1\}^{m \times r}$0 (Miettinen et al., 2020).

Several equivalent viewpoints coexist. A rank-1 Boolean factor corresponds to a biclique in the bipartite graph induced by the 1-entries of the matrix, so a rank-$W \in \{0,1\}^{m \times r}$1 factorization is a superposition of $W \in \{0,1\}^{m \times r}$2 bicliques or Boolean rectangles (Geraci et al., 11 Sep 2025). In formal concept analysis (FCA), a formal concept $W \in \{0,1\}^{m \times r}$3 is a maximal all-1 rectangle, and a collection of such concepts induces factor matrices whose Boolean product recovers the covered entries exactly. Belohlávek–Vychodil universality and optimality results, as summarized in FCA-based work, state that every binary matrix admits a factorization by formal concepts, and that if a factorization of inner dimension $W \in \{0,1\}^{m \times r}$4 exists, then there is a subset of formal concepts of size at most $W \in \{0,1\}^{m \times r}$5 achieving the same exact reconstruction (Nenova et al., 2013).

A further distinction concerns approximation geometry. In “from-below” BMF, the approximation satisfies $W \in \{0,1\}^{m \times r}$6 elementwise, so every 1 in the reconstruction covers only 1s of the data; this eliminates false positives and turns the residual into uncovered 1s only. FCA-based and MDL-guided methods frequently exploit this constraint because it sharpens interpretability and aligns factors with maximal all-1 rectangles (Makhalova et al., 2019).

Quality measures vary with the intended use. The survey literature records four recurring objectives: Hamming reconstruction error, precision/recall on the 1-entries, a two-part minimum-description-length criterion $W \in \{0,1\}^{m \times r}$7, and relative approximation guarantees $W \in \{0,1\}^{m \times r}$8 (Miettinen et al., 2020). This multiplicity is consequential: BMF is not defined by a single optimization criterion, even when the algebraic product is fixed.

2. Boolean rank, uniqueness, and model order

The Boolean rank $W \in \{0,1\}^{m \times r}$9 is the minimum $H \in \{0,1\}^{r \times n}$0 for which $H \in \{0,1\}^{r \times n}$1 or $H \in \{0,1\}^{r \times n}$2 for binary factors of inner dimension $H \in \{0,1\}^{r \times n}$3. When $H \in \{0,1\}^{r \times n}$4 is regarded as the biadjacency matrix of a bipartite graph $H \in \{0,1\}^{r \times n}$5, $H \in \{0,1\}^{r \times n}$6 equals the biclique-cover number $H \in \{0,1\}^{r \times n}$7; this identity is the standard bridge between algebraic factorization and graph covering (Geraci et al., 11 Sep 2025). The same equivalence underlies the survey statement that deciding whether zero reconstruction error is achievable for a given $H \in \{0,1\}^{r \times n}$8 is equivalent to covering a bipartite graph with $H \in \{0,1\}^{r \times n}$9 bicliques (Miettinen et al., 2020).

Model order is central because BMF simultaneously seeks latent patterns and an appropriate rank. Several rank-selection principles appear in the literature. Two-part MDL minimizes the encoding length of the factors and the residual, and is used both as an outer-loop criterion over candidate ranks and as an internal selection mechanism for greedy factor extraction (Miettinen et al., 2020). In MDLGreConD, this criterion is specialized to from-below decompositions, where $$(W \odot H)_{ij} = \bigvee_{k=1}^r (W_{ik} \wedge H_{kj})$$0 is partitioned into a covered part and an uncovered error matrix $$(W \odot H)_{ij} = \bigvee_{k=1}^r (W_{ik} \wedge H_{kj})$$1, and factor concepts are added greedily only when they reduce total description length (Makhalova et al., 2019). TrustPal replaces heuristic rank penalties with false-discovery-rate control: each tile is tested under a Bernoulli noise model by a density bound or a coherence bound, and rank expansion stops when added tiles do not survive the statistical test (Hess et al., 2019). The bfact package scans $$(W \odot H)_{ij} = \bigvee_{k=1}^r (W_{ik} \wedge H_{kj})$$2 and stops when error no longer improves for $$(W \odot H)_{ij} = \bigvee_{k=1}^r (W_{ik} \wedge H_{kj})$$3 consecutive steps, reporting reliable recovery of the true rank in simulated settings (Visscher et al., 7 Sep 2025).

Questions of uniqueness are more specialized but important. “Factorization of Binary Matrices: Rank Relations, Uniqueness and Model Selection of Boolean Decomposition” studies factorizations under standard arithmetic, nonnegative, Boolean, and $$(W \odot H)_{ij} = \bigvee_{k=1}^r (W_{ik} \wedge H_{kj})$$4 operations, examines relationships between the different ranks, and characterizes when a Boolean factorization $$(W \odot H)_{ij} = \bigvee_{k=1}^r (W_{ik} \wedge H_{kj})$$5 has a unique $$(W \odot H)_{ij} = \bigvee_{k=1}^r (W_{ik} \wedge H_{kj})$$6, a unique $$(W \odot H)_{ij} = \bigvee_{k=1}^r (W_{ik} \wedge H_{kj})$$7 for fixed $$(W \odot H)_{ij} = \bigvee_{k=1}^r (W_{ik} \wedge H_{kj})$$8, and when both $$(W \odot H)_{ij} = \bigvee_{k=1}^r (W_{ik} \wedge H_{kj})$$9 and $r$0 are unique, given a rank constraint. The same work introduces BMF$r$1, described as a robust Boolean model-selection method that, on numerical examples, accurately determines the correct number of Boolean latent features and reconstructs the pre-determined factors accurately (DeSantis et al., 2020).

Lower bounds on Boolean rank also matter because they delimit feasible model orders. Beyond the classical induced-matching lower bound, recent commutative-algebra work shows

$r$2

and constructs a second lower bound through the isolation number $r$3, again satisfying $r$4. These bounds can be tight for block-diagonal all-1 structures, but they can also be far from the true rank: for identity-complement matrices $r$5, $r$6 and $r$7, while $r$8 can grow unboundedly with $r$9 (Geraci et al., 11 Sep 2025). A common misconception is therefore that readily computed lower bounds are close surrogates for Boolean rank; the documented gap shows otherwise.

3. Algorithmic paradigms

The algorithmic literature is heterogeneous because no single optimization strategy dominates across ranks, densities, missing-data patterns, and matrix scales.

Paradigm	Representative papers	Core mechanism
FCA and from-below methods	(Nenova et al., 2013, Makhalova et al., 2019, Krajča et al., 14 Mar 2026)	formal concepts, greedy coverage, MDL-guided concept selection
SAT, MaxSAT, and IP exact methods	(Avellaneda et al., 2021, Kovacs et al., 2020, Kolomvakis et al., 2023)	exact encodings, column generation, AO with integer subproblems
Probabilistic and continuous relaxations	(Rukat et al., 2017, Truong et al., 2021, Dalleiger et al., 2023, Dalleiger et al., 2024)	MCMC, auxiliary nonnegative optimization, proximal gradient, federated proximal updates
Large-scale heuristics and hybrids	(Wan et al., 2019, Desouki et al., 2019, Fomin et al., 2022, Visscher et al., 7 Sep 2025)	geometric segmentation, top fibers, constrained clustering, restricted master problem

FCA-based methods exploit the fact that formal concepts are maximal all-1 rectangles. Greedy procedures repeatedly choose the concept maximizing newly covered 1s, remove those entries from the residual matrix, and continue until a target number of factors or a coverage fraction is reached. MDLGreConD modifies this pattern by searching on demand for concepts that reduce a description-length objective, while GreCon3 retains exact-greedy behavior but replaces GreCon2’s memory-intensive bookkeeping with a sparse row-jagged structure and on-the-fly incremental initialization [(Nenova et al., 2013); (Makhalova et al., 2019); (Krajča et al., 14 Mar 2026)].

Exact optimization methods target smaller or medium-scale instances. SAT and MaxSAT encodings introduce factor variables and per-cell error indicators, with the MaxSAT objective $B \in \{0,1\}^{m \times r}$0 equal to Hamming error by construction; for larger matrices, the same work switches to a maximal-biclique edge-cover heuristic (Avellaneda et al., 2021). Column-generation formulations instead treat rank-1 patterns as exponentially many variables and solve a restricted master problem plus a pricing problem over candidate bicliques, thereby providing dual bounds and, in many cases, high-accuracy factorizations without heuristic pattern mining (Kovacs et al., 2020). Alternating-optimization with integer programming decouples one factor update into independent BoolLS subproblems, each solved exactly by an IP; a second IP can then select the best subset of rank-1 factors drawn from multiple runs (Kolomvakis et al., 2023).

Probabilistic and continuous approaches replace the discrete search by relaxations or full generative models. The OrMachine introduces Bernoulli priors on factors, a noisy OR likelihood with global dispersion parameter $B \in \{0,1\}^{m \times r}$1, and Metropolised Gibbs sampling for posterior inference, explicitly supporting completion and uncertainty quantification (Rukat et al., 2017). BANMF transforms exact BMF into a nonnegative auxiliary optimization problem with a matrix $B \in \{0,1\}^{m \times r}$2 constrained to share the Boolean support pattern of $B \in \{0,1\}^{m \times r}$3, then thresholds the nonnegative solution to recover Boolean factors under exact-solution conditions (Truong et al., 2021). Elbmf relaxes BMF to $B \in \{0,1\}^{m \times r}$4-valued factors regularized by an elastic-binary penalty and solves the resulting problem by proximal alternating linearized minimization, gradually increasing the binarizing regularizer to obtain Boolean outputs without a separate threshold-search stage (Dalleiger et al., 2023). In federated settings, Felb keeps local relaxed factors on clients, shares only relaxed component matrices, and aggregates them centrally with a proximal operator while providing convergence and differential-privacy guarantees (Dalleiger et al., 2024).

Heuristic large-scale methods emphasize speed and candidate quality. MEBF permutes the matrix toward an Upper-Triangular-Like arrangement with Simultaneous Consecutive-ones Property heuristics, seeds a dense rectangle from a median row or column, expands it bidirectionally, and subtracts the corresponding rank-1 pattern from the residual (Wan et al., 2019). topFiberM greedily selects the row or column fiber with the largest number of yet-uncovered ones, extends it into a rectangle under a precision threshold, and may replace a weaker stored factor if a later candidate yields better gain (Desouki et al., 2019). PLRBMF translates BMF into constrained clustering inspired by the EPTAS literature and applies Lloyd-style refinement, with recursive handling of larger ranks through sub-rank blocks (Fomin et al., 2022). The bfact system adopts a hybrid restricted-master-problem design: clustering generates a modest candidate set, a warm-started restricted master problem chooses disjoint factors, and a second MIP or heuristic reassignment-and-pruning stage refines the solution (Visscher et al., 7 Sep 2025).

4. Complexity, guarantees, and mathematical structure

Hardness results are foundational. The survey literature states that deciding exact factorization for a given rank is NP-complete, and that, unless ETH-type assumptions fail, there is no algorithm running in time $B \in \{0,1\}^{m \times r}$5 even to distinguish exact from non-exact factorizations. At the same time, parameterized exact solvers with runtime $B \in \{0,1\}^{m \times r}$6, a $B \in \{0,1\}^{m \times r}$7-approximation in time $B \in \{0,1\}^{m \times r}$8, and a 4-pass streaming $B \in \{0,1\}^{m \times r}$9-approximation are all recorded in the theoretical literature (Miettinen et al., 2020).

The EPTAS results of Ban et al. and Fomin et al. established a major approximation-theoretic threshold: for fixed $C \in \{0,1\}^{r \times n}$0 and $C \in \{0,1\}^{r \times n}$1, one can compute a rank-$C \in \{0,1\}^{r \times n}$2 factorization whose Hamming cost is within a $C \in \{0,1\}^{r \times n}$3 factor of optimum. Yet the reported running time,

$C \in \{0,1\}^{r \times n}$4

has double-exponential dependence on $C \in \{0,1\}^{r \times n}$5 and $C \in \{0,1\}^{r \times n}$6, and is described as unrealistic for moderate ranks. This theoretical–practical gap directly motivated constrained-clustering heuristics that preserve some EPTAS ideas while abandoning exact parameter dependence (Fomin et al., 2022).

Optimality guarantees vary by method. When the MaxSAT formulation runs to completion, it yields an optimal solution or a certificate of optimality for the encoded instance, while its large-matrix heuristic carries an $C \in \{0,1\}^{r \times n}$7 competitive ratio in the worst case for biclique cover (Avellaneda et al., 2021). Column generation yields lower bounds from the master LP and valid early-stopping bounds through reduced-cost analysis, and can re-solve a smaller integer model over generated columns for an integer feasible solution with explicit gap information (Kovacs et al., 2020). Continuous methods usually trade such certificates for convergence guarantees: BANMF proves a nonincreasing auxiliary objective; Felb proves convergence of the federated proximal-gradient sequence and Boolean convergence when the regularization parameter grows unboundedly; Elbmf guarantees that every limit point is a critical point of the relaxed objective and uses the regularization schedule to force a Boolean accumulation point (Truong et al., 2021, Dalleiger et al., 2024, Dalleiger et al., 2023).

A second misconception concerns the relation between Boolean and nonnegative factorization. BANMF explicitly frames a nonnegative auxiliary problem whose solution space coincides with exact BMF solutions under the existence of an exact decomposition, but the same work also studies the gap between nonnegative rank and Boolean rank empirically. This should not be conflated with an unconditional equivalence of the two models (Truong et al., 2021).

5. Empirical performance, scalability, and rank recovery

Empirical comparisons are typically organized around reconstruction error, factor sparsity or interpretability, rank recovery, and runtime. The hybrid bfact system reports simulated matrices up to $C \in \{0,1\}^{r \times n}$8, ranks $C \in \{0,1\}^{r \times n}$9, noise $B \odot C$0 up to $B \odot C$1, and densities up to $B \odot C$2. In those experiments, bfact-recon and bfact-MDL recover $B \odot C$3–$B \odot C$4 on signal $B \odot C$5 and estimate $B \odot C$6 with $B \odot C$7 error, outperforming PANDA+, PRIMP, and MDL4BMF; on Human Lung Cell Atlas single-cell RNA-sequencing benchmarks, bfact achieves $B \odot C$8 on $B \odot C$9–$B \circ C$0 using $B \circ C$1–$B \circ C$2 factors where alternatives use $B \circ C$3–$B \circ C$4 (Visscher et al., 7 Sep 2025).

Elbmf reports that on synthetic data it converges quickly, recovers the ground truth precisely, and estimates the simulated rank exactly, while on nine real-world datasets it yields the lowest relative Hamming loss on every dataset. The same study describes its runtime as faster than Primp, comparable to clamped NMF, and often significantly faster than Asso, Grecond, OrM, and Sofa on large sparse matrices (Dalleiger et al., 2023). The OrMachine emphasizes a different performance axis: posterior calibration. On matrix completion tasks with only $B \circ C$5–$B \circ C$6 observations, it recovers more than $B \circ C$7 of missing entries, and on large single-cell data with $B \circ C$8 measurements it runs on commodity hardware using a Metropolised Gibbs sampler with parallel updates (Rukat et al., 2017).

Discrete greedy methods remain competitive because of their speed. topFiberM attained the highest coverage in 16 out of 20 fixed-rank settings and was best in 17 out of 20 fixed-coverage settings on standard Boolean benchmarks; on RDF matrices from Semantic Web Dog Food, it was in average 128 times faster than Asso and 838 times faster than GreConD (Desouki et al., 2019). MEBF reports, on large crime and single-cell benchmarks, higher coverage, lower density factors, and runtimes of 2.9s versus 333s and 10.6s versus 992s against MP on the crime dataset, and 1.8s versus 137s and 5.9s versus 390s on the scRNA-seq dataset (Wan et al., 2019). GreCon3, while preserving the exact greedy factors of GreCon2, is reported to be up to 10× faster on large or dense tables and to succeed on some datasets where GreCon2 runs out of memory (Krajča et al., 14 Mar 2026).

For medium-scale exact or near-exact optimization, AO with integer subproblems and combination IPs delivers strong results. The 2023 integer-programming study reports that multistart AO matches or improves best known errors in 8 of 12 cases within 30 seconds and 10 of 12 cases within 5 minutes, while MS-Comb-AO improves further on several instances (Kolomvakis et al., 2023). Column generation reports strictly the lowest reconstruction error in 15 out of 24 tested instances under a 20-minute budget (Kovacs et al., 2020). SAT/MaxSAT methods report exact optimality on small instances and large-matrix heuristic factorization within 5–10% of the best heuristic in prior work while being 10×–100× faster (Avellaneda et al., 2021).

6. Applications and contemporary directions

Applications span recommendation, security, bioinformatics, language and graph data, hardware design, and privacy-preserving analytics. In collaborative filtering, FCA-based BMF on binarized MovieLens ratings achieves mean average error competitive with SVD, and “it is enough to have a binary-scaled rating data to obtain almost the same quality in terms of MAE by BMF than for the SVD-based algorithm in case of non-scaled data” (Nenova et al., 2013). The broader survey literature also lists multi-label classification, role mining in access control, lifted inference in graphical models, bioinformatics, catalog attribute imputation, and approximate logic synthesis as representative domains (Miettinen et al., 2020).

In computational biology and genomics, BMF is repeatedly used because binary factors map naturally to gene modules, cell states, and marker sets. The OrMachine analyzes single-cell gene expression from 1.3 million mouse brain cells across 11 thousand genes and reports hierarchical gene-set structure under full posterior uncertainty (Rukat et al., 2017). bfact uses Human Lung Cell Atlas single-cell RNA-sequencing data to show high signal recovery with much lower rank, arguing that lower $B \circ C$9 is desirable for biological interpretability via marker-gene sets (Visscher et al., 7 Sep 2025). Elbmf’s TCGA case study reports that one of the 33 recovered components selects 356 genes strongly associated with thyroid hormone generation and thyroid gland development, while a t-SNE embedding of the reconstruction separates 33 cancer types more clearly than the raw binarized input (Dalleiger et al., 2023).

In machine learning and data mining, BMF has been used for topic modeling and imaging. The 2025 integer-programming and heuristic study evaluates tdt2 topic-modeling subsets up to $$\min_{W \in \{0,1\}^{m \times r},\, H \in \{0,1\}^{r \times n}} \|X-(W \odot H)\|_F^2,$$0, reporting lowest errors for IP-AO and MS-Comb-AO and semantically coherent topics such as “Clinton–Lewinsky” and “Nagano Olympics”; on CBCL faces, AO-BMF extracts localized features, while greedy methods produce more global ones (Kolomvakis et al., 3 Dec 2025). Approximate logic synthesis is another distinctive direction: BLASYS factorizes truth tables of circuit subblocks and uses the factors to synthesize less complex subcircuits, reporting up to 63% in power savings while introducing an average relative error of 5% (Hashemi et al., 2018).

Current directions include privacy, distributed optimization, and software engineering for large instances. Felb formulates federated BMF through relaxed local factor updates and central proximal aggregation, proving convergence and differential-privacy guarantees while outperforming federation schemes built from state-of-the-art BMF methods on synthetic and real-world datasets (Dalleiger et al., 2024). The bfact package is a pip-installable Python package depending on Gurobi or CPLEX, scikit-learn, leidenalg, numpy, and scipy, and is described as scalable to matrices up to $$\min_{W \in \{0,1\}^{m \times r},\, H \in \{0,1\}^{r \times n}} \|X-(W \odot H)\|_F^2,$$1 at 10% density (Visscher et al., 7 Sep 2025). GreCon3 emphasizes memory-conscious data structures for FCA-based exact-greedy factorization, while C++ bitset and bitmatrix structures in the 2025 heuristic/IP work report 1 bit per entry storage and up to 5× speedup over Eigen/Armadillo for Boolean multiplication when $$\min_{W \in \{0,1\}^{m \times r},\, H \in \{0,1\}^{r \times n}} \|X-(W \odot H)\|_F^2,$$2 (Krajča et al., 14 Mar 2026, Kolomvakis et al., 3 Dec 2025).

Across these settings, the central tension remains stable: BMF offers binary latent structure that is often more directly interpretable than real-valued factorizations, but the choice among exact, probabilistic, continuous, and greedy methods is driven by rank regime, matrix scale, tolerance for approximation, and the need for guarantees on error, uncertainty, privacy, or model order.