Algoritmo Kahan: Como obter a diferença exata de produtos

Recentemente, retornei à análise de erro de ponto flutuante para refinar alguns dos detalhes na próxima edição do livro Renderização com Base Física . Os números de ponto flutuante são uma área interessante da computação, cheia de surpresas (boas e ruins), além de truques complicados para se livrar de surpresas desagradáveis.

No processo, me deparei com este post no StackOverflow , com o qual aprendi sobre um algoritmo elegante para o cálculo exato $$$a \ times b-c \ vezes d$.

Mas antes de prosseguir com o algoritmo, você precisa entender o que é tão astuto na expressão $$$a \ times b-c \ vezes d$? Tome $$$a = 33962.035$, $$$b = -30438.8$, $$$c = 41563.4$e $$$d = -24871.969$. (Esses são os valores reais que obtive durante o lançamento do pbrt .) Para valores flutuantes de 32 bits, obtemos: $$$a \ times b = -1,03376365 \ times 10 ^ 9$e $$$c \ times d = -1,03376352 \ times 10 ^ 9$. Realizamos a subtração e obtemos $$$-128$. Porém, se você executar os cálculos com precisão dupla e, no final, convertê-los em flutuante, terá $$$-75.1656$. O que aconteceu

O problema é que o valor de cada trabalho pode ir muito além da linha de fundo $$$-1 \ vezes 10 ^ 9$, onde a distância entre valores representáveis ​​de ponto flutuante é muito grande - 64. Ou seja, ao arredondar $$$a \ vezes b$e $$$c \ vezes d$individualmente, para o ponto flutuante representável mais próximo, eles se transformam em números múltiplos de 64. Por sua vez, sua diferença será um múltiplo de 64 e não haverá esperança de que ele se torne $$$-75.1656$mais perto que $$$-64$. No nosso caso, o resultado foi ainda maior devido à forma como os dois trabalhos foram arredondados em $$$-1 \ vezes 10 ^ 9$. Encontraremos diretamente a boa e velha redução catastrófica ¹ .

Aqui está uma solução melhor que ² :

inline float DifferenceOfProducts(float a, float b, float c, float d) { float cd = c * d; float err = std::fma(-c, d, cd); float dop = std::fma(a, b, -cd); return dop + err; } 

DifferenceOfProducts() calcula $$$a \ times b-c \ vezes d$de uma maneira que evite contração catastrófica. Essa técnica foi descrita pela primeira vez pelo lendário William Kahan no artigo Sobre o custo da computação em ponto flutuante sem aritmética extra-precisa . Vale a pena notar que o trabalho de Kahan é interessante de ler como um todo, eles têm muitos comentários sobre o estado atual do mundo dos pontos flutuantes, além de considerações matemáticas e técnicas. Aqui está uma de suas descobertas:

Aqueles de nós que lutaram com as vicissitudes da aritmética de ponto flutuante e mal otimizadas "otimizações" dos compiladores podem justamente se orgulhar da vitória nesta batalha. Mas se passarmos a continuação dessa batalha para as gerações futuras, isso contradirá toda a essência da civilização. Nossa experiência diz que linguagens de programação e sistemas de desenvolvimento são fontes de muito caos com o qual temos que lidar. Muitos erros podem ser dispensados, bem como algumas otimizações atraentes que são seguras para números inteiros, mas às vezes são fatais para números de ponto flutuante.

Depois de prestar homenagem à sua inteligência, voltemos a DifferenceOfProducts() : a base do domínio dessa função é o uso de instruções FMA ^{3 com} multiplicação e adição agregada. Do ponto de vista matemático, as FMA(a,b,c) são $$$a \ vezes b + c$Portanto, a princípio, parece que essa operação é útil apenas como micro otimização: uma instrução em vez de duas. No entanto fma
Tem uma propriedade especial - termina apenas uma vez.

No habitual $$$a \ vezes b + c$primeiro calculado $$$a \ vezes b$e esse valor, que geralmente não pode ser representado no formato de ponto flutuante, é arredondado para o valor mais próximo. Então, a este valor arredondado é adicionado $$$c$, e esse resultado é novamente arredondado para o flutuador mais próximo. A FMA é implementada de forma que o arredondamento seja realizado apenas no final - um valor intermediário $$$a \ vezes b$mantém a precisão suficiente, portanto, depois de adicionar a ele $$$c$o resultado final será o mais próximo do valor real $$$a \ vezes b + c$valor de flutuação.

Tendo lidado com as FMA, retornaremos a DifferenceOfProducts() . Mais uma vez, mostrarei as duas primeiras linhas:

  float cd = c * d; float err = std::fma(-c, d, cd); 

O primeiro calcula o valor arredondado $$$c \ vezes d$e o segundo ... subtrai $$$c \ vezes d$do seu trabalho? Se você não sabe como as FMA funcionam, pode pensar que o err será sempre zero. Porém, ao trabalhar com FMA, a segunda linha realmente extrai o valor do erro de arredondamento no valor calculado $$$c \ vezes d$e salva para err . Depois disso, a saída é muito simples:

  float dop = std::fma(a, b, -cd); return dop + err; 

