move invariant symmetric bilinear form to instance method

[jacksonwalters]

The computation of G-invariant symmetric bilinear forms only depends on the group and representation. Since there doesn't seem to be a generic class for finite group representations, for now we leave this in SymmetricGroupRepresentation_generic_class. Now it can be tested, optimized, and possibly moved in the future.

📝 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.

⌛ Dependencies

#39200

[github-actions]

Documentation preview for this PR (built with commit 7b9cab2; changes) is ready! 🎉
This preview will update shortly after each push to this PR.

[jacksonwalters]

There is probably a better way to get the representation matrix than instantiating a Specht module within SymmetricGroupRepresentation_generic_class, given that representation_matrix is already defined there.

The issue is that if you try to use representation_matrix, then for the unitary case it tries to create a unitary representation using invariant_symmetric_bilinear_form(), creating an infinite loop.

Perhaps checking implementation before computing self._specht, and if it's "unitary" compute the Specht module, otherwise just use representation_matrix.

[DaveWitteMorris]

I think GAP can find invariant bilinear forms. See MTX.InvariantQuadraticForm of The MeatAxe. Sagemath shouldn't reinvent the wheel unless there is a good reason.

Other comments:

• The docstring speaks of the invariant form, but invariant forms are not unique (up to a scalar multiple) unless the representation is absolutely irreducible.
• I haven't taken the time to fully understand your code, but I especially don't understand why you are summing over all group elements. That seems extremely inefficient, because it suffices to find a bilinear form that is invariant under the generators, so a small system of linear equations will suffice.

[jacksonwalters]

I think GAP can find invariant bilinear forms. See MTX.InvariantQuadraticForm of The MeatAxe. Sagemath shouldn't reinvent the wheel unless there is a good reason.

It appears so. I am aware of the MeatAxe algorithm, and I believe this came up in discussion during #39200, #38455, and #38456.

For the extended Cholesky decomposition, we did use the GAP forms package initially via libgap, then decided to rewrite the code in Sage. This gave us a bit more flexibility, and I believe the issue was the forms wasn't a default installed package, so that caused some issues running tests, etc.

If it's possible and better to use GAP, without the issues we ran into before, I am open to that.

Other comments:

• The docstring speaks of the invariant form, but invariant forms are not unique (up to a scalar multiple) unless the representation is absolutely irreducible.

True. I thought the rest of the docs made it clear the solution is a vector space of some dimension. I've changed "the" to "an".

• I haven't taken the time to fully understand your code, but I especially don't understand why you are summing over all group elements. That seems extremely inefficient, because it suffices to find a bilinear form that is invariant under the generators, so a small system of linear equations will suffice.

Great point, thank you. I'll make that change and ensure tests are passing.

[jacksonwalters]

@DaveWitteMorris Replaced G in the sum with G.gens(). Sage gives the n-1 transpositions as generators, and this seems to work. I have updated the code, and it should be faster.

I tried using a minimal set of generators (an n cyclic and a transposition), but it didn't work. It gave a bunch of elementary matrices for one of the computations, which aren't symmetric. Perhaps this has something to do with $S_n$ being a Coxeter group. It shouldn't matter, since if the generators leave the form invariant, then any word in them will.

[DaveWitteMorris]

I don't understand your code (including why you are using a sum). The way I look at it, the equation \rho(g)^T U \rho(g) = U is a system of n^2 linear equations in n^2 variables whose solutions are the forms that are g-invariant. (Actually, since you are only interested in symmetric forms, you could cut the number of equations and variables essentially in half.) This provides a system of equations for each generator of G (and all of the systems of equations use the same variables). I would take the union of these systems of equations, to get one large system of equations, because the solutions to this union are precisely the G-invariant forms. This should work for any generating set, including the 2-element one.

However, it is a basic principle of the sagemath project that we do not want to reinvent the wheel. (One reason for this is to avoid adding to the maintenance burden.) So I think you should have a good reason to code this, instead of just calling GAP.

