Boost Library, cómo obtener un 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 esto inverso, también quiero el determinante. Tengo una buena idea de cómo calcular, multiplicando la diagonal de la matriz U desde la descomposición LU. Hay un problema, puedo calcular el determinante correctamente, excepto por el signo. Dependiendo del pivote, el signo es incorrecto la mitad del tiempo. ¿Alguien tiene alguna sugerencia sobre cómo obtener el letrero correcto todo el tiempo? Gracias por adelantado.

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); 
} 

Fin de 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

matriz La permutación pm contiene la información necesaria para determinar el cambio de signo: usted querrá multiplicar su determinante por el determinante de la matriz de permutación.

Al examinar el archivo de origen lu.hpp encontramos una función llamada swap_rows que explica cómo aplicar una matriz de permutación a una matriz. Es fácilmente modificado para producir el determinante de la matriz de permutación (el signo de la permutación), dado que cada intercambio real contribuye 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 el uso de los métodos lu_factorize y lu_substitute que no haga ningún giro (consulte la fuente, pero básicamente deje caer el pm en las llamadas al lu_factorize y lu_substitute). Ese cambio haría que su cálculo determinante funcione tal como está. Tenga cuidado, sin embargo: al quitar el pivote, el algoritmo será menos estable numéricamente.

Answer 2

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

Sólo cambia determinant *= A(i,i) a determinant *= (pm(i) == i ? 1 : -1) * A(i,i). Lo intenté de esta manera y funciona.

Sé que es en realidad 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).