LoginSignup
8
4

More than 3 years have passed since last update.

@kaityo256(ロボ太)

モンテカルロ法における境界条件の扱いについて

Last updated at Posted at 2019-09-25

はじめに

古典スピン系では、もしクラスターアルゴリズムが使えるならその方が良いが、フラストレート系などではうまく使えない。その場合はメトロポリス法などのSingle Spin Flip Algorithmを使うことになる。この時、周期境界条件をどうやって補正すれば良いかずっと気になっていたので調べてみた。

コードは以下に置いておく。

問題意識

例えば正方格子上の、強磁性イジング模型を考えよう。普通にメトロポリス法でアップデートするとこんな感じのコードになると思う。

  for (int i = 0; i < N; i++) {
    int ns = 0;
    ns += spin[i + 1];
    ns += spin[i - 1];
    ns += spin[i + L];
    ns += spin[i - L];
    ns *= spin[i];
    if (ns < 0 || exp(-2.0 * ns * beta) > ud(mt)) {
      spin[i] = -spin[i];
    }
  }

上記のコードは境界をちゃんと扱っていないため、このままでは変なところを触ってしまう。

で、周期境界条件を貸すなら、こんな感じに補正しなければならない。

  if (ix < 0) ix += L;
  if (ix >= L) ix -= L;
  if (iy < 0) iy += L;
  if (iy >= L) iy -= L;

で、全スピンについて周期境界条件補正をするとこんな感じのコードになる。

int pos2index_if(int ix, int iy) {
  if (ix < 0) ix += L;
  if (ix >= L) ix -= L;
  if (iy < 0) iy += L;
  if (iy >= L) iy -= L;
  return ix + iy * L;
}

void single_flip_calc_if(double beta) {
  for (int i = 0; i < N; i++) {
    int ix = i % L;
    int iy = i / L;
    int ns = 0;
    ns += spin[pos2index_if(ix + 1, iy)];
    ns += spin[pos2index_if(ix - 1, iy)];
    ns += spin[pos2index_if(ix, iy + 1)];
    ns += spin[pos2index_if(ix, iy - 1)];
    ns *= spin[i];
    if (ns < 0 || exp(-2.0 * ns * beta) > ud(mt)) {
      spin[i] = -spin[i];
    }
  }
}

で、これ見ると「遅そうだな」と思えてくる。なのでここからコードを工夫するのだが、どの工夫がどれくらい効果的なのかいつも気になっていた。なので、一度真面目にコードの高速化の効果をまとめておこうと思う。

以下、このコードをcalc+ifと呼び、ベースラインとする。64x64の正方格子上で、全てのスピンを一回ずつ更新したら1モンテカルロステップ(MCs)とし、イジング模型の臨界点直上で100000 MCs計算するのにかかった時間を計る。計算はLinux+Xeonで行った(詳細は一番下に)。で、ベースラインの計算時間は16151 [msec]だった。ここから高速化してみよう。

高速化

calc + mod

まず思うのは、if文が遅いんじゃないか、ということだ。これをmodによる計算に置き換えてみよう。こんな感じになる。

int pos2index_mod(int ix, int iy) {
  ix = (ix + L) % L;
  iy = (iy + L) % L;
  return ix + iy * L;
}

void single_flip_calc_mod(double beta) {
  for (int i = 0; i < N; i++) {
    int ix = i % L;
    int iy = i / L;
    int ns = 0;
    ns += spin[pos2index_mod(ix + 1, iy)];
    ns += spin[pos2index_mod(ix - 1, iy)];
    ns += spin[pos2index_mod(ix, iy + 1)];
    ns += spin[pos2index_mod(ix, iy - 1)];
    ns *= spin[i];
    if (ns < 0 || exp(-2.0 * ns * beta) > ud(mt)) {
      spin[i] = -spin[i];
    }
  }
}

計算時間は 16064 [msec]となり、思ったより早くならなかった。

table + mod

現在、毎回指数関数の計算をしているが、イジング模型で出現する可能性のあるエネルギー差は数種類しかないので、これをテーブルに入れてしまおう。

void make_table(double beta) {
  for (int i = 0; i < 5; i++) {
    int ns = (i - 2) * 4;
    e_table[i] = exp(-ns * beta);
  }
}

