Избегайте численного потери значимости при получении определителя большой матрицы в 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;

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

Как избежать этого недостаточного потока?

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

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

Любая помощь приветствуется.


person Henrik    schedule 11.12.2013    source источник
comment
Вы можете умножить матрицу на скаляр перед вычислением определителя, скажем, 1000 * A. Затем разделить результат на 1000. Или вы можете работать с большими переменными (long double).   -  person lucas92    schedule 12.12.2013
comment
@ lucas92 lucas92 1000 было бы плохим выбором, так как умножение на 1000 обычно неточно. Вместо этого подходящая степень двойки позволит избежать снижения точности вычислений, даже если в этом нет необходимости, если не происходит переполнения.   -  person Pascal Cuoq    schedule 12.12.2013


Ответы (2)


Мне на ум приходят 2 основных варианта:

  1. Произведение собственных значений квадратной матрицы является определителем этой матрицы, поэтому сумма логарифмов каждого собственного значения есть логарифм определителя этой матрицы. Предположим, что 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 сейчас не теряет значимости, но в некоторых случаях оно все еще может также иметь значение.

  2. Выполните умное предварительное кондиционирование; здесь может быть множество вариантов, и вам лучше прочитать о них из некоторых надежных источников, так как это серьезная тема. Самый простой (и, вероятно, самый дешевый) пример предварительной обработки для этой конкретной задачи может состоять в том, чтобы вспомнить эту матрицу det(c * A) = (c ^ n) * det(A), где A — это матрица n на n, и предварительно умножить вашу матрицу на некоторое c, вычислить определитель, а затем разделить его на c ^ n, чтобы получить настоящий.

Обновлять


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

person Alexander Shukaev    schedule 11.12.2013
comment
Первое решение звучит довольно привлекательно. Вы также знаете, как я могу получить сумму значений (учитывая, что у меня есть только их журналы)? - person Henrik; 12.12.2013
comment
Вы имеете в виду, что у вас есть, например, матрицы A и B и соответствующие им логарифмы определителей, то есть log(det(A)) и log(det(B)), и теперь вы действительно хотите знать det(A) + det(B)? Если это то, что вы хотите, то я не вижу никакой проблемы в том, чтобы извлечь, например, det(A) из log(det(A)), выполнив e ^ log(det(A))... В противном случае я неправильно понял ваше намерение с суммой. - person Alexander Shukaev; 12.12.2013
comment
Да, это то, что я намерен. Но проблема по-прежнему в том, что e ^ log(det(A)) равно 0 из-за числового недополнения, поэтому я должен избегать этого. Однако я думаю, что решение с предварительным кондиционированием на самом деле еще проще. Поскольку я не ученый-компьютерщик и не математик, какой надежный источник мог бы проверить приведенную вами формулу? Любая книга по матричной алгебре? - person Henrik; 12.12.2013
comment
В том-то и проблема, что не помню книгу. Но было бы неплохо начать с базовой терминологии из Википедии, поскольку существует множество различных предобуславливателей, которые служат разным целям и в основном предназначены для решения линейных систем. Например, неполная декомпозиция холецкого определенно не для вашего случая. На самом деле, я думаю, вам подойдет самый простой метод, который я уже изложил для вас в ответе. Лучше сначала попробовать, так как копание в предварительном кондиционировании отнимет у вас довольно много времени (это действительно неприятная тема). - person Alexander Shukaev; 12.12.2013
comment
Обновил, посмотри. - person Alexander Shukaev; 12.12.2013
comment
Потрясающе, большое спасибо. Я готов принять, но я подожду, пока я не успешно реализую это (это может занять некоторое время). - person Henrik; 12.12.2013
comment
Нет проблем, рад помочь. - person Alexander Shukaev; 12.12.2013

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

person tmyklebu    schedule 12.12.2013
comment
Нет никакой гарантии, что определитель R не потеряет значимости. Это сильно зависит от основной проблемы (матрицы). - person Alexander Shukaev; 12.12.2013
comment
@Харуган: Правильно. Но становится тривиальным найти логарифм определителя R; это сумма логарифмов диагональных записей. Или вы можете использовать трюк, аналогичный тому, который вы используете для геометрических средних, чтобы получить мантиссу в double и показатель степени в достаточно большом целочисленном типе. - person tmyklebu; 12.12.2013