# Fourth moment of a random-variable with block-tridiagonal structure

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.

• Is $\exp$ missing in the definition of $p$? Sep 25, 2020 at 20:50

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.
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, 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)$$.
• @YaroslavBulatov : You wrote: "if the matrix is block diagonal, one could do this in O(d) and cubic in block size." How can you do that? I see how to do this in $O(d^2b^2)$ steps, where $b$ is the largest block size, for a block-diagonal $J$. Also, in the definition of $p$ do you have $x'Jx$ or $x'J^{-1}x$? If $x'Jx$, then I don't see how the (block-)triagonality of $J$ can help at all. Sep 27, 2020 at 17:45
• You are right, even reading the vector takes $O(d^2)$ steps. I was confused after looking at a setting of v matrix being rank-1, where $O(d)$ is possible by applying Isserlis on a per-block basis Sep 28, 2020 at 4:05