CODEWITHSUNDEEP
pythonnumpymatrixcluster-analysisprecision
kiyopi
kiyopi

Reputation: 143

What drives numerical instability in eigenvalue computations in python?

Let's say I have a data matrix X with num_samples = 1600, dim_data = 2, from which I can build a 1600*1600 similarity matrix S using the rbf kernel. I can normalize each row of the matrix, by multiplying all entries of the row by (1 / sum(entries of the row)). This procedure gives me a (square) right stochastic matrix, which we expect to have an eigenvalue equal to 1 associated to a constant eigenvector full of 1s.

We can easily check that this is indeed an eigenvector by taking its product with the matrix. However, using scipy.linalg.eig the obtained eigenvector associated to eigenvalue 1 is only piecewise constant.

I have tried scipy.linalg.eig on similarly sized matrices with randomly generated data which I transformed into stochastic matrices and consistently obtained a constant eigenvector associated to eigenvalue 1.

My question is then, what factors may cause numerical instabilities when computing eigenvalues of stochastic matrices using scipy.linalg.eig?

Reproducible example:

def kernel(sigma,X):
  """
  param sigma: variance
  param X: (num_samples,data_dim)
  """
  squared_norm = np.expand_dims(np.sum(X**2,axis=1),axis=1) + np.expand_dims(np.sum(X**2,axis=1),axis=0)-2*np.einsum('ni,mi->nm',X,X)
  return np.exp(-0.5*squared_norm/sigma**2)


def normalize(array):
  degrees = []
  M = array.shape[0]
  for i in range(M):
    norm = sum(array[i,:])
    degrees.append(norm)
  degrees_matrix = np.diag(np.array(degrees))
  P = np.matmul(np.linalg.inv(degrees_matrix),array)
  return P

#generate the data
points = np.linspace(0,4*np.pi,1600)
Z = np.zeros((1600,2))
Z[0:800,:] = np.array([2.2*np.cos(points[0:800]),2.2*np.sin(points[0:800])]).T
Z[800:,:] = np.array([4*np.cos(points[0:800]),4*np.sin(points[0:800])]).T

X = np.zeros((1600,2))
X[:,0] = np.where(Z[:,1] >= 0, Z[:,0] + .8 + params[1], Z[:,0] - .8 + params[2])
X[:,1] = Z[:,1] + params[0]

#create the stochastic matrix P
P = normalize(kernel(.05,X))

#inspect the eigenvectors
e,v = scipy.linalg.eig(P)
p = np.flip(np.argsort(e))
e = e[p]
v = v[:,p]

plot_array(v[:,0])

#check on synthetic data:
Y = np.random.normal(size=(1600,2))
P = normalize(kernel(Y))

#inspect the eigenvectors
e,v = scipy.linalg.eig(P)
p = np.flip(np.argsort(e))
e = e[p]
v = v[:,p]

plot_array(v[:,0])

Using the code provided by Ahmed AEK, here are some results on the divergence of the obtained eigenvector from the constant eigenvector.

[-1.36116641e-05 -1.36116641e-05 -1.36116641e-05 ...  5.44472888e-06
  5.44472888e-06  5.44472888e-06]
norm =  0.9999999999999999
max difference =  0.04986484253966891
max difference / element value -3663.3906291852545

UPDATE:

I have observed that a low value of sigma in the construction of the kernel matrix produces a less sharp decay in the (sorted) eigenvalues. In fact, for sigma=0.05, the first 4 eigenvalues produced by scipy.linalg.eig are rounded up to 1. This may be linked to the imprecision in the eigenvectors. When sigma is increased to 0.5, I do obtain a constant eigenvector.

First 5 eigenvectors in the sigma=0.05 case

First 5 eigenvectors in the sigma=0.05 case

First 5 eigenvectors in the sigma=0.5 case

First 5 eigenvectors in the sigma=0.5 case

Upvotes: 2

Views: 511

Answers (1)

Ahmed AEK
Ahmed AEK

Reputation: 17616

the computer has an expected accuracy of 14 digits as of 64 bits of float as shown here, which means that any result will only be accurate up to 14 digits.

using the below code you can check this result:

Y = np.random.normal(size=(1600,2))
P = normalize(kernel(5,Y))
P = P / np.sum(P,axis=1)
#inspect the eigenvectors
e,v = np.linalg.eig(P)
p = np.flip(np.argsort(e))

a = np.isclose(e,1)
e1 = e[a]
v1 = v[:,a]
v11 = v1[:,0]
print(v11)
print('norm = ',np.sum(v11**2))
print('max difference = ',np.amax(np.abs(np.diff(v11))))
print('max difference / element value =',np.amax(np.abs(np.diff(v11)))/v11[0])

result is:

[0.025+0.j 0.025+0.j 0.025+0.j ... 0.025+0.j 0.025+0.j 0.025+0.j]
norm =  (1+0j)
max difference =  1.97758476261356e-16
max difference / element value = (7.91033905045416e-15+0j)

as you can see, the difference is accuate to within 8e-15 which is around 14 digits of precision, the norm will sometimes be 0.99999999999998, which is within 14 digits of precision.

Upvotes: 0

Related Questions