Menemukan vektor eigen terkecil dari matriks jarang besar, lebih dari 100x lebih lambat di SciPy daripada di Octave

12

Saya mencoba untuk menghitung beberapa (5-500) vektor eigen yang sesuai dengan nilai eigen terkecil dari matriks persegi jarang simetris (hingga 30000x30000) dengan kurang dari 0,1% dari nilai yang bukan nol.

Saat ini saya menggunakan scipy.sparse.linalg.eigsh dalam mode shift-invert (sigma = 0,0), yang saya temukan melalui berbagai posting pada topik adalah solusi yang lebih disukai. Namun, dibutuhkan hingga 1 jam untuk menyelesaikan masalah dalam kebanyakan kasus. Di sisi lain fungsinya sangat cepat, jika saya meminta nilai eigen terbesar (sub detik pada sistem saya), yang diharapkan dari dokumentasi.

Karena saya lebih akrab dengan Matlab dari pekerjaan, saya mencoba memecahkan masalah dalam Oktaf, yang memberi saya hasil yang sama menggunakan eigs (sigma = 0) hanya dalam beberapa detik (sub 10s). Karena saya ingin melakukan sapuan parameter dari algoritma termasuk perhitungan vektor eigen, perolehan waktu semacam itu akan bagus untuk dimiliki dalam python juga.

Saya pertama kali mengubah parameter (terutama toleransi), tetapi itu tidak banyak berubah pada rentang waktu.

Saya menggunakan Anaconda di Windows, tetapi mencoba untuk mengganti LAPACK / BLAS yang digunakan oleh Scipy (yang sangat menyakitkan) dari mkl (default Anaconda) ke OpenBlas (digunakan oleh Octave menurut dokumentasi), tetapi tidak dapat melihat perubahan dalam kinerja.

Saya tidak dapat menemukan, apakah ada sesuatu yang berubah tentang ARPACK yang digunakan (dan bagaimana)?

Saya mengunggah testcase untuk kode di bawah ini ke folder dropbox berikut: https://www.dropbox.com/sh/l6aa6izufzyzqr3/AABqij95hZOvRpnnjRaETQmka?dl=0

Dengan 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)

Dalam oktaf:

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

Bantuan apa pun appriciated!

Beberapa opsi tambahan yang saya coba berdasarkan pada komentar dan saran:

Oktaf: eigs(M,6,0)dan eigs(M,6,'sm')beri saya hasil yang sama:

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

saat eigs(M,6,'sa',struct('tol',2))konvergen ke

[1.0423 2.7604 6.1548 11.1310 18.0207 25.3933] 

jauh lebih cepat, tetapi hanya jika nilai-nilai toleransi di atas 2, jika tidak ia tidak menyatu sama sekali dan nilainya sangat berbeda.

Python: eigsh(M,k=6,which='SA')dan eigsh(M,k=6,which='SM')keduanya tidak konvergen (kesalahan ARPACK tanpa konvergensi tercapai). Hanya eigsh(M,k=6,sigma=0.0)memberikan beberapa nilai eigen (setelah hampir satu jam), yang berbeda dengan satu oktaf untuk yang terkecil (bahkan 1 nilai kecil tambahan ditemukan):

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

Jika toleransi cukup tinggi saya juga mendapatkan hasil dari eigsh(M,k=6,which='SA',tol='1'), yang mendekati nilai yang diperoleh lainnya

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

lagi dengan jumlah nilai eigen kecil yang berbeda. Waktu komputasi masih hampir 30 menit. Sementara perbedaan nilai yang sangat kecil mungkin dapat dimengerti, karena mereka mungkin mewakili kelipatan 0, keberagaman yang berbeda membuat saya bingung.

Selain itu, tampaknya ada beberapa perbedaan mendasar dalam SciPy dan Octave, yang saya belum bisa mengetahuinya.

