Effiziente Methode zur Berechnung des geometrischen Mittels vieler Zahlen

Ich muss das geometrische Mittel einer großen Menge von Zahlen berechnen, deren Werte a priori nicht begrenzt sind. Der naive Weg wäre

double geometric_mean(std::vector<double> const&data) // failure
{
  auto product = 1.0;
  for(auto x:data) product *= x;
  return std::pow(product,1.0/data.size());
}

Dies kann jedoch aufgrund von Unterlauf oder Überlauf im Aufstau durchaus scheitern product (Hinweis: long double vermeidet dieses Problem nicht wirklich). Die nächste Option besteht also darin, die Logarithmen zusammenzufassen:

double geometric_mean(std::vector<double> const&data)
{
  auto sumlog = 0.0;
  for(auto x:data) sum_log += std::log(x);
  return std::exp(sum_log/data.size());
}

Das funktioniert, ruft aber an std::log() für jedes Element, das potenziell langsam ist. Kann ich das vermeiden? Zum Beispiel durch Nachverfolgen (des Äquivalents von) des Exponenten und der Mantisse des Akkumulierten product separat?

Die Lösung “Split Exponent und Mantisse”:

double geometric_mean(std::vector<double> const & data)
{
    double m = 1.0;
    long long ex = 0;
    double invN = 1.0 / data.size();

    for (double x : data)
    {
        int i;
        double f1 = std::frexp(x,&i);
        m*=f1;
        ex+=i;
    }

    return std::pow( std::numeric_limits<double>::radix,ex * invN) * std::pow(m,invN);
}

Wenn Sie das befürchten ex könnte überlaufen, können Sie es als Double anstelle von a definieren long longund multipliziere mit invN bei jedem Schritt, aber Sie könnten mit diesem Ansatz viel Präzision verlieren.

BEARBEITEN Für große Eingaben können wir die Berechnung in mehrere Buckets aufteilen:

double geometric_mean(std::vector<double> const & data)
{
    long long ex = 0;
    auto do_bucket = [&data,&ex](int first,int last) -> double
    {
        double ans = 1.0;
        for ( ;first != last;++first)
        {
            int i;
            ans *= std::frexp(data[first],&i);
            ex+=i;
        }
        return ans;
    };

    const int bucket_size = -std::log2( std::numeric_limits<double>::min() );
    std::size_t buckets = data.size() / bucket_size;

    double invN = 1.0 / data.size();
    double m = 1.0;

    for (std::size_t i = 0;i < buckets;++i)
        m *= std::pow( do_bucket(i * bucket_size,(i+1) * bucket_size),invN );

    m*= std::pow( do_bucket( buckets * bucket_size, data.size() ),invN );

    return std::pow( std::numeric_limits<double>::radix,ex * invN ) * m;
}

Ich glaube, ich habe einen Weg gefunden, es zu tun, es kombiniert die beiden Routinen in der Frage, ähnlich wie Peters Idee. Hier ist ein Beispielcode.

double geometric_mean(std::vector<double> const&data)
{
    const double too_large = 1.e64;
    const double too_small = 1.e-64;
    double sum_log = 0.0;
    double product = 1.0;
    for(auto x:data) {
        product *= x;
        if(product > too_large || product < too_small) {
            sum_log+= std::log(product);
            product = 1;      
        }
    }
    return std::exp((sum_log + std::log(product))/data.size());
}

Die schlechte Nachricht ist: Das kommt mit einem Ast. Die gute Nachricht: Der Verzweigungsprädiktor wird dies wahrscheinlich fast immer richtig machen (die Verzweigung sollte nur selten ausgelöst werden).

Die Verzweigung könnte mit Peters Idee einer konstanten Anzahl von Termen im Produkt vermieden werden. Das Problem dabei ist, dass abhängig von den Werten immer noch innerhalb weniger Terme ein Überlauf/Unterlauf auftreten kann.

