simultaneous diagonalization of set matrices

Question

I have a set of integral square invertible symmetric matrices $A_i$ with $A_i^2=I$ (so also $A_i A_i^T=I$). The matrices commute. I'd like to map them simultaneously to a set of diagonal matrices $D_i$ using a matrix $C : C A_i C^{T}= D_i$. The $D_i$'s are the diagonals of the eigenvalues of $A_i$ in some particular order (I know the $D_i$'s already). I know simultaneous diagonalization isn't in general easy, but maybe this special case has a clever solution. Also interested if there's anything in GAP that might help.

In this older post https://mathematica.stackexchange.com/questions/46949/is-there-a-built-in-procedure-for-simultaneous-diagonalization-of-a-set-of-commu the accepted answer says that you can take the eigenvectors of a "random linear combination" of the matrices. My matrices are diagonalizable so the part about geometric/arithmetic multiplicity applies. This approach is not enough for what I'm doing; there's no guarantee that the linear combination happens to be a right one; also the probability of success doesn't seem that high for what I tried.

Answer 1

Simultaneous diagonalization of a set of commuting matrices is quite easy -- since the matrices commute they preserve each others eigenspaces. Thus find a eigenvector basis for the first matrix, and then split each eigenspace using the second matrix and so on. The following GAP code does this:

# Arguments: Field, matrixlist
SimultaneousDiagonalization:=function(F,mats)
local bas,m,nbas,b,c,start,j,ev,eigen,ran,block;
  # basis so far: List of eigenspace bases
  bas:=[IdentityMat(Length(mats[1]),F)];
  for m in mats do
    nbas:=[];
    #rewrite m wrt the eigenspaces so far
    b:=Concatenation(bas);
    c:=b*m/b;
    start:=0;
    for j in bas do
      # now split this eigenspace
      ran:=[start+1..start+Length(j)];
      block:=c{ran}{ran};
      ev:=Eigenvalues(F,block);
      eigen:=List(ev,x->NullspaceMat(block-x*block^0));
      # and store the new spaces to use instead
      Append(nbas,List(eigen,k->k*j));
      start:=start+Length(j);
    od;
    bas:=nbas;
  od;
  return bas;
end;

It returns a list of bases of the common eigenspaces. Thus, if the result is B then with C:=Concatenation(B), you have that $C\cdot A_i\cdot C^{-1}$ is diagonal. To get an orthogonal matrix, you simply run Gram-Schmidt orthonormalization on each of these eigenspace bases.

The result often will be ugly, since GAP does not naturally work with real numbers, thus I'm not attempting to do so.