unitary DFT for symmetric group algebra

[jacksonwalters]

📝 Checklist

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

⌛ Dependencies

• add hermitian_decomposition method #39200: Uses the decomposition.

[jacksonwalters]

#38456

[github-actions]

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

[jacksonwalters]

To handle extensions of \Q, we can work over a number field containing all necessary square roots. Currently there is some redundancy. I compute the diagonalizations to know which roots to take, then compute them again when computing Fourier coefficients.

[jacksonwalters]

getting environment error "pytest is not installed in the venv, skip checking tests that rely on it
Error: Process completed with exit code 1." added ~3 lines to handle case when diagonal contains higher degree algebraic numbers (n=5 and above).

[jacksonwalters]

relevant MSE questions:

• Weyl's unitary trick over number fields
• Roots of a polynomial over $\mathbb{Q}(\sqrt{2},\sqrt{3})$

[dimpase]

I checked that this builds if merged over the latest beta. Could you fix "This branch is out-of-date with the base branch" ?
The interactive rebase feature as provided by GitHub is reasonably easy to understand and work with.

[tscrim]

That isn’t necessary until it is fully reviewed and only as a final check with the most recent changes.

[dimpase]

These rebases unbreak the CI - so that's the main reason to do it.

[jacksonwalters]

It looks to me like the CI tests are passing (one skipped) so far.

[jacksonwalters]

Work in issue 38456 is sufficient to construct this DFT. We compute $G$-invariant symmetric bilinear forms, factor the associated matrix $U=AA^\ast$, and obtain unitary representations $\tilde{\rho}(g)=A^\ast\rho(g)(A^{\ast})^{-1}$. The Fourier coefficients are $\hat{f}(\rho) = \sqrt{\frac{d_\rho}{|G|}}\sum_g f(g)\tilde{\rho}(g)$. However, $DFT.DFT^\ast = S$ where $S$ is a diagonal matrix consisting of signs, +1's and -1's. Factoring $S=RR^*$ with diagonal $R$, $uDFT = R^{-1}.DFT$ is unitary. Note $c=zz^\ast$ has a unique solution $z \in GF(q^2)$ when $c \in GF(q)$.

[tscrim]

So the Forms package is a GAP package that has been included since GAP v4.9, but it is not included by our default installation nor in gap_packages.

However, I am thinking we should just reproduce their algorithm directly to take advantage of Sage's linear algebra packages. The only function really used is the BaseChangeToCanonical, and from looking at the source code, it seems like they are just running Gram-Schmidt by Gaussian elimination.

What do you think?

[jacksonwalters]

So the Forms package is a GAP package that has been included since GAP v4.9, but it is not included by our default installation nor in gap_packages.

However, I am thinking we should just reproduce their algorithm directly to take advantage of Sage's linear algebra packages. The only function really used is the BaseChangeToCanonical, and from looking at the source code, it seems like they are just running Gram-Schmidt by Gaussian elimination.

What do you think?

That's unfortunate that it's not included by default. I don't know what the barriers are to getting it incorporated, or how long that would take. It would certainly make things easy.

I agree. We'd been talking about doing this anyways, and I'm happy to go ahead and just rewrite it in Sage. I'll try to get a basic version working tomorrow. Then there will probably be ways to optimize it.

[tscrim]

I believe it should just be getting the echelon form of the appropriately augmented Gram matrix of the initial sesquilinear form and then taking the appropriate submatrix (possibly transposed). So it should just be a few lines of code (at most).

[jacksonwalters]

I believe it should just be getting the echelon form of the appropriately augmented Gram matrix of the initial sesquilinear form and then taking the appropriate submatrix (possibly transposed). So it should just be a few lines of code (at most).

I agree it should a lot simpler than what is written in GAP. However, a straightforward augmentation and row reduction doesn't yield the correct matrix. For example, when q=11 and la=[2,1,1] with n=4, we have

q = 11
F = GF(q**2)
U=matrix(F,[[1,4,7],[4,1,4],[7,4,1]])
#use libgap.eval for GAP evalutation of BaseChangeToCanonical using `forms` package
def unitary_change_of_basis(U,q):
    if U.nrows() == 1 and U.ncols() == 1:
        return matrix(F,[[factor_scalar(U[0,0])]])
    loaded_forms = libgap.LoadPackage("forms")
    return matrix(F,libgap.BaseChangeToCanonical(libgap([list(row) for row in U]).HermitianFormByMatrix(F))).inverse()

