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 $PR(X, W, \alpha)$ is the polynomial regression loss. Regularization weights $\beta_1$ and $\beta_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 $\mathcal{O}(d^3)$ cost—by employing the spectral radius. For $A \in \mathbb{R}^{d \times d}_{\ge 0}$, acyclicity is guaranteed if $\rho(A) = 0$, since $\mathrm{tr}(A^k) = 0$ for all $k=1, ..., d$ if and only if all eigenvalues vanish.

The spectral radius is approximated efficiently via power iteration: add a small positive perturbation $\epsilon \mathbf{1}$ to $A$, then propagate vectors $v$ and $u$ repeatedly through $A$ and $A^\top$, respectively, normalizing at each step. The leading eigenvalue approximation is $(u^\top \hat{A} v)/(u^\top v)$. The gradient with respect to each $W_{ij}$ is computed using Perron–Frobenius theory: $$\frac{\partial}{\partial W_{ij}} \rho(W \odot W + \epsilon \mathbf{1}) = 2 W_{ij} \frac{v_i u_j}{u^\top v}$$ Both spectral radius evaluation and its gradient scale as $\mathcal{O}(d^2)$, 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 $K$ (set to 3 in experiments): $f_j(x) = \sum_{k=0}^{K} \alpha_{jk} x^k$ The loss function is

$$PR(X, W, \alpha) = \sum_{i=1}^d \sum_{l=1}^n \left[ X_{il} - \sum_{j=1}^d W_{ij} \sum_{k=0}^{K} \alpha_{jk} (X_{jl})^k \right]^2$$

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 $\xi \ge 0$ and penalty $\eta > 0$, each outer iteration alternates between:

1. Inner Unconstrained Minimization over $(W, \alpha)$ using first-order optimizers such as Adam, minimizing: $$F(W, \alpha; \xi, \eta) = \frac{1}{2n} PR(X, W, \alpha) + \beta_1 ( \| W \|_F^2 + \| W \|_1 ) + \beta_2 \| \alpha \|_F^2 + \xi \rho(W \odot W + \epsilon) + \frac{\eta}{2} \rho(W \odot W + \epsilon)^2$$
2. Multiplier Update:

$$\xi \leftarrow \xi + \eta \cdot \rho(W \odot W), \quad \eta \leftarrow 1.1 \cdot \eta$$

Stopping criteria are $\rho(W \odot W) < t_0$ (e.g., $5 \times 10^{-3}$) 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 $\mathcal{O}(nd^2)$ per inner iteration (for $K \ll d$). 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 $>$10× 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 ($d=100,300$, $n=100$ to $2,000$), average precision (AP) strongly favors NO-BEARS above $n \ge 500$, with improvements up to +10 AP points over NO-TEARS and $>20$ 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 $n \gtrsim d$.
• GPU memory constraints limit single-card runs to $\sim12,800$ genes, but multi-GPU or sparse extensions could feasibly scale to $>20,000$ 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 $\mathcal{O}(d^3)$ bottleneck, introducing adaptivity for nonlinear regulation, and leveraging efficient GPU parallelism (Lee et al., 2019).