Boost Library, comment obtenir le déterminant de lu_factorize ()?

Question

J'essaie de calculer un déterminant à l'aide des 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 la façon de 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. En fonction du pivot, le signe est incorrect la moitié du temps. Quelqu'un a-t-il une suggestion sur la manière d'obtenir le bon panneau à 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;

C’est ici que j’ai inséré mon meilleur atout 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 nécessaires pour déterminer le changement de signe: vous souhaitez multiplier votre déterminant par celui de la matrice de permutation.

Parcourir 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 donner 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 solution consisterait à utiliser les méthodes lu_factorize et lu_substitute qui ne pivotent pas (consultez le code source, mais supprimez essentiellement le pm dans les appels à lu_factorize et lu_substitute ). Ce changement rendrait votre calcul déterminant déterminant tel quel. Attention cependant: supprimer le pivotement rend l'algorithme moins stable numériquement.

Answer 2

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

Remplacez simplement déterminant * = A (i, i) par déterminant * = (pm (i) == i? 1: -1) * A (i, i) . J'ai essayé de cette façon et ça marche.

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