Localizando os menores autovetores de matriz esparsa grande, 100x mais lenta no SciPy do que no Octave

Estou tentando calcular alguns (5-500) autovetores correspondentes aos menores autovalores de matrizes esparsas quadradas simétricas grandes (até 30000x30000) com menos de 0,1% dos valores diferentes de zero.

Atualmente, estou usando scipy.sparse.linalg.eigsh no modo shift-invert (sigma = 0,0), que eu descobri através de várias postagens sobre o tópico como a solução preferida. No entanto, leva até 1 hora para resolver o problema na maioria dos casos. Por outro lado, a função é muito rápida, se eu pedir os maiores valores próprios (sub-segundos no meu sistema), o que era esperado na documentação.

Como eu estou mais familiarizado com o Matlab do trabalho, tentei resolver o problema no Octave, o que me deu o mesmo resultado usando eigs (sigma = 0) em apenas alguns segundos (sub 10s). Como eu quero fazer uma varredura de parâmetro do algoritmo, incluindo o cálculo do vetor próprio, esse tipo de ganho de tempo também seria ótimo em python.

Mudei primeiro os parâmetros (especialmente a tolerância), mas isso não mudou muito nas escalas de tempo.

Estou usando o Anaconda no Windows, mas tentei alternar o LAPACK / BLAS usado pelo scipy (que era uma grande dor) de mkl (padrão Anaconda) para OpenBlas (usado pelo Octave de acordo com a documentação), mas não consegui ver uma alteração no desempenho.

Não consegui descobrir se havia algo a mudar sobre o ARPACK usado (e como)?

Fiz upload de um caso de teste para o código abaixo na seguinte pasta dropbox: https://www.dropbox.com/sh/l6aa6izufzyzqr3/AABqij95hZOvRpnnjRaETQmka?dl=0

Em Python

import numpy as np
from scipy.sparse import csr_matrix, csc_matrix, linalg, load_npz   
M = load_npz('M.npz')
evals, evecs = linalg.eigsh(M,k=6,sigma=0.0)

Na oitava:

M=dlmread('M.txt');
M=spconvert(M);
[evecs,evals] = eigs(M,6,0);

Qualquer ajuda é apreciada!

Tentei algumas opções adicionais com base nos comentários e sugestões:

Oitava: eigs(M,6,0)e eigs(M,6,'sm')me dê o mesmo resultado:

[1.8725e-05 1.0189e-05 7.5622e-06 7.5420e-07 -1.2239e-18 -2.5674e-16]

enquanto eigs(M,6,'sa',struct('tol',2))converge para

[1.0423 2.7604 6.1548 11.1310 18.0207 25.3933] 

muito mais rápido, mas somente se os valores de tolerância estiverem acima de 2, caso contrário, eles não convergirão e os valores serão fortemente diferentes.

Python: eigsh(M,k=6,which='SA')e eigsh(M,k=6,which='SM')ambos não convergem (erro ARPACK em nenhuma convergência alcançada). eigsh(M,k=6,sigma=0.0)Fornece apenas alguns autovalores (após quase uma hora), que são diferentes da oitava para os menores (até 1 pequeno valor adicional é encontrado):

[3.82923317e-17 3.32269886e-16 2.78039665e-10 7.54202273e-07 7.56251500e-06 1.01893934e-05]

Se a tolerância é alta o suficiente, também obtenho resultados eigsh(M,k=6,which='SA',tol='1'), que se aproximam dos outros valores obtidos

[4.28732218e-14 7.54194948e-07 7.56220703e-06 1.01889544e-05, 1.87247350e-05 2.02652719e-05]

novamente com um número diferente de pequenos autovalores. O tempo de computação ainda é quase 30min. Embora os diferentes valores muito pequenos possam ser compreensíveis, uma vez que podem representar múltiplos de 0, a multiplicidade diferente me confunde.

Além disso, parece haver algumas diferenças fundamentais em SciPy e Octave, que ainda não consigo descobrir.

Uma conjectura e alguns comentários, já que não tenho Matlab / Octave:

Para encontrar pequenos autovalores de matrizes simétricas com autovalores> = 0, como o seu, o seguinte é muito mais rápido que o shift-invert:

# flip eigenvalues e.g.
# A:     0 0 0 ... 200 463
# Aflip: 0 163 ... 463 463 463
maxeval = eigsh( A, k=1 )[0]  # biggest, fast
Aflip = maxeval * sparse.eye(n) - A
bigevals, evecs = eigsh( Aflip, which="LM", sigma=None ... )  # biggest, near 463
evals = maxeval - bigevals  # flip back, near 463 -> near 0
# evecs are the same

