Multiplying by certain elementary symmetric polynomials

Let me preface this with saying that it's probably not necessary to understand all the math behind this to review my code. Unless you have a lot of spare time or a very strong interest I also wouldn't try too hard to understand the background if you don't already know some of the basics, like what a symmetric polynomial is. In my edit below I provide a description of what the function is actually doing, which is a lot less lofty than it sounds.

This is math code, which as I'm sure you all know has a reputation for being ugly. The language is Python 3. I am implementing a special case of the formula in this paper. Specifically, it is to get a formula (in terms of indexing permutations) for multiplying a Schubert polynomial (which is a polynomial S_p indexed by a permutation p, there is one for each permutation, and any polynomial in any number of variables with integer coefficients can be expressed uniquely as a linear combination of Schubert polynomials with integer coefficients) by a polynomial of the form x_1*x_2*...*x_k. The input is a dictionary with key a tuple of integers (the indexing permutation) and value an int (the coefficient of the corresponding Schubert polynomial) and the result is a dictionary of the same kind.

single_variable calculates the coefficients of Schubert polynomials when multiplied by the monomial x_k. elem_sym_mul takes a similar permutation->coefficient dictionary and a list of values of k and for each k multiplies by x_1, x_2, ..., up to x_k, using single_variable.

Permutations are lists or tuples, which are unfortunately necessarily zero indexed, but the permutation is a permutation of all positive integers where all but finitely many elements are fixed, but if, for example, after index n everything is fixed, we omit everything after index n. For example (2,3,1,4) (when it needs to be hashed) or [2,3,1,4]. permtrim is just a utility function for trimming the permutation if the final elements are redundant. [2,3,1,4,5,6] is considered the same as [2,3,1,4] and permtrim implements this.

I'm looking for notable style issues or, if you see any, optimizations.

Edit: I realized it would probably be helpful to explain exactly what single_variable does. Let us assume the input is

{perm: coeff}


This represents the Schubert polynomial S_perm. The result when multiplying by x_k gives you some Schubert polynomials with a coefficient of coeff, and some Schubert polynomials with a coefficient of -coeff. The ones with a coefficient of coeff are obtained as follows: all permutations perm2 such that perm2 differs from perm by exchanging the element at index k-1 with an element at index j for some j>k-1 such that perm[k-1]<perm[j] and there does not exist any k-1<i<j such that perm[k-1]<perm[i]<perm[j]. In this case, the index j is allowed to go beyond the length of the list/tuple into fixed elements, but it will only ever go one element past the length of the list due to the considerations above, so before the processing one element is added to complete the search space.

The ones with coefficient -coeff are obtained similarly, but instead we are looking for indexes i<k-1 such that perm[i]<perm[k-1] and there does not exist i<j<k-1 such that perm[i]<perm[j]<perm[k-1].

Doing this for the whole dictionary of permutations to coefficients just iterates this for each entry and sums the result.

from typing import Dict, Tuple, List

def permtrim(perm:List[int]):
if len(perm)==1:
return perm
elif perm[len(perm)-1]==len(perm):
return permtrim(perm[:len(perm)-1])
else:
return perm

def single_variable(perm_dict:Dict[Tuple[int],int],k:int):
res_dict = {}
for perm,val in perm_dict.items():
perm2 = (*perm,len(perm)+1)
for i in range(k,len(perm2)):
if perm2[k-1]<perm2[i]:
good = True
for p in range(k,i):
if perm2[k-1]<perm2[p] and perm2[p]<perm2[i]:
good = False
break
if good:
permp = list(perm2)
permo = tuple(permtrim(permp[:k-1]+[permp[i]]+permp[k:i]+[permp[k-1]]+permp[i+1:]))
res_dict[permo] = res_dict.get(permo,0)+val
if res_dict.get(permo,0) == 0:
del res_dict[permo]
for i in range(0,k-1):
if perm2[i]<perm2[k-1]:
good = True
for p in range(i+1,k-1):
if perm2[i]<perm2[p] and perm2[p]<perm2[k-1]:
good = False
break
if good:
permp = list(perm2)
permo = tuple(permtrim(permp[:i]+[permp[k-1]]+permp[i+1:k-1]+[permp[i]]+permp[k:]))
res_dict[permo] = res_dict.get(permo,0)-val
if res_dict.get(permo,0) == 0:
del res_dict[permo]
return res_dict

