Resultados pequenos e imprevisíveis em execuções de um modelo determinístico

Eu tenho um modelo considerável (~ 5000 linhas) escrito em C. É um programa serial, sem geração de números aleatórios em nenhum lugar. Utiliza a biblioteca FFTW para funções que utilizam a FFT - não conheço os detalhes da implementação da FFTW, mas presumo que as funções nela também sejam determinísticas (corrija-me se estiver errado).

O problema que não consigo entender é que estou obtendo pequenas diferenças nos resultados para execuções idênticas na mesma máquina (mesmo compilador, mesmas bibliotecas).

Uso variáveis ​​de precisão dupla e, para gerar o resultado em variável, valuepor exemplo, emito: fprintf(outFID, "%.15e\n", value);ou
fwrite(&value, 1, sizeof(double), outFID);

E eu constantemente
obtinha diferenças como: 2.07843469652206 4 e-16 vs. 2.07843469652206 3 e-16

Passei muito tempo tentando descobrir por que isso acontece. Inicialmente, pensei que um dos meus chips de memória estava com defeito e os encomendei e substituí, sem sucesso. Posteriormente, também tentei executar meu código na máquina Linux de um colega e obtive diferenças da mesma natureza.

O que poderia estar causando isso? É uma questão pequena agora, mas me pergunto se é a "ponta do iceberg" (de um problema sério).

Eu pensei em publicar aqui em vez do StackOverflow, caso alguém que trabalhe com modelos numéricos possa ter encontrado esse problema. Se alguém puder esclarecer isso, eu ficaria muito grato.

Seguimento dos comentários:
Christian Clason e Vikram: primeiro, obrigado por sua atenção à minha pergunta. Os artigos que você vinculou sugerem que: 1. erros de arredondamento limitam a precisão e 2. código diferente (como a introdução de declarações de impressão aparentemente inofensivas) pode afetar os resultados até o epsilon da máquina. Devo esclarecer que não estou comparando os efeitos fwritee fprintffunções. Eu estou usando um OU o outro. Em particular, o mesmo executável é usado para ambas as execuções. Estou simplesmente afirmando que o problema ocorre se eu uso fprintfOU fwrite.

Portanto, o caminho do código (e executável) é o mesmo e o hardware é o mesmo. Com todos esses fatores externos mantidos constantes, de onde vem a aleatoriedade, fundamentalmente? Suspeitei que o bit flip ocorreu devido a uma falha na memória que não estava retendo um pouco corretamente, e foi por isso que substituí os chips de memória, mas esse não parece ser o problema aqui, verifiquei e você indicou. Meu programa gera milhares desses números de precisão dupla em uma única execução, e sempre há um punhado aleatório com inversões aleatórias de bits.

Seguimento do primeiro comentário de Christian Clason: Por que igual a 0 na precisão da máquina? O menor número positivo para um duplo é 2,22e-308, então não deveria ser igual a 0? Meu programa gera milhares de valores na faixa de 10 ^ -16 (variando de 1e-15 a 8e-17) e temos observado variações significativas em nosso projeto de pesquisa, por isso espero que não tenhamos analisado coisas sem sentido. números.$2\cdot 10^{-16}$

Acompanhamento # 2 :
Este é um gráfico da série temporal produzida pelo modelo, para ajudar na discussão dos comentários nos comentários.

Existem aspectos dos sistemas de computação modernos que são inerentemente não determinísticos que podem causar esse tipo de diferença. Desde que as diferenças sejam muito pequenas em comparação com a precisão necessária de suas soluções, provavelmente não há motivo para se preocupar com isso.

Um exemplo do que pode dar errado com base na minha própria experiência. Considere o problema de calcular o produto escalar de dois vetores x e y.

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

$x_{i}y_{i}$

Por exemplo, você pode calcular o produto de dois vetores primeiro como

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

e então como

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

Como isso pode acontecer? Aqui estão duas possibilidades.

1. Cálculos multithread em núcleos paralelos. Os computadores modernos normalmente possuem 2, 4, 8 ou até mais núcleos de processador que podem funcionar em paralelo. Se o seu código estiver usando threads paralelos para calcular um produto pontual em vários processadores, qualquer perturbação aleatória do sistema (por exemplo, o usuário moveu o mouse e um dos núcleos do processador precisa processar esse movimento do mouse antes de retornar ao produto pontual) resultar em uma alteração na ordem das adições.

