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

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?

## 1 Answer

(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.

• Thanks, that's indeed quite relevant! This formula means that my covariance should be written as a sum of two Kronecker products, so I may need to look at kronecker product tricks that apply to inverse of sum of two matrices (it seems one formula is math.stackexchange.com/questions/17776/…) Sep 25 '19 at 18:41
• Thanks for the good starting point, it seems after applying Isserlis' and some Kronecker rearranging, dividing by this covariance matrix reduces to solving Sylvester matrix equation on small matrices: wolframcloud.com/obj/yaroslavvb/newton/isserlis-inversion.nb Sep 26 '19 at 6:05
• Nice, I had missed that Sylvester rearranging. Thanks for sharing this solution, this was an interesting read. Sep 26 '19 at 6:42