Spacekiller23
sumber
2
1 - Saya berasumsi Anda bermaksud menempatkan tanda kurung di sekitar [evals, evecs] dalam kode oktaf? 2 - dapatkah Anda memasukkan contoh kecil untuk M? atau mungkin skrip generator untuk satu jika itu mungkin?
Nick J
1 - Ya persis, saya mengedit posting saya. 2 - Saya mencoba kinerja untuk beberapa sub-matriks data saya dan sepertinya Oktaf selalu lebih cepat, tetapi untuk matriks yang lebih kecil di bawah 5000x5000 hanya faktor 2-5 kali, di atas itu menjadi sangat jelek. Dan karena itu "data nyata" saya tidak bisa memberikan skrip generator. Apakah ada cara standar untuk mengunggah contoh? Karena sparsity, file npz cukup kecil.
Spacekiller23
Saya kira Anda dapat membagikan tautan ke fasilitas penyimpanan cloud apa pun.
Patol75
Terima kasih. Saya menyertakan tautan dropbox di pos asli dan memperbarui kode ke contoh yang berfungsi.
Spacekiller23
1
Juts untuk memperkuat poin Anda, saya diuji pada Matlab R2019b dan mendapat 84 detik vs 36,5 menit dalam Python 3.7, Scipy 1.2.1 (26 kali lebih cepat).
Bill

Jawaban:

1

Dugaan dan beberapa komentar, karena saya tidak punya Matlab / Oktaf:

Untuk menemukan nilai eigen kecil dari matriks simetris dengan nilai eigen> = 0, seperti nilai Anda, berikut ini cara yang lebih cepat daripada 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 )untuk pasangan eigen besar lebih cepat daripada shift-invert untuk kecil, karena A * xlebih cepat daripada solve()shift-invert yang harus dilakukan. Matlab / Oktaf mungkin dapat melakukan ini Aflipsecara otomatis, setelah tes cepat untuk positif-pasti dengan Cholesky.
Bisakah Anda menjalankan eigsh( Aflip )di Matlab / Oktaf?

Faktor lain yang dapat mempengaruhi akurasi / kecepatan:

Default Arpack untuk vektor awal v0adalah vektor acak. Saya menggunakan v0 = np.ones(n), yang mungkin mengerikan untuk beberapa Atetapi dapat direproduksi :)

Ini Amatrix hampir persis sigular, A * ones~ 0.

Multicore: scipy-arpack dengan openblas / Lapack menggunakan ~ 3,9 dari 4 core di iMac saya; apakah Matlab / Octave menggunakan semua core?


Berikut adalah nilai-nilai eigen scipy-Arpack untuk beberapa kdan tol, diambil dari logfile di bawah 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 

Apakah Matlab / Oktaf hampir sama? Jika tidak, semua taruhan dimatikan - periksa dulu kebenarannya, lalu kecepatan.

Mengapa nilai eigen goyah begitu banyak? Tiny <0 untuk matriks yang seharusnya non-negatif-pasti adalah tanda kesalahan pembulatan , tetapi trik biasa dari perubahan kecil A += n * eps * sparse.eye(n), tidak membantu.


Dari mana ini Aberasal, bidang masalah apa? Bisakah Anda menghasilkan yang serupa A, lebih kecil atau lebih jarang?

Semoga ini membantu.

denis
sumber
Terima kasih atas masukan Anda dan maaf atas keterlambatan (sangat). Proyek yang saya gunakan untuk ini sudah selesai, tapi saya masih penasaran, jadi saya memeriksa. Sayangnya, nilai eigen di Ocatve berubah berbeda, untuk k = 10 saya menemukan [-2.5673e-16 -1.2239e-18 7.5420e-07 7.5622e-06 1.0189e-05 1.8725e-05 2.0265e-05 2.268e-05 2.1568e- 05 4.2458e-05 5.1030e-05] yang juga tidak tergantung pada nilai toleransi di kisaran 1e-5 hingga 1e-7. Jadi ada perbedaan lain di sini. Tidakkah menurut Anda aneh bahwa Scipy (termasuk saran Anda) menghasilkan nilai kecil yang berbeda tergantung pada jumlah nilai yang ditanyakan?
Spacekiller23
@ Spacekiller23, ini adalah bug, sekarang diperbaiki di scipy 1.4.1 (lihat scipy / issues / 11198 ); Bisakah Anda memeriksa versinya? Juga tolberantakan untuk nilai eigen kecil - ajukan pertanyaan baru jika Anda suka, beri tahu saya.
denis
1

