How to transform an array using its spectral components

Question

I am trying to transform an array using its spectral components.

Let's consider an 4x4 array arr
It will have an adjacent matrix A, a degree matrix D and Laplacian matrix L = D-A

import numpy as np
from numpy.linalg import eig

Compute its eigenvectors and eigenvalues, and sort it in descending order

eig_val, eig_vec = eig(L)
idx = eig_val.argsort()[::-1]

Sort the eigenvectors accordingly

eig_vec = eig_vec[:,idx]

The product of 2 distincts eigenvectors must be 0. I notice that this is not the case here e.g. the product of the first and second eigenvectors is:

sum(np.multiply(eig_vec[0], eig_vec[1])) = 0.043247527085787975

Is there anything I am missing ?

Compute the spectral components of the input array

spectral = np.matmul(eig_vec.transpose(), arr.flatten())
print(spectral.shape)

Take the first 15 components.

masked = np.zeros(spectral.shape)
m = spectral[:15]
masked[:15] = m

Get the new features

updated_arr = np.matmul(eig_vec, masked)
updated_arr = updated_arr.reshape(4, -1)

The updated array is very different from the original.

Any suggestion or resource to have a look at will be welcome.

Answer 1

To get the eigendecomposition of a symmetric matrix, you should use the eigh function instead of eig. The following code works for me:

import numpy as np
from numpy.linalg import eigh

A = np.asarray([[0., 1., 0., 0., 1., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
                   [1., 0., 1., 0., 1., 1., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
                   [0., 1., 0., 1., 0., 1., 1., 1., 0., 0., 0., 0., 0., 0., 0., 0.],
                   [0., 0., 1., 0., 0., 0., 1., 1., 0., 0., 0., 0., 0., 0., 0., 0.],
                   [1., 1., 0., 0., 0., 1., 0., 0., 1., 1., 0., 0., 0., 0., 0., 0.],
                   [1., 1., 1., 0., 1., 0., 1., 0., 1., 1., 1., 0., 0., 0., 0., 0.],
                   [0., 1., 1., 1., 0., 1., 0., 1., 0., 1., 1., 1., 0., 0., 0., 0.],
                   [0., 0., 1., 1., 0., 0., 1., 0., 0., 0., 1., 1., 0., 0., 0., 0.],
                   [0., 0., 0., 0., 1., 1., 0., 0., 0., 1., 0., 0., 1., 1., 0., 0.],
                   [0., 0., 0., 0., 1., 1., 1., 0., 1., 0., 1., 0., 1., 1., 1., 0.],
                   [0., 0., 0., 0., 0., 1., 1., 1., 0., 1., 0., 1., 0., 1., 1., 1.],
                   [0., 0., 0., 0., 0., 0., 1., 1., 0., 0., 1., 0., 0., 0., 1., 1.],
                   [0., 0., 0., 0., 0., 0., 0., 0., 1., 1., 0., 0., 0., 1., 0., 0.],
                   [0., 0., 0., 0., 0., 0., 0., 0., 1., 1., 1., 0., 1., 0., 1., 0.],
                   [0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 1., 1., 0., 1., 0., 1.],
                   [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 1., 0., 0., 1., 0.]])

degree = A.sum(axis=1)
D = np.diag(degree)

# Laplacian matrix
L = D - A

eig_val, eig_vec = eigh(L)

print(np.mean(np.abs(np.dot(np.dot(eig_vec, np.diag(eig_val)), eig_vec.transpose()) - L)))
# -> 1.5500965941594069e-15

def compare(k):
    cmp_val = eig_val[len(eig_val)-k:]
    cmp_vec = eig_vec[:, len(eig_val)-k:]
    print(np.mean(np.abs(np.dot(np.dot(cmp_vec, np.diag(cmp_val)), cmp_vec.transpose()) - L)))

compare(4)
# -> 0.4905237882047826

compare(8)
# -> 0.3885449860779778

compare(12)
# -> 0.1260071124311158

compare(14)
# -> 0.05127788166023957

compare(15)
# -> 1.5583365306702967e-15

The fact is that your L matrix has repeated eigenvalues (for example the largest eigenvalue has multiplicity of 2) and so your eigenvectors do not have to be necessarily orthogonal (as any linear combination of eigenvectors related to the same eigenvalue is also an eigenvector). And in this case eig does not guarantee returning you orthogonal eigenvectors. You can test it:

eig_val, eig_vec = eig(L)
idx = eig_val.argsort()[::-1]
print(eig_val[idx][:2])
print(np.dot(eig_vec[:, idx][:,0], eig_vec[:, idx][:,1]))
# -> [9.80529155 9.80529155]
# -> 0.017548101393349412

Unfortunately, the documentation of eigh does not promise you orthogonality either, however it mentions that it uses the lapack function _syevd for real matrices of which documentation do promise to return orthogonal eigenvectors. As eig is designed to find the right eigenvectors of general (not necessary symmetric) matrices, it does not even try to orthogonalize them (it uses the lapack routine _geev of which documentation does not mention orthogonalization).

Answer 2

Answering the part about the orthogonal eigenvectors: np.linalg.eig returns the eigenvectors, such that the column eigenvectors[:,i] is the eigenvector corresponding to the eigenvalue eigenvalues[i]. So calculating the product like this

sum(np.multiply(eig_vec[:,0], eig_vec[:,1]))

the sum will be zero.