Seriously, Symmetrize the Input Matrix Before Eigendecomposition
A Quick Review of Eigendecomposition
Given a real symmetric matrix (or an Hermitian matrix in the complex case). The eigendecomposition of is given by , where is an orthogonal matrix (or a unitary matrix in the complex case) consists of the eigenvectors, and is a diagonal matrix whose diagonal elements are the eigenvalues of . 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.
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:
eigfunction 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. This behavior leads to the following two consequences:
- The resulting eigenvalues/eigenvectors may be complex, which may screw up the following computations which are expecting real inputs.
- 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 is an identity matrix. We add a very small perturbation to the second element in the first row of . 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
In : A = np.eye(3) In : A[0,1] = 1e-15 In : [v, E] = np.linalg.eig(A) In : print(v) [ 1. 1. 1.] In : 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 : np.random.seed(192) In : A = np.diag([5., 4., 3., 2., 1., 1.]) In : A += 1e-15 * np.random.randn(*A.shape) In : v, E = np.linalg.eig(A) In : print(v) [ 5. 4. 3. 2. 1. 1.] In : 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 is pretty straightforward. Because we expect our input matrix to be symmetric, we simply symmetrize it before passing it to the
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
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 : np.random.seed(192) In : A = np.diag([5., 4., 3., 2., 1., 1.]) In : A += 1e-15 * np.random.randn(*A.shape) In : v, E = np.linalg.eigh(A) In : print(v) [ 1. 1. 2. 3. 4. 5.] In : print(np.linalg.norm(E @ E.T - np.eye(*E.shape), 'fro')) 1.17359699483e-15
Now the results are correct :)
If you want to learn more about numeric algorithms of eigendecomposition. Check the book Matrix Computations by Gene H. Golub. ↩