title: What drives numerical instability in eigenvalue computations in python?
body: <p>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 <a href="https://en.wikipedia.org/wiki/Radial_basis_function_kernel" rel="nofollow noreferrer">rbf kernel</a>. 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) <a href="https://en.wikipedia.org/wiki/Stochastic_matrix" rel="nofollow noreferrer">right stochastic matrix</a>, which we expect to have an eigenvalue equal to 1 associated to a constant eigenvector full of 1s.</p>
  <p>We can easily check that this is indeed an eigenvector by taking its product with the matrix. However, using <code>scipy.linalg.eig</code> the obtained eigenvector associated to eigenvalue 1 is only piecewise constant.</p>
  <p>I have tried <code>scipy.linalg.eig</code> on similarly sized matrices with randomly generated data which I transformed into stochastic matrices and consistently obtained a constant eigenvector associated to eigenvalue 1.</p>
  <p>My question is then, what factors may cause numerical instabilities when computing eigenvalues of stochastic matrices using <code>scipy.linalg.eig</code>?</p>
  <p>Reproducible example:</p>
  <pre><code>def kernel(sigma,X):
    &quot;&quot;&quot;
    param sigma: variance
    param X: (num_samples,data_dim)
    &quot;&quot;&quot;
    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-&gt;nm',X,X)
    return np.exp(-0.5*squared_norm/sigma**2)
  
  </code></pre>
  <pre><code>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
  </code></pre>
  <pre><code>#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] &gt;= 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])
  </code></pre>
  <p>Using the code provided by <a href="https://stackoverflow.com/users/15649230/ahmed-aek">Ahmed AEK</a>, here are some results on the divergence of the obtained eigenvector from the constant eigenvector.</p>
  <pre><code>[-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
  </code></pre>
  <p>UPDATE:</p>
  <p>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 <code>scipy.linalg.eig</code> 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.</p>
  <p>First 5 eigenvectors in the sigma=0.05 case</p>
  <p><a href="https://i.stack.imgur.com/gOX3z.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/gOX3z.png" alt="First 5 eigenvectors in the sigma=0.05 case" /></a></p>
  <p>First 5 eigenvectors in the sigma=0.5 case</p>
  <p><a href="https://i.stack.imgur.com/jwiAJ.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/jwiAJ.png" alt="First 5 eigenvectors in the sigma=0.5 case" /></a></p>