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

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()))


Thank you in advance for your help


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.

  • $\begingroup$ From what you are saying, I think I need to understand how the eigenvectors are yielded because, I remember doing the proof that they are mutually orthogonal by definition. But maybe it was not in the general case, but in the particular case where there were different eigenvalues. Do you know of any good course that explains how to find eigenvalues ? $\endgroup$ – Tobbey Jun 25 '19 at 13:01
  • $\begingroup$ @Tobbey The classical reference for numerical linear algebra is Golub and van Loan's Matrix Computations. $\endgroup$ – user1551 Jun 25 '19 at 14:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.