Saya tahu ini sudah tua sekarang, tetapi saya memiliki masalah yang sama. Apakah Anda mengulas di sini ( https://docs.scipy.org/doc/scipy/reference/tutorial/arpack.html )?

Sepertinya ketika Anda menetapkan sigma ke angka rendah (0) Anda harus mengatur yang = 'LM', meskipun Anda ingin nilai rendah. Ini karena pengaturan sigma mengubah nilai yang Anda inginkan (rendah dalam kasus ini) menjadi tinggi dan sehingga Anda masih dapat memanfaatkan metode 'LM', yang jauh lebih cepat untuk mendapatkan apa yang Anda inginkan (nilai eigen rendah ).

Anthony Gatti
sumber
Apakah ini benar-benar mengubah kinerja untuk Anda? Itu akan menjadi kejutan bagi saya. Saya tahu tautan yang Anda kirim dan saya juga secara implisit menentukan = 'LM' dalam contoh saya. Karena nilai default untuk unset yaitu 'LM'. Saya masih memeriksa, dan kinerjanya tidak berubah untuk contoh saya.
Spacekiller23
Memang, tampaknya memiliki perbedaan yang sama seperti dari Python ke oktaf untuk Anda. Saya juga punya matriks besar yang saya coba dekomposisi dan akhirnya menggunakan eigsh (matrix, k = 7, yang = 'LM', sigma = 1e-10). Awalnya saya salah menentukan = pemikiran 'SM' yang perlu saya lakukan untuk mendapatkan nilai eigen terkecil dan perlu waktu lama untuk memberi saya jawaban. Kemudian, saya menemukan artikel itu dan menyadari bahwa Anda hanya perlu menentukannya ke 'LM' yang lebih cepat, dan menetapkan k menjadi apa pun yang Anda inginkan dan itu akan mempercepat segalanya. Apakah matriks Anda sebenarnya hermitian?
Anthony Gatti
0

Saya ingin mengatakan pertama bahwa saya tidak tahu mengapa hasil yang Anda dan @Bill laporkan adalah sebagaimana adanya. Saya hanya ingin tahu apakah eigs(M,6,0)dalam Octave sesuai dengan k=6 & sigma=0, atau mungkin itu sesuatu yang lain?

Dengan ceroboh, jika saya tidak memberikan sigma, saya bisa mendapatkan hasil dalam waktu yang layak dengan cara ini.

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)

Saya sepenuhnya tidak yakin apakah ini masuk akal.

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]

Satu-satunya cara saya menemukan menggunakan sigma dan untuk mendapatkan hasil dalam waktu yang layak adalah dengan menyediakan M sebagai LinearOperator. Saya tidak terlalu terbiasa dengan hal ini, tetapi dari apa yang saya pahami implementasi saya mewakili matriks identitas, yang seharusnya menjadi M jika tidak ditentukan dalam panggilan. Alasan untuk ini adalah bahwa alih-alih melakukan penyelesaian langsung (dekomposisi LU), Scipy akan menggunakan pemecah iteratif, yang berpotensi lebih cocok. Sebagai perbandingan, jika Anda memberikan M = np.identity(a.shape[0]), yang harus sama persis, maka eigsh membutuhkan waktu lama untuk membuahkan hasil. Perhatikan bahwa pendekatan ini tidak berfungsi jika sigma=0disediakan. Tetapi saya tidak yakin apakah sigma=0ini benar-benar bermanfaat?

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)))

Sekali lagi, tidak tahu apakah itu benar tetapi pasti berbeda dari sebelumnya. Itu akan bagus untuk mendapatkan masukan dari seseorang dari cipy.

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
sumber
Terima kasih atas masukan dan umpan balik Anda. Saya mencoba beberapa hal untuk memberikan jawaban yang layak untuk poin Anda. 1. Tugas saya saat ini membutuhkan menemukan nilai / vektor eigen terkecil. Karenanya pendekatan yang menggunakan sigma = 0 bahkan diberikan dalam dokumen SciPy: docs.scipy.org/doc/scipy/reference/tutorial/arpack.html 2. Saya mencoba beberapa opsi lagi, yang saya edit ke pertanyaan awal. 3. Ketika saya memahami film dokumenter dari Octave dan SciPy, eigs (M, 6,0) dan k = 6, simga = 0 harus sama.
Spacekiller23
4. Karena matriks saya nyata dan persegi, saya pikir seharusnya tidak ada perbedaan antara SA dan SM sebagai pilihan, tetapi jelas ada, setidaknya dalam perhitungan. Apakah saya di jalur yang salah di sini? Secara keseluruhan itu berarti lebih banyak pertanyaan dan tetapi tidak ada jawaban atau solusi nyata dari saya.
Spacekiller23