標本分散の計算方法による速度差と誤差について

数値計算

(2012-01-27追記) 同僚のご指摘で一部修正 (B会長に感謝!)

後輩のつぶやきがきっかけで標本分散の計算方法について気になって検証してみたのでメモ．

標本分散の計算方法は2つある．以下の定義式:
\frac{1}{N} \sum_i^N (x_i - \bar{x})^2
と，自乗の平均から平均の自乗を引くもの:
\frac{1}{N} \sum_i^N x_i^2 - (\bar{x})^2.

ちゃんとした呼び方を知らないので，本記事では前者を定義式，後者を簡易計算式と呼ぶことにする．計算機で標本分散を計算しようとすると，定義式では一度平均を計算する必要があり2回のループを回す必要があり後者に比べて時間がかかる．1パスで済む簡易計算式の方が計算は速いものの，先に自乗計算を行うため，桁落ちに伴う数値誤差が発生する．ここまではWikipediaにも書いてあること．僕の記憶が確かならば大学の統計の授業でも言われた気がする．

さて後輩のつぶやきにコメントしていて以下の3点が気になった

  • (1) 単精度ならともかく，倍精度なら桁落ち誤差とか気にしなくていいんじゃね?
  • (2) 定義式と簡易計算式では実際に速度差はどれくらいなんだろう?
  • (3) 単精度と倍精度ではどれくらい速度が変わるのか?

思い立ったが吉日ということでさっと実験コードを書いて検証．

予想

(この記事を書いているときには結果を知ってしまっているのだけれど) 実験前には以下の予想を立てた

  • (1) 倍精度ならば気にしなくてよい．
  • (2) 実際あまり変わらない
  • (3) 単精度の方がちょっと速い．

