Resultados pequeños e impredecibles en ejecuciones de un modelo determinista

Tengo un modelo considerable (~ 5000 líneas) escrito en C. Es un programa en serie, sin generación de números aleatorios en ninguna parte. Utiliza la biblioteca FFTW para funciones que usan FFT: no conozco los detalles de la implementación de FFTW, pero supongo que las funciones allí también son deterministas (corríjame si estoy en error).

El problema que no puedo entender es que obtengo pequeñas diferencias en los resultados para ejecuciones idénticas en la misma máquina (mismo compilador, mismas bibliotecas).

Utilizo variables de doble precisión, y para generar el resultado en una variable, valuepor ejemplo, emito: fprintf(outFID, "%.15e\n", value);o
fwrite(&value, 1, sizeof(double), outFID);

Y constantemente obtendría diferencias como:
2.07843469652206 4 e-16 vs. 2.07843469652206 3 e-16

He pasado mucho tiempo tratando de entender por qué es esto. Inicialmente pensé que uno de mis chips de memoria se había estropeado, y los ordené y los reemplacé, sin éxito. Posteriormente también intenté ejecutar mi código en la máquina Linux de un colega, y obtengo diferencias de la misma naturaleza.

¿Qué podría estar causando esto? Es un problema pequeño ahora, pero me pregunto si es la "punta del iceberg" (de un problema grave).

Pensé que publicaría aquí en lugar de StackOverflow en caso de que alguien que trabajara con modelos numéricos pudiera haber encontrado este problema. Si alguien puede arrojar luz sobre esto, estaría muy agradecido.

Seguimiento de comentarios:
Christian Clason y Vikram: primero, gracias por su atención a mi pregunta. Los artículos que vinculó sugieren que: 1. los errores de redondeo limitan la precisión, y 2. códigos diferentes (como la introducción de declaraciones de impresión aparentemente inofensivas) pueden afectar los resultados hasta la máquina epsilon. Debo aclarar que no estoy comparando los efectos fwritey fprintffunciones. Estoy usando uno o el otro. En particular, se utiliza el mismo ejecutable para ambas ejecuciones. Simplemente estoy afirmando que el problema ocurre si uso fprintfOR fwrite.

Entonces, la ruta del código (y el ejecutable) es la misma, y ​​el hardware es el mismo. Con todos estos factores externos mantenidos constantes, ¿de dónde proviene fundamentalmente la aleatoriedad? Sospeché que el cambio de bit ocurrió debido a que la memoria defectuosa no retiene un bit correctamente, es por eso que reemplacé los chips de memoria, pero ese no parece ser el problema aquí, verifiqué e indicó. Mi programa genera miles de estos números de doble precisión en una sola ejecución, y siempre hay un puñado aleatorio que tiene volteos de bits aleatorios.

Seguimiento del primer comentario de Christian Clason: ¿Por qué lo mismo que 0 dentro de la precisión de la máquina? El número positivo más pequeño para un doble es 2.22e-308, entonces ¿no debería ser igual a 0? Mi programa genera miles de valores en el rango de 10 ^ -16 (que van desde 1e-15 a 8e-17) y hemos estado viendo variaciones significativas en nuestro proyecto de investigación, así que espero que no hayamos estado mirando sin sentido números.$2\cdot 10^{-16}$

Seguimiento n. ° 2 :
Esta es una gráfica de la serie de tiempo producida por el modelo, para ayudar en la discusión derivada en los comentarios.

Hay aspectos de los sistemas informáticos modernos que son inherentemente no deterministas que pueden causar este tipo de diferencias. Mientras las diferencias sean muy pequeñas en comparación con la precisión requerida de sus soluciones, probablemente no haya ninguna razón para preocuparse por esto.

Un ejemplo de lo que puede salir mal según mi propia experiencia. Considere el problema de calcular el producto escalar de dos vectores x e y.

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

$x_{i}y_{i}$

Por ejemplo, puede calcular el producto de dos vectores primero como

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

y luego como

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

¿Cómo puede suceder esto? Aquí hay dos posibilidades.

1. Cálculos multiproceso en núcleos paralelos. Las computadoras modernas suelen tener 2, 4, 8 o incluso más núcleos de procesador que pueden funcionar en paralelo. Si su código está usando hilos paralelos para calcular un producto de puntos en múltiples procesadores, entonces cualquier perturbación aleatoria del sistema (por ejemplo, el usuario movió su mouse y uno de los núcleos del procesador tiene que procesar el movimiento del mouse antes de regresar al producto de puntos) resultar en un cambio en el orden de las adiciones.