unitary_change_of_basis(U,q)
[        1         0         0]
[        4  9*z2 + 2         0]
[        7 10*z2 + 1  3*z2 + 3]

which is correct. However,

def base_change_hermitian(U):
    F = U.base_ring()
    augmented_mat = U.augment(identity_matrix(F, U.nrows()))
    echelon_form = augmented_mat.echelon_form()
    return echelon_form[:, U.ncols():]

base_change_hermitian(U)
[7 2 9]
[2 7 2]
[9 2 7]

which is just $U^{-1}$. I'm not sure how to properly augment the matrix. Note also in the linked source code there is the block

if not IsOne(a) then
   # find an b element with norm b*b^t = b^(t+1) equal to a
   b := RootFFE(gf, a, t+1);
    Forms_MultRow(D,i,1/b);
fi; 

which suggests at some point we should be factoring scalars as $a=bb^\ast$.

[tscrim]

Okay, a little bit more care is needed here since it does the RREF. We only want to operate on the strictly lower (or upper) triangle of the matrix as the other part is taken care of by the symmetry. So a short version is to use the inverse of the upper part of the LU decomposition:

sage: U = matrix(F,[[1,4,7],[4,1,4],[7,4,1]])
sage: Up = U.LU()[2]
sage: D = Up.diagonal()
sage: A = ~Up * matrix.diagonal([d.sqrt() for d in D])
sage: A.H * U * A
[ 1  0  0]
[ 0 10  0]
[ 0  0 10]

Of course, this is not very efficient since it is computing far more things than necessary. However, it is easy to code...

[dimpase]

Mathematically it's not echelon form, and not even Gram-Schmid, as you have a twist by the Galois automorphism in the factorisation.

[jacksonwalters]

Mathematically it's not echelon form, and not even Gram-Schmidt, as you have a twist by the Galois automorphism in the factorisation.

For what it's worth, here is a direct translation of the GAP code into Sage which appears to be working now. Perhaps we can reduce it, but I agree with Dima that the twist makes it something slightly different.

