Sparse Convex Biclustering (SpaCoBi)

• Sparse Convex Biclustering (SpaCoBi) is a convex optimization–based method that integrates row/column fusion with group-lasso sparsity to uncover biclusters in high-dimensional data.
• It leverages the Sylvester equation and ADMM for efficient optimization, ensuring global optimality and robustness against noise.
• Empirical results on simulated and transcriptomic datasets show high adjusted Rand Index scores and effective feature selection compared to traditional methods.

Sparse Convex Biclustering (SpaCoBi) is a convex optimization–based method for simultaneous clustering of the rows and columns of high-dimensional data matrices, with integrated feature selection via group-lasso sparsity. SpaCoBi addresses limitations in existing biclustering approaches by directly penalizing noise in features, maintaining global optimality, and employing a stability-based criterion for hyperparameter tuning. Its design yields accurate and robust bicluster recovery in high-dimensional and large-scale applications, as demonstrated on simulated and transcriptomic datasets (Jiang et al., 5 Jan 2026).

1. Mathematical Formulation

Let $X \in \mathbb{R}^{n \times p}$ be a data matrix, where $X_{i\cdot}$ denotes the $i$th row and $x_j$ the $j$th column. SpaCoBi fits a matrix $A \in \mathbb{R}^{n \times p}$, simultaneously biclustering rows and columns while enforcing column-wise sparsity. The method solves the convex program: $$\min_{A\in\mathbb R^{n\times p}} \frac12\sum_{i=1}^n\|X_{i\cdot}-A_{i\cdot}\|_2^2 +\gamma_1\sum_{i<j}w_{ij}\|A_{i\cdot}-A_{j\cdot}\|_2 +\gamma_2\sum_{k<\ell}\tilde w_{k\ell}\|A_{\cdot k}-A_{\cdot\ell}\|_2 +\gamma_3\sum_{j=1}^p u_j\|A_{\cdot j}\|_2$$ where:

• $\frac12\|X-A\|_F^2$ enforces data fidelity,
• Row-fusion term: $\sum_{i<j}w_{ij}\|A_{i\cdot}-A_{j\cdot}\|_2$,
• Column-fusion term: $$\sum_{k<\ell}\tilde w_{k\ell}\|A_{\cdot k}-A_{\cdot\ell}\|_2$$,
• Group-lasso column sparsity term: $\sum_{j=1}^p u_j\|A_{\cdot j}\|_2$.

Introducing auxiliary variables $v_{ij}$ (row pairs), $z_{k\ell}$ (column pairs), and $g_j$ (column groups), SpaCoBi can be written in constrained form with convex objectives and linear constraints linking $A$ and the auxiliary variables.

2. Convexity and Optimization

Each objective term is convex: the quadratic loss is strictly convex in $A$; fusion and group-lasso terms are convex norms. The global solution is unique.

Optimization proceeds via the Alternating Direction Method of Multipliers (ADMM) with the following major steps:

1. $A$-update (Sylvester equation): The matrix $A$ is updated by solving the Sylvester equation $MA + AN = H$, with $M$ and $N$ constructed from row/column graph Laplacian structures and penalty parameters.
2. Proximal updates for $v_{ij}$, $z_{k\ell}$, and $g_j$ enforce respective fusions and sparsity via $\ell_2$-norm proximal operators.
3. Dual variable updates for Lagrange multipliers ensure convergence.

Efficient solution of the Sylvester equation relies on Bartels–Stewart–type or modified Schur methods, exploiting the structure of $M$ and $N$ for computational gains. Multi-block ADMM convergence is guaranteed under mild conditions.

3. Stability-Based Tuning

SpaCoBi hyperparameters $\gamma_1$ (row fusion), $\gamma_2$ (column fusion), and $\gamma_3$ (sparsity) control the tradeoff between fit, cluster granularity, and feature selection. To efficiently select these, SpaCoBi may collapse $(\gamma_1, \gamma_2)$ into a single $\gamma$ and perform a grid search over $(\gamma, \gamma_3)$.