2. Alineación de datos e instrucciones vectoriales. Los procesadores Intel modernos tienen un conjunto especial de instrucciones que pueden funcionar (por ejemplo) para números de coma flotante a la vez. Estas instrucciones vectoriales funcionan mejor si los datos están alineados en límites de 16 bytes. Normalmente, un bucle de producto de punto dividiría los datos en secciones de 16 bytes (4 flotantes a la vez). Si vuelve a ejecutar el código por segunda vez, los datos podrían alinearse de manera diferente con los bloques de memoria de 16 bytes para que las adiciones sean realizado en un orden diferente, lo que resulta en una respuesta diferente.

Puede abordar el punto 1 haciendo que su código se ejecute como un subproceso único y deshabilitando todo el procesamiento paralelo. Puede abordar el punto 2 requiriendo la asignación de memoria para alinear bloques de memoria (normalmente lo haría compilando el código con un interruptor como -align). Si su código todavía está dando resultados que varían, entonces hay otras posibilidades para buscar a.

Esta documentación de Intel trata temas que pueden llevar a la no reproducibilidad de los resultados con la Biblioteca Intel Math Kernel. Otro documento de Intel que analiza los conmutadores de compilador para usar con los compiladores de Intel.

La biblioteca FFTW mencionada podría ejecutarse en un modo no determinista.

Si está usando el modo FFTW_MEASURE o FFTW_PATIENT, los programas verifican en tiempo de ejecución, qué valores de parámetros funcionan más rápido y luego usarán esos parámetros en todo el programa. Debido a que el tiempo de ejecución obviamente fluctuará un poco, los parámetros serán diferentes y el resultado de las transformadas de Fourier será no determinista. Si desea FFTW determinista, use el modo FFTW_ESTIMATE.

Si bien es cierto que los cambios en el orden de evaluación del término de expresión pueden ocurrir muy bien debido a escenarios de procesamiento multi-core / multi-thread, no olvide que puede haber (aunque sea una posibilidad remota) algún tipo de falla de diseño de hardware en el trabajo. ¿Recuerdas el problema Pentium FDIV? (Ver https://en.wikipedia.org/wiki/Pentium_FDIV_bug ). Hace algún tiempo, trabajé en un software de simulación de circuito analógico basado en PC. Parte de nuestra metodología implicaba desarrollar conjuntos de pruebas de regresión, que correríamos contra compilaciones nocturnas del software. Con muchos de los modelos que desarrollamos, métodos iterativos (por ejemplo, Newton-Raphson ( https://en.wikipedia.org/wiki/Newton%27s_method) y Runge-Kutta) se utilizaron ampliamente en los algoritmos de simulación. Con los dispositivos analógicos, a menudo ocurre que los artefactos internos, como voltajes, corrientes, etc., tienen valores numéricos extremadamente pequeños. Estos valores, como parte del proceso de simulación, varían gradualmente a lo largo del tiempo (simulado). La magnitud de estos cambios puede ser extremadamente pequeña, y lo que a menudo observamos fue que las operaciones posteriores de FPU en tales valores delta limitaban con el umbral de "ruido" de la precisión de la FPU (la flotación de 64 bits tiene una mantisa de 53 bits, IIRC). Eso, junto con el hecho de que a menudo tuvimos que introducir el código de registro "PrintF" en los modelos para permitir la depuración (¡ah, los buenos días!), ¡Resultados esporádicos prácticamente garantizados, a diario! Y qué' s todo esto significa? Debe esperar ver diferencias en tales circunstancias, y lo mejor que puede hacer es definir e implementar una forma de decidir (magnitud, frecuencia, tendencia, etc.) cuándo / cómo ignorarlas.

Si bien el redondeo de punto flotante de las operaciones asincrónicas puede ser el problema, sospecho que es algo más banal. El uso de una variable no inicializada que agrega aleatoriedad a su código determinista. Es un problema común que los desarrolladores suelen pasar por alto porque cuando se ejecuta en modo de depuración, todas las variables se inicializan a 0 en la declaración. Cuando no se ejecuta en modo de depuración, la memoria asignada a una variable tiene el valor que tenía la memoria antes de la asignación. La memoria no se pone a cero en la asignación como una optimización. Si esto está sucediendo en su código, será fácil de solucionar, menos en el código de las bibliotecas.