Reconstructing statistics of $x\otimes y$ from E[XX'], E[YY'] and E[XY']

Question

I'm looking at random vectors $z$ of size $d^2$ which can be written as $z=x\otimes y$ where $x$,$y$ are random vectors in $\mathbb{R}^d$ with following second moments known -- $E[XX']$, $E[YY']$ and $E[XY']$. First moments are 0.

Given a new vector $\hat{z}=\hat{x}\otimes \hat{y}$, I need to approximately divide it by covariance matrix of $z$, which needs to be estimated from 3 moments matrices above. If $E[XY']$ is zero, then division is straightforward

$$\text{cov}(z)^{-1}\hat{z} = E[XX']^{-1} \hat{x} \otimes E[YY']^{-1} \hat{y}$$

However, in my case $E[XY']$ is not 0, and this estimate will incur an error. Is there some systematic way to incorporate $E[XY']$ to correct this?

Answer 1

(Converted to an answer from my comments and expanded.)

Basically you need to compute the fourth moments $E[x_i x_j y_k y_l]$ for all $i,j,k,l$, given the second moments. These fourth moments are not universally determined, but they depend on the distribution of your variables, even in the one-dimensional case (that's exactly the reason why kurtosis is a thing). (In particular, if $x$ and $y$ are uncorrelated but not independent, I strongly suspect that even your formula for the case $E[xy']=0$ is wrong since it seems that you derived it using independence.)

Since you are asking for an approximation, though, not everything is lost. You may compute the value of that expression under the assumption that $x$ and $y$ are Gaussian; even if they are not, hopefully it serves as a reasonable approximation of the true covariance if their distribution is not too extreme.

For Gaussian variables, the result is given by Isserlis' theorem: $$ cov(x\otimes y)_{di+k,dj+l} = E[x_ix_jy_ky_l] = E[x_ix_j]E[y_ky_l] + E[x_iy_k]E[x_jy_l] + E[x_iy_l]E[x_jy_k]. $$ If you set $cov(x)=V,cov(y)=W$, then the first summand gives the matrix $V\otimes W$ that you have already found; the second gives the rank-1 matrix $vec(V)vec(W)'$, and I don't see immediately a nice expression for the third. You can compute all the matrix entries easily, but I don't see immediately any linear algebra tricks to perform the inversion by reducing it to the inversion of $n\times n$ matrices, as it happened instead in your special case.