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?

• [X] 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: https://github.com/sagemath/sage/pull/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 $<x,y>=\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=(u0, \ldots, u{n^2})$ as a vector where $n=d_\rho$ , but when I do $\rho(g)U \otimes \lambdag 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 > 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]

[dimpase]

@dimpase I'm not sure I completely understand. We're already demanding that ρ ( g ) U ρ ( g ) ― T = U , so if that condition is satisfied then we have λ g det ( U ) = det ( U ) . We only have λ g ≠ 1 when det ( U ) = 0 , which only occurs for reducible modules.

Putting the fact that you only work with $S_n$ aside, certainly, if you require $\lambda_g=1$ then you'd miss the more general case, which is perfectly possible with other groups.

It would be good to have a general function, for all finite groups. And then you'd need $\lambda_g$.

[jacksonwalters]

@dimpase I don't mind keeping $\lambda_g$ in the code, if it can be not equal to 1 in some cases. The priority at this point is to resolve the actual issues in the code. The linear algebra appears to be fine, we do indeed get solutions $U$ which have the property $\rho(g)U\overline{\rho(g)}^{T} == U$ for all $g \in G$. It's just the way I'm checking the resulting form.

EDIT: It appears that in cases where my symbolic check is failing, $\rho(g)U\overline{\rho(g)}^T == U$ is true, but $\rho(g)^T U\overline{\rho(g)} == U$ is false.

check_form_properties(7,[2,1,1])
False

#NOTE: this is the condition we need for G-invariance of the form
rho(G[6]).transpose()*U_form*conjugate_pos_char(rho(G[6])) == U_form
False

#NOTE: this is the condition we are ensuring with the above linear algebra conditions
rho(G[6])*U_form*conjugate_pos_char(rho(G[6])).transpose() == U_form
True

So I think I actually want $\rho(g)^T U \overline{\rho(g)} == \lambda_g U$, so I've updated that in the code, and now getting verified solutions for every case.

[dimpase]

EDIT: It appears that in cases where my symbolic check is failing, $ρ ( g ) U \overline{ρ ( g )}^\top == U$ is true, but $ρ ( g )^\top U \overline{ρ ( g )} == U$ is false.

It's not a surprise, as these are different actions. One of my students had a more or less the same issue when he was implementing a procedure from my paper, where one has to sum $g^\top \overline{g}$ over a group (in char 0). So mixing it with $\overline{g}^\top g$ gave us quite a headache...

By the way, I don't see how $\lambda_g$ would help you here.

[jacksonwalters]

@dimpase

It's not a surprise, as these are different actions. One of my students had a more or less the same issue when he was implementing a procedure from my paper, where one has to sum $\overline{g}^Tg$ over a group (in char 0). So mixing it with $\overline{g}g^T$ gave us quite a headache...

Yes, in characteristic zero this is essentially Weyl's unitary trick, and it takes you from the $G$-invariant form to the matrix you actually need to do the conjugation. The argument is something like

$\langle x, y \rangle_G = \frac{1}{|G|}\sum_g \langle gx, gy \rangle = \frac{1}{|G|}\sum_g (gx)^Tgy = \frac{1}{|G|}\sum_g x^T (g^Tg)y = x^T \left(\frac{1}{|G|}\sum_g g^Tg \right) y$

So set $P = \left(\frac{1}{|G|}\sum_g g^Tg \right)$, and note this matrix is symmetric, and find a square root $P = Q^2$. Then if you conjugate $Q^{-1}\rho(g)Q$ your matrices will be unitary.

I'm just not sure what to do in this case in characteristic $p$. We now have a $G$-invariant symmetric bilinear form defined by matrix $U$, but I'm not sure how to obtain the matrix I need to conjugate them to a "unitary" representation. I guess the rep'n is already unitary with respect to the form defined by $U$.

[dimpase]

