Optimise instantiation of free nilpotent Lie algebras

[pjt33]

A MathOverflow question asked about inefficiency in the instantiation of free nilpotent Lie algebras. Some basic profiling shows that the bottleneck is the calls to leading_support() when extracting the structural coefficients. A larger refactor to avoid precalculating the table might be worth considering, but simply caching the calls gives a notable improvement, about 50% for LieAlgebra(QQ, 3, 9).

📝 Checklist

• The title is concise and informative.
• The description explains in detail what this PR is about.
• I have linked a relevant issue or discussion.
• I have created tests covering the changes.
• I have updated the documentation and checked the documentation preview.

[tscrim]

Your analysis is good; moreover, about 68% of the time is being spent in line 414 (418 with your PR). However, a relatively simple further optimization that can be done is add to tuples in basis_by_deg with their leading support. Breaking encapsulation a bit, this becomes:

        basis_by_deg = {d: [(ind, elt, max(elt._monomial_coefficients)) for ind, elt in lst]
                        for d, lst in basis_by_deg.items() if d <= s}

        # extract structural coefficients from the free Lie algebra
        s_coeff = {}
        for dx in range(1, s + 1):
            # Brackets are only computed when deg(X) + deg(Y) <= s
            # We also require deg(Y) >= deg(X) by the ordering
            for dy in range(dx, s + 1 - dx):
                if dx == dy:
                    for i, val in enumerate(basis_by_deg[dx]):
                        X_ind, X, X_ls = val
                        for Y_ind, Y, Y_ls in basis_by_deg[dy][i + 1:]:
                            Z = L[X, Y]
                            if not Z.is_zero():
                                s_coeff[(X_ind, Y_ind)] = {W_ind: Z[W_ls]
                                                           for W_ind, W, W_ls in basis_by_deg[dx + dy]}
                else:
                    for X_ind, X, X_ls in basis_by_deg[dx]:
                        for Y_ind, Y, Y_ls in basis_by_deg[dy]:
                            Z = L[X, Y]
                            if not Z.is_zero():
                                s_coeff[(X_ind, Y_ind)] = {W_ind: Z[W_ls]
                                                           for W_ind, W, W_ls in basis_by_deg[dx + dy]}

[tscrim]

From using%lprun, it seems like with the change above, we start running into the sheer size of the computation:

   417   3156069 1047708549.0    332.0     38.1                                  s_coeff[(X_ind, Y_ind)] = {W_ind: Z[W_ls]
   418   1579341  269057779.0    170.4      9.8                                                             for W_ind, W, W_ls in basis_by_deg[dx + dy]}
   419                                           
   420         1     871012.0 871012.0      0.0          names, index_set = standardize_names_index_set(names, index_set)
   421         2  556301141.0    3e+08     20.2          s_coeff = LieAlgebraWithStructureCoefficients._standardize_s_coeff(
   422         1        133.0    133.0      0.0              s_coeff, index_set)
   423                                           
   424         4  671149556.0    2e+08     24.4          NilpotentLieAlgebra_dense.__init__(self, R, s_coeff, names,
   425         1        173.0    173.0      0.0                                             index_set, s,
   426         2        413.0    206.5      0.0                                             category=category, **kwds)

(the columns are: line #, # hits, time, time/hit, % time, code).