  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.