Sie können dies möglicherweise beschleunigen, indem Sie Zahlen wie in Ihrer ursprünglichen Lösung multiplizieren und nur jede bestimmte Anzahl von Multiplikationen in Logarithmen umwandeln (abhängig von der Größe Ihrer Anfangszahlen).

Ein anderer Ansatz, der eine bessere Genauigkeit und Leistung als das Logarithmusverfahren erbringen würde, bestünde darin, außerhalb des Bereichs liegende Exponenten durch einen festen Betrag zu kompensieren, wobei ein exakter Logarithmus des gelöschten Überschusses aufrechterhalten würde. So:

const int EXP = 64; // maximal/minimal exponent
const double BIG = pow(2, EXP); // overflow threshold
const double SMALL = pow(2, -EXP); // underflow threshold

double product = 1;
int excess = 0; // number of times BIG has been divided out of product

for(int i=0; i<n; i++)
{
    product *= A[i];
    while(product > BIG)
    {
        product *= SMALL;
        excess++;
    }
    while(product < SMALL)
    {
        product *= BIG;
        excess--;
    }
}

double mean = pow(product, 1.0/n) * pow(BIG, double(excess)/n);

Alle Multiplikationen mit BIG und SMALL sind genau, und es gibt keine Anrufe zu log (eine transzendente und daher besonders ungenaue Funktion).

Es gibt eine einfache Idee, um die Berechnung zu reduzieren und auch einen Überlauf zu verhindern. Sie können Zahlen gruppieren, sagen wir mindestens zwei gleichzeitig, und ihr Protokoll berechnen und dann ihre Summe auswerten.

log(abcde) = 5*log(K)

log(ab) + log(cde)  = 5*log(k)

Das Summieren von Protokollen zur stabilen Berechnung von Produkten ist vollkommen in Ordnung und ziemlich effizient (wenn dies nicht ausreicht: Es gibt Möglichkeiten, vektorisierte Logarithmen mit ein paar SSE-Operationen zu erhalten – es gibt auch die Vektoroperationen von Intel MKL).

Um einen Überlauf zu vermeiden, besteht eine übliche Technik darin, jede Zahl vorher durch den maximalen oder minimalen Betragseintrag zu dividieren (oder logarithmische Differenzen zu logmax oder logmin zu summieren). Sie können auch Buckets verwenden, wenn die Zahlen stark variieren (z. B. das Protokoll von kleinen Zahlen und großen Zahlen separat summieren). Beachten Sie, dass normalerweise nichts davon benötigt wird, außer für sehr große Mengen, da das Protokoll von a double ist nie riesig (zwischen sagen wir -700 und 700).

Außerdem müssen Sie die Schilder separat im Auge behalten.

Rechnen log x behält in der Regel die gleiche Anzahl signifikanter Stellen wie xausser wenn x liegt in der Nähe 1: Sie verwenden möchten std::log1p wenn du rechnen musst prod(1 + x_n) mit klein x_n.

Schließlich, wenn Sie beim Summieren Rundungsfehler haben, können Sie verwenden Kahan-Summe oder Varianten.

Anstatt den sehr teuren Logarithmus zu verwenden, können Sie die Ergebnisse direkt mit Zweierpotenzen skalieren.

double geometric_mean(std::vector<double> const&data) {
  double huge = scalbn(1,512);
  double tiny = scalbn(1,-512);
  int scale = 0;
  double product = 1.0;
  for(auto x:data) {
    if (x >= huge) {
      x = scalbn(x, -512);
      scale++;
    } else if (x <= tiny) {
      x = scalbn(x, 512);
      scale--;
    }
    product *= x;
    if (product >= huge) {
      product = scalbn(product, -512);
      scale++;
    } else if (product <= tiny) {
      product = scalbn(product, 512);
      scale--;
    }
  }
  return exp2((512.0*scale + log2(product)) / data.size());
}