0
$\begingroup$

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.

Improve this question
$\endgroup$
1
  • 1
    $\begingroup$ Is $\exp$ missing in the definition of $p$? $\endgroup$
    – Iosif Pinelis
    Sep 25, 2020 at 20:50

1 Answer 1

Reset to default
1
$\begingroup$

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)$.

Improve this answer
$\endgroup$
3
  • $\begingroup$ It seems this method destroys block structure. IE, if the matrix is block diagonal, one could do this in O(d) and cubic in block size. I suspect block tridiagonal is cheaper than general case $\endgroup$
    – Yaroslav Bulatov
    Sep 27, 2020 at 17:01
  • $\begingroup$ @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. $\endgroup$
    – Iosif Pinelis
    Sep 27, 2020 at 17:45
  • $\begingroup$ 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 $\endgroup$
    – Yaroslav Bulatov
    Sep 28, 2020 at 4:05

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .