Generating summation formulas for factorized 4th moments

Question

I'm interested in getting summation formulas for the following expression, in Einstein summation notation, with indices ranging from $1$ to $d$ $$c=X_{ik}M_{ijkl}X_{jl}$$

Here $M_{ijkl}$ is specified as a product of terms, $X$ is user-provided $d$-by-$d$ matrix and the goal is to find summation formulas which take less than $d^4$ steps to compute.

Specific example, suppose $$M_{ijkl}=M_{ij}M_{kl}$$ then I can compute this sum faster than $d^4$ operations. With summation notation we have

$$A_{kj}=X_{ik} M_{ij}\\B_{kj} = M_{kl} X_{jl}\\c=A_{kj}B_{kj}$$

A more interesting factorization from Isserlis theorem:

$$M_{ijkl}=M_{ij}M_{kl}+M_{jk}M_{il}+M_{ik}M_{jl}$$

I have an expression to obtain $c$ faster than $d^4$ for this factorization, but I'd like to verify my algebra.

More generally, I want to generate summation formulas for to find $c$ for all relevant factorizations of $M_{ijkl}$, and reuse the ones that that need less than $d^4$ operations in another application. Any suggestions how to approach this in Mathematica?

For more background, my tensor $M_{ijkl}$ represents $E[X_i X_j X_k X_l]$ which I approximate from data by using factorizations obtained from MomentConvert as follows -- convert Moment[{1,1,1,1} to cumulants, set some cumulants to 0, convert back to moments, estimate these moments from data. There are 16 possible cumulants to set to 0, which means $2^{16}$ formulas, but I feel like vast majority of formulas will require $O(d^4)$ steps to compute $c$, so only a few formulas are interesting.

Here's an example of obtaining $il\cdot jk + ik \cdot jl + ij \cdot kl$ factorization from Isserlis theorem example:

(* Prints Isserlis factorization:  *)


(* Convert Moments term to Cumulant term and visa versa *)

convert[a_Moment] := MomentConvert[a, "Cumulant"];
convert[a_Cumulant] := MomentConvert[a, "Moment"];

(* Convert all moment (or cumulant) terms in the expression *)

convertAll[expr_] := (
  MapAt[convert, expr, termPositions[expr]]
  )

(* Get position of every Moment/Cumulant term in the expression *)

termPositions[expr_] := (
   poses0 = Most /@ Position[expr, Moment];
   poses1 = Most /@ Position[expr, Cumulant];
   poses0~Join~poses1
   );

(* drop Cumulant/Moment terms matching given criterion *)

dropTerms[expr_, crit_] := Module[{},
   helper[part_] := 
    If[part[[0]] === Cumulant || part[[0]] === Moment, 
     If[crit[part], 0, part], part];
   MapAt[helper, expr, termPositions[expr]]
   ];

rank[expr_] := Total[expr[[1]]]

moments = Moment[{1, 1, 1, 1}];
cumulants = convertAll@moments;
(* drop rank-4 cumulants *)
cumulants = dropTerms[cumulants, rank[#] > 3 &];
moments = convertAll@cumulants;
(* drop rank-1 cumulants *)
moments = dropTerms[moments, rank[#] == 1 &]

letters[expr_] := 
 StringJoin@Pick[{"i", "j", "k", "l"}, Thread[# == 1]] &[expr[[1]]]
MapAt[letters, moments, termPositions[moments]]


```

Have you seen https://mathematica.stackexchange.com/questions/14654/how-to-define-the-einstein-summation-convention-in-mathematica, https://www.researchgate.net/post/Is_there_any_simple_Mathematica_package_handling_implicit_Einstein_summation, and https://sites.math.washington.edu/~lee/Ricci/ ? — JimB, Oct 04 '19 at 00:42
So what is the desired input and the desired output? How do you want to encode the result? Your code seems to address a different issue. — yarchik, Oct 04 '19 at 09:20
You could encode the einsum computation by using specifying {term1,term2}->{indices-to-sum}. The result can be chained up to more computations. For the example above, you could specify it as {{X[ik],M[ij]}->{i},{M[kl],X[jl]}->{l}}->{kj} — Yaroslav Bulatov, Oct 04 '19 at 16:25
I figured out the most difficult piece after thinking about this problem a bit more, posted as separate question -- https://mathematica.stackexchange.com/questions/207336/generating-an-efficient-way-to-compute-einsum — Yaroslav Bulatov, Oct 04 '19 at 17:16

score 0 · Accepted Answer · answered Oct 04 '19 at 00:10

I must say that I'm not sure what your question is. However, the following code will reproduce your code much quicker (about 90 times faster).

moments = Select[Tuples[{0, 1}, 4], Total[#] == 2 &];
n = Length[moments];
moments = Table[Moment[moments[[i]]] Moment[moments[[n - i + 1]]], {i, n/2}];


letters[expr_] := StringJoin@Pick[{"i", "j", "k", "l"}, Thread[# == 1]] &[expr[[1]]];
MapAt[letters, moments, termPositions[moments]] // Total
(* "il" "jk" + "ik" "jl" + "ij" "kl" *)

So to do what you want in general maybe using Tuples will help speed things up.

Thanks for the tip. Will edit the question, but ultimately I'm asking for an algorithm to break einstein summations like $X_{ik} M_{ij} M_{kl} X_{jl}$ into cheaper summations like the one given in example — Yaroslav Bulatov, Oct 04 '19 at 00:20

Generating summation formulas for factorized 4th moments

1 Answers1