def elem_sym_mul(perm_dict:Dict[Tuple[int],int],ks:List[int]):
dicto = perm_dict
for i in range(1,max(ks)+1):
for k in ks:
if i<=k:
dicto = single_variable(dicto,i)
return dicto


• Avoid recursion. It is expensive. permtrim can be (mechanically) converted into a clean iterative form:

def permtrim(perm):
while len(perm) > 1 and perm[-1] == len(perm):
perm.pop()
return perm

• Avoid naked loops. Every loop implements an algorithm, and deserves a name. For example, he loops for p in range(k, i) and for p in range(i+1, k-1) test that the every element in the range are between than some bounding values. Consider factoring it into a function, like range_is_between,

    if perm2[k-1] < perm2[i]:
if range_is_between (perm2[k:i], perm2[k-1], perm2[i]):
permp = ....
....


Now it should be obvious that the first if essentially is a part of the same logic, and better be delegated to the same function.

     if range_is_between (perm2[k:i], perm2[k-1], perm2[i]):
....


Flat is better than nested.

Also, notice how the good boolean flag disappears. Elimination of a boolean flags is a strong indication that you are going in a right direction.

• DRY. The bodies of two loops are suspiciously similar:

for i in range(0,k-1):
if range_is_between(perm2[i+1:k-1], perm2[i], perm2[k-1]):
permp = list(perm2)
permo = tuple(permtrim(permp[:i] + [permp[k-1]] + permp[i + 1:k-1] + [permp[i]] + permp[k:]))
res_dict[permo] = res_dict.get(permo,0)-val
if res_dict.get(permo,0) == 0:
del res_dict[permo]

for i in range(k,len(perm2)):
if range_is_between(perm2[k:i], perm2[k-1], perm2[i]):
permp = list(perm2)
permo = tuple(permtrim(permp[:k-1] + [permp[i]] + permp[k:i] + [permp[k-1]] + permp[i+1:]))
res_dict[permo] = res_dict.get(permo,0)+val
if res_dict.get(permo,0) == 0:
del res_dict[permo]


This is also an indication that they really want to be functions too:

def single_variable(perm_dict:Dict[Tuple[int],int],k:int):
res_dict = {}
for perm,val in perm_dict.items():
perm2 = (*perm,len(perm)+1)
res_dict = do_important_stuff(0, k, perm2)
res_dict = do_important_stuff(k, len(perm2), perm2)


I have no idea what is the right name for do_important_stuff. I have no domain knowledge. You do.

• I was skeptical about the recursion performance hit being enough to matter since it's generally only one or two layers deep, but I tried your version on a large example that previously took 29 seconds and it took 27 seconds, so definitely an improvement. I'll have to look at the rest later, thanks. – Matt Samuel May 10 '20 at 23:47
• @MattSamuel Don't get me wrong, the recursion point was not so much about performance (I am glad it helped the performance), but about clarity. Just like the other points. – vnp May 11 '20 at 2:31
• Unfortunately that speed gain is lost when I implement your second suggestion, it goes to 30.5 seconds. Apparently function calls in general are expensive. – Matt Samuel May 11 '20 at 10:41
• I ended up naming range_is_between as bruhat_ascent and leaving the two loops separate to save time, but making them look cleaner by naming the indices first_index and last_index. It also saved time to just do a straight swap of the two indices rather than all of that slicing, also clearer. Anyway, thanks. – Matt Samuel May 11 '20 at 11:26
• In the large example I speak of, it spends over 6 seconds in bruhat_ascent, because it is called 33363163 times. Most of the time the examples aren't that large, but if I want to test n=9 I'm probably going to want to use a different language. – Matt Samuel May 11 '20 at 21:10