[jacksonwalters]

I don't understand your code (including why you are using a sum). The way I look at it, the equation \rho(g)^T U \rho(g) = U is a system of n^2 linear equations in n^2 variables whose solutions are the forms that are g-invariant. (Actually, since you are only interested in symmetric forms, you could cut the number of equations and variables essentially in half.) This provides a system of equations for each generator of G (and all of the systems of equations use the same variables). I would take the union of these systems of equations, to get one large system of equations, because the solutions to this union are precisely the G-invariant forms. This should work for any generating set, including the 2-element one.

However, it is a basic principle of the sagemath project that we do not want to reinvent the wheel. (One reason for this is to avoid adding to the maintenance burden.) So I think you should have a good reason to code this, instead of just calling GAP.

@DaveWitteMorris

1. As for the sum, this is list comprehension in Python. The sum is not actually a "sum", but rather used to concatenate lists since [a]+[b] = [a,b] and sum([[g] for g in G],[]) == list(G). Indeed, for each g in G we get a system of equations resulting in a d_rho $\times$ d_rho matrix (let's avoid using n because it is the degree of the symmetric group here) with unknowns $u_{ij}$ as formal variables in a polynomial ring R = GF(q**2)[u_{ij}]. This follows from the requirement $\langle gx, gy \rangle = \langle x,y \rangle$, which means $(gx)^\ast(gy) = x^T \rho(g)^\ast U \rho(g)y == x^T U y$. This holds for each basis vector $e_i$, so the equation $\rho(g)^T U \rho(g) == U$ should be true. Since the rep'ns are valued in the base field $GF(q)$, we can just write transpose rather than conjugate-transpose since conjugation leaves $GF(q)$ fixed. The "sum" operation should probably be called "stack" or something like that (I believe it used to be), and indeed it just gives a larger system of equations, one for each g in G. Note that we start with the empty list [], and then stack on top of that. I wouldn't call it a union (that's like "or", whereas this is "and" for each condition); it's closer to an intersection. That's true that we could reduce the number of variables since $u_{ij} = u_{ji}$. That would be an extra one or two lines, but it means we are assuming the symmetry condition (and presumably not be able to find alternating forms). I agree that it should work for any system of generators, including the 2-element system. That is a bug that could be addressed in a separate issue if this gets merged.

2. As for refactoring this (working) code to use GAP, in this case we'd be re-reinventing the wheel, and there is another PR that depends on this code and it has already been merged. Do you mind writing a minimal example using GAP which calls the libgap interface that computes these computes the forms in the doctests and translates them back into Sage?

[DaveWitteMorris]

That is a bug that could be addressed in a separate issue if this gets merged.

I think the bug needs to be fixed before the PR is merged, not after. I will post more comments later today when I have time.

[jacksonwalters]

That is a bug that could be addressed in a separate issue if this gets merged.

I think the bug needs to be fixed before the PR is merged, not after. I will post more comments later today when I have time.

The point of this PR is to move this computation to a different part of the code, and abstract and streamline the code in the UnitaryRepresentation class. I think optimizing this piece of code (which has already received a big speed improvement by switching from G to G.gens()) is a separate issue. In any case, it shouldn't be too hard; I can tackle it.

[DaveWitteMorris]

I think it's a good idea to have a method to find invariant bilinear forms for general representations, but the developer guide makes it clear that contributions to sagemath need to be high quality. I get the impression that you want to merge a half-baked PR and fix the code later, but that's not the way this project does things.

a. When you wrote "Now it can be tested, optimized, and possibly moved in the future", I thought you meant while it is in draft status. Testing and moving (and obvious optimizations) are to be done before the code is merged. In particular, sagemath's deprecation policy makes it cumbersone to move (or rename) methods after they are part of the codebase.

The algorithm works for (finite or infinite) semigroups, not just groups, so the method may belong somewhere near the top of Representation_abstract.

