W jakiej kolejności należy dodawać liczby zmiennoprzecinkowe, aby uzyskać jak najdokładniejszy wynik?

To było pytanie, które zadałem mi na ostatnim wywiadzie i chcę wiedzieć (tak naprawdę nie pamiętam teorii analizy numerycznej, więc proszę o pomoc :)

Jeśli mamy jakąś funkcję, która gromadzi liczby zmiennoprzecinkowe:

std::accumulate(v.begin(), v.end(), 0.0);

vjest std::vector<float>na przykład.

• Czy lepiej byłoby posortować te liczby przed ich zgromadzeniem?

• Która kolejność daje najbardziej precyzyjną odpowiedź?

Podejrzewam, że sortowania liczb w kolejności rosnącej faktycznie sprawiają, że błąd liczbowy mniej , ale niestety nie mogę tego udowodnić sobie.

PS Zdaję sobie sprawę, że to prawdopodobnie nie ma nic wspólnego z programowaniem w prawdziwym świecie, po prostu jestem ciekawy.

Twój instynkt ma w zasadzie rację, sortowanie w rosnącym porządku (wielkości) zwykle nieco poprawia sytuację. Rozważmy przypadek, w którym dodajemy zmiennoprzecinkowe pojedynczej precyzji (32-bitowe) i mamy 1 miliard wartości równych 1 / (1 miliard) i jedną wartość równą 1. Jeśli 1 pojawi się jako pierwsza, to suma przyjdzie do 1, ponieważ 1 + (1/1 miliarda) to 1 z powodu utraty precyzji. Każdy dodatek nie ma żadnego wpływu na całość.

Jeśli małe wartości pojawią się pierwsze, to przynajmniej zsumują się do czegoś, chociaż nawet wtedy mam ich 2 ^ 30, a po około 2 ^ 25 wracam do sytuacji, w której każda z osobna nie wpływa na sumę już więcej. Więc nadal będę potrzebować więcej sztuczek.

To skrajny przypadek, ale generalnie dodanie dwóch wartości o podobnej wielkości jest dokładniejsze niż dodanie dwóch wartości o bardzo różnych wielkościach, ponieważ w ten sposób „odrzuca się” mniej bitów precyzji z mniejszej wartości. Sortując liczby, grupujesz razem wartości o podobnej wielkości, a dodając je w porządku rosnącym, dajesz małym wartościom „szansę” na skumulowane osiągnięcie wielkości większych liczb.

Jednak jeśli chodzi o liczby ujemne, łatwo jest „przechytrzyć” to podejście. Rozważmy trzy wartości w sumie {1, -1, 1 billionth}. Suma poprawna arytmetycznie to 1 billionth, ale jeśli moje pierwsze dodanie obejmuje małą wartość, wtedy moja suma końcowa wyniesie 0. Z 6 możliwych zamówień tylko 2 są „poprawne” - {1, -1, 1 billionth}i {-1, 1, 1 billionth}. Wszystkie 6 rzędów dają wyniki, które są dokładne w skali największej wartości na wejściu (0,0000001% na zewnątrz), ale dla 4 z nich wynik jest niedokładny w skali rzeczywistego rozwiązania (100% poza). Konkretny problem, który rozwiązujesz, powie ci, czy ten pierwszy jest wystarczająco dobry, czy nie.

