Binary Spike-and-Slab Lasso Biclustering

• Binary Spike-and-Slab Lasso Biclustering (BiSSLB) is a Bayesian method that uses logistic matrix factorization with spike-and-slab priors to recover overlapping biclusters in binary data.
• It automatically infers the number of biclusters through an Indian Buffet Process prior, effectively addressing noise, scalability, and bias challenges in binary matrices.
• The method employs an efficient coordinate-ascent algorithm with proximal updates to ensure computational tractability even for large datasets with variable bicluster sizes.

Binary Spike-and-Slab Lasso Biclustering (BiSSLB) is a Bayesian method for biclustering binary datasets, introduced to recover coherent subsets of rows and columns in matrices such as single-nucleotide polymorphism (SNP) and protein-protein interaction (PPI) data (Fan et al., 19 Mar 2026). It is formulated for settings in which most existing biclustering methods for continuous data are not applicable and in which binary-data alternatives may struggle under noise, scale poorly, or bias results toward biclusters of particular size or character. BiSSLB combines logistic matrix factorization, spike-and-slab lasso priors on latent factors, and an Indian Buffet Process (IBP) prior for automatic determination of the number of biclusters. The resulting model is designed to recover overlapping biclusters of various sizes without prior knowledge of the noise level or bicluster characteristics, while remaining computationally tractable through a coordinate-ascent algorithm with proximal updates (Fan et al., 19 Mar 2026).

1. Problem formulation and conceptual scope

Biclustering seeks simultaneous structure in the row and column dimensions of a data matrix, rather than imposing a single global partition. In binary matrices, this means identifying local submatrices with coherent presence–absence patterns. BiSSLB is explicitly targeted at this setting and is presented as a binary-data biclustering method for problems in which standard continuous-data biclusturing frameworks are inappropriate (Fan et al., 19 Mar 2026).

The method is motivated by three recurrent difficulties in binary biclustering. First, recovery can degrade substantially when the data are noisy. Second, scalability can become limiting when the matrix is large. Third, some algorithms are biased toward biclusters of a particular size or characteristic. BiSSLB addresses these issues by combining a probabilistic binary-data likelihood with sparse latent structure and nonparametric control over model size. A central feature is that biclusters may overlap and may vary in size, rather than being constrained to a fixed partition (Fan et al., 19 Mar 2026).

A useful clarifying point is that BiSSLB is not a fixed-$K$ biclustering method. Although inference is carried out under a truncation level $K^\star$ for practical computation, the effective number of columns is inferred from the data through the IBP stick-breaking construction, and inactive columns are removed during optimization. This distinguishes the practical truncation used for computation from the final active bicluster count $\widehat K$ (Fan et al., 19 Mar 2026).

2. Probabilistic model

Let $Y=(y_{ij}) \in \{0,1\}^{I \times J}$ denote the observed binary matrix. BiSSLB assumes a logistic matrix-factorization model,

$$y_{ij} \mid A,B,\mu \sim \mathrm{Bernoulli}\!\left(\sigma\!\left(\mu_i + (AB^\top)_{ij}\right)\right),$$

with latent score matrix

$$M = \mu\,\mathbf 1^\top + A B^\top, \qquad A \in \mathbb R^{I \times K},\; B \in \mathbb R^{J \times K},\; \mu \in \mathbb R^I.$$

Here $\mu$ is a row-specific intercept, while the factorization $AB^\top$ captures bicluster structure (Fan et al., 19 Mar 2026).

Sparsity is induced through spike-and-slab lasso (SSL) priors on the entries of $A$ and $B$. Using

$K^\star$0

the factor-specific priors are

$K^\star$1

$K^\star$2

with $K^\star$3 and $K^\star$4. The larger “spike” penalties encourage zeros, while the smaller “slab” penalties preserve nonzero factor loadings. In the formulation given for BiSSLB, sparsity in $K^\star$5 and $K^\star$6 is identified directly with bicluster structure (Fan et al., 19 Mar 2026).