b. The documentation of the method is wrong, so needs to be fixed. For example:

The first line of the docstring says the method will return a symmetric form, but there is no guarantee that the form is symmetric.

The method returns a basis of the space of invariant forms, not "a form".

Is it really true that the space of invariant forms will generically be one-dimensional?

Is it really true that "the forms decompose into a sum of symmetric and alternating"? It's not obvious to me in characteristic two.

I also notice that a comment says "#solve the linear system for U for all group elements", even though it's only for the generators, not all group elements.

c. The name invariant_symmetric_bilinear_form of the method is not accurate, so it should be changed. Some classes of groups have an invariant_form method, so that might work. However, this method doesn't return a single form, so invariant_forms might be better. (Some thought should be put into this.)

d. The method sometimes returns a form and sometimes returns a list. It is much better to have a consistent return type. I think it would make sense to return a basis, so the return value would be an empty list when there are no invariant forms (instead of raising an exception).

e. When I suggested getting rid of half of the variables, you responded that it would ignore the alternating forms. Well, since the name of the method is invariant_symmetric_bilinear_form and the documentation says it will return a symmetric form, I think it was reasonable of me to think you were only interested in symmetric forms.

However, now that you bring it up, I certainly agree that the method should find both types. There should be a keyword argument to let the user choose which type they want (and something like None would mean to return all invariant forms). Even in the case of None, I think it would be faster to run the restricted algorithm twice (once for each type) than to use the full set of variables.

f. I don't agree that implementing a libgap interface would be "reinventing the wheel". (It is eliminating code duplication, which is something quite different.) Although I think libgap should have been used from the start, you are right that the calculus has changed now that there is already code in sagemath. (However, I don't think it is accurate to call it "working" code, since there is a bug.) It would be reasonable to merge a corrected version of your code into sagemath, but then I think an issue should be opened to replace it with a libgap wrapper.

g. I do not have time to write a libgap wrapper right now. When I saw the issue, I just thought I should make a few quick comments because I could see that you were headed in some wrong directions. I'm glad you have accepted some of my suggestions, but I have already spent a lot more time on this than I intended.

h. I think additional comments would make the code easier to understand. When I first looked at it, I thought the sum was probably concatenating lists of equations. But then you said that the code didn't work for a certain generating set, but did work for a different generating set. Since the particular generating set does not make any difference to the algorithm I had in mind, I thought you must be doing something else. Instead, we seem to agree that it is a bug in the code.

i. The code needs to work for all (semi)groups, so there will not be any control at all over the generating set. I am glad you agreed to try to fix this. The PR can't get a positive review until the code is working properly.

[jacksonwalters]

I think it's a good idea to have a method to find invariant bilinear forms for general representations, but the developer guide makes it clear that contributions to sagemath need to be high quality. I get the impression that you want to merge a half-baked PR and fix the code later, but that's not the way this project does things.

a. When you wrote "Now it can be tested, optimized, and possibly moved in the future", I thought you meant while it is in draft status. Testing and moving (and obvious optimizations) are to be done before the code is merged. In particular, sagemath's deprecation policy makes it cumbersone to move (or rename) methods after they are part of the codebase.

The algorithm works for (finite or infinite) semigroups, not just groups, so the method may belong somewhere near the top of Representation_abstract.

b. The documentation of the method is wrong, so needs to be fixed. For example:

The first line of the docstring says the method will return a symmetric form, but there is no guarantee that the form is symmetric.

The method returns a basis of the space of invariant forms, not "a form".

Is it really true that the space of invariant forms will generically be one-dimensional?

Is it really true that "the forms decompose into a sum of symmetric and alternating"? It's not obvious to me in characteristic two.

I also notice that a comment says "#solve the linear system for U for all group elements", even though it's only for the generators, not all group elements.

c. The name invariant_symmetric_bilinear_form of the method is not accurate, so it should be changed. Some classes of groups have an invariant_form method, so that might work. However, this method doesn't return a single form, so invariant_forms might be better. (Some thought should be put into this.)

