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.
sumber
Jawaban:
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:
eigsh( Aflip )
untuk pasangan eigen besar lebih cepat daripada shift-invert untuk kecil, karenaA * x
lebih cepat daripadasolve()
shift-invert yang harus dilakukan. Matlab / Oktaf mungkin dapat melakukan iniAflip
secara 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
v0
adalah vektor acak. Saya menggunakanv0 = np.ones(n)
, yang mungkin mengerikan untuk beberapaA
tetapi dapat direproduksi :)Ini
A
matrix 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
k
dantol
, diambil dari logfile di bawah gist.github :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
A
berasal, bidang masalah apa? Bisakah Anda menghasilkan yang serupaA
, lebih kecil atau lebih jarang?Semoga ini membantu.
sumber
tol
berantakan untuk nilai eigen kecil - ajukan pertanyaan baru jika Anda suka, beri tahu saya.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 ).
sumber
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 dengank=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.
Saya sepenuhnya tidak yakin apakah ini masuk akal.
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 jikasigma=0
disediakan. Tetapi saya tidak yakin apakahsigma=0
ini benar-benar bermanfaat?Sekali lagi, tidak tahu apakah itu benar tetapi pasti berbeda dari sebelumnya. Itu akan bagus untuk mendapatkan masukan dari seseorang dari cipy.
sumber