The inclusion indicators themselves follow IBP stick-breaking priors:

$K^\star$7

and similarly

$K^\star$8

The intercept prior is improper uniform, $K^\star$9. For computation, the hierarchy is truncated at $\widehat K$0, with the expectation that the active number of biclusters is pruned to a smaller $\widehat K$1 during optimization (Fan et al., 19 Mar 2026).

3. Posterior optimization and proximal coordinate ascent

BiSSLB is fitted by maximizing the log-posterior

$\widehat K$2

The optimization proceeds by coordinate ascent with proximal updates and momentum, rather than by a slow EM-type scheme (Fan et al., 19 Mar 2026).

For the $\widehat K$3-update, with $\widehat K$4 and $\widehat K$5 fixed, the smooth part of the objective is

$\widehat K$6

with gradient

$\widehat K$7

where

$\widehat K$8

BiSSLB then uses a proximal-gradient step with momentum (FISTA):

$\widehat K$9

Each entry is thresholded through the “refined” SSL proximal operator,

$Y=(y_{ij}) \in \{0,1\}^{I \times J}$0

where the adaptive penalty is

$Y=(y_{ij}) \in \{0,1\}^{I \times J}$1

and

$Y=(y_{ij}) \in \{0,1\}^{I \times J}$2

The $Y=(y_{ij}) \in \{0,1\}^{I \times J}$3-update follows the same recipe with $Y=(y_{ij}) \in \{0,1\}^{I \times J}$4 and $Y=(y_{ij}) \in \{0,1\}^{I \times J}$5 (Fan et al., 19 Mar 2026).

The intercept update is obtained by a minorize–maximize step:

$Y=(y_{ij}) \in \{0,1\}^{I \times J}$6

The sparsity weights are updated by

$Y=(y_{ij}) \in \{0,1\}^{I \times J}$7

and similarly

$Y=(y_{ij}) \in \{0,1\}^{I \times J}$8

After each cycle, columns are reordered by descending $Y=(y_{ij}) \in \{0,1\}^{I \times J}$9 and $$y_{ij} \mid A,B,\mu \sim \mathrm{Bernoulli}\!\left(\sigma\!\left(\mu_i + (AB^\top)_{ij}\right)\right),$$0, zero columns are dropped, and $$y_{ij} \mid A,B,\mu \sim \mathrm{Bernoulli}\!\left(\sigma\!\left(\mu_i + (AB^\top)_{ij}\right)\right),$$1 is reset. This update pattern is central to BiSSLB’s practical ability to infer the active bicluster count while maintaining sparse latent representations (Fan et al., 19 Mar 2026).

4. Computational profile and hyperparameterization

The dominant computational cost arises in the matrix multiplications needed for the factor updates. Each full $$y_{ij} \mid A,B,\mu \sim \mathrm{Bernoulli}\!\left(\sigma\!\left(\mu_i + (AB^\top)_{ij}\right)\right),$$2- or $$y_{ij} \mid A,B,\mu \sim \mathrm{Bernoulli}\!\left(\sigma\!\left(\mu_i + (AB^\top)_{ij}\right)\right),$$3-update requires computing $$y_{ij} \mid A,B,\mu \sim \mathrm{Bernoulli}\!\left(\sigma\!\left(\mu_i + (AB^\top)_{ij}\right)\right),$$4 or $$y_{ij} \mid A,B,\mu \sim \mathrm{Bernoulli}\!\left(\sigma\!\left(\mu_i + (AB^\top)_{ij}\right)\right),$$5, which costs $$y_{ij} \mid A,B,\mu \sim \mathrm{Bernoulli}\!\left(\sigma\!\left(\mu_i + (AB^\top)_{ij}\right)\right),$$6. The $$y_{ij} \mid A,B,\mu \sim \mathrm{Bernoulli}\!\left(\sigma\!\left(\mu_i + (AB^\top)_{ij}\right)\right),$$7-update costs $$y_{ij} \mid A,B,\mu \sim \mathrm{Bernoulli}\!\left(\sigma\!\left(\mu_i + (AB^\top)_{ij}\right)\right),$$8, and the $$y_{ij} \mid A,B,\mu \sim \mathrm{Bernoulli}\!\left(\sigma\!\left(\mu_i + (AB^\top)_{ij}\right)\right),$$9-updates cost $$M = \mu\,\mathbf 1^\top + A B^\top, \qquad A \in \mathbb R^{I \times K},\; B \in \mathbb R^{J \times K},\; \mu \in \mathbb R^I.$$0. Consequently, the overall per-iteration complexity is $$M = \mu\,\mathbf 1^\top + A B^\top, \qquad A \in \mathbb R^{I \times K},\; B \in \mathbb R^{J \times K},\; \mu \in \mathbb R^I.$$1 (Fan et al., 19 Mar 2026).

