Boost Library, ¿cómo obtener determinante de lu_factorize ()?

Question

Estoy tratando de calcular un determinante usando las bibliotecas boost c ++. Encontré el código para la función InvertMatrix () que he copiado a continuación. Cada vez que calculo esta inversa, también quiero el determinante. Tengo una buena idea de cómo calcular, multiplicando la diagonal de la matriz U a partir de la descomposición LU. Hay un problema, puedo calcular el determinante correctamente, excepto por el signo. Dependiendo del pivote, obtengo el signo incorrecto la mitad del tiempo. ¿Alguien tiene una sugerencia sobre cómo hacer la señal correcta cada vez? Gracias de antemano.

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;

Aquí es donde inserté mi mejor oportunidad para calcular el determinante.

 T determinant = 1;

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

Finalizar mi parte del código.

 // 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 matriz de permutación pm contiene la información que necesitará para determinar el cambio de signo: querrá multiplicar su determinante por el determinante de la matriz de permutación.

Examinando el archivo fuente lu.hpp encontramos una función llamada swap_rows que dice cómo aplicar una matriz de permutación a una matriz. Se modifica fácilmente para obtener el determinante de la matriz de permutación (el signo de la permutación), dado que cada intercambio real contribuye con un factor 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;
}

Otra alternativa sería utilizar los métodos lu_factorize y lu_substitute que no hacen ningún pivote (consulte la fuente, pero básicamente elimine el pm en las llamadas a lu_factorize y lu_substitute ). Ese cambio haría que su cálculo determinante funcione tal como está. Sin embargo, tenga cuidado: eliminar el pivote hará que el algoritmo sea menos estable numéricamente.

Answer 2

Según http://qiangsong.wordpress.com / 2011/07/16 / lu-factorisation-in-ublas / :

Simplemente cambie determinant * = A (i, i) a determinant * = (pm (i) == i? 1: -1) * A (i, i) . Lo intenté de esta manera y funciona.

Sé que en realidad es muy similar a la respuesta de Managu y la idea es la misma, pero creo que es más simple (y "2 en 1" si se usa en la función InvertMatrix).