Cara menggunakan CSR_MATRIX pada Python

Here is performance comparison of the three most upvoted answers using Jupyter notebook. The input is a 1M x 100K random sparse matrix with density 0.001, containing 100M non-zero values:

from scipy.sparse import random
matrix = random(1000000, 100000, density=0.001, format='csr')

matrix
<1000000x100000 sparse matrix of type '<type 'numpy.float64'>'
with 100000000 stored elements in Compressed Sparse Row format>

io.mmwrite / io.mmread

from scipy.sparse import io

%time io.mmwrite('test_io.mtx', matrix)
CPU times: user 4min 37s, sys: 2.37 s, total: 4min 39s
Wall time: 4min 39s

%time matrix = io.mmread('test_io.mtx')
CPU times: user 2min 41s, sys: 1.63 s, total: 2min 43s
Wall time: 2min 43s    

matrix
<1000000x100000 sparse matrix of type '<type 'numpy.float64'>'
with 100000000 stored elements in COOrdinate format>    

Filesize: 3.0G.

(note that the format has been changed from csr to coo).

np.savez / np.load

import numpy as np
from scipy.sparse import csr_matrix

def save_sparse_csr(filename, array):
    # note that .npz extension is added automatically
    np.savez(filename, data=array.data, indices=array.indices,
             indptr=array.indptr, shape=array.shape)

def load_sparse_csr(filename):
    # here we need to add .npz extension manually
    loader = np.load(filename + '.npz')
    return csr_matrix((loader['data'], loader['indices'], loader['indptr']),
                      shape=loader['shape'])


%time save_sparse_csr('test_savez', matrix)
CPU times: user 1.26 s, sys: 1.48 s, total: 2.74 s
Wall time: 2.74 s    

%time matrix = load_sparse_csr('test_savez')
CPU times: user 1.18 s, sys: 548 ms, total: 1.73 s
Wall time: 1.73 s

matrix
<1000000x100000 sparse matrix of type '<type 'numpy.float64'>'
with 100000000 stored elements in Compressed Sparse Row format>

Filesize: 1.1G.

cPickle

import cPickle as pickle

def save_pickle(matrix, filename):
    with open(filename, 'wb') as outfile:
        pickle.dump(matrix, outfile, pickle.HIGHEST_PROTOCOL)
def load_pickle(filename):
    with open(filename, 'rb') as infile:
        matrix = pickle.load(infile)    
    return matrix    

%time save_pickle(matrix, 'test_pickle.mtx')
CPU times: user 260 ms, sys: 888 ms, total: 1.15 s
Wall time: 1.15 s    

%time matrix = load_pickle('test_pickle.mtx')
CPU times: user 376 ms, sys: 988 ms, total: 1.36 s
Wall time: 1.37 s    

matrix
<1000000x100000 sparse matrix of type '<type 'numpy.float64'>'
with 100000000 stored elements in Compressed Sparse Row format>

Filesize: 1.1G.

Note: cPickle does not work with very large objects (see this answer). In my experience, it didn't work for a 2.7M x 50k matrix with 270M non-zero values. np.savez solution worked well.

Conclusion

(based on this simple test for CSR matrices) cPickle is the fastest method, but it doesn't work with very large matrices, np.savez is only slightly slower, while io.mmwrite is much slower, produces bigger file and restores to the wrong format. So np.savez is the winner here.

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.