Pensieve header: A full implementation of $gl_n^\epsilon$.

SetDirectory["C:/drorbn/AcademicPensieve/Projects/glneps"];
Once[
<<"KnotTheory`";
<< "Rot.m";
<< "PolyPlot.m";
]

Canonical Forms:

CF[ℰ_] := Expand@Collect[ℰ,Cases[ℰ,(Subscript[(y|a|x|η|β|τ|ξ), __])[_],∞]⋃Cases[ℰ,(y|a|x|η|β|τ|ξ)[_],∞],F]/.F->Factor;
CF[ℰ:𝔼[_][___]]:=CF/@ℰ;

“Define” Code

SetAttributes[Define, HoldAll];
Define[def_,defs__] := (Define[def]; Define[defs];);
Define[op_[is__]=ℰ_] := Module[{SD,ii,jj,kk,isp,nis,nisp,sis},Block[{i,j,k},
ReleaseHold[Hold[
SD[op[nisp,$k_Integer],Block[{i,j,k},op[isp,$k]=ℰ;op[nis,$k]]];
SD[op[isp],op[{is},$k]];
SD[op[sis__],op[{sis}]];
]/. {SD->SetDelayed,
isp->{is}/.{i->i_,j->j_,k->k_},
nis->{is}/.{i->ii,j->jj,k->kk},
nisp->{is}/.{i->ii_,j->jj_,k->kk_}
}]  ]]

(Subscript[utm, n_])[i_,j_->k_]:=(*(Subscript[utm, n])[i,j->k]=*)𝔼[{i,j}->{k}][
B=IdentityMatrix[n];
Do[B=B.MatrixExp[SparseArray[{α,β}->(Subscript[ξ, α,β]),{n,n}]],{β,1,n},{α,1,β}]; (* This specifies the PBW ordering *)
eqns=And@@Thread[Flatten[CF[(Inverse[B/.(Subscript[ξ, αβ__]):>(Subscript[ξ, αβ])[k]]/.(Subscript[ξ, α_,α_])[k]:>(Subscript[ξ, α,α])[i]+(Subscript[ξ, α,α])[j]).(B/.(Subscript[ξ, αβ__]):>(Subscript[ξ, αβ])[i]).(B/.(Subscript[ξ, αβ__]):>(Subscript[ξ, αβ])[j])]]==Flatten[IdentityMatrix[n]]];vars=Union@Cases[eqns,(Subscript[ξ, __])[k],∞];{sol}=Solve[eqns,vars];
CF[(vars/.ξ->x).(vars/.sol)+Sum[(Subscript[x, α,α])[k]((Subscript[ξ, α,α])[i]+(Subscript[ξ, α,α])[j]),{α,n}]]
];

(Subscript[utm, 2])[i,j->k]

Table[First@Timing[n->(Subscript[utm, n])[1,2->3]],{n,10}]