2. Alinhamento de dados e instruções vetoriais. Os processadores Intel modernos têm um conjunto especial de instruções que podem funcionar (por exemplo) para números de ponto flutuante por vez. Essas instruções de vetor funcionam melhor se os dados estiverem alinhados nos limites de 16 bytes. Normalmente, um loop de produto de ponto dividiria os dados em seções de 16 bytes (4 flutuações por vez.) Se você executar novamente o código uma segunda vez, os dados poderão ser alinhados de maneira diferente com os blocos de 16 bytes de memória, para que as adições sejam executado em uma ordem diferente, resultando em uma resposta diferente.

Você pode endereçar o ponto 1, executando seu código como um único thread e desativando todo o processamento paralelo. Você pode endereçar o ponto 2 exigindo a alocação de memória para alinhar os blocos de memória (normalmente, você faria isso compilando o código com uma opção como -align.) às.

Esta documentação da Intel discute problemas que podem levar à não reprodutibilidade dos resultados com a Intel Math Kernel Library. Outro documento da Intel que discute as opções do compilador para usar com os compiladores da Intel.

A biblioteca FFTW mencionada pode ser executada em um modo não determinístico.

Se você estiver usando o modo FFTW_MEASURE ou FFTW_PATIENT, os programas verificarão em tempo de execução, quais valores de parâmetro funcionam mais rapidamente e, em seguida, usarão esses parâmetros em todo o programa. Como o tempo de execução obviamente varia um pouco, os parâmetros serão diferentes e o resultado das transformadas de Fourier será não determinístico. Se você quiser FFTW determinístico, use o modo FFTW_ESTIMATE.

Embora seja verdade que as mudanças na ordem de avaliação do termo de expressão possam muito bem ocorrer devido a cenários de processamento com vários núcleos / vários threads, não se esqueça de que pode haver (embora seja um tiro no escuro) algum tipo de falha de design de hardware em funcionamento. Lembre-se do problema do Pentium FDIV? (Consulte https://en.wikipedia.org/wiki/Pentium_FDIV_bug ). Há algum tempo, trabalhei no software de simulação de circuitos analógicos baseado em PC. Parte de nossa metodologia envolvia o desenvolvimento de conjuntos de testes de regressão, nos quais executaríamos contra versões noturnas do software. Com muitos dos modelos que desenvolvemos, métodos iterativos (por exemplo, Newton-Raphson ( https://en.wikipedia.org/wiki/Newton%27s_method) e Runge-Kutta) foram amplamente utilizados nos algoritmos de simulação. Com dispositivos analógicos, geralmente ocorre que artefatos internos, como tensões, correntes etc., ocorrem com valores numéricos extremamente pequenos. Esses valores, como parte do processo de simulação, variam gradualmente ao longo do tempo (simulado). A magnitude dessas mudanças pode ser extremamente pequena, e o que geralmente observamos foi que as operações subseqüentes da FPU em tais valores delta limitavam o limiar de "ruído" da precisão da FPU (a flutuação de 64 bits possui uma mantissa de 53 bits, IIRC). Isso, juntamente com o fato de que muitas vezes tivemos que introduzir o código de log "PrintF" nos modelos para permitir a depuração (ah, os bons e velhos tempos!), Praticamente garantimos resultados esporádicos diariamente! E daí' tudo isso significa? Você deve esperar ver diferenças nessas circunstâncias, e a melhor coisa a fazer é definir e implementar uma maneira de decidir (magnitude, frequência, tendência, etc.) quando / como ignorá-las.

Embora o arredondamento de ponto flutuante das operações assíncronas possa ser o problema, suspeito que seja algo mais banal. O uso de variável não inicializada que está adicionando aleatoriedade ao seu código determinístico. É um problema comum que geralmente é ignorado pelos desenvolvedores porque, quando você executa no modo de depuração, todas as variáveis ​​são inicializadas como 0 na declaração. Quando não está sendo executada no modo de depuração, a memória atribuída a uma variável tem qualquer valor que a memória tivesse antes da atribuição. A memória não é zerada na atribuição como uma otimização. Se isso estiver acontecendo no seu código, será fácil corrigi-lo, menos no código das bibliotecas.