eigsh( Aflip )para grandes pares próprios é mais rápido que o inverso por turno para pequenos, porque A * xé mais rápido do que o solve()inverso por turno deve fazer. É possível que o Matlab / Octave possa fazer isso Aflipautomaticamente, após um rápido teste de definição positiva com Cholesky.
Você pode executar eigsh( Aflip )no Matlab / Octave?

Outros fatores que podem afetar a precisão / velocidade:

O padrão do Arpack para o vetor inicial v0é um vetor aleatório. Eu uso v0 = np.ones(n), o que pode ser terrível para alguns, Amas é reproduzível :)

Essa Amatriz é quase exatamente sigular, A * ones~ 0.

Multicore: o scipy-arpack com openblas / Lapack usa ~ 3.9 dos 4 núcleos no meu iMac; o Matlab / Octave usa todos os núcleos?

k 10  tol 1e-05:    8 sec  eigvals [0 8.5e-05 0.00043 0.0014 0.0026 0.0047 0.0071 0.0097 0.013 0.018] 
k 10  tol 1e-06:   44 sec  eigvals [0 3.4e-06 2.8e-05 8.1e-05 0.00015 0.00025 0.00044 0.00058 0.00079 0.0011] 
k 10  tol 1e-07:  348 sec  eigvals [0 3e-10 7.5e-07 7.6e-06 1.2e-05 1.9e-05 2.1e-05 4.2e-05 5.7e-05 6.4e-05] 

