Я реализовал алгоритм 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 основных варианта:
Произведение собственных значений квадратной матрицы является определителем этой матрицы, поэтому сумма логарифмов каждого собственного значения является логарифмом определителя этой матрицы. Предполагать 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.
Вы можете использовать удаление Домохозяина, чтобы получить QR-разложение delta0
, Тогда определитель части Q равен +/- 1 (в зависимости от того, делали вы четное или нечетное количество отражений), а определитель части R — произведение диагональных элементов. Оба из них легко вычислить, не впадая в адский поток — и вам может даже не потребоваться первое.