Here is a list of three self-adjoint commuting matrices I would like to simultaneously diagonalize:
SageMath version 8.9, Release Date: 2019-09-29
sage: kk=QQ.algebraic_closure()
sage: L=[[[0,1,0,0,0,0,0,0],[1,4,0,1,1,1,1,0],[0,0,1,1,1,1,1,2],[0,1,1,1,1,1,1,2],[0,1,1,1,1,1,1,2],[0,1,1,1,1,1,1,2],[0,1,1,1,1,1,1,2],[0,0,2,2,2,2,2,3]],[[0,0,1,0,0,0,0,0],[0,0,1,1,1,1,1,2],[1,1,1,1,1,1,1,2],[0,1,1,2,1,1,1,2],[0,1,1,1,2,1,1,2],[0,1,1,1,1,2,1,2],[0,1,1,1,1,1,2,2],[0,2,2,2,2,2,2,3]], [[0,0,0,1,0,0,0,0],[0,1,1,1,1,1,1,2],[0,1,1,2,1,1,1,2],[1,1,2,0,0,2,3,2],[0,1,1,0,2,2,2,2],[0,1,1,2,2,1,1,2],[0,1,1,3,2,1,0,2],[0,2,2,2,2,2,2,4]]]
sage: M=[matrix(kk,l) for l in L]
sage: M
[
[0 1 0 0 0 0 0 0] [0 0 1 0 0 0 0 0] [0 0 0 1 0 0 0 0]
[1 4 0 1 1 1 1 0] [0 0 1 1 1 1 1 2] [0 1 1 1 1 1 1 2]
[0 0 1 1 1 1 1 2] [1 1 1 1 1 1 1 2] [0 1 1 2 1 1 1 2]
[0 1 1 1 1 1 1 2] [0 1 1 2 1 1 1 2] [1 1 2 0 0 2 3 2]
[0 1 1 1 1 1 1 2] [0 1 1 1 2 1 1 2] [0 1 1 0 2 2 2 2]
[0 1 1 1 1 1 1 2] [0 1 1 1 1 2 1 2] [0 1 1 2 2 1 1 2]
[0 1 1 1 1 1 1 2] [0 1 1 1 1 1 2 2] [0 1 1 3 2 1 0 2]
[0 0 2 2 2 2 2 3], [0 2 2 2 2 2 2 3], [0 2 2 2 2 2 2 4]
]
I applied the following process to find the appropriate change of basis matrix p:
sage: p=identity_matrix(kk,8)
sage: q=~p
sage: for m in M:
....: a=q*m*p
....: da, pa = a.jordan_form(transformation=True)
....: p=p*pa
....: q=~p
But the process does not finish (quickly enough, as expected), if I break it (Ctrl+C) after 5 minutes, I get:
sage: p
[1 1 1 1 1 1 1.000000000000000? 1]
[0 0 0 -0.2360679774997897? 4.236067977499789? -1 0 9]
[1 1 1 -1 -1 0 1.000000000000000? 10]
[-3.959901309262505? -0.4551399712853820? -0.3198201782016355? 0 0 1 1.734861458749522? 11]
[-1.747468453907949? 0.1971262975999527? 1.126757059617217? 0 0 1 -2.576414903309221? 11]
[1.441290910412462? 0.5668886275946296? -2.523095590631670? 0 0 1 -0.4850839473754211? 11]
[3.266078852757992? -1.308874953909201? 0.7161587092160885? 0 0 1 0.3266373919351204? 11]
[0 0 0 0.618033988749895? -1.618033988749895? -2 0 18]
[I removed spaces between entries for p being more readable]
and in fact the problem is that SageMath cannot inverte p, the following does not finish (after several hours):
sage: q=~p
whereas the matrix p is invertible:
sage: det(p)
97930.0?
What is the problem? How to solve it?