# Is there a built-in procedure for simultaneous diagonalization of a set of commuting matrices?

Given a set $\{A_1,...,A_m\}$ of $m$ commuting $N\times N$ diagonalizable matrices, it is known that there exists a basis of eigenvectors $\Lambda$ that simultaneously diagonalizes all the $A_i$.

Is there an automated way to compute $\Lambda$ in Mathematica?

For a pair of matrices $\{A_1,A_2\}$ it can be done "by hand" using Eigenvectors to first compute the eigenvectors $\Lambda_1$ of $A_1$, finding similarity transforms $P_k$ that act on the degenerate subspaces of $\Lambda_1$ to make them also eigenvector of $A_2$, and then combining the transformed degenerate subspaces together to form $\Lambda$, but I don't know how to generalize this to more than 2 matrices.

-
I'll mention the related post at Wolfram Community –  Daniel Lichtblau Apr 28 '14 at 16:16
@DanielLichtblau: Thanks, I applied Eigenvectors to the matrix pencil $L(\lambda)=\sum_{j=1}^m\lambda^jA_j$ and it worked for $m=3$ matrices using a random real value for $\lambda$, providing a basis set which simultaneously diagonalized all three. How sensitive is it to floating point errors in the $A_j$ that spoil the commutativity slightly? It seems there is a fair amount of material available online on it, so I'll look into it further. –  DumpsterDoofus Apr 28 '14 at 18:08
I don't know how sensitive it is. I suspect it is fine when the "near commuting" n x n matrices each have n distinct (and reasonably well separated) eigenvalues. At degeneracies it might not be so nice since eigenspaces can start to get pretty cranky. –  Daniel Lichtblau Apr 28 '14 at 18:37
I'll add that use of such matrices shows up in numerically solving polynomial systems via eigensystems. My take is that one should do as much as possible to avoid multiple eigendecompositions, so I've used random linear operators to avoid multiplicity whenever possible. This is, in effect, what you did via random lambdas. How well this generalizes to other applications I cannot say. –  Daniel Lichtblau Apr 28 '14 at 18:40
@DanielLichtblau: Interesting, my 3 test matrices each had several multiple-degenerate subspaces, but it still seemed to work fine without getting cranky. –  DumpsterDoofus Apr 28 '14 at 19:17

That's probably why it worked in my case: all the matrices $A_k$ I had were all either unitary or symmetric, and thus admitted eigendecompositions. Thanks! –  DumpsterDoofus Apr 28 '14 at 23:07
Minor update: Actually, the matrices I used were symmetric. I tried it on unitary matrices (not necessarily symmetric) with root-of-unity eigenvalues and it failed. So some caution and double-checking may be necessary (by verifying that $\Lambda^{-1}\mathbf{A}_k\Lambda$ is in fact diagonal for all $k$). –  DumpsterDoofus Apr 29 '14 at 0:11