#for u in GF(q), we can factor as u=aa^* using gen. z and modular arithmetic
def conj_square_root(u):
    if u == 0:
        return 0  # Special case for 0
    z = F.multiplicative_generator()
    k = u.log(z)  # Compute discrete log of u to the base z
    if k % (q+1) != 0:
        raise ValueError("Unable to factor: u is not in base field GF(q)")
    return z**((k//(q+1))%(q-1))

def base_change_matrix(mat):
    """
    Diagonalizes a Hermitian matrix over a finite field.
    Returns the base change matrix and the rank of the Hermitian form.
    
    Arguments:
        mat: The Gram matrix of a Hermitian form (Sage matrix object)
        gf: The finite field (GF(q))
    
    Returns:
        D: The base change matrix
        r: The number of non-zero rows in D*mat*D^T
    """

    gf = mat.base_ring()
    n = mat.nrows()
    q = gf.order()
    t = sqrt(q)

    A = copy(mat)
    D = identity_matrix(gf, n)
    row = 0

    # Diagonalize A
    while True:
        row += 1

        # Look for a non-zero element on the main diagonal, starting from `row`
        i = row - 1  # Adjust for zero-based indexing in Sage
        while i < n and A[i, i].is_zero():
            i += 1

        if i == row - 1:
            # Do nothing since A[row, row] != 0
            pass
        elif i < n:
            # Swap to ensure A[row, row] != 0
            A.swap_rows(row - 1, i)
            A.swap_columns(row - 1, i)
            D.swap_rows(row - 1, i)
        else:
            # All entries on the main diagonal are zero; look for an off-diagonal element
            i = row - 1
            while i < n - 1:
                k = i + 1
                while k < n and A[i, k].is_zero():
                    k += 1
                if k == n:
                    i += 1
                else:
                    break

            if i == n - 1:
                # All elements are zero; terminate
                row -= 1
                r = row
                break

            # Fetch the non-zero element and place it at A[row, row + 1]
            if i != row - 1:
                A.swap_rows(row - 1, i)
                A.swap_columns(row - 1, i)
                D.swap_rows(row - 1, i)

            A.swap_rows(row, k)
            A.swap_columns(row, k)
            D.swap_rows(row, k)

            b = A[row, row - 1]**(-1)
            A.add_multiple_of_row(row - 1, row, -b**t)
            A.add_multiple_of_column(row, row - 1, -b)
            D.add_multiple_of_row(row - 1, row, -b)

        # Eliminate below-diagonal entries in the current column
        a = -A[row - 1, row - 1]**(-1)
        for i in range(row, n):
            b = A[i, row - 1] * a
            if not b.is_zero():
                A.add_multiple_of_column(i,row - 1, b**t)
                A.add_multiple_of_row(i, row - 1, b)
                D.add_multiple_of_row(i, row - 1, b)

        if row == n - 1:
            break

    # Count how many variables have been used
    if row == n - 1:
        if not A[n - 1, n - 1].is_zero():
            r = n
        else:
            r = n - 1

    # Normalize diagonal elements to 1
    for i in range(r):
        a = A[i, i]
        if not a.is_one():
            # Find an element `b` such that `b*b^t = b^(t+1) = a`
            b = conj_square_root(a)
            D.rescale_row(i, 1 / b)

    return D

[jacksonwalters]

The doctests are currently failing. I'm using output from my laptop, but it doesn't match what's in Sage. For example,

 [       1        0        0]
 [       3 5*z2 + 4        0]
 [       4 3*z2 + 1       z2]

is a factorization of the matrix

[1 3 4]
[3 1 3]
[4 3 1]

when it should be factorization of the matrix associated to the invariant symmetric bilinear form

[1 4 4]
[4 1 4]
[4 4 1]

I don't see any major differences between the code.

EDIT: Confirmed, the Sage code is producing U = [[1 3 4],[3 1 3],[4 3 1]] while a Jupyter notebook is producing U = [[1 4 4],[4 1 4],[4 4 1]] with the exact same code.

Note that [[1 3 4],[3 1 3],[4 3 1]] does not correspond to a S_n-invariant symmetric bilinear form, so there must be something wrong with the way Sage is computing these matrices, despite the code being the same.

[jacksonwalters]

@tscrim The only real difference is how the Specht modules are initiated. In the Sage code, we don't use SGA. It's a little concerning that setting that base_ring doesn't seem to change it:

rho = SymmetricGroupRepresentation(la,ring=GF(7**2))._representation_matrix_uncached; rho(Permutation([2,1,3,4])).base_ring()
Rational Field

It won't even run if you swap this with SGA.specht_module in the notebook version. This corresponds to line 1036 of symmetric_group_representations.py in the Sage code:

self._specht = SymmetricGroupRepresentation(partition, 'specht', ring=parent._ring)

Setting ring=parent._ring will not have the correct behavior, given the above. I'm not even sure why it doesn't just throw the same error when trying to compute the invariant symmetric bilinear form. It throws an error about being over different rings in the notebook.

[tscrim]

That is a bit concerning, but since it is defined over Z, then conputationally it should be the same. So that shouldnt matter.

The invariance of the bilinear form needs to checked against the corresponding representation (and basis).

It should work for whichever set of representation matrices you use (which usually give different matrices). The existence of the unitary form is definitely an isomorphism invariant.

[jacksonwalters]

That is a bit concerning, but since it is defined over Z, then conputationally it should be the same. So that shouldnt matter.

The invariance of the bilinear form needs to checked against the corresponding representation (and basis).

It should work for whichever set of representation matrices you use (which usually give different matrices). The existence of the unitary form is definitely an isomorphism invariant.

It should be the same, but it's the only difference between what I have in the notebook version. The Sage version is definitely computing a different U, so the linear system must be different. I'll take a look at it.

I checked the two different forms, and only [[1 4 4],[4 1 4],[4 4 1]] is invariant while [[1 3 4],[3 1 3],[4 3 1]] fails. Here's the code (of course one must manually set U if it's not being computed above).

https://github.com/jacksonwalters/dft-finite-groups/blob/main/symmetric_group/unitary_repns_over_finite_fields.ipynb

#ensure the resulting form is G-invariant, symmetric, bilinear by symbolic verification
def check_form_properties(q,partition):
    #define the representation matrix corresponding to q, partition
    F = GF(q**2)
    SGA = SymmetricGroupAlgebra(F,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(F, d_rho, var_array='x,y')
    x = matrix([R_xy.gens()[2*i] for i in range(d_rho)]).transpose()
    y = matrix([R_xy.gens()[2*i+1] for i in range(d_rho)]).transpose()
    R_xy_lambda = PolynomialRing(R_xy,'lambda')
    lambda_ = R_xy_lambda.gens()[0]

    #compute the bilinear form matrix. coerce over polynomial ring
    U_mats = invariant_symmetric_bilinear_matrix(SGA,partition)
    if len(U_mats) > 1:
        print("Space of G-invariant symmetric bilinear forms has dimension > 1 for la=",partition)
        print("Dimension of space=",len(U_mats))
    U_mat = U_mats[0]
    U_form = matrix(R_xy_lambda,U_mat)
    
    #check symmetric property
    symmetric = bilinear_form(x,y,U_form) == bilinear_form(y,x,U_form)
    
    #check G-invariance property
    G_invariant = all(bilinear_form(rho(g)*x,rho(g)*y,U_form) == bilinear_form(x,y,U_form) for g in G)
    
    #check bilinear property. ISSUE: lambda_^q is a power of the ring generator, i.e. doesn't simplify.
    first_arg = bilinear_form(lambda_*x,y,U_form) == lambda_*bilinear_form(x,y,U_form)
    second_arg = bilinear_form(x,lambda_*y,U_form) == lambda_*bilinear_form(x,y,U_form) #need to amend for conjugation
    bilinear = first_arg and second_arg

    return symmetric and G_invariant and bilinear

If different matrices give a different form, that's fine, but isn't the space usually one dim'l? So whatever set of matrices is used to generate the linear system, the result should be the linear span of one matrix. In the above examples, they're not scalar multiples of each other.

UPDATE: It appears that using rho = G.algebra(F).specht_module(self._partition).representation_matrix gives the right invariant symmetric bilinear form. Whatever's going on, that isolates the issue.

[tscrim]

That's not how linear algebra works. I just need to change my basis for the rerpreaentation vector space. This changes all of the matrices for the $S_n$ elements too, but I get a completely different matrix representing the same bilinear form.

The way the two implementations construct the representations might be slightly different (i.e., have different matrices). However, I think the problem comes from what you noted, the matrix's base ring is wrong, which makes the conjugate() wrong. That's my guess without testing things.

[jacksonwalters]

The problem I was seeing was that _unitary_change_basis_matrix is outputting [[1 3 4],[3 1 3],[4 3 1]] rather than [[1 4 4],[4 1 4],[4 4 1]]. Yes, one can conjugate to change bases, but when we solve for the invariant symmetric bilinear form, there is a system of linear equations that needs to be solved. That solution space has a dimension, which is usually one. If, for example, the matrix U=[[1 4 4],[4 1 4],[4 4 1]] is in null_space.basis(), then any other solution should be a scalar multiple of it. [[1 3 4],[3 1 3],[4 3 1]] is not a scalar multiple of [[1 4 4],[4 1 4],[4 4 1]], and checking shows it's not invariant. Perhaps I'm missing something.

In any case, you're probably right that it's the conjugation. It's definitely to do with the line

self._specht = SymmetricGroupRepresentation(partition, 'specht', ring=parent._ring)

which I am going to try replacing with

self._specht = Permutations(self._n).algebra(parent._ring).specht_module(self._partition)

temporarily just to see if I can get tests passing, and match what's happening with the code in my notebook.

[jacksonwalters]

Okay, so all tests (other than one of the flaky Meson / Conda ones) are passing.

I would like to do a separate PR addressing the base_ring issues for symmetric_group_representations.py. It'd be nice to have the representations defined over a minimal number field, and address the above concern.

[tscrim]

Moving the conjugate-transpose part into the cached portion.

[tscrim]

Suggested change

	A = self._unitary_change_basis_matrix
	return A.H * rho(permutation) * A.H.inverse()
	A = self._unitary_change_basis_matrix
	return A * rho(permutation) * A.inverse()

[tscrim]

Suggested change

	return U.hermitian_decomposition()
	return U.hermitian_decomposition().H