Boost Library, как получить определитель из lu_factorize ()? - PullRequest
Новая
       28

Boost Library, как получить определитель из lu_factorize ()?

5 голосов
/ 14 сентября 2009

Я пытаюсь вычислить определитель, используя библиотеки boost c ++. Я нашел код для функции InvertMatrix (), который я скопировал ниже. Каждый раз, когда я вычисляю это обратное значение, мне также нужен определитель. У меня есть хорошая идея, как рассчитать, умножив диагональ матрицы U на разложение LU. Есть одна проблема, я могу правильно рассчитать определитель, за исключением знака. В зависимости от поворота я получаю неправильный знак в половине случаев. У кого-нибудь есть предложения о том, как правильно вывесить знак каждый раз? Заранее спасибо. 

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;

Здесь я вставил свой лучший способ расчета определителя. 

 T determinant = 1;

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

Конец моей части кода. 

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

Ответы [ 2 ]

3 голосов
/ 14 сентября 2009

Матрица перестановок pm содержит информацию, необходимую для определения смены знака: вы хотите умножить свой определитель на определитель матрицы перестановок.

Просматривая исходный файл lu.hpp, мы находим функцию с именем swap_rows, которая сообщает, как применить матрицу перестановок к матрице. Его легко изменить, чтобы получить определитель матрицы перестановок (знак перестановки), учитывая, что каждый фактический обмен дает множитель -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;
}

Другой альтернативой может быть использование методов lu_factorize и lu_substitute, которые не делают никакого поворота (обратитесь к источнику, но в основном отбрасываете pm в вызовах на lu_factorize и lu_substitute). Это изменение сделает ваш расчет детерминанта работать как есть. Однако будьте осторожны: устранение поворотов сделает алгоритм менее устойчивым в численном отношении.

1 голос
/ 06 апреля 2013

Согласно http://qiangsong.wordpress.com/2011/07/16/lu-factorisation-in-ublas/:

Просто измените determinant *= A(i,i) на determinant *= (pm(i) == i ? 1 : -1) * A(i,i). Я попробовал этот способ, и он работает.

Я знаю, что это на самом деле очень похоже на ответ Манагу и идея та же, но я считаю, что это проще (и "2 в 1", если используется в функции InvertMatrix).

...