```cpp
void single_flip_table_mod(double beta) {
  for (int i = 0; i < N; i++) {
    int ix = i % L;
    int iy = i / L;
    int ns = 0;
    ns += spin[pos2index_mod(ix + 1, iy)];
    ns += spin[pos2index_mod(ix - 1, iy)];
    ns += spin[pos2index_mod(ix, iy + 1)];
    ns += spin[pos2index_mod(ix, iy - 1)];
    ns *= spin[i];
    if (ns < 0 || e_table[ns / 2 + 2] > ud(mt)) {
      spin[i] = -spin[i];
    }
  }
}

こうすると 7488 [msec]となり、倍近く早くなった。

tableonly + mod

メトロポリス法では、まわりのスピンの和(局所磁場)との内積nsが負の場合、ひっくり返ったほうがエネルギーが下がるので必ずひっくり返し、エネルギーが増える向きには確率的にひっくり返す。

ここで、ns < 0の時には、exp(-2.0 * ns * beta)は常に1より大きい。そこで、if文の中身を

    if (ns < 0 || exp(-2.0 * ns * beta) > ud(mt)) {

から

    if (e_table[ns / 2 + 2] > ud(mt)) {

と簡略化することができる。条件分岐が簡単になるので少し早くなるかと思ったのだが、結果は8009 [msec]と若干遅くなってしまった。

table + list

毎回隣接するスピンを計算で得るのではなく、最初に隣接リストを作ってみよう。

std::vector<std::array<int, 4>> neighbor(N);

void init_neighbors() {
  for (int i = 0; i < N; i++) {
    int ix = i % L;
    int iy = i / L;
    neighbor[i][0] = pos2index_mod(ix + 1, iy);
    neighbor[i][1] = pos2index_mod(ix - 1, iy);
    neighbor[i][2] = pos2index_mod(ix, iy + 1);
    neighbor[i][3] = pos2index_mod(ix, iy - 1);
  }
}

隣接リストを使うと、スピンの更新はこんな感じに書ける。

void single_flip_table_list(double beta) {
  for (int i = 0; i < N; i++) {
    int ns = 0;
    ns += spin[neighbor[i][0]];
    ns += spin[neighbor[i][1]];
    ns += spin[neighbor[i][2]];
    ns += spin[neighbor[i][3]];
    ns *= spin[i];
    if (ns < 0 || e_table[ns / 2 + 2] > ud(mt)) {
      spin[i] = -spin[i];
    }
  }
}

結果は 6815 [msec]となり、若干早くなった。

bulk + border

さて、これが本命である。境界条件補正は端にあるスピンのみに必要であり、バルク、つまり 1 < x < L-1 かつ 1 < y < L-1であるようなスピンに補正は不要である。

したがって、バルクにあるスピンと、端(border)にあるスピンに分けてループをまわすと、不必要な境界補正コードが減る分、早くなることが期待される。というわけでバルクは補正なし、端は隣接リストを使って更新してみる。コードはかなり長くなる。

void single_flip_bulk(double beta) {

  // Bulk
  for (int iy = 1; iy < L - 1; iy++) {
    for (int ix = 1; ix < L - 1; ix++) {
      int i = ix + iy * L;
      int ns = 0;
      ns += spin[i + 1];
      ns += spin[i - 1];
      ns += spin[i + L];
      ns += spin[i - L];
      ns *= spin[i];
      if (ns < 0 || e_table[ns / 2 + 2] > ud(mt)) {
        spin[i] = -spin[i];
      }
    }
  }

  // site_i = (j, 0)
  for (int j = 0; j < L; j++) {
    int i = j;
    int ns = 0;
    ns += spin[neighbor[i][0]];
    ns += spin[neighbor[i][1]];
    ns += spin[neighbor[i][2]];
    ns += spin[neighbor[i][3]];
    ns *= spin[i];
    if (ns < 0 || e_table[ns / 2 + 2] > ud(mt)) {
      spin[i] = -spin[i];
    }
  }

  // site_i = (j, L-1)
  for (int j = 0; j < L; j++) {
    int i = j + (L - 1) * L;
    int ns = 0;
    ns += spin[neighbor[i][0]];
    ns += spin[neighbor[i][1]];
    ns += spin[neighbor[i][2]];
    ns += spin[neighbor[i][3]];
    ns *= spin[i];
    if (ns < 0 || e_table[ns / 2 + 2] > ud(mt)) {
      spin[i] = -spin[i];
    }
  }
  // site_i = (0, j)
  for (int j = 1; j < L - 1; j++) {
    int i = j * L;
    int ns = 0;
    ns += spin[neighbor[i][0]];
    ns += spin[neighbor[i][1]];
    ns += spin[neighbor[i][2]];
    ns += spin[neighbor[i][3]];
    ns *= spin[i];
    if (ns < 0 || e_table[ns / 2 + 2] > ud(mt)) {
      spin[i] = -spin[i];
    }
  }

  // site_i = (L-1, j)
  for (int j = 1; j < L - 1; j++) {
    int i = L - 1 + j * L;
    int ns = 0;
    ns += spin[neighbor[i][0]];
    ns += spin[neighbor[i][1]];
    ns += spin[neighbor[i][2]];
    ns += spin[neighbor[i][3]];
    ns *= spin[i];
    if (ns < 0 || e_table[ns / 2 + 2] > ud(mt)) {
      spin[i] = -spin[i];
    }
  }
}

結果は 6595 [msec]となり、わずかではあるが最も早かった。

結果のまとめ

上記の結果も含め、いくつかの環境、コンパイラで行ったベンチマークの結果を貼っておく。

それぞれ名前、自発磁化、経過時間である。

Linux + GCC

  • g++ (GCC) 7.2.0
  • Intel(R) Xeon(R) Gold 6148 CPU @ 2.40GHz
  • -O3 -std=c++11 -march=native
  • Red Hat Enterprise Linux Server release 7.4 (Maipo)
calc      + if      0.3537  16151[msec]
calc      + mod     0.3537  16064[msec]
table     + mod     0.3537  7488[msec]
tableonly + mod     0.2853  8009[msec]
table     + list    0.3537  6815[msec]
bulk + border       0.5839  6595[msec]

Linux + Intel Compiler

  • icpc (ICC) 18.0.5 20180823
  • Intel(R) Xeon(R) Gold 6148 CPU @ 2.40GHz
  • -O3 -std=c++11 -xHOST
  • Red Hat Enterprise Linux Server release 7.4 (Maipo)
calc      + if      0.3537  33342[msec]
calc      + mod     0.3537  32431[msec]
table     + mod     0.3537  6152[msec]
tableonly + mod     0.2853  6458[msec]
table     + list    0.3537  30158[msec]
bulk + border       0.5839  29927[msec]

Mac + Clang

  • Apple clang version 11.0.0 (clang-1100.0.33.8)
  • Intel Core i5 3.3 GHz
  • macOS Mojave 10.14.6
  • -O3 -std=c++11 -march=native
calc      + if      0.3537  8778[msec]
calc      + mod     0.3537  7553[msec]
table     + mod     0.3537  5622[msec]
tableonly + mod     0.2853  5980[msec]
table     + list    0.3537  5574[msec]
bulk + border       0.5839  5398[msec]

まとめ

周期境界条件補正について、隣接リストを使うのと、バルクと端で分けてループでまわすのとどっちがどれくらい早いか確認した。大昔、おそらくAlpha 21164だったと思うが、隣接リストを使った間接アクセスがかなり遅く、またmod演算もまぁまぁ遅かったので、バルクと端で分けたループを書くのがデフォルトだった記憶があった。しかし時は流れ、x86を使うようになると、間接アクセスやmodが早いので、ループを分けると早くはなるが、コードが複雑にするコストをかけるだけのメリットは得られない。

個人的に以外だったのは、mod演算より隣接リストを使った間接アクセスの方が早かったこと。しかし、インテルコンパイラを使うとなぜか間接アクセスが相当遅くなる。このあたり、アセンブリ見ていないので何が起きてるのかはわからない。まぁ、x86でGCC/Clang使うなら、「テーブルと隣接リストを使う」のがいまのところベストプラクティスではないかな、と思う。

8
4
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
8
4