Boost Library, come ottenere il determinante da lu_factorize()?

Question

Sto provando a calcolare un determinante usando le librerie boost C++. Ho trovato il codice per la funzione InvertMatrix() che ho copiato di seguito. Ogni volta che calcolo questo inverso, voglio anche il determinante. Ho una buona idea di come calcolare, moltiplicando la diagonale della matrice U dalla scomposizione LU. C'è un problema, sono in grado di calcolare correttamente il determinante, ad eccezione del segno. A seconda del pivoting, il segno è sbagliato per metà del tempo. Qualcuno ha un suggerimento su come ottenere il segno giusto ogni volta? Grazie in anticipo.

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; 

Qui è dove ho inserito il mio colpo migliore per calcolare il determinante.

T determinant = 1; 

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

Terminare la mia parte del codice.

// 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

La permutazione matrice pm contiene le informazioni necessarie per determinare il cambiamento di segno: si vorrà moltiplicare il fattore determinante per il determinante della matrice di permutazione.

Per trovare il file sorgente lu.hpp troviamo una funzione denominata swap_rows che indica come applicare una matrice di permutazione a una matrice. È facilmente modificato per produrre il determinante della matrice di permutazione (il segno della permutazione), dato che ogni scambio effettivo contribuisce un fattore di -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; 
} 

Un'altra alternativa sarebbe quella di utilizzare i lu_factorize e lu_substitute metodi non fare alcun pivoting (consultare la fonte, ma in pratica rilasciare il pm nelle chiamate a lu_factorize e lu_substitute). Questo cambiamento renderebbe il tuo calcolo determinante funzionare così com'è. Fai attenzione, tuttavia: la rimozione del pivot renderà l'algoritmo meno numericamente stabile.

Answer 2

Secondo http://qiangsong.wordpress.com/2011/07/16/lu-factorisation-in-ublas/:

Basta cambiare determinant *= A(i,i) a determinant *= (pm(i) == i ? 1 : -1) * A(i,i). Ho provato in questo modo e funziona.

Lo so, che in realtà è molto simile alla risposta di Managu e l'idea è la stessa, ma credo che sia più semplice (e "2 in 1" se usato nella funzione InvertMatrix).