2
$\begingroup$

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!

$\endgroup$
  • 1
    $\begingroup$ 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. $\endgroup$ – Biswajit Banerjee Sep 20 '20 at 7:59
  • $\begingroup$ Yes, if $d$ was not too large, I could ignore the structure and treat it as a regular linear problem $\endgroup$ – Yaroslav Bulatov Sep 20 '20 at 14:12
3
$\begingroup$

I encountered a similar problem in the past and I could find no simple solution either. One of the terms is a Kronecker product, another is a rank-1 modification, but the rest makes the problem more difficult. I don't think there is a closed-form solution; you could try using an iterative method, dropping some terms to get a preconditioner.

But if someone has a better solution, I am all ears!

$\endgroup$
  • $\begingroup$ I've played around with Jacobi algorithm to invert it, and it diverges when AB correlations are large $\endgroup$ – Yaroslav Bulatov Sep 22 '20 at 18:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.