Scalable Parallel Factorizations of SDD Matrices and Efficient Sampling for Gaussian Graphical Models
This work addresses a fundamental problem in computational statistical inference and spectral graph theory, providing an efficient solution for sampling from Gaussian graphical models, which is incremental but with strong performance gains.
The paper tackles the problem of efficiently sampling from Gaussian random fields with SDDM precision matrices by developing a nearly optimal algorithm that factors SDDM matrices into sparse operators, enabling sampling with nearly linear work and minimal randomness. The result achieves the first sampling algorithm requiring only nearly linear work and n i.i.d. Gaussian samples, which is optimal in randomness and nearly optimal in work and parallel complexity.
Motivated by a sampling problem basic to computational statistical inference, we develop a nearly optimal algorithm for a fundamental problem in spectral graph theory and numerical analysis. Given an $n\times n$ SDDM matrix ${\bf \mathbf{M}}$, and a constant $-1 \leq p \leq 1$, our algorithm gives efficient access to a sparse $n\times n$ linear operator $\tilde{\mathbf{C}}$ such that $${\mathbf{M}}^{p} \approx \tilde{\mathbf{C}} \tilde{\mathbf{C}}^\top.$$ The solution is based on factoring ${\bf \mathbf{M}}$ into a product of simple and sparse matrices using squaring and spectral sparsification. For ${\mathbf{M}}$ with $m$ non-zero entries, our algorithm takes work nearly-linear in $m$, and polylogarithmic depth on a parallel machine with $m$ processors. This gives the first sampling algorithm that only requires nearly linear work and $n$ i.i.d. random univariate Gaussian samples to generate i.i.d. random samples for $n$-dimensional Gaussian random fields with SDDM precision matrices. For sampling this natural subclass of Gaussian random fields, it is optimal in the randomness and nearly optimal in the work and parallel complexity. In addition, our sampling algorithm can be directly extended to Gaussian random fields with SDD precision matrices.