Keep order of diagonal elements when diagonalizing a matrix in Python

Question

I have a non-diagonal matrix, but with the non-diagonal elements much smaller than the diagonal ones, so upon diagonalization, I expect that the diagonal elements will change, but just a little. I am currently just using numpy.linalg.eig to get the eigenvalues, however I would like the order of these values to be the same as the (similar) values before diagonalization. So, if the diagonal in the original matrix was (100,200,300), after diagonalization I would like to have as output (101,202,303) and not (202,303,101). I need this in order to keep track of which value got map to which after diagonalization. How can I do this? It seems like there is no clear order to the output of numpy.linalg.eig? Thank you!

Answer 1

When you use numpy.linalg.eig to compute the eigenvalues of a matrix, the eigenvalues returned are not necessarily ordered in any particular way. You can refer to the numpy manual https://numpy.org/doc/stable/reference/generated/numpy.linalg.eig.html

Now to deal with this we can take an example matrix where the off-diagonal elements are non-zero but small and then try to map them to their original order

example_matrix = np.array([[100, 0.2, 0.3], [0.2, 200, 0.5], [0.3, 0.5, 300]])

# Compute eigenvalues and eigenvectors
eigenvalues, eigenvectors = np.linalg.eig(example_matrix)

# Original diagonal elements
original_diagonal = np.diag(example_matrix)

# Sort eigenvalues by how close they are to the original diagonal elements
sorted_indices = np.argsort(np.abs(eigenvalues - original_diagonal))

# Sorted eigenvalues and eigenvectors
sorted_eigenvalues = eigenvalues[sorted_indices]
sorted_eigenvectors = eigenvectors[:, sorted_indices]

# Output the sorted eigenvalues and their corresponding eigenvectors
sorted_eigenvalues, sorted_eigenvectors

Result

(array([ 99.99915299, 199.99789408, 300.00295293]),
 array([[-0.9999969 , -0.001985  ,  0.00150496],
        [ 0.0019925 , -0.9999855 ,  0.00500279],
        [ 0.00149501,  0.00500578,  0.99998635]]))

Basically you need to sort the eigenvalues along with their corresponding eigenvectors based on the proximity to the original diagonal. This way you can track of the permutation that sorts the eigenvalues so you can map them back to their original order.

Answer 2

Frame challenge

I don't think there can be a rigorous solution to this problem, surely not using straight numpy.

From an admittedly very cursory examination of the source code, I believe that to calculate the eigenvalues numpy starts by the hallowed method of calculating the roots of the determinant of the parametric matrix. This is the fastest and most accurate method and is guaranteed to work.

Yet, it cannot preserve eigenvalue order.

For simplicity, consider a two by two matrix:

[[ 10 .1 ] [.1 20]]

this has lambda-det of L^2-30L+199.99

But if you consider this other matrix, with the same values in the opposite order, obtained by swapping the previous diagonal...

[[ 20 .1 ] [.1 10]]

...you end up with the same determinant. So, from that one determinant, you can't tell which matrix you started from, so you can not determine the original diagonal elements' order. Proceeding with the calculations, from here onwards, that information no longer exists. You do get the eigenvalue values, but not their original order.

When you build the determinant, just as if you multiplied x*y, the information about the order of the factors gets lost (in the case of multiplication, it gets completely lost, always).

So the eigenvalues appear to be "shuffled" in respect to the original diagonal.

But what actually happens is not that the original values get "slightly shifted"; they are, each time, all completely lost, and other similar values are recreated in their place. The fact that they "look like" the previous values is a red herring. So, there's no guarantee that a match will be possible.

It's sort of the mathematical version of the "vanishing illusion" puzzles like Teddy and the Lions:

What you can try to do is match the values "a posteriori", i.e. select the permutation of the eigenvalues that most closely approximates the original diagonal.

There are algorithms to do that, using different metrics to determine which permutation is "closer" (usually, least squares).

A quick way would be to sort the initial diagonal values in ascending order, record the permutation used (this is O(n log n)), then sort the eigenvalues in ascending order and apply the permutation in reverse. This gets you a "close match".

