Jaka jest odpowiednia funkcja LAPACK za Matlabem [Q, R, E] = qr (A)?

12

Obecnie próbuje tanio obliczyć szacunkową dobry rangi dla macierzy . Dlatego obliczam kolumnę obrotową dekompostację QR za pomocąA

[Q,R,E]=qr(A)

w Matlabie. Oceniam rangę za pomocąA

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: plot (sort (abs (diag (R)), 1, 'descend'), 'r *')

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 R

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:

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

MK alias Grisu
źródło
Czy możesz sprowokować to samo zachowanie za pomocą mniejszej matrycy, którą możesz udostępnić tutaj?
GertVdE
A
A
Grisu - Chciałbym spojrzeć na twoją matrycę. Jednak link www-e.uni-magdeburg.de/makoehle/A.mtx.gz nie jest już aktywny (w każdym razie w obecnej chwili). Czy masz bieżący link do matrycy? Dzięki, Les Foster
@LeslieFoster - witamy w scicomp!
Aron Ahmadia,

Odpowiedzi:

7

Istnieją tutaj dwie kwestie:

  • A
  • Czy masz taki sam stos oprogramowania jak biblioteki wewnętrzne MATLAB?

Gęsty czy rzadki?

ADGEQP3[Q,R,E] = qr(A)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ą.

A=QRQR

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.

Geoff Oxberry
źródło
Wydaje mi się, że uzyskanie takiego samego stosu oprogramowania jak Matlab jest niemożliwe. Używanie innego środowiska do prototypowania również nie jest pożądane. Problem polega na tym, że kod działa poprawnie w Matlabie, ale nie w C.
MK alias Grisu
@Grisu: Myślę, że bardzo trudno byłoby uzyskać ten sam stos oprogramowania, bez próby łączenia się w ich bibliotekach. To, co mnie myli, to to, skąd wiesz, że wynik w MATLAB jest poprawny, a wynik w C jest nieprawidłowy. Czy to jakaś matryca testowa o znanych właściwościach? Co więcej, AronAhmadia ma rację; poza replikacją architektury i stosu oprogramowania nie można oczekiwać takich samych wyników przy niestabilnym algorytmie. Zasadniczo powiedziano mi to samo na forach MATLAB dwa lata temu.
Geoff Oxberry
Moim zdaniem QR nie jest niestabilny. Nie mogę bezpośrednio sprawdzić, który rozkład QR jest prawidłowy, ale na podstawie rangi i macierzy Q obliczam projektor i tam mogę łatwo sprawdzić, czy dobre lub złe wyniki i ten z Matlaba są dobre. Ale próbuję połączyć się z bibliotekami Matlab.
MK alias Grisu
@Grisu: Istnieje wyraźna różnica między dobrymi wynikami a poprawnymi wynikami. Niedawno zaimplementowałem niepoprawne obliczenia, dzięki którym moje wyniki wyglądają wspaniale. Niemniej jednak obliczenia były nieprawidłowe, a prawidłowe obliczenia sprawiły, że moje wyniki wyglądały mniej imponująco (ale na szczęście pokazuje, że moje wyniki są poprawne). Czy próbujesz obliczyć rzutnik ortogonalny lub rzutnik ukośny? (Pytam, ponieważ znaczna część mojej pracy dotyczy skośnych projektorów.)
Geoff Oxberry
1
@GeoffOxberry: fwiw, na mojej wersji MATLAB, mogę zadzwonić internal.matlab.language.versionPlugins.blasi internal.matlab.language.versionPlugins.lapackuzyskać wersje BLAS i LAPACK
Amro
6

Zobacz stronę Leslie Foster na temat oprogramowania ujawniającego rangę . Zobacz także tę notatkę roboczą LAPACK analizującą awarie ujawniającego rangę QRxGEQP3 .

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

Jed Brown
źródło
Strona oprogramowania ujawniającego rangę nie pomogła. Procedura RRQR opisana jako pierwsza była moim pomysłem, ale daje jeszcze gorsze wyniki niż idea przestawiania kolumny.
MK alias Grisu
2
@Grisu - To powinno ci pomóc. xGEQP3Algorytm 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 jak xGEQPXlub xGEQPY. 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).
Aron Ahmadia
Wiem, że GEQP3 nie ujawnia rangi, ale daje więcej poprawnych wyników niż podprogramy RRQR. Używanie SVD jest zbyt drogie w moim zewnętrznym algorytmie. Porozmawiam również z jednym z autorów LAWN-176 i uważa, że ​​ten błąd nie jest objęty błędem. Próbowałem również DGEQPF / DGEQP3 z LAPACK 3.0.0 z tymi samymi wynikami.
MK alias Grisu