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!

Improve this question
$\endgroup$
2
  • 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
    Commented Sep 20, 2020 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
    Commented Sep 20, 2020 at 14:12

1 Answer 1

Reset to default
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!

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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