I have a system of linear equations (here 9) as below:
$ 2 a_{\{1,2,1,2,2,1,2,1,2,2\}}+2 a_{\{1,1\}}=1, \\ \left(a_{\{1,2,1,2,2,1,2,1,2,2,2\}}+a_{\{1,2,1,2,2,1,2,2,1,2,2\}}\right)+2 a_{\{1,1,2\}}=1, \\ 2 a_{\{1,1,1,2,1,2,2,1,2,1,2,2\}}+2 a_{\{1,1,1,1\}}=2 a_{\{1,1\}}, \\ \left(a_{\{1,2,1,2,2,1,2,1,2,2,2,2\}}+a_{\{1,2,1,2,2,1,2,2,2,1,2,2\}}\right)+2 a_{\{1,1,2,2\}}=1, \\ 2 a_{\{1,2,1,2,2,1,2,2,1,2,2,2\}}+2 a_{\{1,2,1,2\}}=1, \\ \left(a_{\{1,1,1,2,1,2,2,1,2,1,2,2,2\}}+a_{\{1,1,1,2,2,1,2,1,2,2,1,2,2\}}\right)+2 a_{\{1,1,1,1,2\}}=a_{\{1,1\}}+a_{\{1,1,2\}}, \left(a_{\{1,1,2,1,2,1,2,2,1,2,1,2,2\}}+a_{\{1,1,2,1,2,2,1,2,1,2,2,1,2\}}\right)+2 a_{\{1,1,1,1,2\}}=2 a_{\{1,1,2\}}, \left(a_{\{1,2,1,2,2,1,2,1,2,2,2,2,2\}}+a_{\{1,2,1,2,2,1,2,2,2,2,1,2,2\}}\right)+2 a_{\{1,1,2,2,2\}}=1, \\\left(a_{\{1,2,1,2,2,1,2,2,1,2,2,2,2\}}+a_{\{1,2,1,2,2,1,2,2,2,1,2,2,2\}}\right)+2 a_{\{1,2,1,2,2\}}=1 $
and I get solution in Mathematica as below: $\left\{a_{\{1,2,1,2,2,1,2,1,2,2\}}\to \frac{1}{2}-a_{\{1,1\}},a_{\{1,2,1,2,2,1,2,2,1,2,2\}}\to -2 a_{\{1,1,2\}}-a_{\{1,2,1,2,2,1,2,1,2,2,2\}}+1,a_{\{1,1,1,2,1,2,2,1,2,1,2,2\}}\to a_{\{1,1\}}-a_{\{1,1,1,1\}},a_{\{1,2,1,2,2,1,2,2,1,2,2,2\}}\to \frac{1}{2}-a_{\{1,2,1,2\}},a_{\{1,2,1,2,2,1,2,2,2,1,2,2\}}\to -2 a_{\{1,1,2,2\}}-a_{\{1,2,1,2,2,1,2,1,2,2,2,2\}}+1,a_{\{1,1,1,2,2,1,2,1,2,2,1,2,2\}}\to a_{\{1,1\}}+a_{\{1,1,2\}}-2 a_{\{1,1,1,1,2\}}-a_{\{1,1,1,2,1,2,2,1,2,1,2,2,2\}},a_{\{1,1,2,1,2,2,1,2,1,2,2,1,2\}}\to 2 a_{\{1,1,2\}}-2 a_{\{1,1,1,1,2\}}-a_{\{1,1,2,1,2,1,2,2,1,2,1,2,2\}},a_{\{1,2,1,2,2,1,2,2,2,1,2,2,2\}}\to -2 a_{\{1,2,1,2,2\}}-a_{\{1,2,1,2,2,1,2,2,1,2,2,2,2\}}+1,a_{\{1,2,1,2,2,1,2,2,2,2,1,2,2\}}\to -2 a_{\{1,1,2,2,2\}}-a_{\{1,2,1,2,2,1,2,1,2,2,2,2,2\}}+1\right\}$
This takes about 0.08 seconds. I reproduced the same solution in Python using SymPy Solve. It takes about 0.85 seconds.
Why is Mathematica faster than Python and is there a way to improve the timing further?
(** Equivalent forms: trace invariant terms + transpose invariance **)
equivalentForms[nl_]:=If[nl=={},{{}},Join[NestList[RotateLeft,nl,Length[nl]-1],NestList[RotateLeft,Reverse[nl],Length[nl]-1]]];
cForm[nl_]:=First@Sort@equivalentForms[nl];deleteDuplicate[list_]:=DeleteDuplicates[Map[cForm,list]];
bracelets[k_Integer]:=deleteDuplicate@Tuples[{1,2},k];
formal=Subscript[a, #]/.{Subscript[a, {}]->1}&;
process[nl_]:=If[EvenQ[Count[nl,1]],
formal@cForm@Flatten@nl,0];
fh1[1]:={2,1,2,2,1,2,1,2,2};
fh2[1]:={2,2,1,2,1,2,2,1,2};
(** This is the LHS (Left-Hand Side)**)
loopyInteraction[nl_,pos_Integer]:=
If[nl[[pos]]==1,
2*process[nl]+g1( process[Flatten[MapAt[fh1,nl,pos]]]+process[Flatten[MapAt[fh2,nl,pos]]]),
0];
(** This is the RHS (Right-Hand Side)**)
loopyQuad[nl_,pos_Integer]:=
Module[{d=Flatten@DeleteCases[Position[nl,nl[[pos]]],{pos}],td,doubleTr},
If[Length[d] ==0,
0,
If[nl[[pos]]==1,
td=Sort@{pos,#}&/@d;
doubleTr=Map[process,MapAt[Delete[{{1},{-1}}],#,1]&/@(TakeDrop[nl,#]&/@td),{2}];
Total[Times@@@doubleTr]/.{Subscript[a, {}]->1},
0
]
]
];
constr[nl_,pos_Integer]:=loopyInteraction[nl,pos]==loopyQuad[nl,pos];
loop[k_Integer]:=DeleteDuplicates@Flatten@Outer[constr,bracelets[k],Range@k,1];
br=bracelets;
loopAll[k_Integer]:=Flatten[loop/@Range[k]];
brAll[k_Integer]:=formal/@Flatten[br/@Range[k],1];
(** Max degree and run **)
kmax=5;
MatrixForm[DeleteCases[loopAll[kmax](*/.g1\[Rule]1*)/.Subscript[a, {2}]->1/.Subscript[a, {2,2}]->1/.Subscript[a, {2,2,2}]->1/.Subscript[a, {2,2,2,2}]->1/.Subscript[a, {2,2,2,2,2}]->1/.Subscript[a, {2,2,2,2,2,2}]->1,True],TableAlignments->Left];
loopsNumeric[k_Integer,g_]:=loopAll[k]/.g1->g/.Subscript[a, {2}]->1/.Subscript[a, {2,2}]->1/.Subscript[a, {2,2,2}]->1/.Subscript[a, {2,2,2,2}]->1/.Subscript[a, {2,2,2,2,2}]->1/.Subscript[a, {2,2,2,2,2,2}]->1/.Subscript[a, {2,2,2,2,2,2,2}]->1/.Subscript[a, {2,2,2,2,2,2,2,2}]->1/.Subscript[a, {2,2,2,2,2,2,2,2,2}]->1/.Subscript[a, {2,2,2,2,2,2,2,2,2,2}]->1;
sol10=Flatten@Solve[loopsNumeric[5,1]] (* Timing :~ 0.08 sec *)
