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.



Fourth moment of a random-variable with block-tridiagonal structure
Using e.g. the Gauss elimination, we can diagonalize the matrix $(v\_{kl})$, that is, write
$$v\_{kl}=\sum\_{r=1}^d a\_r s\_{rk}t\_{rl}$$
for some real $a\_r,s\_{rk},t\_{rl}$ and all $k,l$; the computational complexity (CC) of this diagonalization is $O(d^3)$; cf. e.g. [this source](https://nholmber.github.io/2018/05/mpi-diagonalization/).


Now we can write
$$u\_{ij}=\sum\_r a\_r\,Ex\_ix\_jX\_rY\_r,$$
where
$$X\_r:=\sum\_k s\_{rk}x\_k,\quad Y\_r:=\sum\_k t\_{rl}x\_l.$$
For each $(i,r)$, the CC of $Ex\_iX\_r$ and $Ex\_iY\_r$, as well as of $EX\_rY\_r$, is $O(d)$, which adds up to $O(d^3)$ for all $(i,r)$.
After that, because $x\_i,x\_j,X\_r,Y\_r$ are zero-mean jointly normal -- see e.g. the [Isserlis formula](https://en.wikipedia.org/wiki/Isserlis%27_theorem), for each $(i,j,r)$ the CC of $Ex\_ix\_jX\_rY\_r$ is $O(1)$ and hence the CC of $u\_{ij}=\sum\_r a\_r\,Ex\_ix\_jX\_rY\_r$ is $O(d)$ for each $(i,j)$.


So, the overall CC is $O(d^3)$.


