I have a Bayesian network where the number of nodes is potentially large. I've conditioned on some of the nodes (observed data) and I'm trying to draw samples from the distribution remaining nodes (often done with markov chain monte carlo). It turns out that the having the covariance matrix, for the remaining nodes is very useful for most MCMC algorithms because it summarizes the scale and orientation of the joint distribution of the nodes.

However, if the number of nodes, $n$, is large then the covariance can be very expensive to work with since it has $n^2$ elements. If you want to estimate it you'll need $O(n^2)$ points, and to multiply a vector with the matrix, you'll need $O(n^2)$ operations.

However, in most Bayesian networks many variables will be conditionally independent from each other because they will not be in each other's Markov Blanket. The covariance matrix should be fully specifiable by the covariances between variables who are in each others markov blankets.

Because of this, it seems like it should be possible to do many operations more cheaply then by finding the full covariance matrix and working with it. Specifically, the operations that usually come up are

- Draw a random vector from a normal distribution with the covariance matrix (or its inverse)
- Do a dot product of the covariance matrix with a vector.

I would expect there to be an algorithm for doing these two operations with $O(n)$ operations as long as the average Markov blanket size isn't increasing as the dimension increases.

This seems like a problem that has probably been solved before, but I haven't been able to find information on it. I am familiar with the MCMC literature, so I am fairly sure that this is not a technique that gets used on Bayesian problems.

**Example**:

I know how to solve this problem for some simple cases, but not how to solve the general problem.

If you have random variables $A,B,C$ and the network $B\leftarrow A \rightarrow C$

(meaning A is a parent of B and C, but B and C are independent of each other conditioning on A)

Then the B and C are only correlated with each other only because they're both correlated with A. Thus, if we know the variance of A, B and C and the variance of A with B and A with C we can calculate the covariance of B and C.

I believe that

$$ Cov(B, C) = \frac{Cov(A,B) \cdot Cov(A,C)}{Cov(A,A)} $$

in this case. More generally, if A is a vector of random variables, I think this should be

$$ Cov(B,C) = Cov(A,B)^T\cdot Cov(A)^{-1} \cdot Cov(A,C) $$

To draw a random sample from the system, you draw a random sample z_A, from A (easy since you have the whole covariance matrix for A) and then compute

$$ z_B = Cov(A)^{-1} \cdot Cov(A,B) \cdot z_A + randnorm(Cov(B) - Cov(A,B)) $$

and likewise for C.