Exploiting self-adjointness when changing basis

Question

I am using Mathematica to analyze a real, self-adjoint matrix $H$ of the size $32 \times 32$, which comes from a physics problem. In the picture there is also a matrix $Q$ which commutes with $H$.

I would like to see what $Q$ looks like in the eigenbasis of $H$, and I run into some technical issues that I was hoping someone might be kind enough to resolve.

To find an eigenbasis of $H$ is not difficult - the functions Eigenvalues and Eigenvectors work well enough. The next step is very simple in theory: combine the eigenvectors into a transition matrix $P$, and compute the matrix of $Q$ relative to the new basis as $Q' = P^{-1}Q P$. In practice, it seems that the computation of the inverse requires a lot of computation, which is time consuming now, and will probably become infeasible for larger $H$. After a moment's thought, I realized that the matters could greatly helped by the fact that $H$ is self-adjoint: it follows that it has an orthonormal eigenbasis, for which computing the inverse amounts to computing the adjoint. This would solve the matters if $H$ had all eigenvalues distinct (then $P$ mentioned above would already be self-adjoint, maybe up to re-scaling the rows), but this is not the case here.

Hence the question: Given a self-adjoint matrix $H$, together with its eigenvectors, assuming that the eigenspaces are more than $1$-dimensional, how can I effectively compute the inverse matrix? (Alternatively: how else can I find the matrix form of $Q$ relative to the eigenbasis of $H$?)

Answer 1

As $P$ is explicitly constructed from eigenvectors of a self-adjoint matrix, it is unitary, i.e $P P^\dagger = I\qquad$ where the $\dagger$ is the conjugate transpose (or Hermitian conjugate, if you prefer). So, calculating the inverse is simply ConjugateTranspose[P] which is much faster than calculating it using Inverse. That said, you have to ensure that the eigenvectors are orthonormal to properly form $P$ as Eigenvectors only guarantees linear independence in degenerate eigenspaces. A simple example of this effect is

Eigensystem[#.DiagonalMatrix[{-1, 1, 1}].ConjugateTranspose[#]]& @
  RotationMatrix[2 Pi/3, {1, 0, 1}]

which gives

{{-1, 1, 1}, {{1/3, Sqrt[2/3], 1}, {-3, 0, 1}, {-Sqrt[6], 1, 0}}}

As you can see, in this case the eigenvectors are not normalized, nor are the eigenvectors for $\lambda = 1$ orthogonal to each other. In this state, $P$ is not unitary:

#.ConjugateTranspose[#]& @ %[[2]]
(* {{16/9, 0, 0}, {0, 10, 3 Sqrt[6]}, {0, 3 Sqrt[6], 7}} *)

This must be corrected to be useful, and it can be done using Orthogonalize. The code is as follows:

P = Orthogonalize@Eigenvectors@H;
Qprime = P.Q.ConjugateTranspose[P];

Note, this uses $P Q P^\dagger$, instead of $P^\dagger Q P$ as Eigenvectors explicitly generates a matrix where the rows are the eigenvectors of $H$.