Obecnie próbuje tanio obliczyć szacunkową dobry rangi dla macierzy . Dlatego obliczam kolumnę obrotową dekompostację QR za pomocą
[Q,R,E]=qr(A)
w Matlabie. Oceniam rangę za pomocą
tol = size(A,n)*eps*norm(A,'fro');
r = sum(abs(diag(R))>tol)
Działa to dobrze, a wykres wszystkich przekątnych wpisów R wygląda następująco:
W przypadku przeniesienia całego algorytmu do C / Fortran I zastępuję [Q, R, E] = qr (A) za pomocą DGEQP3 z LAPACK, który również oblicza kolumnę rozkładającą QR rozkład. Ale jeśli użyję tego samego oszacowania dla rangi, najczęściej coś źle zrozumiem. Wygląda ten sam wykres dla wyprodukowanego z DGEQP3
Matryca wejściowa jest dokładnie taka sama dla obu eksperymentów.
Moje pytanie brzmi teraz: na jakiej funkcji LAPACK opiera się kolumna rozkładająca kod QR z Matlaba?
Dzięki za wszelką pomoc, Grisu
Edycja: DGEQPF daje ten sam zły wynik.
Edycja2:
- Matryca wejściowa jest gęsta i budowana jakoE + s i g n ( E , F )
- jest dostępna tutaj: http://www-e.uni-magdeburg.de/makoehle/A.mtx.gz (Format MatrixMarket)
- Źle : http://www-e.uni-magdeburg.de/makoehle/R_wrong.mtx.gz
- Użyłem LAPACK 3.4.0 z OpenBlas / GotoBLAS (64-bit)
- Matlab 7, 2007b, 2010b Linux 32bit
Edit3: - Za pomocą GDB dowiedziałem się, że Matlab 2010b wywołuje DGEQP3: # 3 0xaa46ce2f w dgeqp3_ () z /usr/ubuntu10.04/matlabr2010b/bin/glnx86/../../bin/glnx86/../. ./bin/glnx86/mllapack.so Dlaczego otrzymuję zły wynik za pomocą LAPACK (3.4.0 zawiera poprawki wymienione w uwadze roboczej 176)?
Odpowiedzi:
Istnieją tutaj dwie kwestie:
Gęsty czy rzadki?
DGEQP3
[Q,R,E] = qr(A)
Czy masz taki sam stos oprogramowania jak biblioteki wewnętrzne MATLAB?
Prawdopodobnie nie, co może być jednym z powodów uzyskania różnych wyników.
Zetknąłem się z tym problemem, gdy jednostka testowała pisaną przeze mnie bibliotekę wykorzystującą faktoryzacje QR. Użyłem MATLAB do prototypowania mojej pracy i uzyskałem inne wyniki niż przy użyciu LAPACK lub NumPy. O ile wiem, ponieważ MathWorks nie ułatwia znalezienia tych informacji, MATLAB używa wersji LAPACK nie wcześniejszej niż wersja 3.1.1 oraz biblioteki MKL BLAS Intela (dla Windows, Intel Mac i Linux) wersja 9.1 lub wyższy (patrz tutaj ). Nie mogłem znaleźć niczego na temat wersji SuiteSparse MATLAB. Przekopując się przez Internet lub przeglądając pliki bibliotek swojego systemu, możesz uzyskać dodatkowe informacje. Możesz spróbować zmienić biblioteki, do których prowadzi MATLAB, aby móc porównać z tymi samymi bibliotekami w różnych pakietach oprogramowania; Eric Chu zapewnia niezły opispokazuje przynajmniej, jak możesz zastąpić bibliotekę BLAS MATLAB-a własną (oczywiście robisz to na własne ryzyko). Sugeruje, że możesz również zrobić to samo z LAPACK. Może być nawet możliwe zastąpienie wersji SuiteSparse, której używa MATLAB, własną wersją.
W końcu użyłem NumPy do prototypowania moich wyników dla faktoryzacji QR, ponieważ używa systemowych bibliotek BLAS i LAPACK. NumPy i SciPy nie zastępują MATLAB-a, ponieważ dwie połączone biblioteki nie mają części funkcjonalności MATLAB, ale dla tego konkretnego zadania algebry liniowej Python + NumPy + SciPy + Matplotlib powinien działać dobrze.
źródło
internal.matlab.language.versionPlugins.blas
iinternal.matlab.language.versionPlugins.lapack
uzyskać wersje BLAS i LAPACKZobacz stronę Leslie Foster na temat oprogramowania ujawniającego rangę . Zobacz także tę notatkę roboczą LAPACK analizującą awarie ujawniającego rangę QR
xGEQP3
.Powinieneś być w stanie dowiedzieć się, jakich procedur używa MATLAB, ustawiając punkty przerwania w debuggerze i sprawdzając stos. Ostatnio patrzyłem, co prawda kilka lat temu, MATLAB używał bibliotek współdzielonych, w którym to przypadku nazwy symboli nie mogą zostać usunięte, więc zobaczysz nazwy funkcji na stosie wywołań (ale nie argumentów, ponieważ zdecydowanie nie przechowuje informacji debugujących).
źródło
xGEQP3
Algorytm nie jest całkowicie bezpieczny dla ujawniając rangę. Jeśli chcesz zagwarantować, że uzyskasz właściwy wynik, powinieneś użyć SVD lub bezpieczniejszego QR, takiego jakxGEQPX
lubxGEQPY
. Nie można oczekiwać, że niestabilny algorytm zwróci ten sam wynik na różnych architekturach lub w różnych implementacjach (MATLAB prawdopodobnie używa starszego LAPACKA).