In practice, $$M = \mu\,\mathbf 1^\top + A B^\top, \qquad A \in \mathbb R^{I \times K},\; B \in \mathbb R^{J \times K},\; \mu \in \mathbb R^I.$$2 is chosen conservatively, with the summary reporting values such as 50–300, and columns are dropped as soon as they go to zero. Under this regime, the algorithm is described as tractable for matrices up to a few thousand by a few thousand on a modern workstation. Further speed-ups can be achieved by parallelizing the matrix–matrix multiplies or by GPU acceleration (Fan et al., 19 Mar 2026).

Hyperparameterization follows a specific template. The slab penalties are fixed at $$M = \mu\,\mathbf 1^\top + A B^\top, \qquad A \in \mathbb R^{I \times K},\; B \in \mathbb R^{J \times K},\; \mu \in \mathbb R^I.$$3. The spike penalties are set equal, $$M = \mu\,\mathbf 1^\top + A B^\top, \qquad A \in \mathbb R^{I \times K},\; B \in \mathbb R^{J \times K},\; \mu \in \mathbb R^I.$$4, and chosen by grid search over

$$M = \mu\,\mathbf 1^\top + A B^\top, \qquad A \in \mathbb R^{I \times K},\; B \in \mathbb R^{J \times K},\; \mu \in \mathbb R^I.$$5

The IBP intensities are fixed at $$M = \mu\,\mathbf 1^\top + A B^\top, \qquad A \in \mathbb R^{I \times K},\; B \in \mathbb R^{J \times K},\; \mu \in \mathbb R^I.$$6, with Beta second-parameters set to 1. The learning rate has default $$M = \mu\,\mathbf 1^\top + A B^\top, \qquad A \in \mathbb R^{I \times K},\; B \in \mathbb R^{J \times K},\; \mu \in \mathbb R^I.$$7 and may be tuned over

$$M = \mu\,\mathbf 1^\top + A B^\top, \qquad A \in \mathbb R^{I \times K},\; B \in \mathbb R^{J \times K},\; \mu \in \mathbb R^I.$$8

The truncation level $$M = \mu\,\mathbf 1^\top + A B^\top, \qquad A \in \mathbb R^{I \times K},\; B \in \mathbb R^{J \times K},\; \mu \in \mathbb R^I.$$9 should overestimate the expected number of biclusters, after which the algorithm prunes to an active $\mu$0 (Fan et al., 19 Mar 2026).

These design choices clarify an important practical point: BiSSLB’s automatic model-complexity selection is not hyperparameter-free. The active bicluster count is inferred from the data, but the quality and efficiency of inference remain influenced by spike strength, learning rate, and truncation. This is consistent with the paper’s explicit acknowledgment that multiple hyperparameters require tuning (Fan et al., 19 Mar 2026).

5. Empirical evaluation