Once you have $U$, you can compute the Cholesky decomposition ($U=L\overline{L}^\top$, with $L$ lower-triangular), or LDLT decomposition ($U=LD\overline{L}^\top$, with $L_{ii}=1$ for all $i$, $L$ lower-triangular, and $D$ a diagonal matrix with entries in the index 2 subfield) of it, which is a purely algebraic, exact, process, and works in any characteristic, as far as I can see.

Cholesky involves taking square roots in the index 2 subfield - but these are always available, as they exist is the quadratic extension (in the case of odd characteristic, or in the field itself in the case of characteristic 2). LDLT doesn't need this, but leaves you with a $D$ for the matrix for the form, rather than $I$. However, $D$ can still be factored by taking square roots, and it works by the same argument.

[jacksonwalters]

@dimpase Sounds like a good idea, however I'm getting some errors when attempting to compute the decomposition, likely because it's over a finite field. From the documentation:

INPUT:

A positive-definite matrix. Generally, the base ring for the entries of the matrix needs to be a subfield of the algebraic numbers (QQbar). Examples include the rational numbers (QQ), some number fields, and real algebraic numbers and the algebraic numbers themselves. Symbolic matrices can also occasionally be factored.

U_mat = invariant_symmetric_sesquilinear_matrix(3,[2,1,1])[0]; U_mat
[0 1 2]
[1 0 1]
[2 1 0]

U_mat.is_hermitian()
True

#compute the Cholesky decomposition of the matrix U associated to the bilinear form
U_mat.cholesky()

---------------------------------------------------------------------------
ZeroDivisionError                         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/matrix/matrix2.pyx:14526, in sage.matrix.matrix2.Matrix._block_ldlt (build/cythonized/sage/matrix/matrix2.c:111667)()
  14525 d.append(one_by_one_space(A_kk))
> 14526 _block_ldlt_pivot1x1(A, k)
  14527 k += 1

