# Generating summation formulas for factorized 4th moments

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]]

$$$$

• – JimB Oct 4 '19 at 0: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 4 '19 at 9: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 4 '19 at 16:25
• I figured out the most difficult piece after thinking about this problem a bit more, posted as separate question -- mathematica.stackexchange.com/questions/207336/… – Yaroslav Bulatov Oct 4 '19 at 17:16

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 4 '19 at 0:20