Как бороться с недостаточным потоком в научных вычислениях?

Я работаю над вероятностными моделями, и, делая вывод о тех моделях, оцененные вероятности могут стать очень маленькими. Чтобы избежать переполнения, я в настоящее время работаю в домене журнала (я сохраняю журнал вероятностей). Умножимые вероятности эквивалентны добавлению, а суммирование производится по формуле:

log(exp(a) + exp(b)) = log(exp(a - m) + exp(b - m)) + m

где m = max(a, b).

Я использую некоторые очень большие матрицы, и для вычисления матрично-векторных умножений я должен принять элементарную экспоненту этих матриц. Этот шаг довольно дорогостоящий, и мне было интересно, существуют ли другие методы борьбы с недостаточным потоком при работе с вероятностями.

Изменить: по соображениям эффективности, я ищу решение с использованием примитивных типов, а не объектов, сохраняющих произвольное представление реальных чисел.

Изменить 2: Я ищу более быстрое решение, чем трюк в лог-домене, а не более точное решение. Я доволен той точностью, которую я сейчас получаю, но мне нужен более быстрый метод. В частности, суммирование происходит при умножении матричных векторов, и я хотел бы иметь возможность использовать эффективные методы BLAS.

Решение: после обсуждения с Джонатаном Дурси я решил разложить каждую матрицу и вектор на свой самый большой элемент и сохранить этот фактор в домене журнала. Умножения просты. Перед добавлением я должен разложить одну из добавленных матриц/векторов на отношение двух факторов. Я обновляю коэффициент каждые десять операций.

Ответ 1

Эта проблема недавно появилась на сайте обмена стеками вычислительных наук, и, несмотря на то, что возникли проблемы с переполнением, проблемы более или менее то же самое.

Преобразование в пространство журналов, безусловно, является одним разумным подходом. Независимо от того, в каком пространстве вы находитесь, для правильного выполнения большого количества сумм, используйте несколько методов, которые вы можете использовать для повышения точности ваших суммирования. Компенсированные подходы суммирования, наиболее известные суммирование Kahan, содержат как сумму, так и то, что эффективно "остаток"; это дает вам некоторые преимущества использования арифметики с более высокой точностью без всякой стоимости (и только с использованием примитивных типов). Остальной термин также дает вам некоторое представление о том, насколько хорошо вы это делаете.

В дополнение к улучшению фактической механики вашего добавления, изменение порядка того, как вы добавляете свои термины, может иметь большое значение. Сортировка ваших терминов, чтобы вы суммировали от самых маленьких до крупнейших, может помочь, так как тогда вы больше не добавляете термины так часто, что очень разные (что может вызвать значительные проблемы округления); в некоторых случаях выполнение log ₂ N повторных попарных сумм также может быть улучшением по сравнению с простолинейной линейной суммой, в зависимости от того, как выглядят ваши условия.

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

Ответ 2

Вариант 1: Commons Math - Библиотека математики Apache Commons

Commons Math - это библиотека с легкими, автономными компонентами математики и статистики, которые касаются наиболее распространенных проблем, а не доступный на языке программирования Java или Commons Lang.

Примечание. API защищает конструкторы для принудительного создания шаблона factory при присвоении имени factory DfpField (а не несколько более интуитивно понятного DfpFac или DfpFactory). Поэтому вы должны использовать

new DfpField(numberOfDigits).newDfp(myNormalNumber)

чтобы создать экземпляр Dfp, вы можете вызвать .multiply или что-то еще. Я думал, что упомянул об этом, потому что это немного запутанно.

Вариант 2: Научная библиотека GNU или Библиотеки Boost С++. В этих случаях вы должны использовать JNI для вызова этих родных библиотек.

Вариант 3:. Если вы можете использовать другие программы и/или языки, вы можете рассмотреть возможность использования программ/языков для численных вычислений, таких как Octave, Scilab и т.д.

Вариант 4: BigDecimal Java.

Ответ 3

Я столкнулся с подобной проблемой много лет назад. Решение заключалось в разработке приближения log (1 + exp (-x)). Диапазон аппроксимации не должен быть таким большим (x от 0 до 40 будет более чем достаточным), и, по крайней мере, в моем случае точность не должна быть особенно высокой.

В вашем случае, похоже, вам нужно вычислить log (1 + exp (-x1) + exp (-x2) +...). Выбросьте эти большие отрицательные значения. Например, предположим, что a, b и c - три вероятности журнала, с 0 > a > b > c. Вы можете игнорировать c, если a-c > 38. Это не будет способствовать вашей общей вероятности регистрации вообще, по крайней мере, если вы работаете с парным разрядом.

Ответ 4

Вместо того, чтобы хранить значения в логарифмической форме, я думаю, вам, вероятно, будет лучше использовать ту же концепцию, что и double s, а именно, представление с плавающей запятой. Например, вы можете сохранить каждое значение как два long s, один для sign-and-mantissa и один для экспоненты. (Реальная плавающая точка имеет тщательно настроенный дизайн, поддерживающий множество крайних случаев и избегающий потери одного бита, но вам, вероятно, не нужно так беспокоиться ни о каком из них, и может сосредоточиться на его разработке таким образом, чтобы простой для реализации.)

Ответ 5

Я не понимаю, почему это работает, но эта формула работает и проще:

c = a + log(1 + exp(b - a))

Где c = log(exp(a)+exp(b))