Let x be a random variable in $\mathbb{R}^d$, $J$ a block tridiagonal $d\times d$ matrix, and probability of $x$ is defined as follows
$$p(x)\propto \exp(-x'Jx)$$
For a fixed $d\times d$ matrix $v$ compute $u$
$$u_{ij}=\sum_{kl}E[x_i x_j x_k x_l] v_{kl}$$
How can this be done efficiently and what is the algorithm to compute this?
It feels like Isserlis theorem should make this computation much faster than $O(d^4)$, but I'm having trouble finding a source that worked through the details.