The reported empirical study comprises two simulation settings and two biological applications, with comparisons against Bibit, Bibit2, Bimax, iBBiG, and GBC in simulation, and against Bibit, Bimax, iBBiG, and GBC in the HapMap SNP analysis (Fan et al., 19 Mar 2026).

Setting	Configuration	Reported result
Simulation I and II	Simulation I: $\mu$1, $\mu$2, $\mu$3, noise 0–20%; Simulation II: logistic generative noise, up to 50%	Under moderate/high noise, BiSSLB consistently had the highest CE/CS and the best recovery of $\mu$4
HapMap SNP	$\mu$5, true $\mu$6	BiSSLB recovered $\mu$7 and achieved highest CE = CS = 0.356
Homo sapiens PPI	$\mu$8 binary matrix	BiSSLB found $\mu$9, with predictive AUC = 0.965 and AUPR = 0.684

In Simulation I and Simulation II, the evaluation metrics were clustering error (CE), consensus score (CS), relevance, recovery, and correctness of $AB^\top$0. Under moderate and high noise, BiSSLB consistently had the highest CE and CS values and the best recovery of the true $AB^\top$1. The summary further reports that other methods either broke down, as in the case of iBBiG, or substantially over- or under-estimated the number of biclusters (Fan et al., 19 Mar 2026).

In the HapMap SNP application, the matrix size was $AB^\top$2 with true $AB^\top$3. BiSSLB recovered $AB^\top$4 and achieved the highest CE and CS, both equal to 0.356. It outperformed Bibit, Bimax, iBBiG, and GBC even when those methods were initialized at $AB^\top$5 (Fan et al., 19 Mar 2026).

In the Homo sapiens PPI application, the data formed a $AB^\top$6 binary matrix. BiSSLB returned $AB^\top$7 and achieved the best predictive AUC, 0.965, and AUPR, 0.684. The result was also described as much more compact than those of competing methods, which returned thousands of biclusters (Fan et al., 19 Mar 2026).

6. Strengths, limitations, and prospective extensions

The strengths attributed to BiSSLB are explicit. It is reported to provide robust recovery under moderate-to-high noise, a property linked to the explicit Bayesian uncertainty model. It automatically infers $AB^\top$8 through the IBP stick-breaking prior. It accommodates overlapping, variable-sized biclusters through sparse factorization. Its coordinate-ascent procedure with proximal and momentum updates is presented as an efficient alternative to slow EM-type convergence (Fan et al., 19 Mar 2026).

The limitations are equally explicit. The per-iteration cost $AB^\top$9 can be high for very large matrices, although early column-dropping mitigates this in practice. The model relies on logistic factorization, and real data may exhibit more complex couplings. Multiple hyperparameters require tuning, specifically including spike strength, learning rate, and truncation. In addition, the posterior is multimodal; although the refined proximal step prunes many local modes, global optimality is not guaranteed (Fan et al., 19 Mar 2026).

These points also clarify what BiSSLB should not be taken to imply. It is not restricted to non-overlapping or fixed-size biclusters, and it does not require prior knowledge of the noise level or bicluster characteristics. At the same time, it does not eliminate the need for practical model-selection decisions, nor does it resolve the general nonconvexity of latent-factor biclustering (Fan et al., 19 Mar 2026).

Several extensions are proposed directly. The Bernoulli likelihood could be replaced with Poisson or negative-binomial models with a log-link for count data such as mutational signatures. Mixed-type biclustering could be obtained by composing several generalized linear model likelihoods for binary, count, and Gaussian observations. Massive matrices could motivate stochastic or mini-batch proximal updates. Side-information, including network or genomic annotations, could be incorporated through structured priors in the IBP intensities (Fan et al., 19 Mar 2026). These proposals indicate that the BiSSLB framework is defined not only by a particular binary-data model, but also by a modular combination of sparse factorization, adaptive shrinkage, and nonparametric control over the number of latent biclusters.