However, it might still be possible to do something akin to what you describe using a different algorithm, e.g. Jacobi's method or other iterative method.

Note that such algorithms needs must preserve the order locally. The value in a place remains there. But it still might change so that you get values in an order different from the one you might expect, e.g. you start with diagonal [2.1, 1.9, 12], the first element shrinks to 1.999, the second grows to 2.001, and you might be led to believe that "the two elements got swapped".

Implementation found around the Interwebs

As I said, my linear algebra days are long past, but I remembered that Jacobi method was iterative and so might fit your bill. I looked for that and stopped at the first result. It looks correct to me. In case it wasn't, a more focused search might still yield a working Python implementation; the theory is proven.

GitHub page

a = array([[2.0,0.1,0.1],[0.1,7.0,0.1],[0.1, 0.1, 11.0]])
lam,x = jacobi(a,tol = 1.0e-9)
print(lam)

# Expected value: close to [2 7 11] in this order
[1.99693421  6.99940092 11.00366487]

A run of several other random values yielded consistent results.

Answer 3

import numpy as np

# Your non-diagonal matrix
A = np.array([[10, 0.1, 0.2],[0.3, 20, 0.4], [0.5, 0.6, 30]])

# Compute eigenvalues and eigenvectors
eigenvalues, eigenvectors = np.linalg.eig(A)

# Sort eigenvalues and corresponding eigenvectors based on magnitudes
sorted_indices = np.argsort(np.abs(eigenvalues))
sorted_eigenvalues = eigenvalues[sorted_indices]
sorted_eigenvectors = eigenvectors[:, sorted_indices]

# Print the sorted eigenvalues and corresponding eigenvectors
print("Original Eigenvalues:", eigenvalues)
print("Sorted Eigenvalues:", sorted_eigenvalues)
print("Corresponding Eigenvectors:", sorted_eigenvectors)

Now, you can use sorted_eigenvalues and sorted_eigenvectors as needed.

Answer 4

Let diag be the list of diagonal elements. Sorting the list with sorted will allow us to find the permutation that takes sorted_diag to diag:

sorted_diag = sorted(diag)
perm = [sorted_diag.index(d) for d in diag]

perm will indeed be the permutation, because [sorted_diag[i] for i in perm] will return diag.

Now, because the set of eigenvalues is "similar" to the set of diagonal elements, using this permutation on the sorted array of eigenvalues should yield the eigenvalues in the order we want. If eigens is the array of eigenvalues, then we can do

sorted_eigens = sorted(eigens)
eigens_good = [sorted_eigens[i] for i in perm]

So the whole thing just goes:

mat = np.array([[203, 0.4, 0.6], [0.3, 301, 0.7], [0.6, 0.1, 101]])

diag = np.diag(mat) # [203, 301, 101]
eigens, _ = np.linalg.ein(mat) # [203.00228619, 100.99612982, 301.00158399]

sorted_diag = sorted(diag) # [101, 203, 301]
perm = [sorted_diag.index(d) for d in diag] # [1, 2, 0]

sorted_eigens = sorted(eigens) # [100.99612982, 203.00228619, 301.00158399]
eigens_good = [sorted_eigens[i] for i in perm] # [203.00228619, 100.99612982, 301.00158399]

And so eigens_good is exactly what we want, as the original order was [203, 101, 301]!

note: it may be possible to optimize this more by using something different than list comprehensions and sorted, but I figure this will likely be good enough.

Same values in diagonal:

The above won't work right away if the diagonal has identical elements, but you can get something similar to work. You just need to generate perm more carefully, for example, by "popping" the elements of a copy of the sorted array.

diag = np.diag(mat)
sorted_diag = sorted(diag)

perm = [None]  * len(diag)
temp = sorted_diag.copy()
for i, d in enumerate(diag):
    index_d = temp.index(d)
    perm[i] = index_d
    temp[index_d] = None

With diag = [203, 101, 404, 101], [sorted_diag[i] for i in perm] correctly results in diag again, so it should be fine to use perm in the same as way before with eigens.