Boost Library, jak uzyskać wyznacznik z lu_factorize()?

Question

Próbuję obliczyć wyznacznik za pomocą bibliotek boost C++. Znalazłem kod funkcji InvertMatrix(), którą skopiowałem poniżej. Za każdym razem, gdy obliczam to odwrotnie, chcę także wyznacznika. Mam dobry pomysł, jak obliczyć, pomnażając przekątną macierzy U z dekompozycji LU. Jest jeden problem, jestem w stanie właściwie obliczyć wyznacznik, za wyjątkiem znaku. W zależności od obrotu otrzymam znak niepoprawny przez połowę czasu. Czy ktoś ma sugestię, jak za każdym razem uzyskać odpowiedni znak? Z góry dziękuję.

template<class T> 
bool InvertMatrix(const ublas::matrix<T>& input, ublas::matrix<T>& inverse) 
{ 
using namespace boost::numeric::ublas; 
typedef permutation_matrix<std::size_t> pmatrix; 
// create a working copy of the input 
matrix<T> A(input); 
// create a permutation matrix for the LU-factorization 
pmatrix pm(A.size1()); 

// perform LU-factorization 
int res = lu_factorize(A,pm); 
if(res != 0) return false; 

Oto gdzie wstawiłem mój najlepszy strzał przy obliczaniu wyznacznika.

T determinant = 1; 

for(int i = 0; i < A.size1(); i++) 
{ 
    determinant *= A(i,i); 
} 

Zakończ moją część kodu.

// create identity matrix of "inverse" 
inverse.assign(ublas::identity_matrix<T>(A.size1())); 

// backsubstitute to get the inverse 
lu_substitute(A, pm, inverse); 

return true; 
} 

Answer 1

Permutacja matryca pm zawiera informacje potrzebne do określenia zmiany znak: będziemy chcieli pomnożyć wyznacznik przez wyznacznik macierzy permutacji.

Przeglądając plik źródłowy lu.hpp znajdujemy funkcję o nazwie swap_rows, która mówi, jak zastosować macierz permutacji do macierzy. Jest łatwo modyfikować w celu wytworzenia wyznacznikiem macierzy permutacji (znak permutacji), biorąc pod uwagę, że każdy rzeczywisty wymiany przyczynia współczynnik 1:

template <typename size_type, typename A> 
int determinant(const permutation_matrix<size_type,A>& pm) 
{ 
    int pm_sign=1; 
    size_type size=pm.size(); 
    for (size_type i = 0; i < size; ++i) 
     if (i != pm(i)) 
      pm_sign* = -1; // swap_rows would swap a pair of rows here, so we change sign 
    return pm_sign; 
} 

Inną alternatywą byłoby użyć metod lu_factorize i lu_substitute które nie wykonuj żadnego obrotu (skonsultuj się ze źródłem, ale po prostu upuść pm w rozmowach na lu_factorize i lu_substitute). Ta zmiana spowodowałaby, że obliczenia wyznaczające działają tak jak są. Bądź jednak ostrożny: usunięcie obrotu spowoduje, że algorytm będzie mniej stabilny liczbowo.

Answer 2

Według http://qiangsong.wordpress.com/2011/07/16/lu-factorisation-in-ublas/:

Wystarczy zmienić determinant *= A(i,i) do determinant *= (pm(i) == i ? 1 : -1) * A(i,i). Próbowałem w ten sposób i to działa.

Wiem, że w rzeczywistości jest bardzo podobny do odpowiedzi Managu i pomysł jest taki sam, ale uważam, że jest prostszy (i "2 w 1", jeśli jest używany w funkcji InvertMatrix).