d. The method sometimes returns a form and sometimes returns a list. It is much better to have a consistent return type. I think it would make sense to return a basis, so the return value would be an empty list when there are no invariant forms (instead of raising an exception).

e. When I suggested getting rid of half of the variables, you responded that it would ignore the alternating forms. Well, since the name of the method is invariant_symmetric_bilinear_form and the documentation says it will return a symmetric form, I think it was reasonable of me to think you were only interested in symmetric forms.

However, now that you bring it up, I certainly agree that the method should find both types. There should be a keyword argument to let the user choose which type they want (and something like None would mean to return all invariant forms). Even in the case of None, I think it would be faster to run the restricted algorithm twice (once for each type) than to use the full set of variables.

f. I don't agree that implementing a libgap interface would be "reinventing the wheel". (It is eliminating code duplication, which is something quite different.) Although I think libgap should have been used from the start, you are right that the calculus has changed now that there is already code in sagemath. (However, I don't think it is accurate to call it "working" code, since there is a bug.) It would be reasonable to merge a corrected version of your code into sagemath, but then I think an issue should be opened to replace it with a libgap wrapper.

g. I do not have time to write a libgap wrapper right now. When I saw the issue, I just thought I should make a few quick comments because I could see that you were headed in some wrong directions. I'm glad you have accepted some of my suggestions, but I have already spent a lot more time on this than I intended.

h. I think additional comments would make the code easier to understand. When I first looked at it, I thought the sum was probably concatenating lists of equations. But then you said that the code didn't work for a certain generating set, but did work for a different generating set. Since the particular generating set does not make any difference to the algorithm I had in mind, I thought you must be doing something else. Instead, we seem to agree that it is a bug in the code.

i. The code needs to work for all (semi)groups, so there will not be any control at all over the generating set. I am glad you agreed to try to fix this. The PR can't get a positive review until the code is working properly.

@DaveWitteMorris Before going on, I need to clarify something: do you understand that this code, aside from being moved to a class function, has already been merged in #38455?

You're saying some conflicting things:

1. this code is "half-baked" and not up to standards which conflicts with
2. the code has been reviewed, and discussed in detail with @tscrim and @dimpase and merged in

Further, you have stated that

1. you don't undestand the code
2. you would like a lot of modifications done which you yourself are not willing to do or open a separate PR for

Apologies, but this PR is simply meant to reorganize this code, not entirely refactor it.

[tscrim]

By moving this code, @jacksonwalters you are making it more general-purpose, so it should work in those cases. I also agree that the return time should be consistent as the user ahead of time will not know the dimension of the invariant forms ahead of time. I agree with most of @DaveWitteMorris's suggestions, although we should be careful to abide by the "the perfect is the enemy of the good" mantra. Let's try to spend a bit of time to optimize this now and fix any bugs we do come across.

Solving $\binom{n}{2}$ equations is much faster than solving $n^2$ equations. So splitting it into alternating and symmetric form computations would be good, with the general method returning the combination of the two (as every matrix splits into a symmetric and skew-symmetric parts).

Note also you do not need to create another representation. You already have $\rho(g)$ from the class you are working with.

[DaveWitteMorris]

@jacksonwalters: In answer to your question: no, I do not agree that this code has already been reviewed and merged. Most of my comments concern problems that arose in converting it to a class function. I do make mistakes when reviewing PR's, but I think all of my comments this time are valid, except perhaps the suggestion to use libgap, which you have already rejected, and the two that I phrased as questions (but now I'm pretty sure that my concern about characteristic two is correct). (However, I think something like "both" would be better than None for requesting all invariant forms, and there are other things that also need more thought.)

Also, I don't think my comments ask you to "entirely refactor" the code (unless you are referring to using libgap, which is not really relevant to my latest comments). They seem to me to be pretty typical (and pretty minor) reviewer comments. The biggest one is probably the suggestion to treat symmetric forms and alternating forms separately, but that should just be a few additional lines of code. (If it turns out to be a big job, it might be ok to leave it to a follow-up PR.)

