Simultaneous diagonalization of commuting Hermitian matrices

Specific Domains Numerics
1

I have two Hermitian matrices, one of which is positive definite, and the matrices commute.
Is there a standard implementation for the simultaneous diagonalization of such two matrices?

To be precise, I have a Hamiltonian H with parity-symmetry and the parity operator P, and I want to find the common eigenbasis. Instead of P, I can consider (2I+P) operator, which is positive definite.

2

I believe a quick and dirty way to do this is either diagonalise a random linear combination or diagonalise A +i *B

I have an implementation of the “Diagonalise then Deflate” approach

But I believe that’s known to be numerically unstable In general.

3

I believe a quick and dirty way to do this is either diagonalise a random linear combination or diagonalise A +i *B

Diagonalizing A+iB is really pretty good. I believe it should be perfectly stable if the similarity is computed from something like

_, Q = schur(A + im*B)

There is a catch if you use eigen for this, however. If you ask for eigenvectors, a general eigenvalue solver will compute the Schur form and solve additional equations for the eigenvectors. For the normal matrix A+iB, you don’t need that and it can mess up the orthogonality. I just did a test of that on A and B that had approximately repeated eigenvalues and the loss of orthogonality was noticeable. Schur is the way to go.

But I believe that’s known to be numerically unstable In general.

Indeed. Stability is generally a lost cause in algorithms that depend in a critical way on identifying clusters of eigenvalues to avoid ill-conditioning, although they can work fine if you have a favorable distribution of eigenvalues. There are also Jacobi methods, which avoid clustering, extend to the case of m matrices, and should be stable, but I think computing the rotations required is not very simple.

1 Like
4

You can just do a projection of H onto smaller matrices corresponding to the degenerate subspaces of P, no?

Suppose H and P are 2m\times 2m. Say that half of the eigenvectors of P are symmetric (even) vectors (Px=x), and half are anti-symmetric (odd) vectors (Px=-x). Let the columns of the 2m \times m matrix E be an orthonormal basis for the even vectors, and the columns of O be an orthonormal basis for the odd vectors. Then to get the even eigenvectors of H you diagonalize the m \times m matrix E^* H E, and multiply the resulting eigenvectors by E. Similarly for O.

(Note that diagonalizing two m \times m matrices should be roughly 4\times faster than diagonalizing one 2m \times 2m matrix, so this is also a performance win. There should also be no issues with numerical stability.)

Moreover, for the parity operator P is so simple that you should be able to find E and O analytically. For example, if you have a grid of 8 points along a line, with P being the mirror-flip operator, then you get

E = \frac{1}{\sqrt{2}} \begin{pmatrix} 1 & & & \\ & 1 & & \\ & & 1 & \\ & & & 1 \\ & & & 1 \\ & & 1 & \\ & 1 & & \\ 1 & & & \end{pmatrix}, \qquad O = \frac{1}{\sqrt{2}} \begin{pmatrix} 1 & & & \\ & 1 & & \\ & & 1 & \\ & & & 1 \\ & & & -1 \\ & & -1 & \\ & -1 & & \\ -1 & & & \end{pmatrix}

(If you are discretizing a PDE, then E^* H E or O^* H O are exactly equivalent to discretizing a computational cell with half the size using a Neumann/even or Dirichlet/odd boundary in the middle, respectively.)

(This is a pretty common procedure in solid-state physics, where you often project the Hamiltonian into smaller subspaces obtained from the crystallographic symmetries. Group representation theory tells you the projection operators, e.g. EE^* and OO^*.)

3 Likes
5

Or if you have a direct description of the symmetry (i.e. orthogonal action and the basis for the linear space) you could try SymbolicWedderburn.jl which will compute:

  • the projections in the group ring (e.g. (1//2*() + 1//2*(1,2) and 1//2*() - 1//2*(1,2))
  • and realize them as orthogonal projections in the basis using the representation.

In the sign action this will reduce to what @stevengj said, but it’s somehow automated if you’re willing to write a few functions. Here’s a direct example for the dihedral group acting on the linear space of polynomials in x and y of degree at most 4 (it is taken/modified from the examples/ex_robinson_form.jl);

julia> using SymbolicWedderburn, DynamicPolynomials, GroupsCore

julia> include("examples/action_polynomials.jl")

julia> include("examples/dihedral.jl")

julia> @polyvar x y
(x, y)

julia> struct DihedralAction <: OnMonomials end

julia> SymbolicWedderburn.coeff_type(::DihedralAction) = Float64

julia> function SymbolicWedderburn.action(
           ::DihedralAction,
           el::DihedralElement,
           mono::AbstractMonomial,
       )
           if iseven(el.reflection + el.id)
               var_x, var_y = x, y
           else
               var_x, var_y = y, x
           end
           sign_x = 1 <= el.id <= 2 ? -1 : 1
           sign_y = 2 <= el.id ? -1 : 1
           return mono([x, y] => [sign_x * var_x, sign_y * var_y])
       end

julia> G = DihedralGroup(4)
DihedralGroup(4)

julia> basis = monomials([x,y], 0:4)
15-element MonomialVector{true}:
 x⁴
 x³y
 x²y²
 xy³
 y⁴
 x³
 x²y
 xy²
 y³
 x²
 xy
 y²
 x
 y
 1

julia> sa_basis = SymbolicWedderburn.symmetry_adapted_basis(Float64, G, DihedralAction(), basis, semisimple=true)
5-element Vector{SymbolicWedderburn.DirectSummand{Float64, SparseArrays.SparseMatrixCSC{Float64, Int64}}}:
 [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 1.0 0.0]
 [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 1.0; -0.7071067811865472 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
 [0.0 -0.7071067811865472 … 0.0 0.0]
 [0.0 0.0 … 0.0 0.0; 0.0 -0.7071067811865472 … 0.0 0.0]
 [-0.7071067811865472 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]

As can be seen, this action has 5 isotypical components (defined by irreducible reps of the dihedral group) and sa_basis defines pairwise orthogonal projections onto their isotypical subspaces (each projection has additionally orthogonal rows); These are the bases for the subspaces in a more human readable form:

julia> sa_basis[1]*basis # action via multiplication by -1
6-element Vector{Polynomial{true, Float64}}:
 x³
 x²y
 xy²
 y³
 x
 y

julia> sa_basis[2]*basis # trivial action = fixed subspaces 
4-element Vector{Polynomial{true, Float64}}:
 x²y²
 1.0
 -0.7071067811865472x⁴ - 0.7071067811865475y⁴
 -0.7071067811865472x² - 0.7071067811865475y²

julia> sa_basis[3]*basis #
1-element Vector{Polynomial{true, Float64}}:
 -0.7071067811865472x³y + 0.7071067811865475xy³

julia> sa_basis[4]*basis
2-element Vector{Polynomial{true, Float64}}:
 xy
 -0.7071067811865472x³y - 0.7071067811865475xy³

julia> sa_basis[5]*basis
2-element Vector{Polynomial{true, Float64}}:
 -0.7071067811865472x⁴ + 0.7071067811865475y⁴
 -0.7071067811865472x² + 0.7071067811865475y²

@Gregstrq I admit that for sign action that might be an overkill but hey, I can showcase my package :smiley:. If you have more complicated symmetries (e.g. [signed] permutations) there’s much more reduction possible with minimal projection system etc. :wink:

2 Likes