W rzeczywistości możesz grać o wiele więcej sztuczek, niż tylko dodawać je w posortowanej kolejności. Jeśli masz wiele bardzo małych wartości, średnią liczbę średnich wartości i niewielką liczbę dużych wartości, najdokładniejsze może być najpierw zsumowanie wszystkich małych, a następnie osobne zsumowanie średnich i dodanie tych dwóch sum razem, a następnie dodaj duże. Znalezienie najdokładniejszej kombinacji dodawań zmiennoprzecinkowych nie jest wcale trywialne, ale aby poradzić sobie z naprawdę złymi przypadkami, możesz zachować cały szereg bieżących sum o różnych wielkościach, dodać każdą nową wartość do sumy, która najlepiej pasuje do jej wielkości, a kiedy suma bieżąca zacznie być zbyt duża dla swojej wielkości, dodaj ją do następnej sumy i rozpocznij nową. Doprowadzony do jego logicznego ekstremum, proces ten jest równoważny wykonaniu sumy w typie z dowolną precyzją (więc zrób to). Biorąc jednak pod uwagę uproszczony wybór dodawania rosnącego lub malejącego rzędu wielkości, zwiększanie jest lepszym rozwiązaniem.

Ma to pewien związek z programowaniem w świecie rzeczywistym, ponieważ istnieją przypadki, w których obliczenia mogą pójść bardzo źle, jeśli przypadkowo odetniesz „ciężki” ogon składający się z dużej liczby wartości, z których każda jest zbyt mała, aby mieć na nią indywidualny wpływ suma lub jeśli odrzucisz zbyt dużą precyzję z wielu małych wartości, które indywidualnie wpływają tylko na kilka ostatnich bitów sumy. W przypadkach, gdy ogon i tak jest znikomy, prawdopodobnie nie obchodzi cię to. Na przykład, jeśli na początku dodajesz tylko niewielką liczbę wartości i używasz tylko kilku znaczących cyfr z sumy.

Istnieje również algorytm przeznaczony do tego rodzaju operacji akumulacji, zwany sumowaniem Kahana , o którym prawdopodobnie powinieneś wiedzieć.

Według Wikipedii

Algorytm sumowanie Kahan (znany także jako wyrównaną sumowania ) znacząco zmniejsza błąd liczbowe w całości otrzymanego przez dodanie sekwencji liczb zmiennoprzecinkowych skończonej precyzji w stosunku do oczywistym podejściem. Odbywa się to poprzez utrzymywanie oddzielnej bieżącej kompensacji (zmienna do kumulowania małych błędów).

W pseudokodzie algorytm to:

function kahanSum(input)
 var sum = input[1]
 var c = 0.0          //A running compensation for lost low-order bits.
 for i = 2 to input.length
  y = input[i] - c    //So far, so good: c is zero.
  t = sum + y         //Alas, sum is big, y small, so low-order digits of y are lost.
  c = (t - sum) - y   //(t - sum) recovers the high-order part of y; subtracting y recovers -(low part of y)
  sum = t             //Algebraically, c should always be zero. Beware eagerly optimising compilers!
 next i               //Next time around, the lost low part will be added to y in a fresh attempt.
return sum

Wypróbowałem skrajny przykład w odpowiedzi udzielonej przez Steve'a Jessopa.

#include <iostream>
#include <iomanip>
#include <cmath>

int main()
{
    long billion = 1000000000;
    double big = 1.0;
    double small = 1e-9;
    double expected = 2.0;

    double sum = big;
    for (long i = 0; i < billion; ++i)
        sum += small;
    std::cout << std::scientific << std::setprecision(1) << big << " + " << billion << " * " << small << " = " <<
        std::fixed << std::setprecision(15) << sum <<
        "    (difference = " << std::fabs(expected - sum) << ")" << std::endl;

    sum = 0;
    for (long i = 0; i < billion; ++i)
        sum += small;
    sum += big;
    std::cout  << std::scientific << std::setprecision(1) << billion << " * " << small << " + " << big << " = " <<
        std::fixed << std::setprecision(15) << sum <<
        "    (difference = " << std::fabs(expected - sum) << ")" << std::endl;

    return 0;
}

Otrzymałem następujący wynik:

1.0e+00 + 1000000000 * 1.0e-09 = 2.000000082740371    (difference = 0.000000082740371)
1000000000 * 1.0e-09 + 1.0e+00 = 1.999999992539933    (difference = 0.000000007460067)

