Solution for system of linear equations for unknown variables in Mathematica and Python

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}}+1right}$

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 *)