Próbowałem różnych rzadkich rozwiązań dostępnych w Pythonie 3 i porównując wydajność między nimi, a także przeciwko Octave i Matlab. Wybrałem zarówno bezpośrednie, jak i iteracyjne podejścia, wyjaśniam to bardziej szczegółowo poniżej.

Aby wygenerować odpowiednią rzadką matrycę, ze strukturą pasmową, problem Poissona rozwiązuje się przy użyciu elementów skończonych z kwadratowymi siatkami N = 250, N = 500 i N = 1000. Powoduje to wymiarach matrycy A = N ^ 2XN ^ 2 i wektora b = n ^ 2x1, tj. Największy NXN jest milion. Jeśli ktoś jest zainteresowany replikowaniem moich wyników, przesłałem matryce A i wektory B w poniższym ogniu (wygaśnie ona 30 dni) Pobierz systemy używane tutaj. Matryce są przechowywane w Trouble I, J, V, tj. Pierwsze dwie kolumny są odpowiednio wskaźnikami wierszy i kolumn, a trzecia kolumna są wartościami odpowiadającymi takim indeksom. Obserwuj, że istnieją pewne wartości w V, które są prawie zerowe, są celowe. Mimo to, że taśmowa struktura jest zachowana po poleceniu matrycy "szpiegowskiej" w Matlabie i Pythonie.

Dla porównania, użyłem następujących solvers:

Matlab i Octave, Solver Direct: The Canonical x=A\b.