File /private/var/tmp/sage-10.4-current/local/var/lib/sage/venv-python3.12.4/lib/python3.12/site-packages/sage/matrix/matrix2.pyx:19039, in sage.matrix.matrix2._block_ldlt_pivot1x1 (build/cythonized/sage/matrix/matrix2.c:145160)()
  19038              (A.get_unsafe(k+1+i, k+1+j) -
> 19039               A.get_unsafe(k+1+i, k)*A.get_unsafe(k, k+1+j)/pivot))
  19040 A.set_unsafe(k+1+j,

File /private/var/tmp/sage-10.4-current/local/var/lib/sage/venv-python3.12.4/lib/python3.12/site-packages/sage/structure/element.pyx:1734, in sage.structure.element.Element.__truediv__ (build/cythonized/sage/structure/element.c:21248)()
   1733 if HAVE_SAME_PARENT(cl):
-> 1734     return (<Element>left)._div_(right)
   1735 if BOTH_ARE_ELEMENT(cl):

File /private/var/tmp/sage-10.4-current/local/var/lib/sage/venv-python3.12.4/lib/python3.12/site-packages/sage/rings/finite_rings/element_givaro.pyx:1132, in sage.rings.finite_rings.element_givaro.FiniteField_givaroElement._div_ (build/cythonized/sage/rings/finite_rings/element_givaro.cpp:18029)()
   1131 if (<FiniteField_givaroElement>right).element == 0:
-> 1132     raise ZeroDivisionError('division by zero in finite field')
   1133 self._cache.objectptr.div(r, self.element,

ZeroDivisionError: division by zero in finite field

During handling of the above exception, another exception occurred:

ValueError                                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/matrix/matrix2.pyx:13297, in sage.matrix.matrix2.Matrix.cholesky (build/cythonized/sage/matrix/matrix2.c:105779)()
  13296 try:
> 13297     _, L, d = self._block_ldlt(True)
  13298 except ValueError:

File /private/var/tmp/sage-10.4-current/local/var/lib/sage/venv-python3.12.4/lib/python3.12/site-packages/sage/matrix/matrix2.pyx:14530, in sage.matrix.matrix2.Matrix._block_ldlt (build/cythonized/sage/matrix/matrix2.c:111727)()
  14529 except ZeroDivisionError:
> 14530     raise ValueError("matrix has no classical LDL^T factorization")
  14531 

ValueError: matrix has no classical LDL^T factorization

During handling of the above exception, another exception occurred:

ValueError                                Traceback (most recent call last)
Cell In[36], line 2
      1 #compute the Cholesky decomposition of the matrix U associated to the bilinear form
----> 2 U_mat.cholesky()

File /private/var/tmp/sage-10.4-current/local/var/lib/sage/venv-python3.12.4/lib/python3.12/site-packages/sage/matrix/matrix2.pyx:13301, in sage.matrix.matrix2.Matrix.cholesky (build/cythonized/sage/matrix/matrix2.c:105861)()
  13299     # If the matrix was positive-definite, that would
  13300     # have worked.
> 13301     raise ValueError("matrix is not positive definite")
  13302 
  13303 F = L.base_ring()  # field really

ValueError: matrix is not positive definite

U_mat.block_ldlt()
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[37], line 1
----> 1 U_mat.block_ldlt()

File /private/var/tmp/sage-10.4-current/local/var/lib/sage/venv-python3.12.4/lib/python3.12/site-packages/sage/matrix/matrix2.pyx:14981, in sage.matrix.matrix2.Matrix.block_ldlt (build/cythonized/sage/matrix/matrix2.c:113251)()
  14979 cdef Matrix P, L, D   # output matrices
  14980 
> 14981 p, L, d = self._block_ldlt(classical)
  14982 MS = L.matrix_space()
  14983 P = MS.matrix(lambda i, j: p[j] == i)

File /private/var/tmp/sage-10.4-current/local/var/lib/sage/venv-python3.12.4/lib/python3.12/site-packages/sage/matrix/matrix2.pyx:14542, in sage.matrix.matrix2.Matrix._block_ldlt (build/cythonized/sage/matrix/matrix2.c:111812)()
  14540 omega_1 = 0
  14541 for i in range(k+1, n):
> 14542     a_ik_abs = A.get_unsafe(i, k).abs()
  14543     if a_ik_abs > omega_1:
  14544         omega_1 = a_ik_abs

File /private/var/tmp/sage-10.4-current/local/var/lib/sage/venv-python3.12.4/lib/python3.12/site-packages/sage/structure/element.pyx:2898, in sage.structure.element.RingElement.abs (build/cythonized/sage/structure/element.c:30371)()
   2896         ArithmeticError: absolute value not defined on integers modulo n.
   2897     """
-> 2898     return abs(self)
   2899 
   2900 def is_prime(self):

TypeError: bad operand type for abs(): 'sage.rings.finite_rings.element_givaro.FiniteField_givaroElement'

[dimpase]

Check that it works and fix it! I never claimed that Sage's Cholesky works over GF(q), perhaps you just need to get a better version of abs() - which does not make sense over GF(q), but you don't need it.

Or just implement one yourself, it is an easy algorithm.

[jacksonwalters]

@dimpase I went ahead and implemented my own algorithm, but learned that it fails in certain cases, too. Namely, whenever the updated diagonal element is zero. It turns out Sage's Cholesky does work over GF(q^2), just not in all cases:

def cholesky_finite_field(A):
    """
    Perform Cholesky decomposition of a symmetric matrix A over a finite field,
    assuming A is decomposable and the field allows square roots of all required values.

    Parameters:
        A (Matrix): A symmetric matrix over a finite field (e.g., GF(q^2) where all elements have square roots).

    Returns:
        L (Matrix): A lower triangular matrix such that A = L * L.transpose(), if decomposition is possible.
    """
    n = A.nrows()
    L = Matrix(A.base_ring(), n, n, 0)  # Initialize an n x n matrix of zeros in the same field as A

    for i in range(n):
        # Calculate the diagonal element L[i, i]
        diag_sum = sum(L[i, k]^2 for k in range(i))
        diag_val = A[i, i] - diag_sum

        # Check if diag_val has a square root in the field
        if not diag_val.is_square():
            raise ValueError(f"No square root exists for the value at diagonal {i}: {diag_val}")

        L[i, i] = diag_val.sqrt()

        # Calculate the off-diagonal elements L[j, i] for j > i
        for j in range(i + 1, n):
            off_diag_sum = sum(L[j, k] * L[i, k] for k in range(i))
            L[j, i] = (A[j, i] - off_diag_sum) / L[i, i]

    return L

Now we have both working:

#the G-invariant, symmetric, bilinear form is not the standard form. likely a different basis
U_mat = invariant_symmetric_sesquilinear_matrix(11,[3,1,1])[0]; U_mat
[1 4 7 4 7 0]
[4 1 4 4 0 7]
[7 4 1 0 4 7]
[4 4 0 1 4 4]
[7 0 4 4 1 4]
[0 7 7 4 4 1]

#compute the Cholesky decomposition of the matrix U associated to the bilinear form
try:
    print("matrix is Hermition:",U_mat.is_hermitian())
    print(U_mat.cholesky())
except ValueError as e:
    print(e)

matrix is Hermition: True
[        1         0         0         0         0         0]
[        4  3*z2 + 5         0         0         0         0]
[        7  7*z2 + 8  9*z2 + 4         0         0         0]
[        4  9*z2 + 4         0  7*z2 + 8         0         0]
[        7 10*z2 + 2  5*z2 + 1  9*z2 + 4 10*z2 + 2         0]
[        0  3*z2 + 5 6*z2 + 10  9*z2 + 4  7*z2 + 8         5]

cholesky_finite_field(U_mat)
[        1         0         0         0         0         0]
[        4  3*z2 + 5         0         0         0         0]
[        7  7*z2 + 8  9*z2 + 4         0         0         0]
[        4  9*z2 + 4         0  7*z2 + 8         0         0]
[        7 10*z2 + 2  5*z2 + 1  9*z2 + 4 10*z2 + 2         0]
[        0  3*z2 + 5 6*z2 + 10  9*z2 + 4  7*z2 + 8         5]

Over GF(3^2), the computation fails:

#the G-invariant, symmetric, bilinear form is not the standard form. likely a different basis
U_mat = invariant_symmetric_sesquilinear_matrix(3,[3,1,1])[0]; U_mat
[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]

#compute the Cholesky decomposition of the matrix U associated to the bilinear form
try:
    print("matrix is Hermition:",U_mat.is_hermitian())
    print(U_mat.cholesky())
except ValueError as e:
    print(e)
matrix is Hermition: True
matrix is not positive definite

try:
    print(cholesky_finite_field(U_mat))
except ZeroDivisionError as e:
    print(e)
division by zero in finite field

[dimpase]

ah, OK. You need to "do a pivot", i.e. conjugate with a permutation matrix to get a non-zero diagonal element. The resulting decomposition will be less pretty, but that's OK. Probably in Sage code, rather than using abs(), you can have your own version of abs(), sending 0 to 0 and everything else to 1.

[jacksonwalters]

@dimpase Is there a canonical choice of pivot matrix?

Do you mean in cholesky() in SageMath changing abs()?

In any case, it seems for large enough $q$ we can do the Cholesky decomposition for every partition. $p > n$ seems to work for the cases I tried ($n=3,4,5$, so $p=5,7$). Perhaps we can leave the case $p <= n$ until a bit later.

Once we have the Cholesky decomposition $U = L L^\ast$, what are we doing with it? Is this the change-of-basis matrix, so that our unitary matrices are $L \rho(g) L^\ast$?

[dimpase]

Diagonalising a Hermitean form $x^\top U\overline{x}$ means applying a linear transformation to the variables: $x^\top U\overline{x}=(Lx)^\top\overline{Lx}$. This is how to understand what's going on, and should suffice to figure out how $g$ leaving $U$ should be transformed by $L$.

[dimpase]

As long as $\det U\neq 0$, there is no 0 diagonal element to pivot on, in characteristic 0 (this is from positive definiteness). https://en.m.wikipedia.org/Cholesky_decomposition (see Cholesky in outer form there.)

In general one can conjugate by a permutation matrix to get nonzero diagonal element in the right place, as long as these are available. In positive characteristic one can reach a stage where $A^{(i)}$ has all 0s in the lower right block.

E.g. $f(x,y)=y\overline{\alpha x}+\alpha x\overline{y}$ is a Hermitean form, with nonzero eigenvalues of the matrix. Diagonalisible in odd characteristic looks possible, but not possible (?) in even characteristic.

[jacksonwalters]

Diagonalising a Hermitan form $x^TU\overline{x}$ means applying a linear transformation to the variables

$x^TU\overline{x} = (Lx)^T \overline{L}\overline{x}$

This is how to understand what's going on, and should suffice to figure out how $g$ leaving $U$ should be transformed by $L$.

On the RHS, this becomes $x^\top L^\top\overline{L} \overline{x}$, but $U = L\overline{L}^\top \ne L^\top\overline{L}$.

Perhaps if we apply $L^\top$ to $x$ then we get $(L^\top x)^\top\overline{L}^\top\overline{x} = x^\top L L^*\overline{x} = x^\top U \overline{x}$.

Note the SageMath built-in Cholesky is returning $L$ such that $U=LL^T$. I'm trying to modify mine to return the appropriate $U=LL^*$.

This paper touches on positive-definiteness: https://arxiv.org/pdf/2202.04012

[jacksonwalters]

For a unitary change-of-basis, we must have a matrix A s.t.:

(A\rho(g)A.inverse())^* A \rho(g)A.inverse() = A.inverse()^* \rho(g)^*A^*A\rho(g)A.inverse()

then we must have A^*A = U = LL^* ==> A = L^*

= (L^*).inverse()^* \rho(g)^* L^{**}L^* \rho(g) L^*.inverse() #substitution
= L.inverse() \rho(g)^T LL^* \rho(g) L^*.inverse() #assuming \rho is defined over the base field F_q fixed by x |--> x^q
= L.inverse() \rho(g)^T U \rho(g) L^*.inverse() 
= L.inverse() U L^*.inverse() #by G-invariance
= Id

Thus \tilde(\rho}(g) = L^* \rho(g) L^*.inverse() should be a unitary matrix. What remains is to find the correct Cholesky decomposition of U over F_{q^2}.

[dimpase]

one of concrete examples of forms for which you cannot compute a base change to get the canonical form $\sum_k x_k\overline{x_k}$ using your procedure is in https://github.com/sagemath/sage/issues/38456#issuecomment-2447868158, right?

(sorry for a slow response, we were on the move UK->USA (Northwestern U.), still lots of move-related things to do...)

[jacksonwalters]

In that comment that was before I realized the code only does the decomposition so that $U = LL^T$. It does not give $U = LL^*$.

At this point I’m interested in whether or not the Cholesky decomposition exists in general. The above proof shows if so, then you can basis-change to get unitary rep’ns with respect to the standard form. If not, it’s unclear what happens.

To that end, I think it’s worth just brute forcing finding $L$ for small $q$, $n$. That’ll give some valuable information, as I don’t see how to adapt the SageMath algorithm to get factorization with $x \mapsto x^q$ conjugation.

No worries! Hope the move went well. Quite a trip.

[dimpase]

There is up to a linear, i.e. conjunction in $L(n,p^{2k})$ (?) isomorphism only one $U(n,p^k)$ for all $n$ and primes $p$. See e.g. the Altas of Finite groups by Conway et al. (if you need a e-copy, let me know).

This would imply the existence of a base change.

[jacksonwalters]

@dimpase Interesting. Funny story - I checked that book out of the library for a presentation on the Monster group to show people the character table, and it is absolutely massive. Yes, if you could send to jacksonwalters@gmail.com that would be great. Thanks a ton!

[dimpase]

@dimpase Interesting. Funny story - I checked that book out of the library for a presentation on the Monster group to show people the character table, and it is absolutely massive. Yes, if you could send to jacksonwalters@gmail.com that would be great. Thanks a ton!

Shared on googl drive. (among other books, some relevant maybe)

[dimpase]

brute forcing finding L for small q , n .

here is a complete algorithm for $n=2$. Let $H=\begin{pmatrix} \alpha & \omega \\ \overline{\omega} & \beta \end{pmatrix}$, with $\alpha=\overline{\alpha}$, $\beta=\overline{\beta}$, look for $L=\begin{pmatrix} 1 & 0 \\ a & 1 \end{pmatrix}$ so that $LHL^*$ is diagonal. If $\omega=0$ we are done ($H$ is diagonal, and we can factor $\alpha=\tau\overline{\tau}$, $\beta=\mu\overline{\mu}$ to get $H=I$.)

We have $LHL^*=\begin{pmatrix} \alpha & \alpha\overline{a}+\omega \\ a\alpha + \overline{\omega} & a\overline{a}\alpha+a\omega+\overline{a\omega}+\beta \end{pmatrix}.$

If $\alpha\neq 0$ we can set $a:=-\overline{\omega}/\alpha$ and get the result diagonal. If $\beta\neq 0$, we can pre-conjugate $H$ by the permutation matrix $\Pi$ of $(1,2)$ and get nonzero $\alpha$, and proceed as above. So in this case the conjugation is by $L\Pi $.

If $\alpha=\beta=0$, in $LHL^*$ we get new $\beta:=a\omega +\overline{a\omega}$, so we apply conjugation by $\Pi$ to get nonzero $\alpha$, and proceed as above to find a good $L'$. So in this case the conjugation is by $L'\Pi L$.

[jacksonwalters]

Ah, okay, that's fantastic. So that could be coded up. I did do a brief brute force search for a couple of matrices $U$ in the $n=3$ case and they turned up empty, suggesting this is simply not possible in general. The brute force search does work for $n=2$ in the case I tried ($U_0$). Please double check:

#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 cholesky_brute_force(U,q):
    """
        Compute the Cholesky decomposition U = LL^* by brute force, where * denotes conjugate transpose and we are using the conjugation x |--> x^q.

        This will only work for small values of q, n. We can look at representations of dimension d_rho = 2, 3, so that d_rho^2 = 4, 9. 

        We work over GF(q^2) where we have square roots, for q=2, 3. If n=3,4, then p|n is problematic, but we can try q=5. 

        For example, for q=3, q^2=9, n=4, d_rho = 3, d_rho^2 = 9, we would have 9^9 = 387420489 possibilities.
    """
    from itertools import product

    # Define the finite field GF(3^2)
    F = GF(q^2)
    elements = list(F) # List of all elements in GF(3^2)

    # Define the size of the matrix
    n = U.nrows()
    num_entries = n + (n * (n - 1)) // 2
    num_iters = 0

    # Iterate over all possible combinations of elements
    for a in product(elements, repeat=num_entries):
        num_iters += 1
        # Construct the lower triangular matrix
        if n == 2:
            L = matrix(F,[[a[0],0],[a[1],a[2]]])
        if n ==3:
            L = matrix(F, [[a[0], 0, 0],[a[1], a[2], 0],[a[3], a[4], a[5]]])

        if U == L * conjugate_pos_char(L).transpose():
            return L

        if num_iters % 10^4 == 0:
            print(float(num_iters / (q^2)^6))

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

cholesky_brute_force(U_0,q)
[  z2 + 1        0]
[2*z2 + 2        0]

[dimpase]

I'm sure it can be generalized for all n. Maybe it's even publishable :-)