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:


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

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


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.

  • $\begingroup$ 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 $\endgroup$ – Yaroslav Bulatov Oct 4 '19 at 0:20

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.