Небольшие, непредсказуемые результаты в прогонах детерминированной модели

У меня есть значительная модель (~ 5000 строк), написанная на C. Это последовательная программа, нигде не генерирующая случайные числа. Она использует библиотеку FFTW для функций, использующих FFT - я не знаю деталей реализации FFTW, но я предполагаю, что функции в ней также являются детерминированными (поправьте меня, если я ошибаюсь).

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

Я использую переменные двойной точности, и для вывода результата в переменную, valueнапример, я выдаю: fprintf(outFID, "%.15e\n", value);или
fwrite(&value, 1, sizeof(double), outFID);

И я бы постоянно получал различия, такие как:
2.07843469652206 4 e-16 против 2.07843469652206 3 e-16

Я потратил много времени, пытаясь понять, почему это так. Сначала я думал, что один из моих чипов памяти испортился, и я заказал и заменил их, но безрезультатно. Впоследствии я также попытался запустить свой код на компьютере коллеги с Linux, и у меня возникли различия той же природы.

Что может быть причиной этого? Сейчас это небольшая проблема, но мне интересно, является ли это «верхушкой айсберга» (серьезной проблемы).

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

Продолжение комментариев:
Кристиан Класон и Викрам: во-первых, спасибо за внимание к моему вопросу. Статьи, на которые вы ссылаетесь, предполагают, что: 1. ошибки округления ограничивают точность, и 2. другой код (например, введение, казалось бы, безвредных операторов печати) может повлиять на результаты вплоть до epsilon машины. Я должен уточнить, что я не сравниваю эффекты fwriteи fprintfфункции. Я использую один или другой. В частности, один и тот же исполняемый файл используется для обоих запусков. Я просто заявляю, что проблема возникает, использую ли я fprintfИЛИ fwrite.

Таким образом, путь к коду (и исполняемый файл) одинаков, а оборудование одинаково. Когда все эти внешние факторы остаются постоянными, откуда берется случайность? Я подозревал, что переворот произошел из-за того, что неисправная память не сохранила немного корректно, поэтому я заменил микросхемы памяти, но здесь это не проблема, я это проверил, и вы указали. Моя программа выводит тысячи этих чисел двойной точности за один прогон, и всегда есть случайная кучка, которая имеет случайные перевороты битов.

$2\cdot 10^{-16}$

Продолжение № 2 :
Это график временных рядов, выводимых моделью, чтобы помочь в обсуждении ответвления в комментариях.

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

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

$d=\sum_{i=1}^{n} x_{i}y_{i}$

$x_{i}y_{i}$

Например, вы можете сначала вычислить произведение двух векторов как

$d=((x_{1}y_{1})+(x_{2}y_{2}))+(x_{3}y_{3})$

а потом как

$d=(x_{1}y_{1})+((x_{2}y_{2})+(x_{3}y_{3}))$

Как это могло случиться? Вот две возможности.

1. Многопоточные вычисления на параллельных ядрах. Современные компьютеры обычно имеют 2, 4, 8 или даже больше процессорных ядер, которые могут работать параллельно. Если ваш код использует параллельные потоки для вычисления точечного продукта на нескольких процессорах, то любое случайное возмущение системы (например, пользователь переместил свою мышь и одно из ядер процессора должно обработать это движение мыши, прежде чем вернуться к точечному продукту), может привести к изменению порядка дополнений.

2. Выравнивание данных и векторных инструкций. Современные процессоры Intel имеют специальный набор инструкций, которые могут работать (например) для чисел с плавающей запятой за раз. Эти векторные инструкции работают лучше всего, если данные выровнены по 16-байтовым границам. Как правило, точечный цикл произведения разбивает данные на 16-байтовые секции (по 4 с плавающей точкой за раз). Если вы повторно запустите код, данные могут быть выровнены по-разному с 16-байтовыми блоками памяти, так что сложения будут выполняется в другом порядке, что приводит к другому ответу.

Вы можете обратиться к пункту 1, запустив ваш код как один поток и отключив всю параллельную обработку. Вы можете обратиться к пункту 2, требуя выделения памяти для выравнивания блоков памяти (как правило, это можно сделать, скомпилировав код с помощью параметра, например -align.) Если ваш код все еще дает результаты, которые отличаются, то есть другие возможности для просмотра в.

В этой документации Intel обсуждаются проблемы, которые могут привести к невоспроизводимости результатов с библиотекой Intel Math Kernel. Другой документ от Intel, в котором обсуждаются переключатели компиляторов для использования с компиляторами Intel.

Упомянутая библиотека FFTW может работать в недетерминированном режиме.

Если вы используете режим FFTW_MEASURE или FFTW_PATIENT, программы проверяют во время выполнения, какие значения параметров работают быстрее всего, а затем будут использовать эти параметры во всей программе. Поскольку время выполнения, очевидно, будет немного колебаться, параметры будут другими, и результат преобразований Фурье будет недетерминированным. Если вам нужен детерминированный FFTW, используйте режим FFTW_ESTIMATE.

Хотя верно и то, что изменения порядка вычисления терминов выражения могут очень хорошо происходить из-за сценариев многоядерной / многопоточной обработки, не забывайте, что может быть (даже если это слишком далеко) некоторый недостаток проектирования аппаратного обеспечения в работе. Помните проблему с Pentium FDIV? (См. Https://en.wikipedia.org/wiki/Pentium_FDIV_bug ). Некоторое время назад я работал над программным обеспечением для моделирования аналоговых цепей на ПК. Часть нашей методологии включала разработку наборов регрессионных тестов, которые мы будем использовать против ночных сборок программного обеспечения. Со многими из моделей, которые мы разработали, итерационные методы (например, Ньютон-Рафсон ( https://en.wikipedia.org/wiki/Newton%27s_method)) и Рунге-Кутта) широко использовались в алгоритмах моделирования. В аналоговых устройствах часто случается, что внутренние артефакты, такие как напряжения, токи и т. Д., Имеют очень малые числовые значения. Эти значения, как часть процесса моделирования, постепенно изменяются в течение (смоделированного) времени. Величина этих изменений может быть чрезвычайно мала, и мы часто наблюдали, что последующие операции FPU с такими значениями дельты граничат с порогом «шума» точности FPU (64-битное плавающее имеет 53-битную мантиссу, IIRC). Это, в сочетании с тем фактом, что нам часто приходилось вводить код журналирования «PrintF» в модели для обеспечения возможности отладки (ах, добрых старых дней!), Практически гарантировало спорадические результаты ежедневно! И что' Это все значит? Вы должны ожидать увидеть различия при таких обстоятельствах, и лучшее, что нужно сделать, - это определить и реализовать способ решить (величину, частоту, тренд и т. Д.), Когда / как их игнорировать.

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