# Numerical stability: cannot unitarily diagonalize normal matrices

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)


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.