Enumerating all matrices in Jordan normal form

Introduction

For a course in noncommutative geometry / representation theory we discussed two examples of classifying orbits of matrices, more specifically Examples 1.1 and 1.2 in Noncommutative Geometry and Cayley-smooth Orders. I won't discuss the details of the mathematics behind them, the reference should suffice.

I was intrigued by the computational aspects of the problem (the computer scientist in me pops up from time to time), and I recently gained access to a small supercomputer with a Magma installation. Not knowing the complete classification, which is given in Section 2.8 of Noncommutative Geometry and Cayley-smooth Orders under the monicker of the Gerstenhaber-Hesselink theorem (I always read the Gerstenlink-Haber theorem for some reason), I have implemented a general algorithm capable of giving you the details for all matrix rings if you are patient enough.

There are three parts involved in the classification:

1. studying the structure of the variety of singularities
2. enumerating all representatives of orbits, which are determined by Jordan normal forms
3. determining the dimensions of the orbits

The first part is computationally absolutely horrendous, but that could be my fault. I will discuss it if I feel compelled, but at the moment I am not happy with it. Degrees two and three are worked out in detail in the book, degree four is reasonably fast in my implementation while degree five already takes almost an hour... The results are not that uninteresting, so maybe I'll write about them, maybe not.

Jordan normal form enumeration

The second part is what this post is about. Recall that any (complex) matrix can be conjugated with an invertible matrix to a block diagonal matrix, whose blocks are matrices with a single eigenvalue of the original on the diagonal, ones on the superdiagonal and zeroes everywhere else. For details, check Wikipedia's Jordan normal form/Complex matrices. When trying to reproduce the diagram of Example 1.2 containing the orbit dimensions we obviously need to enumerate all possible Jordan normal forms. This is actually quite easy:

1. consider all partitions of $n$, where the numbers are given in decreasing order giving all possible configurations of the eigenvalues (these correspond to the columns in the diagram)
2. fill in all legal superdiagonals, only considering decreasing numbers of ones for eigenvalues with the same multiplicity

The eigenvalues of the representatives are taken to be unknowns in a function field, for maximal flexibility. So the following code

n := 5;
field := FunctionField(RationalField(), n);
variables := ["x", "y", "z", "u", "v", "w"];
AssignNames(~field, variables[1 .. n]);
for partition in Partitions(n) do
  representative := DiagonalFromPartition(field, partition);
  "Enumerating the orbits with representative", representative;
  for M in EnumerateSuperdiagonals(representative, partition) do
    M, " ";
  end for;
end for;

DiagonalFromPartition := function(ring, partition)
  M := Matrix(ring, Weight(partition), Weight(partition), []);
  i := 1; j := 1;
  for multiplicity in partition do
    for k in [1 .. multiplicity] do M[j, j] := ring.i; j := j + 1; end for;
    i := i + 1;
  end for;

  return M;
end function;

while EnumerateSuperDiagonals is

EnumerateSuperdiagonals := function(matrix, partition)
  matrices := {};

  configurations := CartesianProduct([[0 .. multiplicity - 1] : multiplicity in partition]);
  for configuration in configurations do
    if IsNondecreasing(partition, configuration) then
      matrices := matrices join {ConstructJordanNormalForm(matrix, partition, configuration)};
    end if;
  end for;

  return matrices;
end function;

The filter IsNondecreasing is implemented by

IsNondecreasing := function(partition, configuration)
  maximum := Infinity();
  length := partition[1];

  for i in [1 .. #partition] do
    if partition[i] lt length then
      maximum := Infinity();
      length := partition[i];
    elif configuration[i] gt maximum then
      return false;
    else
      maximum := configuration[i];
    end if;
  end for;

  return true;
end function;

while the actual construction as implemented in ConstructJordanNormalForm is given by

ConstructJordanNormalForm := function(matrix, partition, configuration)
  position := 1;
  for i in [1 .. #configuration] do
    for j in [0 .. configuration[i] - 1] do
      matrix[position + j, position + j + 1] := 1;
    end for;
    position := position + (partition[i]);
  end for;

  return matrix;
end function;

The output for this code (where $n=5$) is given in this pastebin.

The calculation of the dimensions

The third part is a bit of (symbolic) linear algebra, where you have to determine the dimension of the orbit by calculating the dimension of the stabilizer group, which is a Zariski dense subset of the centralizer of a matrix, hence suffices for giving the correct dimension. Mathematical details can be found in the text, the actual implementation is given by

OrbitDimension := function(matrix)
  field := BaseRing(matrix);
  n := Rank(field);
  ring := PolynomialRing(field, n^2);

  unknowns := Matrix(ring, n, n, [Name(ring, i) : i in [1 .. n^2]]);
  system := Matrix(ring, n^2, n^2, []);

  for i in [1 .. n] do
    for j in [1 .. n] do
      for k in [1 .. n^2] do
        system[(i-1)*n + j, k] := Coefficient((matrix*unknowns)[i, j], k, 1) - Coefficient((unknowns*matrix)[i, j], k, 1);
      end for;
    end for;
  end for;

  return n^2 - Rank(Kernel(system));
end function;

If we slightly change the iteration of the first code example, calling OrbitDimension for each matrix, the output will be as given in this pastebin in case of $n=5$.

What I have done is not the most relevant contribution to the world of mathematics, but in case you ever feel the urge to enumerate Jordan normal forms, I hope this will help you.