La forma más rápida de calcular k valores autovalores más grandes y vectores propios correspondientes con numpy

Question

Tengo una gran matriz simétrica densa NxN y quiero los autovectores correspondientes a los k autovalores más grandes. ¿Cuál es la mejor manera de encontrarlos (preferiblemente usando numpy pero quizás en general usando blas/atlas/lapack si ese es el único camino a seguir)? En general, N es mucho más grande que k (digamos N> 5000, k < 10).

Parece que Numpy solo tiene funciones para encontrar los k autovalores más grandes si mi matriz inicial es escasa.

Answer 1

En SciPy, puede usar la función linalg.eigh, con el parámetro eigvals.

eigvals: tupla (lo, hi) índices de los valores propios más pequeño y más grande (en orden ascendente) y los vectores propios correspondientes sean devuelto: 0 < = Mín < hi < = M-1. Si se omite, se devuelven todos los valores propios y autovectores.

En su caso, debe establecerse en (N-k,N-1).

Answer 2

En realidad las rutinas escasos trabajan para las matrices numpy densas también, creo que utilizan algún tipo de Krylov iteración subespacio, por lo tanto necesitan para calcular varios matriz-vector productos, lo que significa que si su k < < N, el las rutinas dispersas podrían ser (¿marginalmente?) más rápido.

Mira la documentación http://docs.scipy.org/doc/scipy/reference/tutorial/arpack.html

y el siguiente código (van a tomar un buen café con los amigos hasta que termina)

import numpy as np 
from time import clock 
from scipy.linalg import eigh as largest_eigh 
from scipy.sparse.linalg.eigen.arpack import eigsh as largest_eigsh 

np.set_printoptions(suppress=True) 
np.random.seed(0) 
N=5000 
k=10 
X = np.random.random((N,N)) - 0.5 
X = np.dot(X, X.T) #create a symmetric matrix 

# Benchmark the dense routine 
start = clock() 
evals_large, evecs_large = largest_eigh(X, eigvals=(N-k,N-1)) 
elapsed = (clock() - start) 
print "eigh elapsed time: ", elapsed 

# Benchmark the sparse routine 
start = clock() 
evals_large_sparse, evecs_large_sparse = largest_eigsh(X, k, which='LM') 
elapsed = (clock() - start) 
print "eigsh elapsed time: ", elapsed