PaperlayerView exact revision on arXiv

Source-linked AI summary

Massive parallelization of serial inference algorithms for a complex generalized linear model

Marc A. Suchard, Shawn E. Simpson, Ivan Zorych, Patrick Ryan, David Madigan

arXiv:1208.0945v1stat.COmath.OC

TL;DR

Massive observational databases make drug–adverse-event association estimates with covariate adjustment computationally challenging. The paper parallelizes cyclic coordinate descent for a Bayesian conditioned GLM; on its largest dataset, sparse GPU operations cut MAP estimation from 51 hours to 29 seconds, making cross-validation and bootstrap practical.

  • Problem

    Massive observational databases make estimating drug–adverse-event associations while controlling for covariates computationally difficult, and naive MAP fitting nearly stalls with millions of outcomes and thousands of predictors.

    Source: §1 · pp. 1–2

  • Method

    The paper optimizes and massively parallelizes cyclic coordinate descent for a Bayesian conditioned GLM, vectorizing computation across GPU data.

    Source: §2.4 · p. 6, §4 · p. 16

  • Results

    51 hours to 29 seconds: sparse GPU operations reduce MAP estimation time on the largest dataset, making cross-validation and bootstrap practical.

    Source: §3 · p. 14, §4 · p. 15

  • Takeaways & Limitations

    Cross-validation and bootstrapping can estimate model hyperparameters and report uncertainty around point estimates.

    Source: §4 · p. 15

    Coefficient uncertainty relies on a simple nonparametric bootstrap, while computationally efficient, theoretically well-supported standard-error procedures for L1 or L2 regularization remain out of reach.

    Source: §2.5 · p. 12

Abstract

from arXiv · show

Following a series of high-profile drug safety disasters in recent years, many countries are redoubling their efforts to ensure the safety of licensed medical products. Large-scale observational databases such as claims databases or electronic health record systems are attracting particular attention in this regard, but present significant methodological and computational concerns. In this paper we show how high-performance statistical computation, including graphics processing units, relatively inexpensive highly parallel computing devices, can enable complex methods in large databases. We focus on optimization and massive parallelization of cyclic coordinate descent approaches to fit a conditioned generalized linear model involving tens of millions of observations and thousands of predictors in a Bayesian context. We find orders-of-magnitude improvement in overall run-time. Coordinate descent approaches are ubiquitous in high-dimensional statistics and the algorithms we propose open up exciting new methodological possibilities with the potential to significantly improve drug safety.

1 Motivation and Background

Large observational databases offer important data for post-marketing drug safety studies, but their scale makes estimating drug–adverse-event associations while adjusting for covariates computationally difficult. The paper explores optimizing and parallelizing this fitting task with GPUs.

1 Motivation and Background

  • Large claims databases can contain over 50 million lives, up to 10 years of data per life, and exposure to thousands of drugs.These observational records are used to associate adverse events with drug use after products enter widespread clinical use.Source: §1 · p. 1
  • Estimating each drug–adverse-event association while adjusting for concomitant drugs, demographics, and disease is computationally challenging at this scale.Source: §1 · pp. 1–2
  • Even maximum a posteriori fitting with regularized GLMs can nearly stall in standard packages at millions of outcomes and thousands of predictors, while coefficient-uncertainty estimates become out of reach.Source: §1 · p. 2
  • Distributed clusters often cannot offset communication latency and thread-management costs during iterative GLM fitting, including for MAP estimates and standard errors.Source: §1 · p. 2
  • The paper optimizes cyclic coordinate descent for a Bayesian self-controlled case series model by exploiting sparse design matrices and GPU parallelization.The authors also report rough prior-hyperparameter estimates and limited coefficient-uncertainty measures.Source: §1 · p. 2

2 Methods

The methods develop a Bayesian self-controlled case series model and adapt cyclic coordinate descent for efficient, parallel fitting on large observational databases. They describe GPU-based sparse updates with hardware-tailored storage and precision, cross-validation for selecting σ², and bootstrap estimation of coefficient uncertainty.

2.1 Bayesian Self-Controlled Case Series Model

  • The BSCCS models era-specific adverse-event counts as Poisson outcomes modulated by drug exposures, then conditions on each subject’s total events so only cases contribute.Source: §2.1 · pp. 3–4
  • Normal and Laplacian priors shrink drug effects toward zero, yielding posterior modes analogous to ridge (L2) and lasso (L1) conditioned-Poisson regression.Source: §2.1 · p. 4

2.2 Maximum A Posteriori Estimation using Cyclic Coordinate Descent

  • Cyclic coordinate descent updates one coefficient at a time without inverting the full Hessian, using an adaptive trust-region bound and a convergence threshold of 0.0005.Source: §2.2 · pp. 4–5

2.3 Computational Work

  • A sparse K×J era-by-drug design and elementwise operations express the likelihood and expose parallelizable calculations for one-dimensional gradients and Hessians.Source: §2.3 · pp. 5–6

