Bibliothèque Boost, comment obtenir un déterminant de lu_factorize()?

Question

J'essaie de calculer un déterminant en utilisant les bibliothèques boost C++. J'ai trouvé le code pour la fonction InvertMatrix() que j'ai copié ci-dessous. Chaque fois que je calcule cet inverse, je veux aussi le déterminant. J'ai une bonne idée de comment calculer, en multipliant la diagonale de la matrice U de la décomposition LU. Il y a un problème, je suis capable de calculer le déterminant correctement, sauf pour le signe. Selon le pivotement, je reçois le signe incorrect la moitié du temps. Quelqu'un at-il une suggestion sur la façon d'obtenir le bon signe à chaque fois? Merci d'avance.

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; 

Voici où j'ai inséré ma meilleure photo pour calculer le déterminant.

T determinant = 1; 

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

Terminez ma partie du code.

// 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 matrice de permutation pm contient les informations dont vous aurez besoin pour déterminer le changement de signe: vous voulez multiplier votre déterminant par le déterminant de la matrice de permutation. En examinant le fichier source lu.hpp, nous trouvons une fonction appelée swap_rows qui explique comment appliquer une matrice de permutation à une matrice. Il est facilement modifié pour obtenir le déterminant de la matrice de permutation (le signe de la permutation), étant donné que chaque échange réel contribue un facteur de -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; 
} 

Une autre alternative serait d'utiliser les méthodes lu_factorize et lu_substitute qui ne pas faire pivoter (consulter la source, mais en gros laisser tomber le pm dans les appels à lu_factorize et). Ce changement ferait en sorte que votre calcul déterminant fonctionne tel quel. Attention toutefois: supprimer le pivotement rendra l'algorithme moins stable numériquement.

Answer 2

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

il suffit de changer determinant *= A(i,i) à determinant *= (pm(i) == i ? 1 : -1) * A(i,i). J'ai essayé de cette façon et cela fonctionne.

Je sais, c'est en fait très similaire à la réponse de Managu et l'idée est la même, mais je crois que c'est plus simple (et "2 en 1" si utilisé dans la fonction InvertMatrix).