# Simultaneous diagonalization of commuting Hermitian matrices

**URL:** <https://discourse.julialang.org/t/simultaneous-diagonalization-of-commuting-hermitian-matrices/95495>\
**Category:** Numerics\
**Created:** [March 3, 2023, 12:45pm UTC](https://discourse.julialang.org/t/simultaneous-diagonalization-of-commuting-hermitian-matrices/95495 "2023-03-03T12:45:55Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![Gregstrq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gregstrq/32/20620_2.png) [@Gregstrq](https://discourse.julialang.org/u/Gregstrq)\
**Post date:** [March 3, 2023, 12:45pm UTC](https://discourse.julialang.org/t/simultaneous-diagonalization-of-commuting-hermitian-matrices/95495/1 "2023-03-03T12:45:55Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [March 4, 2023, 10:08am UTC](https://discourse.julialang.org/t/simultaneous-diagonalization-of-commuting-hermitian-matrices/95495/2 "2023-03-04T10:08:43Z")

</div>

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

> <https://github.com/dlfivefifty/NumericalRepresentationTheory.jl/blob/5cd8f76c6b5624597386ef5767add3168339eeb9/src/NumericalRepresentationTheory.jl#L327>

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

---

<div class="post-metadata">

**Author:** ![mstewart](https://avatars.discourse-cdn.com/v4/letter/m/b5a626/32.png) [@mstewart](https://discourse.julialang.org/u/mstewart)\
**Post date:** [March 4, 2023, 2:04pm UTC](https://discourse.julialang.org/t/simultaneous-diagonalization-of-commuting-hermitian-matrices/95495/3 "2023-03-04T14:04:23Z")

</div>

> 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

```julia
_, 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.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 4, 2023, 2:52pm UTC](https://discourse.julialang.org/t/simultaneous-diagonalization-of-commuting-hermitian-matrices/95495/4 "2023-03-04T14:52:54Z")

</div>

> [@Gregstrq](#):
>
> 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.

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^\*.)

---

<div class="post-metadata">

**Author:** ![abulak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abulak/32/28314_2.png) [@abulak](https://discourse.julialang.org/u/abulak)\
**Post date:** [March 7, 2023, 3:50pm UTC](https://discourse.julialang.org/t/simultaneous-diagonalization-of-commuting-hermitian-matrices/95495/5 "2023-03-07T15:50:25Z")

</div>

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
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
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 😃. If you have more complicated symmetries (e.g. [signed] permutations) there’s much more reduction possible with minimal projection system etc. 😉