c++ - relacion - Manera eficiente de calcular la media geométrica de muchos números

promedio geometrico formula

Question

Necesito calcular la media geométrica de un gran conjunto de números, cuyos valores no están limitados a priori. La manera ingenua sería

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());
}

Sin embargo, esto puede fallar debido a un desbordamiento o desbordamiento en el product acumulado (nota: el long double realmente no evita este problema). Entonces, la siguiente opción es resumir los logaritmos:

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());
}

Esto funciona, pero llama a std::log() para cada elemento, que es potencialmente lento. ¿Puedo evitar eso? ¿Por ejemplo, haciendo un seguimiento de (el equivalente de) el exponente y la mantisa del product acumulado por separado?

Answer 1

Creo que encontré una manera de hacerlo, combinando las dos rutinas en la pregunta, similar a la idea de Peter. Aquí hay un código de ejemplo.

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());
}

La mala noticia es: esto viene con una rama. La buena noticia: es probable que el predictor de rama consiga esto casi siempre correcto (la rama solo debe activarse rara vez).

La rama podría evitarse utilizando la idea de Peter de un número constante de términos en el producto. El problema con esto es que el desbordamiento / subdesbordamiento aún puede ocurrir dentro de unos pocos términos, dependiendo de los valores.

Answer 2

En lugar de utilizar logaritmos, que son muy caros, puede escalar directamente los resultados por potencias de dos. En pseudocódigo:

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

Answer 3

Es posible que pueda acelerar esto multiplicando números como en su solución original y solo convirtiendo a logaritmos cada cierto número de multiplicaciones (dependiendo del tamaño de sus números iniciales).

Answer 4

Hay una idea simple para reducir la computación y también para evitar el desbordamiento. Puede agrupar números, por lo menos dos a la vez, calcular su registro y luego evaluar su suma.

log(abcde) = 5*log(K)

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

Answer 5

La solución "exponente dividido y mantisa":

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);
}

Si le preocupa que ex puede desbordarse, puede definirlo como un doble en lugar de un long long , y multiplicar por invN en cada paso, pero podría perder mucha precisión con este enfoque.

EDITAR Para entradas grandes, podemos dividir el cálculo en varios grupos:

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;
}

Answer 6

Sumar registros para calcular productos de manera estable es perfectamente correcto y bastante eficiente (si no es suficiente: hay formas de obtener logaritmos vectorizados con algunas operaciones SSE, también existen las operaciones vectoriales de Intel MKL).

Para evitar el desbordamiento, una técnica común es dividir cada número por la entrada de magnitud máxima o mínima de antemano (o sumar las diferencias de log al log max o log min). También puede usar depósitos si los números varían mucho (por ejemplo, sume el registro de números pequeños y números grandes por separado). Tenga en cuenta que, por lo general, no se necesita nada de esto, excepto para conjuntos muy grandes, ya que el registro de un double nunca es enorme (entre, digamos, -700 y 700).

Además, es necesario realizar un seguimiento de los signos por separado.

Normalmente, el log x cálculo log x mantiene el mismo número de dígitos significativos que x , excepto cuando x está cerca de 1 : desea usar std::log1p si necesita calcular prod(1 + x_n) con un pequeño x_n .

Finalmente, si tiene problemas de error de redondeo al sumar, puede usar Kahan suma o variantes.

Answer 7

Un enfoque diferente que proporcionaría una mayor precisión y rendimiento que el método logarítmico sería compensar los exponentes fuera del rango en una cantidad fija, manteniendo un logaritmo exacto del exceso cancelado. Al igual que:

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);

Todas las multiplicaciones por BIG y SMALL son exactas, y no hay llamadas para log (una función trascendental, y por lo tanto particularmente imprecisa).