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.