A segunda FMA calcula a diferença entre os trabalhos usando a FMA, realizando o arredondamento apenas no final. Consequentemente, é resistente à redução catastrófica, mas precisa trabalhar com um valor arredondado. $$$c \ vezes d$. A return "corrige" esse problema, adicionando o erro destacado na segunda linha. Em um artigo de Jeannenrod et al. é mostrado que o resultado é verdadeiro até 1,5 ulps (medidas de precisão da unidade), o que é ótimo: FMA e operações simples de ponto flutuante são precisas até 0,5 ulps, portanto o algoritmo é quase perfeito.

Nós usamos um novo martelo

Quando você começa a procurar maneiras de usar o DifferenceOfProducts() , ele é surpreendentemente útil. Calculando o discriminante da equação quadrática? DifferenceOfProducts(b, b, 4 * a, c) chamada de produtos DifferenceOfProducts(b, b, 4 * a, c) ⁴ . Calculando o determinante de uma matriz 2x2? O algoritmo resolverá esse problema. Na próxima versão do pbrt, encontrei cerca de 80 usos para ele. De todos eles, a função do produto vetorial é mais amada. Sempre foi uma fonte de problemas, por causa dos quais você teve que levantar as mãos e usar o dobro na implementação para evitar uma redução catastrófica:

 inline Vector3f Cross(const Vector3f &v1, const Vector3f &v2) { double v1x = v1.x, v1y = v1.y, v1z = v1.z; double v2x = v2.x, v2y = v2.y, v2z = v2.z; return Vector3f(v1y * v2z - v1z * v2y, v1z * v2x - v1x * v2z, v1x * v2y - v1y * v2x); } 

E agora podemos continuar trabalhando com float e usar DifferenceOfProducts() .

 inline Vector3f Cross(const Vector3f &v1, const Vector3f &v2) { return Vector3f(DifferenceOfProducts(v1.y, v2.z, v1.z, v2.y), DifferenceOfProducts(v1.z, v2.x, v1.x, v2.z), DifferenceOfProducts(v1.x, v2.y, v1.y, v2.x)); } 

Esse exemplo astuto do início do post é realmente parte de um trabalho vetorial. Em um certo estágio, o código pbrt precisa calcular o produto vetorial dos vetores $$$(33962.035, 41563.4, 7706.415)$e $$$(- 24871.969, -30438,8, -5643.727)$. Quando calculado usando float, obteríamos um vetor $$$(1552, -1248, -128)$. (Regra geral: se nos cálculos de ponto flutuante em que grandes números estão envolvidos, você obtém valores inteiros não tão grandes, esse é um sinal quase certo de que ocorreu uma redução catastrófica.)

Com precisão dupla, o produto vetorial é $$$(1556.0276, -1257.5151, -75.1656)$. Vemos isso com float $$$x$parece normal $$$y$já é muito ruim também $$$z$torna-se o desastre que se tornou a motivação para encontrar uma solução. E que resultados obteremos com os valores DifferenceOfProducts() e float? $$$(1556.0276, -1257.5153, -75.1656)$. Valores $$$x$e $$$z$correspondem a dupla precisão e $$$y$ligeiramente compensado - daí a ulp extra veio.

E a velocidade? DifferenceOfProducts() executa dois FMAs, além de multiplicação e adição. Um algoritmo ingênuo pode ser implementado com uma FMA e uma multiplicação, que, ao que parece, deve levar metade do tempo. Mas, na prática, verifica-se que, depois de obter os valores dos registros, isso não importa: no benchmark sintético do meu laptop, o DifferenceOfProducts() apenas 1,09 vezes mais caro que o algoritmo ingênuo. A operação de precisão dupla foi 2,98 vezes mais lenta.

Assim que você aprende sobre a possibilidade de uma redução catastrófica, todos os tipos de expressões de aparência inocente no código começam a parecer suspeitas. DifferenceOfProducts() parece ser uma boa cura para a maioria deles. É fácil de usar e não temos nenhuma razão específica para não usá-lo.

Anotações

1. A redução catastrófica não é um problema ao subtrair quantidades com sinais diferentes ou ao adicionar valores com o mesmo sinal. No entanto, pode se tornar um problema ao adicionar valores com sinais diferentes. Ou seja, os valores devem ser considerados com a mesma suspeita que as diferenças.
2. Como exercício para o leitor, sugiro escrever a função SumOfProducts() , que fornece proteção contra contrações catastróficas. Se você deseja complicar a tarefa, explique por que, em DifferenceOfProducts() , dop + err == dop , porque os sinais a*b c*d são diferentes.
3. A instrução FMA está disponível na GPU há mais de uma década e na maioria das CPUs por pelo menos cinco anos. Pode ser necessário adicionar sinalizadores de compilador à CPU para gerá-los diretamente ao usar std::fma() ; em gcc e clang, -march=native trabalhos -march=native .
4. No formato de ponto flutuante IEEE, a multiplicação por potências de dois é realizada exatamente, portanto 4 * a não causa nenhum erro de arredondamento, a menos que ocorra um estouro.