I apologize if "half-baked" was a poor choice of words, but it seems clear to me that the code is not ready to be merged, for the reasons I have listed, and my understanding was that you thought it was fine to have problems like this (even known bugs!) in code that is merged. It's not, and reviewers are supposed to point out these problems (and possible problems).

[DaveWitteMorris]

I should also have said that it's fine if you disagree with some specific reviewer comments. Just give a good reason why.

[jacksonwalters]

By moving this code, @jacksonwalters you are making it more general-purpose, so it should work in those cases. I also agree that the return time should be consistent as the user ahead of time will not know the dimension of the invariant forms ahead of time. I agree with most of @DaveWitteMorris's suggestions, although we should be careful to abide by the "the perfect is the enemy of the good" mantra. Let's try to spend a bit of time to optimize this now and fix any bugs we do come across.

@tscrim Yes, the goal of this PR is to make it more general purpose. My idea was to do it in stages: 1) move this out of UnitaryRepresentations to general SymmetricGroupRepresentation so it at least has doctests 2) actually allow for the case of a higher dim'l solution space (we were only returning nullspace.basis()[0] before) 3) find out what needs to happen to move this all the way out to finite group representations generally (which there doesn't seem to be an obvious general class for.

Solving ( n 2 ) equations is much faster than solving n 2 equations. So splitting it into alternating and symmetric form computations would be good, with the general method returning the combination of the two (as every matrix splits into a symmetric and skew-symmetric parts).

With respect to bugs, this only happens if you replace G with a 2-generating set of $S_n$. I'm not sure I would call that a bug, since we're not doing that. I would call it an obstacle to a further optimization. In any case, I'm happy to look into this and try to fix it and get it working, and yes it may reveal an actual bug or something with the code as is.

Reducing to $\binom{n}{2}$ equations using the symmetry condition is a good optimization. Perhaps this will need an option like form_type=symmetric/form_type=alternating. This will require some further doctests to find forms that are in fact alternating, as I don't think we're currently finding any.

Note also you do not need to create another representation. You already have ρ ( g ) from the class you are working with.

If I try to call representation_matrix, then in the unitary case it tries to call symmetric_invariant_bilinear_form() to construct it, which then would call representation_matrix again and creates an infinite loop. My current thought is to keep track of the form type in say self._form, and then if self._form == unitary, set self._specht to a Specht representation if it's not already set, and otherwise use the representation_matrix of the class (which could be specht, seminormal, or `orthogonal). Let me know what you think.

In answer to your question: no, I do not agree that this code has already been reviewed and merged. Most of my comments concern problems that arose in converting it to a class function. I do make mistakes when reviewing PR's, but I think all of my comments this time are valid, except perhaps the suggestion to use libgap, which you have already rejected, and the two that I phrased as questions (but now I'm pretty sure that my concern about characteristic two is correct). (However, I think something like "both" would be better than None for requesting all invariant forms, and there are other things that also need more thought.)

@DaveWitteMorris I think we are (hopefully) in agreement - this current PR of course has not been merged, and I'll agree that more work is necessary if it is to be separated from UnitaryRepresentation of the symmetric group. What I meant was that this is not the first time this code has been presented to SageMath, and everything except the doctests was merged in two other PRs #38455 , #39200, #38456.

Also, I don't think my comments ask you to "entirely refactor" the code (unless you are referring to using libgap, which is not really relevant to my latest comments). They seem to me to be pretty typical (and pretty minor) reviewer comments. The biggest one is probably the suggestion to treat symmetric forms and alternating forms separately, but that should just be a few additional lines of code. (If it turns out to be a big job, it might be ok to leave it to a follow-up PR.)