それぞれそう思った理由を述べてみる．倍精度 (符号付きdouble型) は32bit OSでも64bit OSでも8バイトで52bitの仮数部が存在する．(参考: http://www.drk7.jp/MT/archives/000851.html) そのため，10進数でいうところの約15桁程度を表現することができる．2乗することで桁落ちする可能性があるが，半分の7桁程度が保証される…のではなかろうか．数値計算の素人なので，ここらへん適当な予測．

単精度 (符号付きfloat型) の場合には仮数部は23bitなので保証される桁数は6-7桁といったところ．これを自乗してしまうと3桁程度になってしまうのではないか．これは少し影響がありそうな気がする．

(2)は，確かに2回ループが行われているが，中で行っている演算を考えるとそんなに変わらない気がする．N個データがあった場合には，

  • 定義式では
    • 平均の計算にN回の足し算
    • 各データと平均の差をN回計算 (= N回の足し算)
    • 上記差の自乗を計算 (= N回のかけ算) 上記差の自乗と和を計算 (= N回のかけ算 + N回の足し算) (2012-01-27修正)
    • よって2N回の足し算とN回のかけ算 よって3N回の足し算とN回のかけ算
  • 簡易計算式では
    • 自乗の平均の計算にN回のかけ算と足し算
    • 平均の計算にN回の足し算
    • よって2N回の足し算とN回のかけ算

当たり前かもしれないけれど演算回数はほぼ同じの気がする．N回の足し算分だけ演算コストの差が発生している (2012-01-27修正) なお，計算機的には引き算は足し算と同じはずなので (符号を反転して足し算)，足し算と解釈している．

なので，あとは定義式は2回ループを行っているのでカウンタ変数のインクリメントやデータをメモリからフェッチするオーバーヘッド分くらい遅くなるけれど，他のコストに比べて大したことない気がする．

(3)についてはCPU演算的には多分変わらないが，float型だとdouble型に比べて2倍の数だけキャッシュに乗るので，これによって差が出るのではないかと予想．

という予想を立てたうえで3つの疑問を解決するために実験を行った (実験コードは付録に記載)


実験環境は以下のとおり

結果

  • {double型, float型} X {定義式，簡易計算式} の組み合わせで検証
  • 100,000,000個の[0,1]乱数を発生させて標本分散の値と計算にかかった時間を計測
% g++ -O2 calc_var.cpp
% ./a.out
# double型 定義式
1/N \sum_i^N (x_i - \mu)^2:
var: 0.0833248
elapsed time: 0.291043 sec.

# double型 簡易計算式
1/N \sum_i^N x_i^2 - (1/N \sum_i^N x_i)^2:
var: 0.0833248
elapsed time: 0.151974 sec.

# float型 定義式
1/N \sum_i^N (x_i - \mu)^2:
var: 0.0833327
elapsed time: 0.218493 sec.

# float型 簡易計算式
1/N \sum_i^N x_i^2 - (1/N \sum_i^N x_i)^2:
var: 0.0833327
elapsed time: 0.111026 sec.

あれ，float型でも誤差が出ていない．


なお，-O2オプションを付け忘れた場合の結果は以下のとおり

% g++ calc_var.cpp
% ./a.out
# double型 定義式
1/N \sum_i^N (x_i - \mu)^2:
var: 0.0833248
elapsed time: 1.06947 sec.

# double型 簡易計算式
1/N \sum_i^N x_i^2 - (1/N \sum_i^N x_i)^2:
var: 0.0833248
elapsed time: 1.25177 sec.

# float型 定義式
1/N \sum_i^N (x_i - \mu)^2:
var: 0.149424
elapsed time: 1.04607 sec.

# float型 簡易計算式
1/N \sum_i^N x_i^2 - (1/N \sum_i^N x_i)^2:
var: 0.139625
elapsed time: 1.01841 sec.

O2オプションをつけないと，つけた場合に比べて

  • 定義式の方が速い!
  • float型において大きな誤差が発生する!

という2つの違いが起こった．特に後者についてはO2オプション付きとは挙動がかなり違う．ただ，当初の予想では単精度においてもっと誤差が出るはずだったので，むしろこちらの方が「ふつうの単精度」の挙動であって，O2オプション付きの場合はfloatも機械語に置き換えるときにdouble型扱いしているのではないか．

けれど計算速度も速くなっているからなぁ．．

どこかで聞いたこともあるような気がするのだけれど，正確なソースを見つけることができなかったのでとりあえず下記のページなどを読んでみた

後述の付録2にあるように環境によっても変わるらしい．うーん，ここらへんはきちんと仕組みを理解しないと理由がわからなさそう．今回はスルー

答え合わせ?

というわけで予想は大体合っているという結果を得た．

  • ○ (1) 倍精度ならば気にしなくてよい．
  • △ (2) 実際あまり変わらない
  • ○ (3) 単精度の方がちょっと速い．

(3)についてはキャッシュミスの差なのかということは面倒なので検証していない．O2オプションを付けない際の挙動があるので，別の理由があるのかもしれない．

眠くなってきたので特に考察はしない (投げやり)

結論

というわけで今回の実験から結論をまとめてみる

  • 簡易計算式(「自乗の平均 - 平均の自乗」)の方が速い．ただしO2オプションを忘れずに
  • 倍精度 (double型) を使えば簡易計算式でも誤差は気にせずともよい
  • 単精度 (float型) の場合には誤差を気を付ける．O2オプションつければ気にしなくてよい
  • 倍精度よりも単精度の方が速い．それなりに．

今回はデータ数が1億でこの程度の差なので，

  • データ数 x 繰り返し数のオーダーが1億くらいからようやく秒単位の差が出始める

ということを頭の片隅においておいて，無視できるほど小さい場合にはコードを書く時間が短い方を選択するという結論に至った．

付録: 実験に利用したコード

#include <iostream>
#include <cstdlib>
#include <vector>

#include <time.h>
#include <sys/time.h>

#define RAND ((double)rand()/(double)RAND_MAX)

double
gettimeofday_sec ()
{
  struct timeval tv;
  gettimeofday( &tv, NULL );
  return tv.tv_sec + tv.tv_usec * 1e-6;
}


template <class T> void
test (int num)
{
  std::vector<T> vec;
  for (int i = 0; i < num; i++) {
    vec.push_back( RAND );
  }

  // 計算1
  double t1 = gettimeofday_sec();

  T mean = 0.0;
  for (int i = 0; i < num; i++) {
    mean += vec[ i ];
  }
  mean /= (T)num;

  T var = 0.0;
  for (int i = 0; i < num; i++) {
    T v = mean - vec[ i ];
    var += v * v;
  }
  var /= (T)num;
  double t2 = gettimeofday_sec();

  std::cout << "1/N \\sum_i^N (x_i - \\mu)^2: " << std::endl;
  std::cout << "var: " << var << std::endl;
  std::cout << "elapsed time: " << t2 - t1 << " sec." << std::endl;
  std::cout << std::endl;

  // 計算2
  double t3 = gettimeofday_sec();

  T sum    = 0.0;
  T sq_sum = 0.0;
  for (int i = 0; i < num; i++) {
    sum    += vec[ i ];
    sq_sum += vec[ i ] * vec[ i ];
  }

  T m    = (sum / (T)num);
  T var2 = (sq_sum / (T)num) - (m * m);

  double t4 = gettimeofday_sec();

  std::cout << "1/N \\sum_i^N x_i^2 - (1/N \\sum_i^N x_i)^2:" << std::endl;
  std::cout << "var: " << var2 << std::endl;
  std::cout << "elapsed time: " << t4 - t3 << " sec." << std::endl;
  std::cout << std::endl;
}


int
main (int argc, char *argv[])
{
  test<double>( 100000000 );
  test<float>( 100000000 );

  return 0;
}

付録2

別の環境で実験した際にはO2オプションを付けてもfloat型のvar誤差が変わらなかったので結果をメモ．うーん，gccの影響なのか．

  • O2オプションなし
% g++ calc_var.cpp
rx7:/dev/shm% ./a.out
1/N \sum_i^N (x_i - \mu)^2:
var: 0.0833248
elapsed time: 1.10016 sec.

1/N \sum_i^N x_i^2 - (1/N \sum_i^N x_i)^2:
var: 0.0833248
elapsed time: 1.04039 sec.

1/N \sum_i^N (x_i - \mu)^2:
var: 0.149424
elapsed time: 1.06327 sec.

1/N \sum_i^N x_i^2 - (1/N \sum_i^N x_i)^2:
var: 0.139625
elapsed time: 1.03972 sec.
  • O2オプション有
% g++ -O2 calc_var.cpp
% ./a.out
1/N \sum_i^N (x_i - \mu)^2:
var: 0.0833248
elapsed time: 0.20865 sec.

1/N \sum_i^N x_i^2 - (1/N \sum_i^N x_i)^2:
var: 0.0833248
elapsed time: 0.105695 sec.

1/N \sum_i^N (x_i - \mu)^2:
var: 0.149424
elapsed time: 0.194086 sec.

1/N \sum_i^N x_i^2 - (1/N \sum_i^N x_i)^2:
var: 0.139625
elapsed time: 0.097558 sec.