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

12

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.

Spacekiller23
fonte
2
1 - Suponho que você pretendia colocar colchetes em torno de [evals, evecs] no código da oitava? 2 - você pode incluir um pequeno exemplo para M? ou talvez um script gerador para um, se isso for possível?
Nick J
1 - Sim, exatamente, editei minha postagem. 2 - Experimentei o desempenho para algumas sub-matrizes dos meus dados e parece que o Octave é sempre mais rápido, mas para matrizes menores abaixo de 5000x5000 é apenas um fator de 2-5 vezes, acima do que fica realmente feio. E desde os seus "dados reais" não posso dar um script gerador. Existe uma maneira padrão de enviar um exemplo de alguma forma? Devido à escassez, um arquivo npz é razoavelmente pequeno.
Spacekiller23
Eu acho que você pode compartilhar um link para qualquer instalação de armazenamento em nuvem.
Patol75
Valeu. Incluí um link de caixa de depósito na postagem original e atualizei o código para um exemplo de trabalho.
Spacekiller23
11
Para reforçar seu argumento, testei no Matlab R2019b e obtive 84 segundos vs 36,5 minutos no Python 3.7, Scipy 1.2.1 (26 vezes mais rápido).
Bill

Respostas:

1

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?


Aqui estão os autovalores do scipy-Arpack para vários ke tol, recebidos de arquivos de log em gist.github :

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.


De onde isso Avem, qual área problemática? Você pode gerar similar A, menor ou mais escasso?

Espero que isto ajude.

denis
fonte
Obrigado pela sua contribuição e desculpe pela resposta (muito) tardia. O projeto para o qual eu usei já está concluído, mas ainda estou curioso, então verifiquei. Infelizmente, os autovalores no Ocatve são diferentes, pois k = 10 acho [-2.5673e-16 -1.2239e-18 7.5420e-07 7.5622e-06 1.0189e-05 1.8725e-05 2.0265e-05 2.1568e- 05 4.2458e-05 5.1030e-05], que também é independente do valor de tolerância no intervalo 1e-5 a 1e-7. Portanto, há outra diferença aqui. Você não acha estranho que o scipy (incluindo sua sugestão) produza diferentes valores pequenos, dependendo do número de valores consultados?
Spacekiller23
@ Spacekiller23, este foi um bug, agora corrigido no scipy 1.4.1 (veja scipy / issues / 11198 ); você pode verificar sua versão? Também tolé complicado para pequenos autovalores - faça uma nova pergunta sobre isso, se você quiser, me avise.
denis
1

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 )

Anthony Gatti
fonte
Isso realmente mudou o desempenho para você? Isso seria uma surpresa para mim. Eu sabia o link que você postou e também especifiquei implicitamente qual = 'LM' no meu exemplo. Porque o valor padrão para um unset que é 'LM'. Ainda verifiquei, porém, e o desempenho não mudou no meu exemplo.
Spacekiller23
De fato, parece ter uma diferença semelhante à do Python para a oitava para você. Eu também tinha uma matriz grande que estava tentando decompor e acabei usando o eigsh (matriz, k = 7, que = 'LM', sigma = 1e-10). Originalmente, estava incorretamente especificando qual = 'SM' pensei que precisava fazer isso para obter os menores autovalores e que estava demorando uma eternidade para me dar uma resposta. Então, eu encontrei esse artigo e percebi que você só precisava especificá-lo para o 'LM' mais rápido, e definir k para o que você quisesse e isso aceleraria as coisas. Sua matriz é realmente eremita?
Anthony Gatti
0

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]
Patol75
fonte
Obrigado pela sua contribuição e feedback. Eu tentei algumas coisas para dar uma resposta decente aos seus pontos. 1. Minha tarefa em questão exige encontrar os k menores autovalores / vetores. Portanto, a abordagem usando sigma = 0 é dada nos documentos SciPy: docs.scipy.org/doc/scipy/reference/tutorial/arpack.html 2. Tentei mais algumas opções, que editei na pergunta original. 3. Pelo que entendi os documentários de Octave e SciPy, eigs (M, 6,0) ek = 6, simga = 0 deve ser o mesmo.
Spacekiller23
4. Como minha matriz é real e quadrada, pensei que não deveria haver diferença entre SA e SM como uma opção, mas obviamente existe, pelo menos na computação. Estou no caminho errado aqui? No geral, isso significa mais perguntas, mas nenhuma resposta ou solução real minha.
Spacekiller23