Matlab i oktawę, Solver PCG: Wstępny koniugowany gradient, Solver PCG pcg(A,b,1e-5,size(b,1)) (nie jest używany (nie wstępny).

Scipy (Python), Solver Direct: linalg.spsolve(A, b) gdzie a jest wcześniej sformatowany w formacie csr_matrix.

Scipy (Python), Solver PCG: sp.linalg.cg(A, b, x0=None, tol=1e-05)

Scipy (Python), Umfpack Solver: spsolve(A, b) za pomocą from scikits.umfpack import spsolve. Ten solver jest najwyraźniej dostępny (tylko?) W Linux, ponieważ wykorzystuje libsuitesparse [Timothy Davis, Texas A & AMP; M]. W Ubuntu musi to najpierw zainstalować jako sudo apt-get install libsuitesparse-dev.

Ponadto wyżej wymienione solvers Pythona są testowane w:

  1. Okna.
  2. Linux.
  3. System operacyjny Mac.

Warunki:

  • Czas odbywa się tuż przed i po rozwiązaniu systemów. To, że nad głową czytanie macierzy nie jest brane pod uwagę.
  • Czas odbywa się dziesięć razy dla każdego systemu i obliczane jest średnia i odchylenie standardowe.

Sprzęt komputerowy:

  • Windows i Linux: Dell Intel (R) Core (TM) I7-8850H CPU @ 2,6 GB 2.59GHZ, 32 GB RAM DDR4.
  • Mac OS: MacBook Pro Retina Mid 2014 Intel (R) Quad-Core (TM) I7 2,2 GHz 16 GB RAM DDR3.

Wyniki: Czas rozwiązania rozmiar problemu

Obserwacje:

  • MATLAB A B jest najszybszy pomimo bycia w starszym komputerze.
  • Istnieją znaczne różnice między wersjami Linux i Windows. Patrz na przykład Solver Direct w NXN = 1E6. Jest to pomimo systemu Linux działa pod Windows (WSL).
  • Można mieć ogromny rozproszenie w Scipy Solvers. To jest, jeśli to samo rozwiązanie jest uruchomione kilka razy, jeden z czasów może po prostu zwiększyć więcej niż dwa razy.
  • Najszybszą opcją w Pythonie może być prawie cztery razy wolniej niż Matlab działa w bardziej ograniczonym sprzęcie. Naprawdę?

Jeśli chcesz odtworzyć testy, pozostawię tutaj bardzo proste skrypty. Dla Matlab / Octave:

IJS=load('KbN1M.txt');
b=load('FbN1M.txt');

I=IJS(:,1);
J=IJS(:,2);
S=IJS(:,3);

Neval=10;
tsparse=zeros(Neval,1);
tsolve_direct=zeros(Neval,1);
tsolve_sparse=zeros(Neval,1);
tsolve_pcg=zeros(Neval,1);
for i=1:Neval
    tic
    A=sparse(I,J,S);
    tsparse(i)=toc;
    tic
    x=A\b;
    tsolve_direct(i)=toc;        
    tic
    x2=pcg(A,b,1e-5,size(b,1));
    tsolve_pcg(i)=toc;
end

save -ascii octave_n1M_tsparse.txt tsparse
save -ascii octave_n1M_tsolvedirect.txt tsolve_direct
save -ascii octave_n1M_tsolvepcg.txt tsolve_pcg

Dla Pythona:

import time
from scipy import sparse as sp
from scipy.sparse import linalg
import numpy as np
from scikits.umfpack import spsolve, splu #NEEDS LINUX


b=np.loadtxt('FbN1M.txt')
triplets=np.loadtxt('KbN1M.txt')

I=triplets[:,0]-1
J=triplets[:,1]-1
V=triplets[:,2]

I=I.astype(int)
J=J.astype(int)
NN=int(b.shape[0])

Neval=10
time_sparse=np.zeros((Neval,1))
time_direct=np.zeros((Neval,1))
time_conj=np.zeros((Neval,1))
time_umfpack=np.zeros((Neval,1))
for i in range(Neval):
    t = time.time()
    A=sp.coo_matrix((V, (I, J)), shape=(NN, NN))
    A=sp.csr_matrix(A)
    time_sparse[i,0]=time.time()-t
    t = time.time()
    x=linalg.spsolve(A, b)
    time_direct[i,0] = time.time() - t
    t = time.time()
    x2=sp.linalg.cg(A, b, x0=None, tol=1e-05)
    time_conj[i,0] = time.time() - t
    t = time.time()
    x3 = spsolve(A, b) #ONLY IN LINUX
    time_umfpack[i,0] = time.time() - t

np.savetxt('pythonlinux_n1M_tsparse.txt',time_sparse,fmt='%.18f')
np.savetxt('pythonlinux_n1M_tsolvedirect.txt',time_direct,fmt='%.18f')
np.savetxt('pythonlinux_n1M_tsolvepcg.txt',time_conj,fmt='%.18f')
np.savetxt('pythonlinux_n1M_tsolveumfpack.txt',time_umfpack,fmt='%.18f')

Czy istnieje sposób na dalsze poprawę rzadkich czasów rozwiązań za pomocą Pythona? lub przynajmniej być w podobnym porządku wydajności jako matlab? Jestem otwarty na sugestie za pomocą C / C ++ lub Fortran i opakowania dla Pythona, ale wierzę, że nie będzie już znacznie lepszy niż wybór UMFPack. Sugestie są bardzo mile widziane.

Str.s. Jestem świadomy poprzednich postów, np. Scipy Powolne rzadki Matrix Solver problemy z systemem Scipipy.sparse.Linalg System liniowy Solvers . . Jak korzystać z NUMBA, aby przyspieszyć rzadkie solvery liniowe w Pythonie, które są dostarczane w Scipy.sparse.Linalg? Ale myślę, że żaden nie jest tak wyczerpujący jak ten, podkreślając jeszcze więcej problemów między systemami operacyjnymi podczas korzystania z bibliotek Python.

Edytuj_1: Dodam nową działkę z wynikami za pomocą solvera QR z Intel MKL za pomocą opakowania Pythona, jak sugeruje w komentarzach. Jest to jednak nadal za występem Matlab. Aby to zrobić, trzeba dodać:

from sparse_dot_mkl import sparse_qr_solve_mkl

I

sparse_qr_solve_mkl(A.astype(np.float32), b.astype(np.float32))

Do skryptów podanych w oryginalnym postie. Można pominąć ".Syntype (np.float32)", a wydajność pobiera Grorelty (około 10%) dla tego systemu. Dodano MKL Solver, Umfpack Podświetlony

5
uom0 17 październik 2020, 13:37

1 odpowiedź

Najlepsza odpowiedź

Spróbuję odpowiedzieć na siebie. Aby zapewnić odpowiedź, próbowałem jeszcze bardziej wymagającego przykładu, z matrycą wielkości (N, N) około pół miliona o pół miliona i odpowiedniego wektora (N, 1). Jest to jednak znacznie mniej rzadkie (bardziej gęste) niż jeden podany w pytaniu. Matryca ta przechowywana w ASCII wynosi około 1,7 GB, w porównaniu z jednym z przykładów, który ma około 0,25 GB (pomimo "rozmiar" jest większa). Zobaczyć jego kształt tutaj,

enter image description here

Następnie próbowałem rozwiązać AX = B, używając ponownie Matlab, Octave i Pythona przy użyciu wyżej wymienionych solidów z Scipy, Intel MKL Wrapper, Umfpack z Tim Davis. Moim pierwszym zaskoczeniu jest to, że zarówno Matlab, jak i Octave mogą rozwiązać systemy za pomocą A B, który nie jest dla pewności, że jest to bezpośredni solver, ponieważ wybiera najlepszy solver na podstawie właściwości macierzy, patrz Matlab's X = A b. Jednakże, Python's linalg.spsolve, opakowanie MKL i UMFPack wyrzucały błędy z pamięci w systemie Windows i Linux. W Mac, linalg.spsolve był w jakiś sposób obliczający rozwiązanie, a Alghouth był z bardzo słabą wydajnością, nigdy przez błędy pamięci. Zastanawiam się, czy pamięć jest obsługiwana inaczej w zależności od systemu operacyjnego. Dla mnie wydaje się, że mac zamienił pamięć do dysku twardego, a nie używa go z pamięci RAM. Wydajność solvera CG w Pythonie była raczej słaba, w porównaniu z Matlab. Jednak w celu poprawy wydajności w Solver CG w Pythonie można uzyskać ogromną poprawę wydajności, jeśli obliczany jest A = 0,5 (A + A ") jest obliczany (jeśli jest oczywiście, mieć system symetryczny). Korzystanie z warunku wstępnego w Pythonie nie pomogło. Próbowałem użyć metody sp.linalg.spilu wraz z sp.linalg.LinearOperator, aby obliczyć wstępny wstępny, ale wydajność była raczej słaba. W MATLAB można użyć niekompletnej rozkładu Choliky.

W przypadku problemu poza pamięcią roztwór miał zastosowanie rozkładu LU i rozwiązać dwa zagnieżdżone układy, takie jak AX = B, A = LL ', Y = L i X = Y L' '.

Włożyłem tu min. Czasy rozwiązania,

Matlab mac, A\b = 294 s.
Matlab mac, PCG (without conditioner)= 17.9 s.
Matlab mac, PCG (with incomplete Cholesky conditioner) = 9.8 s.
Scipy mac, direct = 4797 s.
Octave, A\b = 302 s.
Octave, PCG (without conditioner)= 28.6 s.
Octave, PCG (with incomplete Cholesky conditioner) = 11.4 s.
Scipy, PCG (without A=0.5(A+A'))= 119 s.
Scipy, PCG (with A=0.5(A+A'))= 12.7 s.
Scipy, LU decomposition using UMFPACK (Linux) = 3.7 s total.

Odpowiedź brzmi więc tak, istnieją sposoby na poprawę czasów rozwiązań w Scipy. Wykorzystanie opakowań do Umfpack (Linux) lub Intel MKL QR Solver jest wysoce zalecane, jeśli mimigłowa stacji roboczej pozwala mu. W przeciwnym razie wykonywanie A = 0,5 (A + A ") przed użyciem koniugatu Solver Gradient może mieć pozytywny wpływ na wydajność roztworu, jeśli ktoś zajmuje się systemami symetrycznymi. Daj mi znać, jeśli ktoś byłby zainteresowany posiadaniem tego nowego systemu, więc mogę go przesłać.

4
uom0 20 październik 2020, 12:56