여러 숫자의 기하평균을 계산하는 효율적인 방법

여러 숫자의 기하평균을 계산하는 효율적인 방법

값이 선험적으로 제한되지 않는 큰 숫자 집합의 기하평균을 계산해야 합니다.순진한 방법은

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

는 누적된 로 인해 할 수 .product 고:long double실제로 이 문제를 피하지는 않습니다.)다음 옵션은 로그를 요약하는 것입니다.

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

은 효과가 을 라고 부릅니다.std::log()잠재적으로 느린 모든 요소에 대해.피하면 되나요?예를 들어, 누적된 지수와 가수를 추적(동등치)함으로써product따로?

"분할 지수와 가수" 솔루션:

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

이 그 을 한다면.ex할 수 있습니다.다 두 할 수 . 당신은 그것을 대신에 더블로 정의할 수 있습니다.long long , invN모든 단계에서, 하지만 이 접근법으로 많은 정밀도 말이죠.

EDIT 대용량 입력의 경우 계산을 여러 버킷으로 분할할 수 있습니다.

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

방법을 찾은 것 같아요 피터의 생각과 비슷하게 질문에 나오는 두 가지 루틴을 결합했어요여기 코드 예시가 있습니다.

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

나쁜 소식은 이것은 나뭇가지와 함께 온다는 것입니다.좋은 소식은 분기 예측 변수가 거의 항상 이 문제를 정확하게 파악할 가능성이 높다는 것입니다(분기는 거의 트리거되지 않아야 함).

제품에 항수가 일정하다는 피터의 아이디어를 이용하면 지점을 피할 수 있었습니다.문제는 오버플로/언더플로가 값에 따라 몇 개의 항 안에서 여전히 발생할 수 있다는 점입니다.

원래 솔루션에서와 같이 숫자를 곱하고 특정 곱셈 횟수마다 로그로 변환(초기 숫자의 크기에 따라 다름)함으로써 이를 가속화할 수 있습니다.

로그 방법보다 더 나은 정확성과 성능을 제공하는 다른 접근 방식은 범위를 벗어난 지수를 고정된 양으로 보상하고 취소된 초과의 정확한 로그를 유지하는 것입니다.이와 같습니다.

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

의 BIG그리고.SMALL하고,다에 .log(초월 함수이며, 따라서 특히 부정확한 함수).

계산량을 줄이고 오버플로우를 방지하는 간단한 방법이 있습니다.한 번에 두 개 이상의 숫자를 그룹화하여 로그를 계산한 다음 합계를 평가할 수 있습니다.

log(abcde) = 5*log(K)

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

안정적으로 제품을 계산하기 위해 로그를 합산하는 것은 완벽하게 괜찮고 효율적입니다(이 정도로 충분하지 않다면 몇 개의 SSE 연산으로 벡터화된 로그를 얻을 수 있는 방법이 있습니다). Intel MKL의 벡터 연산도 있습니다.

오버플로우를 방지하기 위해 일반적인 기술은 모든 숫자를 미리 최대 또는 최소 크기 항목으로 나누는 것입니다(또는 로그 최대 또는 로그 최소에 대한 합계 로그 차이).숫자가 많이 달라지는 경우(예: 작은 숫자와 큰 숫자의 로그를 별도로 합산) 버킷을 사용할 수도 있습니다.일반적으로 로그 이후 매우 큰 집합을 제외하고는 이 둘 중 어느 것도 필요하지 않습니다.double결코 크지 않습니다(예: -700에서 700 사이).

그리고 표시사항을 따로 파악해야 합니다.

컴퓨팅log x는 일반적으로 다음과 같은 유효 자릿수를 유지합니다.x, 때를 제외하고는x에 가깝습니다.1: 당신이 사용하고싶은std::log1p계산이 필요한 경우prod(1 + x_n)소소하게x_n.

마지막으로 합할 때 반올림 오류 문제가 있는 경우 Kahan 합이나 변형을 사용할 수 있습니다.

값이 매우 비싼 로그를 사용하는 대신 2의 거듭제곱으로 결과를 직접 축척할 수 있습니다.

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

언급URL : https://stackoverflow.com/questions/19980319/efficient-way-to-compute-geometric-mean-of-many-numbers