k 20  tol 1e-05:   18 sec  eigvals [0 5.1e-06 4.5e-05 0.00014 0.00023 0.00042 0.00056 0.00079 0.0011 0.0015 0.0017 0.0021 0.0026 0.003 0.0037 0.0042 0.0047 0.0054 0.006
k 20  tol 1e-06:   73 sec  eigvals [0 5.5e-07 7.4e-06 2e-05 3.5e-05 5.1e-05 6.8e-05 0.00011 0.00014 0.00016 0.0002 0.00025 0.00027 0.0004 0.00045 0.00051 0.00057 0.00066
k 20  tol 1e-07:  267 sec  eigvals [-4.8e-11 0 7.5e-07 7.6e-06 1e-05 1.9e-05 2e-05 2.2e-05 4.2e-05 5.1e-05 5.8e-05 6.4e-05 6.9e-05 8.3e-05 0.00011 0.00012 0.00013 0.00015

k 50  tol 1e-05:   82 sec  eigvals [-4e-13 9.7e-07 1e-05 2.8e-05 5.9e-05 0.00011 0.00015 0.00019 0.00026 0.00039 ... 0.0079 0.0083 0.0087 0.0092 0.0096 0.01 0.011 0.011 0.012
k 50  tol 1e-06:  432 sec  eigvals [-1.4e-11 -4e-13 7.5e-07 7.6e-06 1e-05 1.9e-05 2e-05 2.2e-05 4.2e-05 5.1e-05 ... 0.00081 0.00087 0.00089 0.00096 0.001 0.001 0.0011 0.0011
k 50  tol 1e-07: 3711 sec  eigvals [-5.2e-10 -4e-13 7.5e-07 7.6e-06 1e-05 1.9e-05 2e-05 2.2e-05 4.2e-05 5.1e-05 ... 0.00058 0.0006 0.00063 0.00066 0.00069 0.00071 0.00075

versions: numpy 1.18.1  scipy 1.4.1  umfpack 0.3.2  python 3.7.6  mac 10.10.5 

Matlab / Octave são a mesma coisa? Caso contrário, todas as apostas estão fora - primeiro verifique a exatidão e depois a velocidade.

Por que os autovalores oscilam tanto? Minúsculo <0 para uma matriz definida supostamente não negativa é um sinal de erro de arredondamento , mas o truque usual de uma pequena mudança, A += n * eps * sparse.eye(n)não ajuda.

Espero que isto ajude.

Sei que agora é antigo, mas tive o mesmo problema. Você revisou aqui ( https://docs.scipy.org/doc/scipy/reference/tutorial/arpack.html )?

Parece que quando você define sigma para um número baixo (0), deve definir qual = 'LM', mesmo que você queira valores baixos. Isso ocorre porque a configuração sigma transforma os valores que você deseja (baixo neste caso) para parecer alto e, portanto, ainda é possível tirar proveito dos métodos 'LM', que são muito mais rápidos para obter o que você deseja (os baixos valores próprios )

Quero dizer primeiro que não tenho idéia de por que os resultados que você e a @Bill relataram são do jeito que são. Eu simplesmente me pergunto se eigs(M,6,0)em Octave corresponde a k=6 & sigma=0, ou talvez seja outra coisa?

Com o scipy, se eu não fornecer sigma, posso obter um resultado em um tempo decente dessa maneira.

import numpy as np
from scipy.sparse import csr_matrix
from scipy.sparse.linalg import eigsh
from time import perf_counter
M = np.load('M.npz')
a = csr_matrix((M['data'], M['indices'], M['indptr']), shape=M['shape'])
t = perf_counter()
b, c = eigsh(a, k=50, which='SA', tol=1e1)
print(perf_counter() - t)
print(b)

Não tenho muita certeza se isso faz sentido.

0.4332823531003669
[4.99011753e-03 3.32467891e-02 8.81752215e-02 1.70463893e-01
 2.80811313e-01 4.14752072e-01 5.71103821e-01 7.53593653e-01
 9.79938915e-01 1.14003837e+00 1.40442848e+00 1.66899183e+00
 1.96461415e+00 2.29252666e+00 2.63050114e+00 2.98443218e+00
 3.38439528e+00 3.81181747e+00 4.26309942e+00 4.69832271e+00
 5.22864462e+00 5.74498014e+00 6.22743988e+00 6.83904055e+00
 7.42379697e+00 7.97206446e+00 8.62281827e+00 9.26615266e+00
 9.85483434e+00 1.05915030e+01 1.11986296e+01 1.18934953e+01
 1.26811461e+01 1.33727614e+01 1.41794599e+01 1.47585155e+01
 1.55702295e+01 1.63066947e+01 1.71564622e+01 1.78260727e+01
 1.85693454e+01 1.95125277e+01 2.01847294e+01 2.09302671e+01
 2.18860389e+01 2.25424795e+01 2.32907153e+01 2.37425085e+01
 2.50784800e+01 2.55119112e+01]

A única maneira que encontrei de usar o sigma e obter um resultado em um tempo decente é fornecer M como um LinearOperator. Não estou muito familiarizado com isso, mas pelo que entendi minha implementação representa uma matriz de identidade, que é o que M deveria ser se não for especificado na chamada. A razão para isso é que, em vez de executar uma resolução direta (decomposição da LU), o scipy usará um solucionador iterativo, potencialmente mais adequado. Como comparação, se você fornecer o M = np.identity(a.shape[0])que deve ser exatamente o mesmo, o eigsh levará uma eternidade para produzir um resultado. Observe que essa abordagem não funciona se sigma=0for fornecida. Mas não tenho certeza se isso sigma=0é realmente útil?

import numpy as np
from scipy.sparse import csr_matrix
from scipy.sparse.linalg import eigs, eigsh, LinearOperator
from time import perf_counter


def mv(v):
    return v


M = np.load('M.npz')
a = csr_matrix((M['data'], M['indices'], M['indptr']), shape=M['shape'])
t = perf_counter()
b, c = eigsh(a, M=LinearOperator(shape=a.shape, matvec=mv, dtype=np.float64),
             sigma=5, k=50, which='SA', tol=1e1, mode='cayley')
print(perf_counter() - t)
print(np.sort(-5 * (1 + b) / (1 - b)))

Novamente, não faço ideia se está correto, mas definitivamente diferente de antes. Seria ótimo ter a contribuição de alguém do scipy.

1.4079377939924598
[3.34420263 3.47938816 3.53019328 3.57981026 3.60457277 3.63996294
 3.66791416 3.68391585 3.69223712 3.7082205  3.7496456  3.76170023
 3.76923989 3.80811939 3.81337342 3.82848729 3.84137264 3.85648208
 3.88110869 3.91286153 3.9271108  3.94444577 3.97580798 3.98868207
 4.01677424 4.04341426 4.05915855 4.08910692 4.12238969 4.15283192
 4.16871081 4.1990492  4.21792125 4.24509036 4.26892806 4.29603036
 4.32282475 4.35839271 4.37934257 4.40343219 4.42782208 4.4477206
 4.47635849 4.51594603 4.54294049 4.56689989 4.58804775 4.59919363
 4.63700551 4.66638214]