NO-BEARS: GRN Inference Algorithm

• NO-BEARS is a continuous optimization algorithm that infers causal gene regulatory networks by extending NO-TEARS with nonlinear polynomial regression.
• It employs a spectral radius-based acyclicity constraint with O(d²) complexity, enabling near real-time inference on large-scale transcriptomic datasets.
• Experimental results demonstrate that NO-BEARS outperforms NO-TEARS and GENIE3 in speed and precision while effectively modeling higher-order gene interactions.

NO-BEARS (NEw Optimizer for B serà უდReAring Structures) is a continuous-optimization algorithm for estimating directed acyclic gene regulatory networks (GRNs) from transcriptomic data. Developed as an extension of the NO-TEARS framework, NO-BEARS introduces a computationally efficient acyclicity constraint based on the spectral radius and incorporates polynomial regression to model nonlinear gene interactions. Implementation on modern GPU hardware allows for near real-time inference at previously intractable scales, facilitating the recovery of causal transcriptome networks involving thousands of genes (Lee et al., 2019).

1. Optimization Problem Formulation

Let $X \in \mathbb{R}^{n \times d}$ represent the gene expression matrix for $n$ samples and $d$ genes, with $W \in \mathbb{R}^{d \times d}$ as the weighted adjacency matrix to be inferred. NO-BEARS parameterizes per-gene nonlinearities with polynomials of degree $K$ via coefficients $\alpha \in \mathbb{R}^{d \times (K+1)}$.

The algorithm solves the constrained optimization problem: $$\min_{W,\alpha} \frac{1}{2n} PR(X, W, \alpha) + \beta_{1}\left( \left\| W \right\|_F^2 + \left\| W \right\|_{1} \right) + \beta_{2}\left\| \alpha \right\|_F^2$$ subject to

$\rho(W \odot W) = 0$

where $\odot$ denotes elementwise multiplication, $\rho(\cdot)$ is the spectral radius, and $n$0 is the polynomial regression loss. Regularization weights $n$1 and $n$2 control model complexity.

2. Acyclicity Enforcement via the Spectral Radius

NO-BEARS improves upon NO-TEARS’ acyclicity constraint—which relies on the matrix exponential’s trace and incurs $n$3 cost—by employing the spectral radius. For $n$4, acyclicity is guaranteed if $n$5, since $n$6 for all $n$7 if and only if all eigenvalues vanish.

The spectral radius is approximated efficiently via power iteration: add a small positive perturbation $n$8 to $n$9, then propagate vectors $d$0 and $d$1 repeatedly through $d$2 and $d$3, respectively, normalizing at each step. The leading eigenvalue approximation is $d$4. The gradient with respect to each $d$5 is computed using Perron–Frobenius theory: $d$6 Both spectral radius evaluation and its gradient scale as $d$7, enabling scalability.

3. Polynomial Regression Loss for Nonlinear Gene Interactions

To model nonlinear gene dependency, NO-BEARS replaces the linear structural equation model with a per-gene polynomial of degree $d$8 (set to 3 in experiments): $d$9 The loss function is

$W \in \mathbb{R}^{d \times d}$0

This construction enables recovery of higher-order regulatory relationships, which are pervasive in biological networks.

4. Augmented Lagrangian Optimization and Algorithmic Workflow

The constrained optimization is addressed using an augmented Lagrangian/penalty multiplier approach. For a current Lagrange multiplier $W \in \mathbb{R}^{d \times d}$1 and penalty $W \in \mathbb{R}^{d \times d}$2, each outer iteration alternates between:

1. Inner Unconstrained Minimization over $W \in \mathbb{R}^{d \times d}$3 using first-order optimizers such as Adam, minimizing: $W \in \mathbb{R}^{d \times d}$4
2. Multiplier Update:

$W \in \mathbb{R}^{d \times d}$5

Stopping criteria are $W \in \mathbb{R}^{d \times d}$6 (e.g., $W \in \mathbb{R}^{d \times d}$7) and regression-loss rebound (details in Section 5.4 of the paper).

The entire workflow parallelizes efficiently on the GPU, with all major matrix operations dispatched to hardware-accelerated libraries.

5. Computational Complexity and GPU Implementation

The primary computational bottlenecks are polynomial regression loss evaluation and acyclicity constraint checking, both scaling as $W \in \mathbb{R}^{d \times d}$8 per inner iteration (for $W \in \mathbb{R}^{d \times d}$9). GPU implementation in TensorFlow enables rapid runtime:

# Genes	NO-BEARS (s)	NO-TEARS (s)	GENIE3 (s)
100	0.96	1.05	1,200
3,200	1.14	36.97	89,431.2
12,800	9.70	OOM	OOM

OOM: out of memory. NO-BEARS achieves speedups up to $K$010× over NO-TEARS and several orders of magnitude over CPU-based GENIE3, enabling tractable inference for thousands of genes.

6. Experimental Results: Accuracy, Precision, and Resource Use

Extensive benchmarking on synthetic (SynTReN) and realistic (GeneNetWeaver) networks demonstrates the empirical advantages of NO-BEARS over NO-TEARS and GENIE3:

• On SynTReN ($K$1, $K$2 to $K$3), average precision (AP) strongly favors NO-BEARS above $K$4, with improvements up to +10 AP points over NO-TEARS and $K$5 over GENIE3.
• On GeneNetWeaver:
  • Ecoli-1565: AP = 36.1% (NO-BEARS), 35.2% (NO-TEARS), 12.2% (GENIE3); ROC = 78.7%, 83.0%, 88.0%; time to solution = 154.9 s (NO-BEARS), 352.5 s (NO-TEARS), 5,538.5 s (GENIE3).
  • Yeast-4441: AP = 64.5% (NO-BEARS), 56.9% (NO-TEARS), 9.8% (GENIE3); ROC = 95.6%, 88.9%, 92.9%; time = 2,484.1 s, 4,842.8 s, 89,431.2 s, respectively.

NO-BEARS also delivers consistent AP improvement through acyclicity enforcement over initial unconstrained fits, especially in imbalanced graph settings.

7. Limitations and Future Directions

NO-BEARS’ key limitations are dictated by data and hardware scale:

• Reliable edge recovery requires $K$6.
• GPU memory constraints limit single-card runs to $K$7 genes, but multi-GPU or sparse extensions could feasibly scale to $K$8 genes.
• Real biological data introduces confounders, latent variables, and measurement noise, necessitating future developments in robustness (e.g., latent-variable models) and richer nonlinearities (splines, neural nets).
• Prospective extensions include prior-informed regularization (e.g., known transcription factor-target relationships) and streaming inference for time-series or single-cell datasets.

The algorithm establishes a computational framework for full-transcriptome causal network discovery by removing the previous $K$9 bottleneck, introducing adaptivity for nonlinear regulation, and leveraging efficient GPU parallelism (Lee et al., 2019).