A Quick Review of Eigendecomposition

Given a real symmetric matrix $\mathbf{A}$ (or an Hermitian matrix in the complex case). The eigendecomposition of $\mathbf{A}$ is given by $\mathbf{E}\mathbf{\Lambda}\mathbf{E}^{-1}$, where $\mathbf{E}$ is an orthogonal matrix (or a unitary matrix in the complex case) consists of the eigenvectors, and $\mathbf{\Lambda}$ is a diagonal matrix whose diagonal elements are the eigenvalues of $\mathbf{A}$. The existence of eigendecomposition of real symmetric matrices (or Hermitian matrices) is backed by the spectral decomposition theorem. Eigendecomposition has many applications. In array signal processing, the subspace-based direction-finding algorithms are based on the eigendecomposition of the sample covariance matrix. In spectral graph theory, the eigendecomposition of the graph Laplacian is utilized to analyze the characteristics of graphs. We can also use eigendecomposition to compute matrix functions such as the matrix exponential as follows:

which finds important applications in linear differential equation systems.

The Issue

As the title suggests, this article is about avoiding a potential pitfall when computing eigendecompositions using MATLAB, Python, etc. I first encountered this issue when implementing the MUSIC algorithm for DOA estimation. Recently, I encountered this issue again so I think it is good idea to write an article about it. Here is the issue:

The eig function in MATLAB works for both symmetric (Hermitian) and non-symmetric matrices. If your input matrix is theoretically symmetric but in fact non-symmetric due to floating-point errors, MATLAB will use the eigendecomposition algorithm for non-symmetric matrices^[1]. This behavior leads to the following two consequences:

1. The resulting eigenvalues/eigenvectors may be complex, which may screw up the following computations which are expecting real inputs.
2. The resulting eigenvectors may not be necessarily orthogonal, which may screw up the following computations requiring orthogonality.

The same applies to Python.

To demonstrate this issue, let us first consider a very simple example where $\mathbf{A}$ is an identity matrix. We add a very small perturbation to the second element in the first row of $\mathbf{A}$. Let us first try it in MATLAB:

>> A = eye(3)
A =
     1     0     0
     0     1     0
     0     0     1

>> A(1,2) = 1e-15
A =
    1.0000    0.0000         0
         0    1.0000         0
         0         0    1.0000

>> [E, V] = eig(A)
E =
    1.0000   -0.9762         0
         0    0.2168         0
         0         0    1.0000

V =
     1     0     0
     0     1     0
     0     0     1

Oops. The eigenvalues are correct. However, the eigenvectors are definitely not correct. They are not orthogonal.

Next, let us try it in Python 3.5 with numpy 1.13.3:

In [48]: A = np.eye(3)

In [49]: A[0,1] = 1e-15

In [50]: [v, E] = np.linalg.eig(A)

In [52]: print(v)
[ 1.  1.  1.]

In [53]: print(E)
[[ 1.         -0.97622376  0.        ]
 [ 0.          0.21676522  0.        ]
 [ 0.          0.          1.        ]]

Oops again. Actually, this is expected because MATLAB and numpy should be using similar LAPACK routines.

Now let us consider a more complicated example. This time, we use a larger diagonal matrix that is not an identity matrix. We add very small Gaussian noise to each element to simulate floating-point errors. Let us try MATLAB first:

>> rng(192)
>> A = diag([5 4 3 2 1 1]);
>> A = A + 1e-15 * randn(size(A)); % Adds very small perturbations.
>> [E, V] = eig(A);
>> diag(V)
ans =
   5.0000 + 0.0000i
   4.0000 + 0.0000i
   3.0000 + 0.0000i
   2.0000 + 0.0000i
   1.0000 + 0.0000i
   1.0000 - 0.0000i

>> norm(E*E' - eye(size(E)), 'fro') % Checks orthogonality.
ans =
    0.7857

Oops. The eigenvalues now have very tiny imaginary parts and the eigenvectors are again not orthogonal. Now let us try Python:

In [72]: np.random.seed(192)

In [73]: A = np.diag([5., 4., 3., 2., 1., 1.])

In [74]: A += 1e-15 * np.random.randn(*A.shape)

In [75]: v, E = np.linalg.eig(A)

In [76]: print(v)
[ 5.  4.  3.  2.  1.  1.]

In [77]: print(np.linalg.norm(E @ E.conj().T - np.eye(*E.shape), 'fro'))
0.271374583947

Oops again. Although the perturbation is very small, it has significant effect on the computation of eigendecompositions.

The Fix

The fix is pretty straightforward. Because we expect our input matrix to be symmetric, we simply symmetrize it before passing it to the eig function:

• In MATLAB, use eig(0.5 * (A + A')) instead of eig(A), works for both real and complex cases.

• In Python, we do not need to explicitly symmetrize the input. Use np.linalg.eigh(A) instead of np.linalg.eig(A).

Not convinced? Let us run the above examples with our fix again:

>> rng(192)
>> A = diag([5 4 3 2 1 1]);
>> A = A + 1e-15 * randn(size(A));
>> [E, V] = eig(0.5 * (A + A'));
>> diag(V)
ans =
    1.0000
    1.0000
    2.0000
    3.0000
    4.0000
    5.0000

>> norm(E*E' - eye(size(E)), 'fro')
ans =
   8.9011e-16

In [82]: np.random.seed(192)

In [83]: A = np.diag([5., 4., 3., 2., 1., 1.])

In [84]: A += 1e-15 * np.random.randn(*A.shape)

In [85]: v, E = np.linalg.eigh(A)

In [86]: print(v)
[ 1.  1.  2.  3.  4.  5.]

In [87]: print(np.linalg.norm(E @ E.T - np.eye(*E.shape), 'fro'))
1.17359699483e-15

Now the results are correct :)