Błąd w pierwszej linii jest ponad dziesięciokrotnie większy w drugiej.

Jeśli zmienię doubles na floatsw powyższym kodzie, otrzymam:

1.0e+00 + 1000000000 * 1.0e-09 = 1.000000000000000    (difference = 1.000000000000000)
1000000000 * 1.0e-09 + 1.0e+00 = 1.031250000000000    (difference = 0.968750000000000)

Żadna z odpowiedzi nie jest nawet bliska 2,0 (ale druga jest nieco bliżej).

Używając sumowania Kahana (ze doubles), jak opisał Daniel Pryden:

#include <iostream>
#include <iomanip>
#include <cmath>

int main()
{
    long billion = 1000000000;
    double big = 1.0;
    double small = 1e-9;
    double expected = 2.0;

    double sum = big;
    double c = 0.0;
    for (long i = 0; i < billion; ++i) {
        double y = small - c;
        double t = sum + y;
        c = (t - sum) - y;
        sum = t;
    }

    std::cout << "Kahan sum  = " << std::fixed << std::setprecision(15) << sum <<
        "    (difference = " << std::fabs(expected - sum) << ")" << std::endl;

    return 0;
}

Dostaję dokładnie 2,0:

Kahan sum  = 2.000000000000000    (difference = 0.000000000000000)

I nawet jeśli zmienię doubles na floatsw powyższym kodzie, otrzymam:

Kahan sum  = 2.000000000000000    (difference = 0.000000000000000)

Wydawałoby się, że Kahan to najlepsza droga!

Istnieje pewna klasa algorytmów, które rozwiązują dokładnie ten problem, bez konieczności sortowania lub zmiany kolejności danych .

Innymi słowy, sumowanie można wykonać jednym przejściem po danych. To sprawia, że ​​takie algorytmy mają zastosowanie w sytuacjach, gdy zbiór danych nie jest z góry znany, np. Gdy dane docierają w czasie rzeczywistym i trzeba zachować sumę bieżącą.

Oto streszczenie ostatniego artykułu:

Przedstawiamy nowatorski algorytm online do dokładnego sumowania strumienia liczb zmiennoprzecinkowych. Przez „online” rozumiemy, że algorytm musi widzieć tylko jedno wejście naraz i może przyjąć strumień wejściowy o dowolnej długości z takich danych wejściowych, wymagając tylko stałej pamięci. Przez „dokładne” rozumiemy, że suma wewnętrznej tablicy naszego algorytmu jest dokładnie równa sumie wszystkich danych wejściowych, a zwracany wynik jest sumą prawidłowo zaokrągloną. Dowód poprawności jest ważny dla wszystkich danych wejściowych (w tym liczb nieznormalizowanych, ale przepełnienia pośredniego modulo) i jest niezależny od liczby sumy lub numeru warunku sumy. Algorytm asymptotycznie potrzebuje tylko 5 FLOPów na szczyt, a ze względu na paralelizm na poziomie instrukcji działa tylko około 2-3 razy wolniej niż oczywiste, szybka, ale głupia „zwykła rekurencyjna pętla sumowania”, gdy liczba sumowań jest większa niż 10 000. Zatem, zgodnie z naszą wiedzą, jest to najszybszy, najdokładniejszy i najbardziej wydajny w pamięci spośród znanych algorytmów. Rzeczywiście, trudno jest zrozumieć, w jaki sposób szybszy algorytm lub taki, który wymaga znacznie mniejszej liczby FLOP-ów, mógłby istnieć bez ulepszeń sprzętowych. Dostarczany jest wniosek o wiele zjazdów.

Źródło: Algorytm 908: Dokładne sumowanie online strumieni zmiennoprzecinkowych .

Opierając się na odpowiedzi Steve'a, aby najpierw posortować liczby w porządku rosnącym, przedstawiłbym jeszcze dwa pomysły:

1. Zdecyduj się na różnicę w wykładniku dwóch liczb, powyżej której możesz zdecydować, że stracisz zbyt dużą precyzję.

2. Następnie dodaj liczby w kolejności, aż wykładnik akumulatora będzie zbyt duży dla następnej liczby, a następnie umieść akumulator w tymczasowej kolejce i rozpocznij akumulator z następną liczbą. Kontynuuj, aż wyczerpiesz oryginalną listę.

Powtarzasz ten proces z tymczasową kolejką (po posortowaniu) i prawdopodobnie większą różnicą wykładnika.

Myślę, że będzie to dość powolne, jeśli będziesz musiał przez cały czas obliczać wykładniki.

Szybko przeszedłem z programem i wynik był 1,99903

Myślę, że możesz zrobić coś lepszego niż sortowanie liczb, zanim je zbierzesz, ponieważ podczas procesu akumulacji akumulator staje się coraz większy. Jeśli masz dużą liczbę podobnych liczb, szybko zaczniesz tracić precyzję. Oto, co proponuję zamiast tego:

while the list has multiple elements
    remove the two smallest elements from the list
    add them and put the result back in
the single element in the list is the result

Oczywiście ten algorytm będzie najbardziej efektywny z kolejką priorytetową zamiast listy. Kod C ++:

template <typename Queue>
void reduce(Queue& queue)
{
    typedef typename Queue::value_type vt;
    while (queue.size() > 1)
    {
        vt x = queue.top();
        queue.pop();
        vt y = queue.top();
        queue.pop();
        queue.push(x + y);
    }
}

kierowca:

#include <iterator>
#include <queue>

template <typename Iterator>
typename std::iterator_traits<Iterator>::value_type
reduce(Iterator begin, Iterator end)
{
    typedef typename std::iterator_traits<Iterator>::value_type vt;
    std::priority_queue<vt> positive_queue;
    positive_queue.push(0);
    std::priority_queue<vt> negative_queue;
    negative_queue.push(0);
    for (; begin != end; ++begin)
    {
        vt x = *begin;
        if (x < 0)
        {
            negative_queue.push(x);
        }
        else
        {
            positive_queue.push(-x);
        }
    }
    reduce(positive_queue);
    reduce(negative_queue);
    return negative_queue.top() - positive_queue.top();
}

Liczby w kolejce są ujemne, ponieważ topdaje największą liczbę, ale chcemy najmniejszą . Mogłem podać więcej argumentów szablonu do kolejki, ale takie podejście wydaje się prostsze.

To nie do końca odpowiada na twoje pytanie, ale mądrą rzeczą jest dwukrotne obliczenie sumy, raz w trybie zaokrąglania „zaokrąglenie w górę” i raz za pomocą „zaokrąglenia w dół”. Porównaj te dwie odpowiedzi, a wiesz, / jak / niedokładne są wyniki, i jeśli w związku z tym musisz użyć sprytniejszej strategii sumowania. Niestety, większość języków nie sprawia, że ​​zmiana trybu zaokrąglania zmiennoprzecinkowego jest tak łatwa, jak powinna, ponieważ ludzie nie wiedzą, że jest on faktycznie przydatny w codziennych obliczeniach.

Spójrz na arytmetykę interwałową, w której wykonujesz wszystkie obliczenia matematyczne w ten sposób, zachowując najwyższe i najniższe wartości. Prowadzi to do interesujących wyników i optymalizacji.

Najprostszym sortowaniem poprawiającym dokładność jest sortowanie według rosnącej wartości bezwzględnej. Dzięki temu najmniejsze wartości wielkości mają szansę na akumulację lub anulowanie przed interakcją z większymi wartościami wielkości, które spowodowałyby utratę precyzji.

To powiedziawszy, możesz zrobić to lepiej, śledząc wiele nienakładających się sum częściowych. Oto artykuł opisujący technikę i przedstawiający dowód dokładności: www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps

Ten algorytm i inne podejścia do dokładnego sumowania zmiennoprzecinkowego są zaimplementowane w prostym Pythonie pod adresem : http://code.activestate.com/recipes/393090/ Przynajmniej dwa z nich można w trywialny sposób przekonwertować na C ++.

W przypadku numerów IEEE 754 o pojedynczej lub podwójnej precyzji lub w znanym formacie, inną alternatywą jest użycie tablicy liczb (przekazywanych przez wywołującego lub w klasie dla C ++) indeksowanych przez wykładnik. Podczas dodawania liczb do tablicy dodawane są tylko liczby z tym samym wykładnikiem (do momentu znalezienia pustego pola i zapisania liczby). Gdy żądana jest suma, tablica jest sumowana od najmniejszej do największej, aby zminimalizować obcięcie. Przykład pojedynczej precyzji:

/* clear array */
void clearsum(float asum[256])
{
size_t i;
    for(i = 0; i < 256; i++)
        asum[i] = 0.f;
}

/* add a number into array */
void addtosum(float f, float asum[256])
{
size_t i;
    while(1){
        /* i = exponent of f */
        i = ((size_t)((*(unsigned int *)&f)>>23))&0xff;
        if(i == 0xff){          /* max exponent, could be overflow */
            asum[i] += f;
            return;
        }
        if(asum[i] == 0.f){     /* if empty slot store f */
            asum[i] = f;
            return;
        }
        f += asum[i];           /* else add slot to f, clear slot */
        asum[i] = 0.f;          /* and continue until empty slot */
    }
}

/* return sum from array */
float returnsum(float asum[256])
{
float sum = 0.f;
size_t i;
    for(i = 0; i < 256; i++)
        sum += asum[i];
    return sum;
}

przykład podwójnej precyzji:

/* clear array */
void clearsum(double asum[2048])
{
size_t i;
    for(i = 0; i < 2048; i++)
        asum[i] = 0.;
}

/* add a number into array */
void addtosum(double d, double asum[2048])
{
size_t i;
    while(1){
        /* i = exponent of d */
        i = ((size_t)((*(unsigned long long *)&d)>>52))&0x7ff;
        if(i == 0x7ff){         /* max exponent, could be overflow */
            asum[i] += d;
            return;
        }
        if(asum[i] == 0.){      /* if empty slot store d */
            asum[i] = d;
            return;
        }
        d += asum[i];           /* else add slot to d, clear slot */
        asum[i] = 0.;           /* and continue until empty slot */
    }
}

/* return sum from array */
double returnsum(double asum[2048])
{
double sum = 0.;
size_t i;
    for(i = 0; i < 2048; i++)
        sum += asum[i];
    return sum;
}

Twoje pływaki powinny być dodawane z podwójną precyzją. Zapewni to większą precyzję niż jakakolwiek inna technika. Aby uzyskać nieco większą precyzję i znacznie większą prędkość, możesz utworzyć powiedzmy cztery sumy i zsumować je na końcu.

Jeśli dodajesz liczby podwójnej precyzji, użyj long double jako sumy - jednak będzie to miało pozytywny wpływ tylko w implementacjach, w których long double faktycznie ma większą precyzję niż double (zazwyczaj x86, PowerPC w zależności od ustawień kompilatora).

Jeśli chodzi o sortowanie, wydaje mi się, że jeśli spodziewasz się anulowania, liczby powinny być dodawane w porządku malejącym , a nie rosnącym. Na przykład:

((-1 + 1) + 1e-20) da 1e-20

ale

((1e-20 + 1) - 1) da 0

W pierwszym równaniu dwie duże liczby są anulowane, podczas gdy w drugim człon 1e-20 zostaje utracony po dodaniu do 1, ponieważ nie ma wystarczającej precyzji, aby go zachować.

Ponadto sumowanie parami jest całkiem przyzwoite do sumowania wielu liczb.