There are two separate issues here.
- How to use efficient solvers for $Ax=b$ in order to apply $A^{1/2}b$.
- How to compute the determinant.
The short answers are 1) use rational matrix function approximations, and 2) you don't, but you don't need to anyways. I address both of these issues below.
Matrix square root approximations
The idea here is to convert a rational function approximation for scalar functions into a rational function approximation for matrix functions.
We know that there exist rational functions that can approximate the square root function extremely well,
$$\sqrt{x} \approx r(x) := \frac{a_1}{x+b_1} + \frac{a_2}{x+b_2} + \dots + \frac{a_N}{x+b_N},$$
for positive $b_i$. Indeed, to get high accuracy on the interval $[m,M]$, you need $O(\log \frac{M}{m})$ terms in the series. To get the appropriate weights ($a_i$) and poles ($-b_i$), just look up rational function approximation online or in a book.
Now consider applying this rational function to your matrix:
$$r(A) = a_1(A + b_1 I)^{-1} + a_2(A + b_2 I)^{-1} + \dots + a_N(A + b_N I)^{-1}.$$
Owing to the symmetry of $A$, we have
\begin{align}
||A^{1/2} - r(A)||_2 &= ||U\left(\Sigma^{1/2} - r(\Sigma)\right)U^*||_2, \\
&= \max_i |\sqrt{\sigma_i} - r(\sigma_i)|
\end{align}
where $A = U \Sigma U^*$ is the singular value decomposition (SVD) of $A$. So, the quality of the rational matrix approximation is equivalent to the quality of the rational function approximation at the location of the eigenvalues.
Denoting the condition number of $A$ by $\kappa$, we can apply $A^{1/2}b$ to any desired tolerance by performing $O(\log \kappa)$ positively shifted graph Laplacian solutions of the form,
$$(A + bI)x=b.$$
These solutions can be done with your favorite graph Laplacian solver—I prefer multigrid type techniques, but the one in the paper you cite should be fine, too. The extra $bI$ only helps the convergence of the solver.
For an excellent paper discussing this, as well as more general complex analysis techniques that apply to nonsymmetric matrices, see Computing $A^α$, $\log(A)$, and related matrix functions by contour integrals, by Hale, Higham, and Trefethen (2008).
Determinant "computation"
The determinant is harder to compute. As far as I know, the best way is to compute the Schur decomposition $A = Q U Q^*$ using the QR algorithm, then read off the eigenvalues from the diagonal of the upper-triangular matrix $U$. This takes $O(n^3)$ time, where $n$ is the number of nodes in the graph.
However, calculating determinants is an inherently ill-conditioned problem, so if you ever read a paper that relies on computing determinants of a large matrix, you should be very skeptical of the method.
Luckily, you probably don't actually need the determinant. For example,
- To draw samples from a single Gaussian distribution $N(0,A^{-1})$, the normalization constant is the same at all points so you never need to compute it.
- If your Laplacian matrix $A = A_x$ represents the inverse covariance of a local Gaussian approximation at point $x$ to a non-Gaussian distribution, then the determinant does indeed change from point to point. However, in every effective sampling scheme I know (including Markov chain Monte Carlo, importance sampling, etc.) what you really need is the determinant ratio,
$$\det(A_{x_0}^{-1}A_{x_p}),$$
where $x_0$ is the current point, and $x_p$ is the proposed next sample.
We can view $A_{x_0}^{-1}A_{x_p}$ as a low-rank update to the identity,
$$A_{x_0}^{-1}A_{x_p} = I + Q D Q^*,$$
where the effective numerical rank, $r$, of the low-rank update is a local measure of how non-Gaussian the true distribution is; typically this is much lower than the full rank of the matrix. Indeed, if $r$ is large, then the true distribution is locally so non-Gaussian that one ought to question the whole strategy of trying to sample this distribution using local Gaussian approximations.
The low-rank factors $Q$ and $D$ can be found with randomized SVD or Lanczos by applying the matrix
$$A_{x_0}^{-1}A_{x_p} -I$$
to $O(r)$ different vectors, each application of which requires one graph Laplacian solution. Thus the overall work for getting these low rank factors is $O(r \max(n,E))$.
Knowing $D = \text{diag}(d_1,d_2,\dots,d_r)$, the determinant ratio is then
$$\det(A_{x_0}^{-1}A_{x_p}) = \det(I + Q D Q^*) = \exp\left(\sum_{i=1}^r \log d_i\right).$$
These low rank determinant ration calculation techniques can be found in A Stochastic Newton MCMC Method for Large-Scale Statistical Inverse Problems with Application to Seismic Inversion, by Martin, et al. (2012). In this paper it is applied to continuum problems, so the "graph" is a grid in 3D space and the graph Laplacian is the actual Laplacian matrix. However, all the techniques apply to general graph Laplacians. There are probably other papers applying this technique to general graphs by now (the extension is trivial and basically what I just wrote).