We had very lengthy discussions about using libgap to do parts of this, and after implementing all of that code, we translated everything back to Sage code. Again, while Im not opposed to making this code as good as can be, I'm hesitant to just start rewriting a bunch of code which works for its current intended purpose (getting change-of-basis matrices for unitary representations). I can try to write some minimal examples in GAP, and see how it goes. It's not that I want to reject this, it's that I think we tried it, and because forms wasn't a default package, we didn't want to use it.

I apologize if "half-baked" was a poor choice of words, but it seems clear to me that the code is not ready to be merged, for the reasons I have listed, and my understanding was that you thought it was fine to have problems like this (even known bugs!) in code that is merged. It's not, and reviewers are supposed to point out these problems (and possible problems).

No worries. I also apologize if I was being curt. I also want this code to be as good as it can possibly, and I'm willing to streamline and improve it. I agree that the code in this PR is not where it should be, but I maintain that the code we did merge in was thoroughly tested and as good as it could be for the cases we were interested in. Now that we're generalizing, I think it can be a lot better. I'll try to address the issue you mentioned and I concur with Travis that they're good observations.

[jacksonwalters]

I think it makes more sense to move this to the specht representation. As written, it's circular because it depends on the type of representation for the group to act. In other words, moving it up to SymmetricGroupRepresentation_generic_class is too general.

[jacksonwalters]

@DaveWitteMorris I've resolved the bug with the small generating set. It was just that [1,2] and [1,2,3,4,5] both correspond to the identity permutation, since the map is from the index of the list to the element. I should've used [2,1] and [2,3,4,5,1]. There are different ways to construct permutations, but everything seems to be working fine now.

It makes sense, too, because the basis was just every elementary matrix since there were no conditions to be satisfied.

[jacksonwalters]

@DaveWitteMorris @tscrim Here is an updated version using a symmetric=True to cut the number of variables from $n^2$ to $\binom{n+1}{2}$. It can be adapted pretty easily for alternating forms. Seems to be working.

