Избегайте численного занижения при получении определителя большой матрицы в Eigen

Question

Избегайте численного занижения при получении определителя большой матрицы в Eigen

Я реализовал алгоритм MCMC в C ++ с использованием библиотеки Eigen. Основная часть алгоритма представляет собой цикл, в котором сначала выполняются некоторые матричные вычисления, после чего определитель полученной матрицы получается и добавляется к выходным данным. Например.:

MatrixXd delta0;
NumericVector out(3);

out[0] = 0;
out[1] = 0;

for (int i = 0; i < s; i++) {
...
delta0 = V*(A.cast<double>()-(A+B).cast<double>()*theta.asDiagonal());
...
I = delta0.determinant()
out[1] += I;
out[2] += std::sqrt(I);
}
return out;

Теперь, к сожалению, в некоторых матрицах наблюдается уменьшение числового значения, так что определитель выводится как ноль (а на самом деле это не так).

Как я могу избежать этого недостатка?

Одним из решений было бы получить вместо определителя лог определителя. Тем не мение,

Я не знаю, как это сделать;
как я мог тогда сложить эти журналы?

Любая помощь с благодарностью.

2

c++determinants eigen floating-point underflow

Решение

Другие решения

Вы можете использовать удаление Домохозяина, чтобы получить QR-разложение delta0, Тогда определитель части Q равен +/- 1 (в зависимости от того, делали вы четное или нечетное количество отражений), а определитель части R — произведение диагональных элементов. Оба из них легко вычислить, не впадая в адский поток — и вам может даже не потребоваться первое.

1

Источник

Accepted Answer

На мой взгляд, есть 2 основных варианта:

Произведение собственных значений квадратной матрицы является определителем этой матрицы, поэтому сумма логарифмов каждого собственного значения является логарифмом определителя этой матрицы. Предполагать det(A) = a а также det(B) = b для компактной записи. После применения вышеупомянутых 2 матриц A а также Bмы заканчиваем с log(a) а также log(b)тогда на самом деле верно следующее:
```
log(a + b) = log(a) + log(1 + e ^ (log(b) - log(a)))
```
Да, мы получаем логарифм суммы. Что бы вы сделали с этим дальше? Я не знаю, зависит от того, что вы должны. Если вы должны удалить логарифм по e ^ log(a + b) = a + bтогда вам может повезти, что значение a + b в настоящее время не теряется, но в некоторых случаях оно также может снижаться.
Выполните умную предварительную подготовку; здесь может быть множество вариантов, и вам лучше прочитать о них из надежных источников, поскольку это серьезная тема. Самым простым (и, возможно, самым дешевым) примером предварительной подготовки для этой конкретной проблемы может быть напоминание о том, что det(c * A) = (c ^ n) * det(A), где A является n от n матрица, и для предварительного умножения вашей матрицы с некоторыми c, вычислить определитель, а затем разделить его на c ^ n чтобы получить актуальный.

Обновить

Я подумал еще об одном варианте. Если на последних этапах # 1 или # 2 вы по-прежнему испытываете слишком частое понижение, то было бы неплохо повысить точность специально для этих последних операций, например, используя GNU MPFR.

3