unitary DFT for the symmetric group algebra

[jacksonwalters]

Problem Description

At present the only two DFT's available for the symmetric group algebra are the seminormal basis and modular DFT. Neither of these are unitary transformations.

Proposed Solution

In Beals ['97], Quantum computation of Fourier transforms over symmetric groups he remarks that the DFT can be made unitary with two modifications:

1. include a factor of \sqrt{d_\rho/|G|} in front of each Fourier coefficient at \rho
2. the rep'ns being summed over must be themselves unitary

Alternatives Considered

The seminormal form is orthogonal in that $AA^T$ is diagonal, but the columns are not orthonormal. Young's orthogonal form has the property $AA^T = Id$, but not $AA^* = Id$.

Additional Information

No response

Is there an existing issue for this?

• I have searched the existing issues for a bug report that matches the one I want to file, without success.

[jacksonwalters]

Weyl's unitary trick effectively solves this for the complex numbers and number fields.

For finite fields, the situation is more delicate.

There is: "orthogonal" – for Young’s orthogonal representation for sage.combinat.symmetric_group_representations.SymmetricGroupRepresentation.

However, when ring=GF(7) is set, we obtain TypeErrors for partition=[2,1] corresponding to n=3. What appears to be happening is sqrt(3) is not able to be converted to an element of the finite field, nor are extensions being created. Perhaps it is possible to work over an extension so that these matrices become "orthogonal".

TypeError: self must be a numeric expression

[jacksonwalters]

PR: #38455

The error with square roots has been resolved by first using self._ring on beta before passing to 2x2 matrix element in symmetric_group_representations.py on line 662.

Regarding unitary representations of finite groups over finite fields, see:

https://mathoverflow.net/questions/327823/unitary-representations-of-finite-groups-over-finite-fields

One must check that the action of the q^th power of Frobenius on the Brauer character (which is just the complex character in the non-modular case) is complex conjugation. Even then, it is not immediately clear to me that Young's orthogonal representations will be unitary, or which sesquilinear form they will be unitary with respect to, exactly.

[jacksonwalters]

With the fix implemented locally, we do get orthogonal representations over GF(q^2). However, these representations are not unitary w.r.t. the form <x,y> = \sum_i x_i*y_i^q since x |--> x^q does not fix elements outside of GF(q). It may be that the matrices are unitary w.r.t. a different sesquilinear form.

[dimpase]

Perhaps the invariant form can be obtained by summing $g\overline{g}^\top$ over the group. (At least this works in characteristic 0 case, here it's a question of luck)

[jacksonwalters]

@dimpase I tried this after doing the characteristic zero case (complex numbers, number fields) and it didn't seem to work. This led to this question:

https://math.stackexchange.com/questions/4952613/is-there-a-version-of-weyls-unitary-trick-for-positive-characteristic

It does seem natural to try to just sum. One similar thought I haven't tried is averaging the sesquilinear form $&lt;x,y&gt;=\sum_i x_i y_i^q$ to get a $G$-invariant form.

Updates: representations of the symmetric group have Frobenius-Schur indicator +1, so Lemma 4.4.1 in the linked MO question doesn't apply. However, this question seems of relevance:

https://mathoverflow.net/questions/271932/formula-for-the-frobenius-schur-indicator-of-a-finite-group

This appears to imply (please check) that since $\chi$ is real, and irreducible rep'ns of the symmetric group are self-dual (theorem 1.5 in Gordon James' book), and we know the F-S indicator is +1, then as long as $d(\chi,\phi)$ is odd (it is often +1, see tables in Appendix of James' book), then there should exist a $G$-invariant symmetric bilinear form.

I just don't know how to find it. Perhaps I am overlooking something obvious.

[dimpase]

an obvious way to find invariant forms is by linear algebra in the representation $\rho(g)$. Take a matrix of unknowns $U$, and for each generator $g$ of $G$, write $\rho (g)U=\lambda_g U\overline{\rho( g)}^{-1^\top}$, where $\lambda_g:=\det \rho (g) \det \overline{\rho( g)}^{\top}$.
In the solution space pick a nonzero element.

EDIT: nontrivial $\lambda_g$ arise when one has sesquilinear forms. (Corrected above). Also note that a convenient way to form the system of equations is to take Kronecker products $\rho (g)U\otimes \lambda_g U\overline{\rho( g)}^{-1^\top}$. Then $U$ can be recovered from a common, for all the generators, eigenvector with eigenvalue 1, of these Kronecker products. such an eigenvector gives you an embedding into the full unitary group of the corresponding form - I don't recall off top of my head whether it will be unique for the irreducible representation $\rho$ in positive characteristic. @jacksonwalters

[jacksonwalters]

@dimpase Ah, okay, great. Thanks very much! This is very concrete. I will attempt to implement this, ideally by the end of the week, and see what I can say.

[dimpase]

There is also a GAP package forms which does this in GAP, for matrix groups over finite fields. We probably should add it to the list of GAP packages we can easily install. For now, you can install gap_packages SPKG and then use InstallPackage to install forms.

[jacksonwalters]

I tried installing the GAP package forms via extracting the .tar to the appropriate directory ($SAGE_ROOT/local/share/gap/pkg/). It appears to load the package successfully, then I get all kinds of errors and claims that functions don't exist:

sage: gap('LoadPackage("forms");')
true
sage: gap('Packages();')
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
File /private/var/tmp/sage-10.4-current/local/var/lib/sage/venv-python3.12.4/lib/python3.12/site-packages/sage/interfaces/gap.py:700, in Gap_generic._eval_line(self, line, allow_use_file, wait_for_prompt, restart_if_needed)
    699     error = error.replace('\r', '')
--> 700     raise RuntimeError("%s produced error output\n%s\n   executing %s" % (self, error, line))
    701 if not normal:

RuntimeError: Gap produced error output
Error, Variable: 'Packages' must have a value

   executing \$sage2:=Packages();;;

During handling of the above exception, another exception occurred:

RuntimeError                              Traceback (most recent call last)
File /private/var/tmp/sage-10.4-current/local/var/lib/sage/venv-python3.12.4/lib/python3.12/site-packages/sage/interfaces/expect.py:1519, in ExpectElement.__init__(self, parent, value, is_name, name)
   1518 try:
-> 1519     self._name = parent._create(value, name=name)
   1520 # Convert ValueError and RuntimeError to TypeError for
   1521 # coercion to work properly.

File /private/var/tmp/sage-10.4-current/local/var/lib/sage/venv-python3.12.4/lib/python3.12/site-packages/sage/interfaces/interface.py:517, in Interface._create(self, value, name)
    516 name = self._next_var_name() if name is None else name
--> 517 self.set(name, value)
    518 return name

File /private/var/tmp/sage-10.4-current/local/var/lib/sage/venv-python3.12.4/lib/python3.12/site-packages/sage/interfaces/gap.py:1354, in Gap.set(self, var, value)
   1353 cmd = ('%s:=%s;;' % (var, value)).replace('\n', '')
-> 1354 self._eval_line(cmd, allow_use_file=True)

File /private/var/tmp/sage-10.4-current/local/var/lib/sage/venv-python3.12.4/lib/python3.12/site-packages/sage/interfaces/gap.py:734, in Gap_generic._eval_line(self, line, allow_use_file, wait_for_prompt, restart_if_needed)
    733     else:
--> 734         raise RuntimeError(exc)
    736 except KeyboardInterrupt:

RuntimeError: Gap produced error output
Error, Variable: 'Packages' must have a value

   executing \$sage2:=Packages();;;

During handling of the above exception, another exception occurred:

TypeError                                 Traceback (most recent call last)
Cell In[2], line 1
----> 1 gap('Packages();')

File /private/var/tmp/sage-10.4-current/local/var/lib/sage/venv-python3.12.4/lib/python3.12/site-packages/sage/interfaces/interface.py:299, in Interface.__call__(self, x, name)
    296         pass
    298 if isinstance(x, str):
--> 299     return cls(self, x, name=name)
    300 try:
    301     # Special methods do not and should not have an option to
    302     # set the name directly, as the identifier assigned by the
    303     # interface should stay consistent. An identifier with a
    304     # user-assigned name might change its value, so we return a
    305     # new element.
    306     result = self._coerce_from_special_method(x)

File /private/var/tmp/sage-10.4-current/local/var/lib/sage/venv-python3.12.4/lib/python3.12/site-packages/sage/interfaces/expect.py:1524, in ExpectElement.__init__(self, parent, value, is_name, name)
   1522 except (RuntimeError, ValueError) as x:
   1523     self._session_number = -1
-> 1524     raise TypeError(*x.args)
   1525 except BaseException:
   1526     self._session_number = -1

TypeError: Gap produced error output
Error, Variable: 'Packages' must have a value

   executing \$sage2:=Packages();;;

[dimpase]

Does it work in a separate GAP session?
(start Sage shell: ./sage -sh and then type gap)

[jacksonwalters]

@dimpase Ah, that fixed it! Thank you.

gap> LoadPackage("forms");
---------------------------------------------------------------------
Loading 'Forms' 1.2.12 (30/08/2024)
by John Bamberg (http://school.maths.uwa.edu.au/~bamberg/)
   Jan De Beule (http://www.debeule.eu)
For help, type: ?Forms 
---------------------------------------------------------------------
true
gap> IsPackageLoaded("forms");
true
gap> gf := GF(8);
GF(2^3)
gap> r := PolynomialRing(gf,3);
GF(2^3)[x_1,x_2,x_3]
gap> poly := r.1^2 + r.2 * r.3;
x_1^2+x_2*x_3
gap> form := QuadraticFormByPolynomial(poly,r);
< quadratic form >

I'm still parsing the linear algebra approach. I generated some (non-functional) code but I need to tweak it. The main difficulty right now is just which rings to define these matrices over, in particular U which is symbolic. I want to solve these systems over GF(q^2) but I'm not quite sure how to do symbolic linear algebra in Sage.

TypeError: positive characteristic not allowed in symbolic computations

q = 7
n = 3
#partitions should be such that decomposition numbers d(\chi,\phi) is odd
partition = [2,1]

# Define the group G and its representation as a Specht module
SGA = SymmetricGroupAlgebra(GF(q^2), n)
SM = SGA.specht_module(partition)

# Get representation from Specht module
rho = SM.representation_matrix

# Define the matrix of unknowns U (2x2 matrix in this case)
n = 2  # Dimension of U
# Create symbolic variables in the symbolic ring (SR)
u_11, u_12, u_21, u_22 = var('u_11 u_12 u_21 u_22')
U = Matrix(SR, n, n, [[u_11, u_12], [u_21, u_22]])  # Use SR for symbolic ring

# Create a list to hold the Kronecker products
kron_products = []

# Iterate over each generator g in G
for g in G.gens():
    # Compute rho(g) - you need to replace this with your actual representation logic
    rho_g = rho(Permutation(g))
    rho_g_inv_T = rho_g.inverse().transpose()

    # Compute lambda_g
    det_rho_g = det(rho_g)  # This will be an element in GF(q)
    lambda_g = det_rho_g * (det_rho_g ** q)  # Correct conjugation operation

    # Convert lambda_g to a matrix for the Kronecker product
    lambda_g_matrix = lambda_g * identity_matrix(GF(q^2), n)

    # Form the Kronecker product
    # Note: Use rho_g and its inverse correctly in the tensor product
    kron_prod = (rho_g * U).tensor_product(lambda_g_matrix * U * rho_g_inv_T)
    kron_products.append(kron_prod)

# Now, we need to find the common eigenvector with eigenvalue 1
common_eigenvectors = []
for k in kron_products:
    eigenspaces = k.eigenvectors_right(1)  # Find eigenvectors for eigenvalue 1
    common_eigenvectors.extend(eigenspaces)

# Check if there are non-zero solutions
non_zero_solutions = [vec for eigens in common_eigenvectors for vec in eigens if vec[0] != 0]

# Display results
if non_zero_solutions:
    print("Found non-zero common eigenvectors:")
    for vec in non_zero_solutions:
        print(vec)
else:
    print("No non-zero common eigenvectors found.")

[dimpase]

you don't need to make U a variable!
Just write down these Kronecker products, stack them up in a matrix, etc

anyway, the symbolic ring is the last resort, and it doesn't work in positive char. So all you have are polynomial rings.

[dimpase]

Computing an eigenvector $v$ of $A$ with a known eigenvalue $\lambda$ is the same as finding a vector in the nullspace of $A-\lambda I$, and the nullspace (or kernel) are methods of Sage matrices.
So you can stack up these matrices $A-\lambda I$ and call the (right?) nullspace for this non-square matrix.
This way avoids dealing with explicit variables etc.

[jacksonwalters]

@dimpase Yes, one can find eigenvalues by computing $\text{null}(A-\lambda I)$ and we just want the common $\lambda=1$ eigenspace for each $g \in G$, so we can stack the $A-\lambda I$ where $A=\rho(g)U\otimes \lambda_g U\overline{\rho(g)}^{-1 \top}$.

Can you just write in code how you would instantiate $U$? You seem to be saying two contradictory things:

1. $U$ is a matrix of unknowns
2. $U$ is not a variable

So I don't know how to create $U$, therefore any of these matrices because U=matrix(GF(q^2),d_rho,d_rho) is just the zero matrix. I mean, it's clear that this results in a system of linear equations, but I don't see how it wouldn't be in the variables $u_{ij}$, so I'm confused.

Also, you write $\det(\overline{\rho(g)}^{\top})$. Isn't the determinant of the transpose just the determinant?

EDIT: Okay, beginning to see. We are just setting $\lambda=1$, so this system can be arranged so that $U$ is the unknown vector of a linear system which is non-square, so we just need to find the matrix itself. $U$ won't appear explicitly, e.g. when solving $Ax = b$ we use the augmented system $[A|b]$ with no mention of $x$. I think I do need to work over a polynomial ring $GF(q^2)[u0,...,u_{d_\rho^2}]$ though.

[dimpase]

yeah, right, matrix of unknowns was there for the purpose of explanation. By the way, I am not sure that $\lambda=1$ always holds in sesquilinear case. There perhaps could be different $\lambda$ for different groups generators.

In the end you can fill in the entries of $U$ with the solution entries of the system, no need for variables

[jacksonwalters]

@dimpase Since $U$ occurs in both factors of the Kroenecker product, don't we get terms quadratic in the $u_{ij}$?

[dimpase]

no, you "vectorise" $U$ and get its $n^2$ entries as a vector. So you need to fill them back in as a matrix.

[jacksonwalters]

Yes, we can write $U=(u_0, \ldots, u_{n^2})$ as a vector where $n=d_\rho$ , but when I do $\rho(g)U \otimes \lambda_g U \rho(g)$ I am forming the block matrix with first block $(\rho(g)U)_{11}\lambda_g U \rho(g)$ so I will have terms like $u_0^2$.

[dimpase]

Yes, we can write U = ( u 0 , … , u n 2 ) as a vector where n = d ρ , but when I do ρ ( g ) U ⊗ λ g U ρ ( g ) I am forming the block matrix with first block ( ρ ( g ) U ) 11 λ g U ρ ( g ) so I will have terms like u 0 2 .

Why do you need to do ρ ( g ) U ⊗ λ g U ρ ( g ) ? All I was saying is that the action of the group on the space of $n\times n$ matrices defining sesquilinear forms is the tensor product of the suitable $n$-dimensional representations.
So for fixed $A$, $B$ the system of linear equations $AU=UB$ in $(u_{ij})=U$ can be written as something like $(A\otimes B) u=0$, with $u=(u_{11},u_{12},\dots, u_{nn})$. You solve the latter system to obtain a solution of the former. $u_{ij}$ as explicit variables don't appear anywhere.

[jacksonwalters]

Here is a version without using the Kroenecker products which should work (at least in cases where the decomposition number is odd) for n=3, q=3, partition = [2,1].

EDIT: Corrected swapped indices in conjugate matrix.

# Setup
q = 3
n = 3
partition = [2, 1]

# Define the group G and its representation as a Specht module
SGA = SymmetricGroupAlgebra(GF(q^2), n)
SM = SGA.specht_module(partition)
G = SGA.group()

# Get representation from Specht module
rho = SM.representation_matrix
d_rho = SM.dimension()

# Initialize U as a matrix of variables over GF(q^2)
R = PolynomialRing(GF(q^2), 'u', d_rho^2)
U_vars = R.gens()  # List of variable generators for U
U = Matrix(R, d_rho, d_rho, U_vars)  # U is a d_rho x d_rho matrix of variables

#define conjugation as x |--> x**q, an order two automorphism of F_q^2. note x**q == x for x \in F_q.
def conjugate_pos_char(A):
    assert A.nrows() == A.ncols()
    field_size = A.base_ring().order()
    q = sqrt(field_size) if field_size.is_square() else field_size
    return matrix(GF(q**2),[[A[i][j]**q for j in range(A.nrows())] for i in range(A.nrows())])

# for each generator of G, form the augmented system 
def augmented_matrix(g):

    rho_g = rho(Permutation(g))
    rho_g_conj_inv_T = conjugate_pos_char(rho_g.inverse().transpose())

    # Compute lambda_g
    det_rho_g = det(rho_g)
    lambda_g = det_rho_g * (det_rho_g ** q)

    # Form the matrix equation
    equation_matrix = rho_g * U - lambda_g * U * rho_g_conj_inv_T

    # Initialize a list to hold rows of the augmented system
    augmented_system = []

    # Extract coefficients for each linear equation in the matrix
    for i in range(d_rho):
        for j in range(d_rho):
            # Get the (i, j) entry of the equation matrix, which is a linear combination of the u variables
            linear_expression = equation_matrix[i, j]
        
            # Extract the coefficients of each u_k in the linear expression
            row = [linear_expression.coefficient(u) for u in U_vars]
        
            # Append the row to the augmented system
            augmented_system.append(row)

    # Convert the augmented system to a matrix
    A = Matrix(GF(q^2), augmented_system)

    return A

total_system = matrix(GF(q^2),0,d_rho^2)
for g in G:
    total_system = total_system.stack(augmented_matrix(g))

null_space = total_system.right_kernel()

print(null_space)

Vector space of degree 4 and dimension 0 over Finite Field in z2 of size 3^2
Basis matrix:
[]

[dimpase]

You made a typo in conjugate_pos_char(A). Its last line should be

return matrix(GF(q**2),[[A[i][j]**q for j in range(A.nrows())] for i in range(A.nrows())])

i.e the inner loop is over j.
With this change your code prints

Vector space of degree 4 and dimension 1 over Finite Field in z2 of size 3^2
Basis matrix:
[1 1 1 1]

[jacksonwalters]

Ah, great! Thanks so much for spotting that. Yes, getting lots of interesting forms for different choices of n and partition, e.g. q=3, n=5, [3,1,1]:

Vector space of degree 36 and dimension 1 over Finite Field in z2 of size 3^2
Basis matrix:
[0 1 2 1 2 0 1 0 1 1 0 2 2 1 0 0 1 2 1 1 0 0 1 1 2 0 1 1 0 1 0 2 2 1 1 0]

I'll do some testing and verification to ensure these really define G-invariant sesquilinear forms tomorrow, and for which cases are we getting solutions. If for a given q, n there are solutions for every partition, it may be possible to construct a unitary DFT.

[jacksonwalters]

The code to generate the matrices, as @dimpase described, appears to be working for all cases.

I wrote some code to check the desired properties of the associated form. Since the linear algebra was set up to guarantee the three conditions (symmetric, G-invariant, bi/sesquilinear), it's surprising that the conditions are sometimes met and sometimes not met. It could be a typo or bug on my end. Code attached (updated):

[(q,la,check_form_properties(2,la)) for la in Partitions(3)]
[(3, [3], True), (3, [2, 1], True), (3, [1, 1, 1], True)]

[(q,la,check_form_properties(2,la)) for la in Partitions(4)]
[(3, [4], True),(3, [3, 1], False),(3, [2, 2], True),(3, [2, 1, 1], False),(3, [1, 1, 1, 1], True)]

[(q,la,check_form_properties(2,la)) for la in Partitions(5)]
[(3, [5], True),(3, [4, 1], True),(3, [3, 2], False),(3, [3, 1, 1], True),(3, [2, 2, 1], False),(3, [2, 1, 1, 1], True),(3, [1, 1, 1, 1, 1], True)]

[(q,la,check_form_properties(3,la)) for la in Partitions(3)]
[(3, [3], True), (3, [2, 1], False), (3, [1, 1, 1], True)]

[(q,la,check_form_properties(3,la)) for la in Partitions(4)]
[(3, [4], True),(3, [3, 1], False),(3, [2, 2], False),(3, [2, 1, 1], False),(3, [1, 1, 1, 1], True)]

[(q,la,check_form_properties(3,la)) for la in Partitions(5)]
[(3, [5], True),(3, [4, 1], True),(3, [3, 2], False),(3, [3, 1, 1], True),(3, [2, 2, 1], False),(3, [2, 1, 1, 1], True),(3, [1, 1, 1, 1, 1], True)]

Made a few adjustments from the previous version:

"""

NOTE: the theorem (https://mathoverflow.net/questions/271932/formula-for-the-frobenius-schur-indicator-of-a-finite-group) only guarantees a bilinear form, not a sesqilinear form.

""";

#define conjugation as x |--> x**q, an order two automorphism of F_q^2. note x**q == x for x \in F_q.
def conjugate_pos_char(A):
    assert A.nrows() == A.ncols()
    field_size = A.base_ring().order()
    q = sqrt(field_size) if field_size.is_square() else field_size
    return matrix(GF(q**2),[[A[i][j]**q for j in range(A.nrows())] for i in range(A.nrows())])

def invariant_symmetric_sesquilinear_matrix(q,partition):
    """
    Computes the matrix of a S_n-invariant symmetric sesquilinear 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)U = \lambda_g*U*\rho(g)^{-1 T} where \lambda_g = \det(\rho(g))\overline{\det(\rho(g))}^T.

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

    # Define the group G and its rep'n as a Specht module, dimension
    n = sum(partition)
    SGA = SymmetricGroupAlgebra(GF(q^2), n)
    SM = SGA.specht_module(partition)
    G = SGA.group()
    rho = SM.representation_matrix
    d_rho = SM.dimension()
    
    # Initialize U as a matrix of variables over GF(q^2)
    R = PolynomialRing(GF(q^2), 'u', d_rho^2)
    U_vars = R.gens()  # List of variable generators for U
    U = Matrix(R, d_rho, d_rho, U_vars)  # U is a d_rho x d_rho matrix of variables
    
    # for each generator of G, form the augmented system 
    def augmented_matrix(g):
    
        rho_g = rho(Permutation(g))
        rho_g_conj_inv_T = conjugate_pos_char(rho_g.inverse().transpose())
    
        # Compute lambda_g
        det_rho_g = det(rho_g)
        lambda_g = det_rho_g * (det_rho_g ** q)
    
        # Form the matrix equation
        equation_matrix = rho_g * U - lambda_g * U * rho_g_conj_inv_T
    
        # Initialize a list to hold rows of the augmented system
        augmented_system = []
    
        # Extract coefficients for each linear equation in the matrix
        for i in range(d_rho):
            for j in range(d_rho):
                # Get the (i, j) entry of the equation matrix, which is a linear combination of the u variables
                linear_expression = equation_matrix[i, j]
            
                # Extract the coefficients of each u_k in the linear expression
                row = [linear_expression.coefficient(u) for u in U_vars]
            
                # Append the row to the augmented system
                augmented_system.append(row)
    
        # Convert the augmented system to a matrix
        return Matrix(GF(q^2), augmented_system)
    
    total_system = matrix(GF(q^2),0,d_rho^2)
    for g in G:
        total_system = total_system.stack(augmented_matrix(g))
    
    #compute the null space of the overall matrix
    null_space = total_system.right_kernel()
    
    #return a d_rho x d_rho matrix over GF(q^2) from the 1 dim'l null space given as vector
    return matrix(GF(q^2),d_rho,d_rho,null_space.basis()[0])

def sesquilinear_form(x,y,U):
    return x*U*y.transpose()

#ensure the resulting form is G-invariant, symmetric, bi/sesquilinear
def check_form_properties(q,partition):
    #define the representation matrix corresponding to q, partition
    SGA = SymmetricGroupAlgebra(GF(q^2),sum(partition))
    SM = SGA.specht_module(partition)
    rho = SM.representation_matrix
    d_rho = SM.dimension()
    G = SGA.group()

    #define variables as polynomial generators
    R_xy = PolynomialRing(GF(q^2), d_rho, var_array='x,y')
    x = matrix([R_xy.gens()[2*i] for i in range(d_rho)])
    y = matrix([R_xy.gens()[2*i+1] for i in range(d_rho)])
    R_xy_lambda = PolynomialRing(R_xy,'lambda')
    lambda_ = R_xy_lambda.gens()[0]

    #compute the bilinear form matrix. cast over ring
    U_mat = invariant_symmetric_sesquilinear_matrix(q,partition)
    U_form = matrix(R_xy_lambda,U_mat); U_form
    
    #check symmetric property
    symmetric = sesquilinear_form(x,y,U_form) == sesquilinear_form(y,x,U_form)
    #check G-invariance property
    G_invariant = all(sesquilinear_form((rho(g)*x.transpose()).transpose(),(rho(g)*y.transpose()).transpose(),U_form) == sesquilinear_form(x,y,U_form) for g in G)
    #check sesquilinear property. ISSUE: lambda_^q is a power of the ring generator, i.e. doesn't simplify.
    first_arg = sesquilinear_form(lambda_*x,y,U_form) == lambda_*sesquilinear_form(x,y,U_form)
    second_arg = sesquilinear_form(x,lambda_*y,U_form) == lambda_*sesquilinear_form(x,y,U_form) #need to amend for conjugation

    return symmetric and G_invariant and first_arg and second_arg

[dimpase]

you'd use assert statements to make things fail clearly, without having to guess which of the 3 properties failed

[jacksonwalters]

@dimpase Yes, sorry, I deleted the print statements I was using. G-invariance appears to be the one failing in cases where it fails overall.

I didn’t include the \lambda_g in the checks. Can you explain where that term came from? I don’t see why it’s necessary for G-invariance:

[jacksonwalters]

This question/answer seems to imply if a module V over a finite field supports both G-invariant symmetric and skew-symmetric forms, then the number of orthogonal constituents and the number of symplectic constituents in its decomposition must both be zero, so that V consists solely of pairs of non-self-dual components:

https://math.stackexchange.com/questions/832173/structure-of-g-invariant-bilinear-forms-over-finite-fields

[dimpase]

@dimpase Yes, sorry, I deleted the print statements I was using. G-invariance appears to be the one failing in cases where it fails overall.

I didn’t include the $\lambda_g$ in the checks. Can you explain where that term came from? I don’t see why it’s necessary for $G$-invariance:

I guess I didn't take into account that your $G$ is a symmetric group. If you are putting in Sage a general function to compute such forms for a finite group, then $\lambda_g$ should be present.
Indeed, suppose you want $gU\overline{g}^\top=U$ for all $g\in \rho(G)$, where $\rho$ is our representation. Take determinants on both sides, and cancel $\det U$ --- the latter is nonzero for irreducible modules - however, in finite group theory one sometimes works with reducible modules, e.g. subgroups of $GO(2m+1, q)$. This implies $\det g \det (\overline{g}^\top)=\det g \det \overline{g}=1$, for all $g$. It's an extra restriction on $\rho$, but do you want it? To drop this restriction, use $\lambda_g:=\det g \det \overline{g}$ - you can still have examples where $\lambda_g\neq 1$, which make perfect sense geometrically, as in the underlying projective space you'd still have a $G$-invariant structure.

The representations you work here with are $GF(q^2)$-representations, and all you know is that $\lambda_g\in GF(q)^*$, so for $q=2$ no $\lambda_g$ is needed, but otherwise it's unclear if it can always be dropped. The determinant is a group homomorphism (and a 1-dimensional representation of the group). This seems to imply that the only possible values of $\det g$ are $\pm 1$, so indeed, $\lambda_g=1$ in your case.

[jacksonwalters]

@dimpase I'm not sure I completely understand. We're already demanding that $\rho(g)U\overline{\rho(g)}^T = U$, so if that condition is satisfied then we have $\lambda_g \det(U) = \det(U)$. We only have $\lambda_g \ne 1$ when $\det(U)=0$, which only occurs for reducible modules.

I am working with the symmetric group here, and in particular with Specht modules. These are not always simple, when say $p|n!$ over $\mathbb{F}_p$. The simple modules in this case are $D^\lambda = S^\lambda / (S^\lambda \cap S^{\lambda \bot})$. To start with, I am fine with just considering $p \nmid n!$, and coming back to $p|n$ later.

If one were to include $\lambda_g$, it would have to be as an extra linear variable, right? If you just compute $\lambda_g$ as I am in the code, you just get $\lambda_g=1$ in every case. Note that $GF(q)*$ is a cyclic group of order $q-1$, so it's not clear to me that the determinant homomorphism just maps to ${\pm 1}$.

Regardless, I don't think this condition is needed here. I'm not interested in projective representations.

---

The bigger issue is explaining what's going on with the current calculation, and why the checks are failing even though we are guaranteeing in the initial equations the properties we expect to see. It does depend on the characteristic, but I don't think the issue is simply about p|n! and simple modules. We can work over a slightly bigger field where $p \nmid n!$, and we still see the conditions failing sometimes:

p=2
[([3], True), ([2, 1], True), ([1, 1, 1], True)]
[([4], True), ([3, 1], False), ([2, 2], True), ([2, 1, 1], False), ([1, 1, 1, 1], True)]
[([5], True), ([4, 1], True), ([3, 2], False), ([3, 1, 1], True), ([2, 2, 1], False), ([2, 1, 1, 1], True), ([1, 1, 1, 1, 1], True)]
---------
p=3
[([3], True), ([2, 1], False), ([1, 1, 1], True)]
[([4], True), ([3, 1], False), ([2, 2], False), ([2, 1, 1], False), ([1, 1, 1, 1], True)]
[([5], True), ([4, 1], True), ([3, 2], False), ([3, 1, 1], True), ([2, 2, 1], False), ([2, 1, 1, 1], True), ([1, 1, 1, 1, 1], True)]
---------
p=5
[([3], True), ([2, 1], False), ([1, 1, 1], True)]
[([4], True), ([3, 1], True), ([2, 2], False), ([2, 1, 1], True), ([1, 1, 1, 1], True)]
[([5], True), ([4, 1], False), ([3, 2], False), ([3, 1, 1], False), ([2, 2, 1], False), ([2, 1, 1, 1], False), ([1, 1, 1, 1, 1], True)]
---------
p=7
[([3], True), ([2, 1], False), ([1, 1, 1], True)]
[([4], True), ([3, 1], False), ([2, 2], False), ([2, 1, 1], False), ([1, 1, 1, 1], True)]
[([5], True), ([4, 1], False), ([3, 2], False), ([3, 1, 1], False), ([2, 2, 1], False), ([2, 1, 1, 1], False), ([1, 1, 1, 1, 1], True)]
---------
p=11
[([3], True), ([2, 1], False), ([1, 1, 1], True)]
[([4], True), ([3, 1], False), ([2, 2], False), ([2, 1, 1], False), ([1, 1, 1, 1], True)]
[([5], True), ([4, 1], False), ([3, 2], False), ([3, 1, 1], False), ([2, 2, 1], False), ([2, 1, 1, 1], False), ([1, 1, 1, 1, 1], True)]
---------
p=13
[([3], True), ([2, 1], False), ([1, 1, 1], True)]
[([4], True), ([3, 1], False), ([2, 2], False), ([2, 1, 1], False), ([1, 1, 1, 1], True)]
[([5], True), ([4, 1], False), ([3, 2], False), ([3, 1, 1], False), ([2, 2, 1], False), ([2, 1, 1, 1], False), ([1, 1, 1, 1, 1], True)]
---------
p=17
[([3], True), ([2, 1], False), ([1, 1, 1], True)]
[([4], True), ([3, 1], False), ([2, 2], False), ([2, 1, 1], False), ([1, 1, 1, 1], True)]
[([5], True), ([4, 1], False), ([3, 2], False), ([3, 1, 1], False), ([2, 2, 1], False), ([2, 1, 1, 1], False), ([1, 1, 1, 1, 1], True)]
---------
p=19
[([3], True), ([2, 1], False), ([1, 1, 1], True)]
[([4], True), ([3, 1], False), ([2, 2], False), ([2, 1, 1], False), ([1, 1, 1, 1], True)]
[([5], True), ([4, 1], False), ([3, 2], False), ([3, 1, 1], False), ([2, 2, 1], False), ([2, 1, 1, 1], False), ([1, 1, 1, 1, 1], True)]
---------

In fact, it appears when $p \nmid n!$, i.e. $p &gt; n$ we are getting a pretty clear pattern: it's only the "extreme" partitions which are coming up true, i.e. the trivial representation and sign representation.

---

Note, there is one case (so far) where we are getting a 2 dim'l nullspace. I thought there were some theoretical reasons why this space should be 1 dim'l. It occurs when $p=2$:

dimension of space of G-invariant symmetric bilinear forms > 1
[3, 1, 1]
Vector space of degree 36 and dimension 2 over Finite Field in z2 of size 2^2
Basis matrix:
[1 1 1 1 1 0 1 1 1 1 0 1 1 1 1 0 1 1 1 1 0 1 1 1 1 0 1 1 1 1 0 1 1 1 1 1]
[0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0]

[jacksonwalters]

@dimpase Great, please feel free to email it when you've written it down.

I could be interested in those questions. I think Young's orthogonal form handles the orthogonal case for finite fields if I'm not mistaken. I don't know much about symplectic embeddings, but I assume you just mean conjugating the representations into the symplectic group over a finite field as described here:

https://mathoverflow.net/questions/230885/symplectic-group-over-integers-and-finite-fields

[dimpase]

gap> LoadPackage("forms");

forms has BaseChangeToCanonical, which according to the docs, does the base change you'd like to see.
E.g.

$ gap
 ┌───────┐   GAP 4.13.1 of 2024-06-11
 │  GAP  │   https://www.gap-system.org
 └───────┘   Architecture: amd64
 Configuration:  gmp 6.3.0, GASMAN, readline
 Loading the library and packages ...
 Packages:   Alnuth 3.2.1, AtlasRep 2.1.8, AutPGrp 1.11, CRISP 1.4.6, CTblLib 1.3.9, 
             FactInt 1.6.3, FGA 1.5.0, GAPDoc 1.6.7, IO 4.8.2, IRREDSOL 1.4.4, LAGUNA 3.9.6, 
             Polenta 1.3.10, Polycyclic 2.16, PrimGrp 3.4.4, RadiRoot 2.9, ResClasses 4.7.3, 
             SmallGrp 1.5.3, Sophus 1.27, TomLib 1.2.11, TransGrp 3.6.5, utils 0.84
 Try '??help' for help. See also '?copyright', '?cite' and '?authors'
gap> LoadPackage("forms");
---------------------------------------------------------------------
Loading 'Forms' 1.2.12 (30/08/2024)
by John Bamberg (http://school.maths.uwa.edu.au/~bamberg/)
   Jan De Beule (http://www.debeule.eu)
For help, type: ?Forms 
---------------------------------------------------------------------
true
gap> u2:=Z(3)^0*[[0,1,2],[1,0,1],[2,1,0]];
[ [ 0*Z(3), Z(3)^0, Z(3) ], [ Z(3)^0, 0*Z(3), Z(3)^0 ], [ Z(3), Z(3)^0, 0*Z(3) ] ]
gap> form := HermitianFormByMatrix( u2, GF(3^2) );
< hermitian form >
gap> b:=BaseChangeToCanonical(form);
[ [ Z(3)^0, Z(3^2), 0*Z(3) ], [ Z(3^2)^3, Z(3^2)^2, 0*Z(3) ], [ Z(3^2)^3, Z(3^2)^7, Z(3^2)^7 ] 
 ]
gap> bstar:=List(TransposedMat(b),x->List(x,y->y^3));
[ [ Z(3)^0, Z(3^2), Z(3^2) ], [ Z(3^2)^3, Z(3^2)^6, Z(3^2)^5 ], [ 0*Z(3), 0*Z(3), Z(3^2)^5 ] ]
gap> Display(b*u2*bstar);
 1 . .
 . 1 .
 . . 1
gap> 

So the (science) question is how it's implemented there, can we add anything new?

[jacksonwalters]

@dimpase Looking good, so here the matrix associated to the form is $u_2$, and we're able to obtain $bu_2b^\ast = Id$.

Then we can write $u_2 = b^{-1}b^{-1 \ast}$. This makes sense, because we're getting a decomposition but $b^{-1}$ is not lower triangular.

b^-1;
[ [ Z(3^2)^3, Z(3^2)^6, 0*Z(3) ], [ Z(3)^0, Z(3^2), 0*Z(3) ], [ Z(3^2)^5, Z(3^2)^7, Z(3^2) ] ]

This requirement was never strictly necessary, and it explains why my brute force search failed. Now we are free to write

$\tilde{\rho}(g) = b^{-1 \ast} \rho(g) b^\ast$

to obtain unitary representations of $S_n$ over $F_q$. Awesome!

As for anything new, my hope was to use these unitary rep'ns to construct a DFT which is unitary w.r.t. the standard form. It would be interesting to study those eigenvalues. The cyclic group case is easy (fourth roots of unity), but I'm not sure what happens for $S_n$. It seems like if you don't have a unitary DFT, then the eigenvalues would be "scattered", but normalizing in this way might (speculatively) clean them up.

[dimpase]

Beals doesn't say much about $\hat{f}(\rho)$ for an irreducible representation $\rho$ of $G$. These are some kinds of linear combinations of the characters of $G$, i.e. $\hat{f}(\rho)$ is constant on any conjugacy class, right?

For the symmetric group, let's look at Vershik-Okounkov approach. At least in characteristic 0 it allows very quick algorithms.

[jacksonwalters]

Beals just says to get a unitary DFT, one should use unitary representations for each $\rho$. See:

Unitary DFT for Finite Groups

$\hat{f}(\rho) = \sqrt{\frac{d_\rho}{|G|}}\sum_g f(g)\rho(g)$

and then just transform each basis vector $e_g$ to get the entire DFT. In order to get the unitary rep'ns in characteristic zero, one can use Weyl's unitary trick. See:

38455

This works for number fields and over $\mathbb{C}$. Approaching the eigenvalues over number fields, I compute the Galois group of the characteristic polynomial for $n=3$ in the preprint

The Modular DFT of the Symmetric Group

I'm happy to read more about the Vershik-Okounkov approach. I think in this case though we're already pretty much there with your method using GAP. I'll translate that this week, and see if the above formula for a unitary DFT works over finite fields. At the very least I should now be able to get representations that are unitary in the $x \mapsto x^q$ conjugation sense.

[jacksonwalters]

@dimpase @tscrim I need to automate transitioning the basis-change computation from GAP to Sage, but it appears to be working in that we are getting unitary rep'ns of the symmetric group over $F_{q^2}$ w.r.t. the standard form.

#define the representation matrix corresponding to q, partition
q = 3; la = [3,1]
SGA = SymmetricGroupAlgebra(GF(q^2),sum(la))
SM = SGA.specht_module(la)
rho = SM.representation_matrix
d_rho = SM.dimension()
G = SGA.group()

Need to automate this part, just checking $q=3$, $la=[3,1]$ for now:

F = GF(3^2)
z = F.gen()

#U_1 is the S_n-invariant symmetric bilinear form associated to the rep'n la = [3, 1], q = 3 over GF(q^2)
#U_1 = matrix(GF(3^2),[[1,2,2],[2,1,2],[2,2,1]])
U_1 = invariant_symmetric_sesquilinear_matrix(q,la)[0]

#GAP version of the base change matrix b1 associated to u1
#b1 := [[Z(3)^0, 0*Z(3), 0*Z(3)],[Z(3^2)^2,Z(3)^0, Z(3^2)], [Z(3),Z(3^2)^3,Z(3^2)^2]]
# Define the matrix using the elements from GAP
# n.b. Z(3) is the multiplicative generator of Z/3Z which is 2 == -1
b1 = matrix(F, [
    [1, 0, 0],           # First row: Z(3)^0, 0*Z(3), 0*Z(3)
    [z^2, 1, z],         # Second row: Z(3^2)^2, Z(3)^0, Z(3^2)
    [2, z^3, z^2]        # Third row: Z(3), Z(3^2)^3, Z(3^2)^2
])

#we use the matrix factorization U_1 = AA^* where A is not necessarily lower triangular
#set A = b1.inverse() for convenience
A = b1.inverse()
A_star = conjugate_pos_char(A).transpose()
U_1 == A*A_star
True

Checking the unitary condition:

#construct the unitary representations as \tilde{\rho}(g) = A^*\rho(g)A^*.inverse()
all_unitary = True
for g in G:
    tilde_rho_g = A_star*rho(g)*A_star.inverse()
    unitary_check = conjugate_pos_char(tilde_rho_g).transpose()*tilde_rho_g == identity_matrix(F,d_rho)
    all_unitary = all_unitary and unitary_check
print(all_unitary)
True

[dimpase]

Regarding #38455, all these square roots, are they really needed? (well, they are in the original formulation, but is it need to get the same or better efficiency of algorithms?) E.g. that normalisation by $1/\sqrt{n!}$ should not be carried in the computations (having it there will blow everything up for quite small $n$'s). Then, do we need the canonical Hermitian forms, or it's enough to have $\sum_{k} d_k x_k\overline{x_k}$, with all $d_k=\overline{d_k}$ ? If not, can we still carry all the computations in a small cyclotomic field?

[jacksonwalters]

@dimpase @tscrim Good question. In 38455 the square roots I take are the diagonal elements (I need $P=Q^2$, so I take a matrix square root), then just the normalization factors $\frac{d_\rho}{|G|}$. How do we avoid taking $\sqrt{n!}$ while still getting a unitary representation?

One thought I had was to note that if these representations are real (if we don't take square roots of negative numbers), then rather than doing all this work over number fields, we can just use the orthogonal forms, right? The drawback is I think you get more information here, in that you can compute, say, the splitting field of the characteristic polynomial explicitly (even if it is of degree 192 over $\mathbb{Q}$ even for $n=3$).

On a separate note, I have automated the conversion from GAP to Sage to get the base-change matrices in the finite field case:

U_1 = matrix(GF(3^2),[[1,2,2],[2,1,2],[2,2,1]])
U_2 = matrix(GF(3^2),[[0,1,2],[1,0,1],[2,1,0]])

#define conjugation as x |--> x**q, an order two automorphism of F_q^2. note x**q == x for x \in F_q.
def conjugate_pos_char(A):
    assert A.nrows() == A.ncols()
    field_size = A.base_ring().order()
    q = sqrt(field_size) if field_size.is_square() else field_size
    return matrix(A.base_ring(),[[A[i][j]**q for j in range(A.nrows())] for i in range(A.nrows())])

#define finite field F_9 we are working over
q = 3
F = GF(q^2)
z = F.gen()
w = F.subfield(1).multiplicative_generator()

import re
def convert_to_sage_gens(expression, q):
    
    # Convert Z(q) to w and Z(q^2) to z
    expression = re.sub(r'Z\((\d+)\)', str('w'), expression)  # Z(q)
    expression = re.sub(r'Z\((\d+)\^2\)', str('z'), expression)  # Z(q^2)
    expression = re.sub(r'\^', '**', expression) # Convert all occurrences of ^ to **
    
    return expression

#just use gap.eval to do gap evaluations in sequence
def base_change_matrix(U,q):
    gap.eval('LoadPackage("forms");')
    U_str = str([list(row) for row in U])
    gap.eval(f"u := Z({q})^0*{U_str};")
    gap.eval(f"form := HermitianFormByMatrix(u, GF({q}^2));")
    b_gap = gap.eval("b := BaseChangeToCanonical(form);")
    converted_expr = convert_to_sage_gens(b_gap, q)
    b = matrix(F,eval(converted_expr))
    return b

#verify base change for U_1
b1 = base_change_matrix(U_1,q)
A = b1.inverse()
A_star = conjugate_pos_char(A).transpose()
U_1 == A*A_star

#verify base change for U_2
b2 = base_change_matrix(U_2,q)
A = b2.inverse()
A_star = conjugate_pos_char(A).transpose()
U_2 == A*A_star

[dimpase]

I don't quite understand that "automatic conversion". I think we are already able to convert between GAP matrices over cyclotomics (resp. GF(q)) and Sage's matrices over cyclotomics (resp. GF(q)).

No need to do any string parsing.

[jacksonwalters]

@dimpase Thanks for the tip, here's the updated version:

#use gap.eval for GAP evalutation of BaseChangeToCanonical using `forms` package
def base_change_matrix(U,q):
    gap.eval('LoadPackage("forms");')
    U_str = str([list(row) for row in U])
    gap.eval(f"u := Z({q})^0*{U_str};")
    gap.eval(f"form := HermitianFormByMatrix(u, GF({q}^2));")
    b_str = gap.eval("b := BaseChangeToCanonical(form);")
    return matrix(F,gap(b_str))

[dimpase]

There is seldom any need to call eval.
And don't use gap. - it is deprecated. Use libgap..

E.g. loading GAP package should be done like this:

libgap.LoadPackage("forms")

[dimpase]

And we do have a direct way to pass matrices to (lib)GAP and get them back.
Please check out docs and code in e.g. src/sage/gap/matrix_gps/finitely_generated_gap.py

[dimpase]

basically you can create handles to GAP objects
like

foo=FooBar(...)

where FooBar is a GAP function, and convert them into Sage objects like foo.sage().

You can call GAP's functions Baz(arg1,...) which take foo as the 1st argument as methods: foo.Baz(...). Etc.

[dimpase]

regarding square roots: for making representations unitary they are an overkill.

Indeed, you only need to represent the $d_k$ in $\sum_k
d_k x_k\overline{x_k}$ as $d_k=t_k\overline{t_k}$, i.e. as sums of two squares, not as squares, to make them all 1.

[jacksonwalters]

Ah, didn't know gap was deprecated. Switched to libgap:

#use libgap.eval for GAP evalutation of BaseChangeToCanonical using `forms` package
def base_change_matrix(U,q):
    libgap.LoadPackage("forms")
    U_str = str([list(row) for row in U])
    libgap.eval(f"u := Z({q})^0*{U_str};")
    libgap.eval(f"form := HermitianFormByMatrix(u, GF({q}^2));")
    b_str = libgap.eval("b := BaseChangeToCanonical(form);")
    return matrix(F,libgap(b_str))

And we do have a direct way to pass matrices to (lib)GAP and get them back.
Please check out docs and code in e.g. src/sage/gap/matrix_gps/finitely_generated_gap.py

I'm looking in src/sage/ and I don't see gap/. I just see game_theory/, games/, geometry/, etc.

[dimpase]

does

u_gap=libgap([list(row) for row in U])

work?

[dimpase]

low level libgap interface is in src/sage/libs/

interfaces to GAP's group theory stuff are in
src/sage/groups/

[tscrim]

Dima meant sage/groups/matrix_gps I think

[jacksonwalters]

@dimpase Apologies for the delay. Yes, thanks that cleans things up a bit:

#use libgap.eval for GAP evalutation of BaseChangeToCanonical using `forms` package
def base_change_matrix(U,q):
    libgap.LoadPackage("forms")
    u_gap = libgap([list(row) for row in U])
    libgap.eval("u := " + str(u_gap))
    libgap.eval(f"form := HermitianFormByMatrix(u, GF({q}^2));")
    b_str = libgap.eval("b := BaseChangeToCanonical(form);")
    return matrix(F,libgap(b_str)).inverse()

Not sure if there is a way to remove the str call completely.

@tscrim Great, thank you. I also notice unitary.py in there as well.

[dimpase]

u_gap.HermitianFormByMatrix(GF(q^2))
probably works. Else, it's probably libgap.GF rather than just GF.

[jacksonwalters]

@dimpase Nice, here it is as a one-liner:

#use libgap.eval for GAP evalutation of BaseChangeToCanonical using `forms` package
def base_change_matrix(U,q):
    libgap.LoadPackage("forms")
    return matrix(GF(q^2),libgap.BaseChangeToCanonical(libgap([list(row) for row in U]).HermitianFormByMatrix(GF(q^2)))).inverse()

[jacksonwalters]

@dimpase The 1x1 matrix edge case throws an error:

libgap.LoadPackage("forms")
libgap.eval("m := [ [ 1 ] ];")
libgap.eval("f := HermitianFormByMatrix(GF(3^2), m);")

---------------------------------------------------------------------------
GAPError                                  Traceback (most recent call last)
Cell In[63], line 3
      1 libgap.LoadPackage("forms")
      2 libgap.eval("m := [ [ 1 ] ];")
----> 3 libgap.eval("f := HermitianFormByMatrix(GF(3^2), m);")

File /private/var/tmp/sage-10.4-current/local/var/lib/sage/venv-python3.12.4/lib/python3.12/site-packages/sage/libs/gap/libgap.pyx:406, in sage.libs.gap.libgap.Gap.eval (build/cythonized/sage/libs/gap/libgap.c:7707)()
    404 
    405         initialize()
--> 406         elem = make_any_gap_element(self, gap_eval(gap_command))
    407 
    408         # If the element is NULL just return None instead

File /private/var/tmp/sage-10.4-current/local/var/lib/sage/venv-python3.12.4/lib/python3.12/site-packages/sage/libs/gap/util.pyx:364, in sage.libs.gap.util.gap_eval (build/cythonized/sage/libs/gap/util.c:8232)()
    362 try:
    363     GAP_Enter()
--> 364     result = GAP_EvalString(cmd)
    365     # We can assume that the result object is a GAP PList (plain list)
    366     # and we should use functions for PLists directly for now; see

GAPError: Error, no method found! Error, no 1st choice method found for `HermitianFormByMatrix' on 2 arguments

Here is a quick fix. In the 1x1 case, we can use a generator and factor.

#handle 1x1 case by factoring u=aa^* using generator z of GF(q^2). u=z^k, a=z^l, aa^*=z^{(q+1)l} ==> k == (q+1)l (mod q^2-1)
    if U.nrows() == 1 and U.ncols() == 1:
        u = U[0, 0]
        if u == 0:
            return matrix(GF(q**2), [[0]])  # Special case for 0
        F = GF(q**2)
        z = F.multiplicative_generator()
        k = u.log(z)  # Compute discrete log of u to the base z
        if k % (q + 1) != 0:
            raise ValueError("Matrix entry is not a unitary form (k is not divisible by q+1).")
        m = k // (q+1)  # Quotient
        l = m % (q-1)  # Reduce modulo (q-1)
        a = z**l
        return matrix(F, [[a]])

[dimpase]

it's a forms GAP package bug I suppose. Report it there

[jacksonwalters]

@dimpase I'll report it there. So now we are almost done. Here is putting everything together in a DFT which should be unitary (using unitary rep'ns, using \sqrt{d_\rho/|G|}) to normalize. It's slightly off for n=3 in that I get both +1's and -1's on the diagonal:

n = 3; q = 5
SGA = SymmetricGroupAlgebra(GF(q^2), n)

#define conjugation as x |--> x**q, an order two automorphism of F_q^2. note x**q == x for x \in F_q.
def conjugate_pos_char(A):
    assert A.nrows() == A.ncols()
    field_size = A.base_ring().order()
    q = sqrt(field_size) if field_size.is_square() else field_size
    return matrix(GF(q**2),[[A[i][j]**q for j in range(A.nrows())] for i in range(A.nrows())])

def invariant_symmetric_bilinear_matrix(q,partition):
    """
    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.
    
    """

    # Define the group G and its rep'n as a Specht module, dimension
    n = sum(partition)
    SGA = SymmetricGroupAlgebra(GF(q^2), n)
    SM = SGA.specht_module(partition)
    G = SGA.group()
    rho = SM.representation_matrix
    d_rho = SM.dimension()
    
    # Initialize U as a matrix of variables over GF(q^2)
    R = PolynomialRing(GF(q^2), 'u', d_rho^2)
    U_vars = R.gens()  # List of variable generators for U
    U = Matrix(R, d_rho, d_rho, U_vars)  # U is a d_rho x d_rho matrix of variables
    
    # for each generator of G, form the augmented system 
    def augmented_matrix(g):

        #compute \rho(g), transpose, conjugate
        rho_g = rho(Permutation(g))
        rho_g_T = rho_g.transpose()
        rho_g_conj = conjugate_pos_char(rho_g)
    
        # Compute lambda_g
        det_rho_g = det(rho_g)
        lambda_g = det_rho_g * (det_rho_g ** q)
    
        # Form the matrix equation \rho(g)^T*U*\overline{\rho(g)} = \lambda_g * U
        equation_matrix = rho_g_T*U*rho_g_conj - lambda_g * U
    
        # Initialize a list to hold rows of the augmented system
        augmented_system = []
    
        # Extract coefficients for each linear equation in the matrix
        for i in range(d_rho):
            for j in range(d_rho):
                # Get the (i, j) entry of the equation matrix, which is a linear combination of the u variables
                linear_expression = equation_matrix[i, j]
            
                # Extract the coefficients of each u_k in the linear expression
                row = [linear_expression.coefficient(u) for u in U_vars]
            
                # Append the row to the augmented system
                augmented_system.append(row)
    
        # Convert the augmented system to a matrix
        return Matrix(GF(q^2), augmented_system)

    #stack linear systems for each g in G
    total_system = matrix(GF(q^2),0,d_rho^2)
    for g in G:
        total_system = total_system.stack(augmented_matrix(g))
    
    #compute the null space of the overall matrix
    null_space = total_system.right_kernel()
    
    #return a d_rho x d_rho matrix over GF(q^2) from the 1 dim'l null space given as vector
    U_mats = [matrix(GF(q^2),d_rho,d_rho,b) for b in null_space.basis()]

    #verify that a solution to the linear system satisfies the G-invariance property
    assert all(rho(g).transpose()*U_mats[0]*conjugate_pos_char(rho(g)) == U_mats[0] for g in G)
    
    return U_mats

#use libgap.eval for GAP evalutation of BaseChangeToCanonical using `forms` package
def unitary_change_of_basis(U,q):
    libgap.LoadPackage("forms")
    if U.nrows() == 1 and U.ncols() == 1:
        u = U[0, 0]
        if u == 0:
            return matrix(GF(q**2), [[0]])  # Special case for 0
        F = GF(q**2)
        z = F.multiplicative_generator()
        k = u.log(z)  # Compute discrete log of u to the base z
        if k % (q + 1) != 0:
            raise ValueError("Matrix entry is not a unitary form (k is not divisible by q+1).")
        m = k // (q+1)  # Quotient
        l = m % (q-1)  # Reduce modulo (q-1)
        a = z**l
        return matrix(F, [[a]])
    return matrix(GF(q^2),libgap.BaseChangeToCanonical(libgap([list(row) for row in U]).HermitianFormByMatrix(GF(q^2)))).inverse()

#define the Fourier coefficient at the rep'n specht_module corresponding to partition
def hat(f,partition,SGA,unitary=False):
    specht_module = SGA.specht_module(partition)
    rho = specht_module.representation_matrix
    if unitary:
        U = invariant_symmetric_bilinear_matrix(q,partition)[0]
        A = unitary_change_of_basis(U,q)
        A_star = conjugate_pos_char(A).transpose()
        unitary_factor = specht_module.dimension()/SGA.group().cardinality()
        sqrt_unitary_factor = sqrt(GF(q**2)(unitary_factor))
        return sqrt_unitary_factor*sum(f(g)*A_star*rho(g)*A_star.inverse() for g in SGA.group())
    else:
        return sum(f(g)*rho(g) for g in SGA.group())

#for each basis element g \in G compute the Fourier coefficients \hat{\delta_g}(partition) for all partitions
from sage.misc.flatten import flatten
delta = lambda s: lambda t: 1 if t == s else 0 #delta function \delta_s(t)
def dft(SGA,unitary=False):
    fourier_transform = [flatten([hat(delta(g),partition,SGA,unitary).list() for partition in Partitions(SGA.group().degree())]) for g in SGA.group()]
    if unitary:
        return matrix(GF(q**2),fourier_transform).transpose()
    else:
        return matrix(fourier_transform).transpose()

unitary_dft = dft(SGA,unitary=True); unitary_dft
[       1        1        1        1        1        1]
[4*z2 + 3 2*z2 + 4   z2 + 2 3*z2 + 1 3*z2 + 1 2*z2 + 4]
[       0 2*z2 + 2        0 3*z2 + 3 2*z2 + 2 3*z2 + 3]
[       0 2*z2 + 1        0 2*z2 + 1 3*z2 + 4 3*z2 + 4]
[4*z2 + 3 3*z2 + 1 4*z2 + 3 3*z2 + 1 3*z2 + 1 3*z2 + 1]
[       1        4        4        1        1        4]

unitary_dft*conjugate_pos_char(unitary_dft).transpose()
[1 0 0 0 0 0]
[0 4 0 0 0 0]
[0 0 4 0 0 0]
[0 0 0 4 0 0]
[0 0 0 0 4 0]
[0 0 0 0 0 1]

[dimpase]

Good question. In 38455 the square roots I take are the diagonal elements (I need P = Q 2 , so I take a matrix square root), then just the normalization factors d ρ | G | . How do we avoid taking n ! while still getting a unitary representation?

Obviously you can normalise the form by any integer square, it does not affect anything. If you venture out to $\mathbb{Z}[i]$, then you can use that about half the primes are sums of 2 squares, and so for the cost of going complex you can get rid of more non-1 diagonal elements $d_k$ by writing them as $d_k=u\overline{u}=(a+ib)(a-ib)$. Similarly, you can keep going all the way to Hamiltonian quaternions: as every positive integer is a sum of 4 squares, you can write any $d_k=u\overline{u}$, with $u=a+ib+jc+kd$ a Hamiltonian quaternion.

[jacksonwalters]

it's a forms GAP package bug I suppose. Report it there

gap-packages/forms#73

[jacksonwalters]

@dimpase @tscrim I will try to parse your comment re: the number field case.

For the finite field case, it looks like we get a mix of -1's and +1's on the diagonal, e.g. for n=4, q=11 we have

$DFT.DFT^* = diag([10,10,10,10,10,10,10,10,10,10,1, 1,1,1,10,10, 10,10, 10,10,10,10,10,10])$

[SGA.specht_module(la).dimension()**2 for la in Partitions(4)]
[1, 9, 4, 9, 1]

Perhaps the normalizing factor for each representation

sqrt_unitary_factor = sqrt(GF(q**2)(specht_module.dimension()/SGA.group().cardinality()))

is important, as there are two distinct opposite possibilities. I'm not sure which one to choose. We may even need to use a fourth root of unity.

[jacksonwalters]

@dimpase What we ended up doing is just factoring each sign on the diagonal $c_i$ as $c_i = z_i z_i^$, so if $DFT.DFT^ = S$ then we write $S = RR^*$. Then $uDFT = R^{-1}.DFT$ is unitary. This more or less closes the issue, modulo comments about ways to improve the existing code in the number field case. I am going to begin updating the PR tomorrow.

[jacksonwalters]

The PR #38455 has been updated. The techniques here for obtaining unitary representations over finite fields via $S_n$-invariant symmetric bilinear forms and Hermitian decomposition/factoring are used.

For now, we're just using the orthogonal rep'ns in the number field case since they're also unitary. What remains is to ensure we're working over a proper number field rather than the symbolic ring. That shouldn't be too difficult.

Good question. In 38455 the square roots I take are the diagonal elements (I need P = Q 2 , so I take a matrix square root), then just the normalization factors d ρ | G | . How do we avoid taking n ! while still getting a unitary representation?

Obviously you can normalise the form by any integer square, it does not affect anything. If you venture out to Z [ i ] , then you can use that about half the primes are sums of 2 squares, and so for the cost of going complex you can get rid of more non-1 diagonal elements d k by writing them as d k = u u ― = ( a + i b ) ( a − i b ) . Similarly, you can keep going all the way to Hamiltonian quaternions: as every positive integer is a sum of 4 squares, you can write any d k = u u ― , with u = a + i b + j c + k d a Hamiltonian quaternion.

The PR #38455 has been updated. The techniques here for obtaining unitary representations over $F_{q^2} $via $S_n$-invariant symmetric bilinear forms and Hermitian decomposition/factoring are used.

For now, we're just using the orthogonal rep'ns in the number field case since they're also unitary. What remains is to ensure we're working over a proper number field rather than the symbolic ring. That shouldn't be too difficult.