2.4 Targets for Parallelization

  • Even for the smallest database studied, over 99.5% of runtime is spent computing coordinate gradients and Hessians rather than serial coefficient updates.Source: §2.4 · p. 6
  • Sparse updates extend beyond Xβ to subject-level gradient and Hessian components, while GPU parallelization targets both dense exponentiation-heavy updates and sparse fused reductions.Source: §2.4 · p. 7
  • Fusing transformations with simultaneous reductions avoids materializing W and repeatedly moving it through memory, reducing memory-bandwidth demands on CPUs and GPUs.Source: §2.4 · p. 8
  • GPU off-loading reduces each cycle step’s transfer from 2 × N floating-point numbers to a single floating-point ∆βj, while atomic updates avoid race conditions and fuse sparse update steps.The sparse CUDA kernel uses one thread per non-zero entry of Xj.Source: §2.4 · p. 10
  • Sparse storage is matched to access patterns: CCS for columns of X, CRS for dense-X multiplication by M, and COO with structure-of-arrays storage for precomputed MXj when X is sparse.Because the non-zero entries of X and N are all one, they need not be stored.Source: §2.4 · p. 10
  • Single precision suffices for massive BSCCS fitting because the MAP problem is convex and the computations avoid subtracting nearly equal or multiplying small quantities.The paper notes that hj(β) < 0 and every element of [W × (1 − W)] is non-negative.Source: §2.4 · p. 12

2.5 Hyperparameter Selection and Measures of Coefficient Uncertainty

  • A 10-fold cross-validation scheme selects σ^2 by maximizing held-out predictive log-likelihood across a log-scale grid, using warm-starts as prior variance increases.The authors note that the predictive log-likelihood surface is relatively flat near its maximum, so precise hyperparameter estimation is unnecessary.Source: §2.5 · p. 12
  • For coefficient uncertainty, the authors use a pragmatic non-parametric bootstrap: 200 samples provide 95% quantile confidence intervals and, under a Laplacian prior, observed non-zero MAP proportions as a poor-man’s posterior probability proxy.They note that computationally efficient and theoretically well-supported uncertainty procedures for L1- or L2-regularized estimates remain out of reach.Source: §2.5 · p. 12

3 Demonstration

The demonstrations benchmark BSCCS fitting across large observational databases and adverse events, with datasets spanning millions of cases-only patients and tens of millions of exposure eras. Sparse GPU computation sharply reduces MAP-estimation time, while cross-validation and bootstrap analysis illustrate hyperparameter selection and uncertainty estimation.

3 Demonstration

  • The benchmarks covered two claims databases and four adverse events, with 115K–3.6M cases-only patients, 1,224–1,428 drugs, and 3.8M–75M exposure eras per analysis.Benchmarking used an Intel Xeon X5570 CPU and one NVIDIA Tesla C2050 GPU.Source: §3 · p. 13
  • Up to 181-fold speed-up from sparsity and up to another 37-fold from GPU computation; the largest dataset’s MAP fit fell from over 51 hours to 29 seconds.The comparison covers dense CPU, sparse CPU, and sparse GPU implementations.Source: Figure 2 · p. 13, §3 · p. 14
  • Cross-validation selected variances of 0.05–0.15 for L1 and 0.02–0.13 for L2; in the angioedema example, Drotrecogin alfa’s fourth-largest effect estimate had a 95% interval covering zero.The bootstrap analysis plotted intervals for 441 drug effects with nonzero coefficients in more than half of replicates, with shading indicating the proportion of nonzero estimates.Source: §3 · p. 14, Figure 3 · p. 14

4 Discussion

The discussion attributes the speed gains to sparse-model optimization and GPU vectorization within serial coordinate descent. The reduced run-time makes cross-validation and bootstrapping feasible, while the approach may extend to richer priors and other optimization problems.

4 Discussion

  • The BSCCS achieves an approximate 30-fold speed-up on a single GPU versus a single CPU core, while sparse design-matrix calculations yield over a 100-fold speed-up versus the most advanced CCD algorithms considered.The over-100-fold comparison concerns sparse subject-specific gradient and Hessian contributions; the approximate 30-fold comparison is for the BSCCS model overall.Source: §4 · p. 15
  • Hierarchical priors could borrow strength across drugs grouped by mode of action or therapeutic target, or across related outcomes, with run-times expected to remain quick.The paper considered only independent and identically distributed priors over drug effect sizes; prior gradient and Hessian computations remain lightweight.Source: §4 · p. 15
  • Orders-of-magnitude lower run-time makes nested cross-validation and bootstrapping feasible for estimating model hyperparameters and measures of uncertainty around point estimates.A preliminary R implementation required around 5.3 hours for one MAP estimate on a dataset with N = 7460.Source: §4 · p. 15
  • The serial CCD approach decomposes high-dimensional optimization into low-dimensional updates; one-dimensional updates are preferred because higher-dimensional operations can incur costly data dependencies.The authors envision extensions beyond GLMs, while noting related GPU applications to matrix factorization and multidimensional scaling.Source: §4 · p. 15
  • GPU vectorization performs branch-free operations across large data arrays, with speed-up scaling most directly with data quantity rather than the device’s critical data dimension.Unlike distributing independent tasks, this strategy avoids thread divergence but does not scale with the critical dimension as parallel devices grow.Source: §4 · p. 16
Loading 1208.0945v1…