def invariant_symmetric_bilinear_matrix(SGA,partition,symmetric=False):
    """
    Computes the matrix of a S_n-invariant symmetric bilinear form.
    Sets up and solves system of linear equations based on writing U as an unknown in polynomial ring generators. 
    The equations are \rho(g)^T*U*\overline{\rho(g)} = \lambda_g*U where \lambda_g = \det(\rho(g))\overline{\det(\rho(g))}.
    The variables for U can be extracted to yield a matrix over GF(q^2) for each g.
    These are stacked to get the overall system, and we find the one dim'l null space to get a solution vector, and format as a matrix.
    NOTE: one could also form the Kroenecker products \rho(g) \otimes \rho(g)^{-1 T} to explicitly obtain the system.
    NOTE: the solution space is typically 1 dim'l, except in "modular" cases where p divides |G|, and there is multiplicity in the decomposition factors.
    """
    n = sum(partition)
    G = Permutations(n)
    G_gens = [G([2, 1]), G([i%n+1 for i in range(1,n+1)])] #small gen set {(1 2), (1 2 ... n)}. can use G.gens() generally
    F = SGA.base_ring()
    specht_module = SGA.specht_module(partition)
    rho = specht_module.representation_matrix
    d_rho = specht_module.dimension()
    if symmetric:
        R = PolynomialRing(F, 'u', d_rho * (d_rho + 1) // 2)
    else:
        R = PolynomialRing(F, 'u', d_rho**2)
    U_vars = R.gens()
    if symmetric:
        U_temp = matrix(R, d_rho, d_rho, lambda i, j: U_vars[i * d_rho - (i * (i + 1)) // 2 + j] if i <= j else 0)
        U_temp += U_temp.transpose() - diagonal_matrix(U_temp.diagonal())
    else:
        U_temp = matrix(R, d_rho, d_rho, U_vars)

    def augmented_matrix(g):
        rho_g = rho(g)
        equation_matrix = rho_g.transpose() * U_temp * rho_g.conjugate() - U_temp
        augmented_system = []
        for i in range(d_rho):
            for j in range(d_rho):
                linear_expression = equation_matrix[i, j]
                row = [linear_expression.coefficient(u) for u in U_vars]
                augmented_system.append(row)
        return augmented_system

    total_system = sum((augmented_matrix(g) for g in G_gens), [])
    null_space = matrix(F, total_system).right_kernel()
    U_list = []
    if symmetric:
        for B in null_space.basis():
            upper_triangle = matrix(F, d_rho, d_rho, lambda i, j: B[i * d_rho - (i * (i + 1)) // 2 + j] if i <= j else 0)
            U_list.append(upper_triangle + upper_triangle.transpose() - diagonal_matrix(upper_triangle.diagonal()))
    else:
        U_list = [matrix(F, d_rho, d_rho, B) for B in null_space.basis()]
    return U_list

[dimpase]

the solution space for reducible representations can easily be as big as the number of irreducible subrepresentations, or even bigger.

so, "typically 1-dimensional" isn't really correct without extra assumptions

[jacksonwalters]

the solution space for reducible representations can easily be as big as the number of irreducible subrepresentations, or even bigger.

so, "typically 1-dimensional" isn't really correct without extra assumptions

Specht modules are irreducible in characteristic $p > n$. I note the modular case, where the Specht modules are indeed reducible.

Is it true that there are no alternating forms for the symmetric group? Soon there'll be code for this, but it would be good to know some cases to expect.

[dimpase]

Yes, Specht modules are irreducible, in char 0 at least, but what does this have to do with the general functionality this PR is implementing?

Yes, all the representations of the symmetric groups are real, so no alt. forms will arise in this case.

[jacksonwalters]

Yes, Specht modules are irreducible, in char 0 at least, but what does this have to do with the general functionality this PR is implementing?

Yes, all the representations of the symmetric groups are real, so no alt. forms will arise in this case.

When you say representations of the symmetric group are real, you mean defined over $GF(q)$ rather than $GF(q^2)$?

I'm viewing this PR as a half-step on the way to generalizing this for any representation of a finite group. For now, I'd just like this to be a method with appropriate doctests (noting that without this PR, it doesn't ever return more than one linearly independent matrix). This PR will include at least two optimizations now (using a 2-element generating set, cutting the number of variables in half by imposing the symmetry condition).

Currently moving the code to Specht is causing issues.

[dimpase]

I was talking about char 0 case. I don't know much about the modular case. @darijgr is an expert on all things $S_n$ and could say much more.

[darijgr]

Sorry for not reading the whole thread, but here are a few things I do know about symmetric groups:

For any Specht module S^D over any commutative ring (with D being any diagram -- not necessarily a partition shape), there is a canonical S_n-invariant symmetric bilinear form on S^D obtained by restricting the obvious form on the Young module M^D (the one which makes all tabloids orthonormal) to the submodule S^D. When the ring is a field of characteristic 0, this form is nondegenerate. To see this, it suffices to show it when the field is QQ (after all, this is a question of a determinant being nonzero), and there we can use the positive definiteness of the form on M^D and the fact that restricting a positive definite form to a subspace yields a positive definite form again. On the other hand, if the field has characteristic p between 1 and n, then this form can have a nontrivial kernel, even when D is a partition shape. This kernel is used to construct the actual irreps (see §5 in Mark Wildon's https://www.ma.rhul.ac.uk/~uvah099/Maths/Sym/SymGroup2014.pdf ).

Over a field of characteristic 0, the actual partition-shaped Specht modules are irreducible, so the space of S_n-invariant bilinear forms is 1-dimensional, and thus all such forms are symmetric because the canonical one is. Over other fields or for other shapes, this might not be true. I definitely see how it's false for non-partition shapes (e.g., because of irreps occurring with multiplicities), but I'm less definite about the partition-shaped char-p case. Might be a good question. Corollary 13.17 in James's "The Representation Theory of the Symmetric Groups" might be useful.