Numerical stability: cannot unitarily diagonalize normal matrices

Question

I am trying to setup small numerical experiments to see if a unitary matrix can be unitarily diagonalized thanks to the spectral theorem: https://en.wikipedia.org/wiki/Spectral_theorem (see normal matrices).

However, in my example, although I can check that the matrix is unitary up to some good absolute precision, I cannot obtain the diagonalization with unitary matrix with an acceptable precision.

In my opinion, there is a mathematical assumption that I am doing that is clearly wrong, maybe related to the behaviour of the eigen decomposition, but I cannot see which one it is.

Here is my example in python using numpy and eigen decomposition:

If you try to run it, you will notice that the check for unitarity of the matrix of eigenvectors fails, and that the result of the product of v with v.T.conj() is wrong by more than 10e-1 in some matrix coefficients.

I really wonder what is the reason underlying this pretty large discrepancy, although the condition number of both G and V should in theory be close to the smallest existing one (1).

import numpy as np

N=4
x = np.linspace(0, N-1, N) 
y = np.linspace(0, N-1, N)
xg, yg = np.meshgrid(x, y)
F = np.exp(2*np.pi*1j*xg*yg/N)
G = F/np.sqrt(N)

def check_unitary_matrices(M):
    # Now check invert
    assert np.allclose(np.dot(M.T.conj(), M), np.identity(N, dtype=np.complex64), atol=1e-5)
    assert np.allclose(np.dot(M, M.T.conj()), np.identity(N, dtype=np.complex64), atol=1e-5)

def check_unitarily_diagonalizable(M):
    w,v = np.linalg.eig(G)
    print(np.dot(v, v.T.conj()))
    check_unitary_matrices(v)

check_unitary_matrices(G)
check_unitarily_diagonalizable(G)

Thank you in advance for your help

Answer 1

When $N=4$, your unitary matrix $G$ has four eigenvalues $1,1,-1,i$. Since $1$ is a repeated eigenvalue, numpy.linalg.eig will produce two linearly independent eigenvectors in the corresponding eigenspace, but these two eigenvectors are not guaranteed to be mutually orthogonal.

In general, given a normal matrix $G$, its Schur triangulation $G=URU^\ast$ (apparently there is a SciPy implementation) will automatically be a unitary diagonalisation. However, to clean up possible rounding errors in the calculation of the strictly upper triangular part of $R$, you should either zero out the strictly upper triangular part or extract the diagonal part of $R$ in practice.