Stability selection is employed: two bootstrap samples of rows yield biclustering solutions $(\psi_1, \psi_2)$ at a given $(\gamma, \gamma_3)$. The clustering distance $d_F$ is computed: $$d_F(\psi_1, \psi_2) = \mathbb E_{x, y \sim F} \left| I\{\psi_1(x) = \psi_1(y)\} - I\{\psi_2(x) = \psi_2(y)\} \right|$$ Estimated from resampled pairs, the $(\gamma, \gamma_3)$ minimizing $d_F$ is selected for maximal stability.

4. Computational Complexity and Scaling

The per-iteration cost is dominated by the Sylvester solve in the $A$-update. Naively, this requires $O(n^3 + p^3)$ time, but this is mitigated by:

• The Laplacian structure of $M$ and $N$,
• Fast generalized Schur/Bartels–Stewart algorithms.

Proximal updates for $v$, $z$, and $g$ scale with the size of row and column edge sets. Practical implementation uses $m$-nearest neighbor graphs, rendering $|\mathcal E_1| = O(mn)$ and $|\mathcal E_2| = O(mp)$ for small $m$.

A warm-starting strategy—using solutions from nearby parameter grid points as initializations—accelerates tuning by $20\%$–$100\%$ on large-scale problems.

5. Empirical Results and Benchmarking

Simulation studies using synthetic checkerboard biclusters and known informative columns ($p_\mathrm{true}$) demonstrate:

• Adjusted Rand Index (ARI): SpaCoBi achieves mean ARI in $[0.75, 0.96]$ versus Bi-ADMM $L_2$-norm $[0.12, 0.80]$; COBRA approaches zero in high noise.
• Feature-selection: False Negative Rate $=0$–$0.07$, False Positive Rate $=0.04$–$0.27$, AUC $=0.76$–$0.90$. Bi-ADMM (without sparsity) FPR $=1$.

In mouse olfactory bulb (MOB) single-cell RNA-seq data ($n=305$, $p=1250$) with known three-class structure:

• SpaCoBi recovers clusters perfectly (ARI $= 1.0$) and selects marker genes such as Pbxip1, Pdlim2, and Isg15.
• Bi-ADMM $L_2$-norm yields ARI $= 0.12$ and cannot suppress noise features.

6. Implementation and Recommended Practices

Selection of fusion weights and group-lasso factors is critical:

• Row/Column weights: $m$-nearest neighbor Gaussian kernel,

$$w_{ij} = \mathbf{1}\{j \in \mathrm{NN}_m(i)\} \exp(-\phi \|X_{i\cdot} - X_{j\cdot}\|_2^2)$$

with $m=5$, $\phi=0.5$.

• Group-lasso weights: Adaptive $u_j = 1 / \|a_j^{(0)}\|_2$, using the solution with $\gamma_3=0$, to penalize uninformative features.
• Rescaling: Normalize $\{w\}$, $\{\tilde w\}$, $\{u\}$ so that sums of parameters are $1/\sqrt{p}$, $1/\sqrt{n}$, $1/\sqrt{n}$, ensuring comparable tuning parameter magnitudes.

Limitations arise when both $n$ and $p$ are very large, as the Sylvester solve becomes a bottleneck; approximate methods or Block-ADMM provide possible remedies. Extensions to other sparsity norms (e.g., $\ell_1$) or overlapping-group penalties can be considered, provided convexity and ADMM compatibility are retained.

7. Context and Directions

SpaCoBi unifies row/column fusion and group-lasso feature selection within a convex optimization paradigm, ensuring global optimality and robust bicluster detection. Its empirical superiority over non-sparse convex biclustering is particularly pronounced in high-dimensional, noisy settings. Potential extensions include more general sparsity-inducing penalties and scalable iterative solvers, motivating ongoing research for truly massive omics applications (Jiang et al., 5 Jan 2026).