# Solving MX=N where M is structured as a Gaussian 4th-moment tensor

I'm looking to solve numerically the following equation for $$(d,d)$$ variable $$X$$, in Einstein summation notation

$$M_{ijkl}X_{kl}=N_{ij}$$

Where $$M$$ is a $$(d,d,d,d)$$ 4th-moment tensor of random variable $$(a,b)$$ and admits representation in terms of 2nd and 1st moment tensors based on the Isserlis theorem

$$M_{ijkl}=E[a_i a_j b_k b_k]=AA_{ij}BB_{kl}+AB_{ik}AB_{jl}+AB_{il}AB_{jk}-2A_iA_jB_kB_l$$

relying on shorthand $$AA=E[aa'],BB=E[bb'],AB=E[ab'],A=E[a],B=E[b]$$

If my factorization only had the first term, an explicit solution is formed by inverting $$AA$$ and $$BB$$ matrices

$$X=(BB)^{-1}N(AA)^{-1}$$

However, the presence of $$AB$$ terms is complicating things, probably an iterative method is required. Is there a standard procedure for this kind of structured linear equation system? Any literature pointers appreciated!

• If $d$ is not too large, you can represent $M$ by a $d^2 \times d^2$ matrix and arrange $X$ and $N$ into $d^2 \times 1$ matrices. Some care is needed to make sure all terms have been accounted for. Sep 20 '20 at 7:59
• Yes, if $d$ was not too large, I could ignore the structure and treat it as a regular linear problem Sep 20 '20 at 14:12