周期カーネルに対するFeature Map

ガウス過程回帰の課題のひとつとしてスケーラビリティがあります。
データ数をNとしてfitting の時間計算量 O(N^3)、predict の時間計算量が O(N^2)であることから大量のデータを扱いたいときは何らかの近似が必要です。

そのような近似手法のひとつとして、カーネルを有限次元の特徴ベクトルの内積で近似してベイズ線形回帰に帰着するものがあります。
特徴量ベクトルの次元をDとして、ベイズ線形回帰のfitting の時間計算量は O(D^3) 、predict の時間計算量は O(D)なので特徴量ベクトルの次元がデータ数よりも小さければより高速に計算ができることになります。

与えられたカーネルに対応する特徴量ベクトルを作る方法としてRandom Fourier Features
があります。*1
Gaussian カーネルやLaplacian カーネルについてはRandom Feature Map で特徴量ベクトルを得られますが、周期カーネルにこの方法を適用するのは難しそうでした。

そこでいろいろ調べたところ、周期カーネルを近似する論文*2を見つけたので概要と実装を備忘録としてまとめておきます。

(追記: 一部の数式が上手く表示されないようなので画像で置き換えました。何もわからない。)

Random Fourier Features

そもそもRandom Fourier Features がどのようなアルゴリズムなのか簡単に述べておきます。
 \mathbb{R}^{D}上のstationaryなカーネル k(\mathbf{x}, \mathbf{y}) = k(\mathbf{x} - \mathbf{y})がpositive definite であることと、ある有限非負測度のフーリエ変換で書けることが同値である(Bochner の定理)ことを使います。

Bocher の定理より、適当に定数倍したpositive definite なカーネルはある確率密度関数 p(\mathbf{\omega})フーリエ変換で書けます。
\begin{align}
k(\mathbf{x} - \mathbf{y})
&= \int_{\mathbb{R}^{D}} p(\mathbf{\omega}) e^{i (\mathbf{x} - \mathbf{y})^{T} \mathbf{\omega}} d\mathbf{\omega} \\
&= \mathbb{E}_{\mathbf{\omega}, b} \left[ z_{\mathbf{\omega}}(\mathbf{x}) z_{\mathbf{\omega}}(\mathbf{y})\right]
\quad (\mbox{where} \quad z_{\mathbf{\omega}}(\mathbf{x}) := \sqrt{2} \cos(\mathbf{x}^T \mathbf{\omega} + b))
\end{align}
最後の式において、 \mathbf{\omega} p(\mathbf{\omega})に従い、bは [0, 2\pi]の一様分布に従います。

 \mathbf{\omega}とbをそれぞれD個ずつサンプリングして
\begin{align}
\mathbf{z}(\mathbf{x}) := \sqrt{\frac{2}{D}} \left( \cos(\mathbf{x}^T \mathbf{\omega}_{1} + b_{1}), \cdots, \cos(\mathbf{x}^T \mathbf{\omega}_{D} + b_{D})\right)^{T}
\end{align}
とすると、カーネルモンテカルロ積分の要領で近似できます。
\begin{align}
k(\mathbf{x} - \mathbf{y})
&\simeq \frac{1}{D} \sum_{j = 1}^{D} z_{\mathbf{\omega}_{j}}(\mathbf{x}) z_{\mathbf{\omega}_{j}}(\mathbf{y}) \\
&= \mathbf{z}(\mathbf{x})^{T} \mathbf{z}(\mathbf{y})
\end{align}

例えばGaussian カーネル正規分布フーリエ変換で書けるのでこの方法で対応する特徴量ベクトルを得ることができます。
これは既にsklearn に実装されています。*3

Fourier Features for Periodic Kernels

次に周期カーネルを考えます。
具体的には次のようなカーネルです。*4
\begin{align}
k_{perSE}(t, t') = k_{perSE}(t, t')
&= \exp \left( -2 \left( \frac{\sin (\omega_{0} (t - t'))}{l} \right)^{2} \right) \\
&= \exp \left( \frac{\cos \omega_{0} (t - t') - 1}{l^{2}} \right)
\end{align}
ここで\omega_{0} = \frac{2\pi}{T}で、このカーネルは周期Tの関数を説明したい時などに使われます。

Random Fourier Features で特徴量ベクトルを作るなら、周期カーネルを逆フーリエ変換して得られる確率密度関数に従って周波数をサンプリングすることになりますが、どうみても簡単にサンプリングができる確率分布を得られそうにありません。
そこで、周期カーネルフーリエ級数展開して近似することを考えます。

周期カーネルフーリエ級数展開を次のように書きます。
f:id:lan496:20180716013606p:plain

k_{perSE}(\tau)は偶関数であることと、第一種変形ベッセル関数の積分表示I_{n}(x) = \frac{1}{\pi} \int_{0}^{\pi} e^{x \cos \theta} \cos n \theta\, d\theta \, (n \in \mathbb{Z})
*5を使うと、 c_{k}は次のように計算できます。
f:id:lan496:20180716013619p:plain

フーリエ級数の第K項以降を切り捨てて近似します。
f:id:lan496:20180716013630p:plain
よって、周期カーネルを2K次元の特徴量ベクトル \mathbf{\phi}(t)で近似できることがわかりました。

実装とテスト

論文の公式実装が見当たらなかったので自分で書きました。*6

論文ではランダムに2000個の入力に対して、周期T=2の周期カーネルでグラム行列 \mathbf{K}を計算したものと、第C項までを使って近似で得られる行列 \tilde{\mathbf{K}}のフロベニウスノルムでの誤差 \frac{||  \tilde{\mathbf{K}} -
 \mathbf{K}||_{F}}{|| \mathbf{K} ||_{F}} を測っていたので、自分の実装で検証してみました。

特徴量ベクトルの次元とハイパーパラメターlength_scale(l)を変えた時の誤差は下図のとおりです。
図の左上部分の誤差が元論文と異なりますが、傾向としては同じ結果が得られました。
f:id:lan496:20180716011413p:plain

下図はlength_scale(l)=0.2の場合にカーネルの値をプロットしたものです。
フーリエ級数展開の16次まで取り込めば、元の周期カーネルをほぼ再現できていることがわかります。
f:id:lan496:20180716012425p:plain

化合物の構造を表現する記述子

無機化合物の化学式とユニットセルの情報及び原子の座標から物性値(結合エネルギー・バンドギャップ)を予測するKaggleコンペに参加していた。

Nomad2018 Predicting Transparent Conductors | Kaggle

feature engineering の過程で分子および無機結晶の構造から特徴量抽出する手法に関する論文をいろいろ読んだので、以下にまとめておく。 基本的に記述子の作り方の部分しか読んでいないので、詳しくは論文を参照。

分子に対する手法

Coulomb Matrix (CM)

ここで Z_{I}は原子Iの形式電荷 r_{IJ}は原子Iと原子Jの距離。

  • 分子を構成する原子Iと原子Jの間の反発力を表している。 対角成分については、原子の自己エネルギーを多項式でフィッティングしたもの。

  • Coulomb Matrix  M固有値を絶対値の大きい順で並べ、適当に0でpadding して次元を揃えたものを記述子として使う(分子ごとに原子数が異なるので記述子の次元を揃える必要がある)。

  • 分子の回転・並進移動の元で記述子は不変。

Bag of Bonds (BoB)

Bonds, angles, machine learning (BAML)

Extended Connectivity Fingerprints (ECFP)

Molecular atomic radial angular distribution (MARAD)

Histogram of distances, angles and dihedral angles (HDAD)

Fourier series of radial distance distributions

Molecular Graph (MG)

無機結晶に対する手法

分子に対する手法ほど進展していない印象。 分子の場合は並進移動と回転に対する対称性を考えればよかったが、結晶の場合は更にユニットセルの選び方に関して不変であるべきなのでより難しそう。 crystal net *2を一意に表す方法があればいいのだが、そもそもグラフ同型判定問題はNPに属しているという厳しさがある。

Partial radial distribution function (PRDF)

Atom centered symmetry functions

 f_{c}(R) = \begin{cases} \frac{1}{2} \left( 1 + \cos \left( \frac{\pi R}{R_{cutoff}} \right) \right) \quad (R < R_{cutoff}) \\  0 \quad (R > R_{cutoff}) \end{cases}

  • 各原子 iに対して次のような動径分布関数を考える。  \eta, R_{s}をいろいろ変えて記述子とする。

 G_{i}^{rad} (\eta, R_{s})= \sum_{j} e^{-\eta (R_{ij} - R_{s})^{2} } f_{c} (R_{ij})

  • 各原子 iに対して次のような角度分布関数を考える。  \zeta, \etaをいろいろ変えて記述子とする。

 G_{i}^{ang} (\zeta, \eta) = \frac{1}{2^{\zeta - 1}} \sum_{j, k \neq i} \left( 1 \pm \cos \theta_{ijk} \right)^{\zeta} \, e^{-\eta (R_{ij}^{2} + R_{jk}^{2} + R_{ki}^{2})} f_{c} (R_{ij}) f_{c} (R_{jk}) f_{c} (R_{ki})

Bond orientational order parameter (BOP, BOOP)

Bispectrum

F-Fingerprint

 F_{AB}(R) = \sum_{A_{i} \in unit cell} \sum_{B_{j} s.t. R_{ij} \lt R_{max}} \frac{\delta (R - R_{ij}) }{4 \pi R_{ij}^{2} \frac{N_{A} N_{B}}{V} \Delta}  - 1

  • ここで、 R_{ij}は原子 A_{i}と原子 B_{j}の距離、 N_{A}, N_{B}はユニットセル内の原子A,Bの個数、 Vはユニットセルの体積、 \Deltaはこの関数をヒストグラムにするビン幅、 \delta (\cdot )はgaussian-smeared delta function

Ewald Sum Matrix

  • [1503.07406] Crystal Structure Representations for Machine Learning Models of Formation Energies
  • Coulomb Matrix における M_{IJ} をサイトIを占めるすべての原子とサイトJを占めるすべての原子との間の静電ポテンシャルの和に修正。
  • Ewald法を用いて静電ポテンシャルの総和を計算する。
  • pymatgen にEwald法を用いて結晶の静電ポテンシャルの総和を計算する関数があるが、計算する式が少し違う。ただEwald Sum Matrix を実際に実装するときに大いに参考にはなる。

Extended Coulomb-like Matrix

Sine Matrix

Crystal Graph (CG)

マルチカノニカル法で個数推定

はじめに

去年の11,12月にマルコフ連鎖モンテカルロ法(MCMC)の集中講義を受けてきました。 講義の中で紹介されたマルチカノニカル法を応用してレアイベントをサンプリングする手法が特に面白かったので、復習も兼ねて実装してみました。

マルチカノニカル法

まずレアイベントサンプリング用のマルチカノニカル法を説明します。

状態xの確率分布が\pi(x)が与えられているときに、統計量\xi (x)の分布が知りたいとします。 単純にxをサンプリングしてもいいのですが、統計量の周辺分布P(\xi)に極めて確率の小さい区間があるとその部分をうまくサンプリングすることができません。

この問題を回避するために適当な重みQ(\xi )を用意して、\pi (x)の代わりに重みQ(\xi (x)) \pi(x) に従ってxをサンプリングすることを考えます。 これは普通のメトロポリスヘイスティングス法で実現できます。 すなわち、状態xから状態x' に遷移する確率を

 \mathrm{min}(1, \frac{Q(\xi (x')) \pi(x') p(x \to x')}{Q(\xi (x)) \pi(x)
p(x' \to x)})

とすればいいわけです。

重みQ(\xi )として、求めたい周辺分布P(\xi)に近い分布\tilde{P}(\xi)の逆数を用いるのがマルチカノニカル法です。

重み \tilde{P}(\xi (x))^{-1} \pi (x)に従って xをサンプリングしたのときの統計量 \xi の周辺分布 P'(\xi )

 P'(\xi ) \propto P(\xi)  \tilde{P}(\xi )^{-1}

となって、この分布はほぼ平坦であることがわかります。 平坦な周辺分布 P'(\xi )からサンプリングすることで、普通にサンプリングすると確率0と推定されてしまうようなレアイベントにもMCMCで訪れられるのがポイントです。

Wang-Landau 法

残る問題は、求めたい周辺分布P(\xi)に近い分布\tilde{P}(\xi)をどのように得るかです。

周辺分布P(\xi)を大まかに求める方法としてWang-Landau 法というものがあります。 詳細はwiki擬似コードを読むのが一番分かりやすいと思います(https://en.wikipedia.org/wiki/Wang_and_Landau_algorithm)。

大まかに言うと、 \xiヒストグラムが平坦になるように出現回数の低いところの重みを逐次大きくしていき、ある程度平坦になったらより細かく重みを調整していく感じです。

演習問題

具体例として、次のような数え上げ問題を考えます。

L×Lの正方格子において、角と向かい側の角とをつなぐ自己交差のない経路の個数を求めよ。

文章だと分かりづらいですが、要は↓の動画で扱われている問題です。 www.youtube.com

一辺の長さがLのとき、辺は全部で2L(L+1)本あるので、各辺について使うか否かを考えると辺の配置は全部で2^{2L(L+1)}通りです。探索する配置の数が指数関数で増大するので、愚直に全通りを試していてはすぐに計算が終わらなくなってしまいます。

ただ、この問題には効率的なアルゴリズムがあり、L=26までの厳密解が求まっているようです(http://oeis.org/A007764?language=japanese)。 なので個数推定をする意義はあまりないのですが、演習問題を解いているということでそこら辺は気にせずいきます。

この問題をマルチカノニカル法で扱うために、正方格子の各辺について使うか否か決めたものC (以下これを「状態C 」と呼ぶ)に対してエネルギーを次式で定めます。

 \displaystyle {
E(C) = \sum_{i : \mbox{左上と右下以外の頂点}} |\mbox{(iの次数)} - 2| + \sum_{i : \mbox{左上と右下の頂点}} |\mbox{(iの次数)} - 1|
+ (\mbox{2個以上の頂点からなる連結要素の数)}
}

ここで、E(C) = 0であることと状態Cが左上と右下の頂点をつなぐ自己交差のない経路であることが同値になっています。 先のマルチカノニカル法を用いてエネルギーの分布P(E)を求めて、P(E=0)に状態の総数2^{2L(L+1)}をかければ今求めたい経路の数を推定することができます。

具体的な実装はこちら。 github.com

結果

以上の方法で自己交差のない経路の数を100回推定し、その平均と不偏分散の平方根を真値と比較したものが下図です。 有効数字2桁くらいの精度が得ることができました。

ちなみにL=6が手元のMacBook Air で4時間くらいかかっています。 L=6の状態の数が 2^{2 \times 6 \times 7} \sim 3 \times 10^{50}で、その中で自己交差のない経路になっているものはおよそ 6.4 \times 10^{42}分の1しかないことを考えると、これほどレアな事象をそれなりの精度でサンプリングできるのはすごいと思えます。

L 推定値 誤差 真値
3 185.072 2.6228 184
4 8241.03 275.432 8512
5 1.2242e+06 75237.9 1262816
6 5.58901e+08 6.03206e+07 575780564

参考文献

http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0125062 https://sites.google.com/site/iwanamidatascience/vol2/number-place http://www.springer.com/us/book/9783642031625

AOJ 1320 City Merger

問題文(http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=1320)

AOJ2022(http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=2022)とほとんど同じ問題。他の都市名の部分文字列になっている都市名を予め取り除いてからbitDPをする。

#include <bits/stdc++.h>
using namespace std;

#define REP(i,n) for(int i=0;i<(int)(n);i++)
#define ALL(x) (x).begin(), (x).end()

typedef long long ll;
typedef long double ld;

const int INF = 1e9;
const ld EPS = 1e-8;

bool is_contained(string s, string t) {
  for(int i = 0; i + s.length() - 1 < t.length(); ++i) {
    if(s == t.substr(i, s.length())) return true;
  }
  return false;
}

int d[16][16];
int memo[16][1<<16];

int dp(int n, int k, int S) {
  //cout << k << " " << S << endl;
  if(memo[k][S] != INF) return memo[k][S];
  if(S == (1<<k)) return d[k][k];
  int res = INF;
  REP(i, n){
    if(i == k) continue;
    if(S & (1<<i)){
      res = min(res, d[k][i] + dp(n, i, S - (1<<k)));
    }
  }
  //cout << k << " " << S << ": " << res << endl;
  return memo[k][S] = res;
}

int main(){
  int n;
  while(cin >> n && n) {
    vector<string> cities(n);
    REP(i,n) cin >> cities[i];
    vector<string> tmp;
    REP(i,n){
      bool flag = true;
      REP(j,n){
        if(i == j) continue;
        if(is_contained(cities[i], cities[j])){
          flag = false;
          break;
        }
      }
      if(flag) tmp.push_back(cities[i]);
    }
    cities = tmp;
    n = cities.size();

    REP(i,16)REP(j,16) d[i][j] = INF;
    REP(i,n)REP(j,n){
      string s = cities[i];
      string t = cities[j];
      bool flag = true;
      for(int k = max(0, (int)(s.length() - t.length())); k < s.length(); ++k){
        if(s.substr(k) == t.substr(0, s.length() - k)){
          d[i][j] = k;
          flag = false;
          break;
        }
      }
      if(flag) d[i][j] = s.length();
    }
    REP(i,n) d[i][i] = cities[i].length();
    
    int res = INF;
    REP(i,16)REP(j,1<<16) memo[i][j] = INF;
    REP(i,n){
      int tmp = dp(n, i, (1<<n) - 1);
      res = min(res, tmp);
    }
    cout << res << endl;
  }
  return 0;
}

AOJ 2534 Dictionary

問題文(http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=2534)

まず、与えられた文字列の組(string_i, string_j)(i < j)が辞書式順序を満たすようにアルファベットの間に順序を入れていく。アルファベットx, yに順序が入っていることは有向グラフに(x, y)の辺を張ることで表現する。文字列の長さが違うときは'@'(='a' - 1)で埋めておいた。

このグラフの構築段階で(i,j)に辺があるのに(j,i)にも辺が張られているば矛盾しているのでfalseを返す。また、'@'は最小であるはずなのでそれが満たされていない時もfalseを返す。

'@'を始点として辺の張られているアルファベットに深さ優先探索で移動していくとき、もしも矛盾がなければ高々深さは26にしかならないから、それ以上の深さを探索しようとしていたらfalseを返す。

#include <bits/stdc++.h>
using namespace std;

#define REP(i,n) for(int i=0;i<(int)(n);i++)

bool dfs(int d[][27], int k, int cnt) {
  if(cnt > 30) return false;
  bool flag = true;
  REP(j,27) if(d[k][j] == 1) flag = false;
  if(flag) return true;

  bool res = true;
  REP(j,27){
    if(d[k][j] == 1){
      res = res && dfs(d, j, cnt + 1);
    }
  }
  return res;
}

bool solve(vector<string> &s) {
  int n = s.size();
  int d[27][27] = {0};
  REP(i,n)for(int j = i + 1; j < n; ++j) {
    int pos = 0;
    while(pos < s[i].length() && s[i][pos] == s[j][pos]) pos++;
    if(pos == s[i].length()) continue;
    if(d[s[j][pos]-'a'+1][s[i][pos]-'a'+1] == 1) return false;
    d[s[i][pos]-'a'+1][s[j][pos]-'a'+1] = 1;
  }
  REP(i,27) if(d[i][0] == 1) return false;
  
  return dfs(d, 0, 0);
}

int main(){
  int n;
  while(cin >> n && n) {
    vector<string> s(n);
    REP(i,n) cin >> s[i];
    int lmx = 0;
    REP(i,n) lmx = max(lmx, (int)s[i].length());
    REP(i,n)REP(j,lmx-s[i].length()) s[i] += 'a' - 1;
    if(solve(s)) cout << "yes" << endl;
    else cout << "no" << endl;
  }
  return 0;
}

Codeforces Round #341 (Div. 2) D. Rat Kwesh and Cheese

問題文はこちら(Problem - D - Codeforces)

まず\( (x^y)^z = (x^z)^y \)であることから\(a_4, a_8, a_{12}\)が候補から外れる。
単純に与えられた式を計算するとオーバーフローしてしまうため、対数をとって比較していく。ただし、\(x < 1\)のとき\( \mathrm{log} x < 0 \)より\( \mathrm{log} \mathrm{log} x\)は(実数上で)存在しない。よってx, y, zの中に1より小さいものが含まれているかどうかで場合分けをする。

(i) x,y,zが全て1より大きい時

\( \{\mathrm{log} \mathrm{log} a_i\}\)の中で最大のものを探す。ただし、最大値を取るものが複数あるときはindexの最も小さいものを選ぶ。

(ii) x, y, zのいずれかが1よりも大きい時

1よりも小さい数が底ならば、今の場合指数は必ず正であるから最終的な値は1よりも小さくなる。逆に1以上の数が底ならば、最終的な値は1以上になる。従って、底が1よりも大きいもののみを(i)と同じように比較すればよい(底が1より小さいものは適当に-INFとかで埋めておく)。

(iii) x,y,zが全て1以下の時

オーバーフローする心配が無いため、\( \mathrm{log} a_i\) の中で最大のものを探せばよい。


コンテストでは\( \mathrm{log} \mathrm{log} x\)が死にうることに気づかず提出してしまった。

#include <bits/stdc++.h>
using namespace std;

#define REP(i,n) for(int i=0;i<(int)(n);i++)
#define ALL(x) (x).begin(), (x).end()

typedef long long ll;
typedef long double ld;

const int INF = 1e9;
const ld EPS = 1e-8;

bool comp(pair<ld,int> p1, pair<ld,int> p2) {
  if(p1.first != p2.first) return p1.first < p2.first;
  return p1.second > p2.second;
}

int main(){
  double x, y, z;
  cin >> x >> y >> z;
  string ans[] = {"x^y^z", "x^z^y", "(x^y)^z", "y^x^z", "y^z^x", "(y^x)^z", "z^x^y", "z^y^x", "(z^x)^y"};
  vector<pair<ld,int>> val(9);

  if(x > 1.0 || y > 1.0 || z > 1.0) {
    if(x > 1.0) {
      val[0] = make_pair(z * log(y) + log(log(x)), 0);
      val[1] = make_pair(y * log(z) + log(log(x)), 1);
      val[2] = make_pair(log(y * z) + log(log(x)) , 2);
    }else{
      val[0] = make_pair(-INF, 0);
      val[1] = make_pair(-INF, 1);
      val[2] = make_pair(-INF, 2);
    }
    if(y > 1.0) {
      val[3] = make_pair(z * log(x) + log(log(y)), 3);
      val[4] = make_pair(x * log(z) + log(log(y)), 4);
      val[5] = make_pair(log(z * x) + log(log(y)) , 5);
    }else{
      val[3] = make_pair(-INF, 3);
      val[4] = make_pair(-INF, 4);
      val[5] = make_pair(-INF, 5);
    }
    if(z > 1.0) {
      val[6] = make_pair(y * log(x) + log(log(z)), 6);
      val[7] = make_pair(x * log(y) + log(log(z)), 7);
      val[8] = make_pair(log(y * x) + log(log(z)) , 8);
    }else{
      val[6] = make_pair(-INF, 6);
      val[7] = make_pair(-INF, 7);
      val[8] = make_pair(-INF, 8);
    }
  }else{
    val[0] = make_pair(pow(y, z) * log(x), 0);
    val[1] = make_pair(pow(z, y) * log(x), 1);
    val[2] = make_pair(y * z * log(x), 2);
    val[3] = make_pair(pow(x, z) * log(y), 3);
    val[4] = make_pair(pow(z, x) * log(y), 4);
    val[5] = make_pair(x * z * log(y), 5);
    val[6] = make_pair(pow(x, y) * log(z), 6);
    val[7] = make_pair(pow(y, x) * log(z), 7);
    val[8] = make_pair(x * y * log(z), 8);
  }

  sort(ALL(val), comp);
  cout << ans[val[8].second] << endl;
  return 0;
}

焼きなまし法でThomson問題の解を探索する

Thomson問題

Thomson問題とは、三次元単位球の表面にN個の電荷の等しい粒子を配置するとき、クーロンポテンシャル\( U = \sum_{i < j} \frac{1}{|\rm{r}_{i} - \rm{r}_{j}|}\)が最小となる配置となるを求める問題です。詳しくはwikiThomson problem - Wikipedia, the free encyclopedia。N=5のときにて三方両錐形が解になったり、N=8で立方体が解にならないなどいろいろと面白い問題です。

去年の12月頃、そもそもこの問題設定に名前が付いていることをある方から教えていただき、N=5の場合で実際に初期条件を与えて時間発展させることで本当に三方両錐形になるのか実験してみました。

一応それらしい配置は得られましたが、Uの値が初期値を変える度に大きく変化するのが目につきました。後日、Uの値はNに関して指数関数的に極小点をもつようになり、一般にN粒子のThomson問題はNP困難であることを教えてもらいました。
(参考)http://estamine.net/fc0809/JJFerreira-fisica_computacional.pdf

今回、焼きなまし法で再びこの問題について考えてみました。

焼きなまし法

焼きなまし法自体については以下を主に参照しました。
焼きなまし法 - Wikipedia
shindannin.hatenadiary.com

N個の粒子それぞれの座標を極座標表示したときの偏角の組を状態とします。
近傍状態としては、この2N個の値の中からランダムに一つ選び、平均値が元の値で標準偏差が\( \sqrt{2T}\) (Tは温度)であるような正規分布に従ってランダムに変化させたものを取ります。
また、焼きなましスケジュールは温度が毎回0.9倍になるようにします。

実際のコードは以下の通り

#include <iostream>
#include <vector>
#include <cmath>
#include <climits>
#include <stdexcept>
#include <ctime>
#include <fstream>
using namespace std;

unsigned int randxor(){
  static unsigned int x = 123456789, y = 362436069, z = 521288629, w = 88675123;
  unsigned int t;
  t = (x ^ (x << 11)); x = y; y = z; z = w;
  return (w = (w ^ (w >> 19)) ^ (t ^ (t >> 8)));
}

double uniform(){
  return (double)(randxor()) / UINT_MAX;
}

double rand_normal(double mu, double sigma){
  double z = sqrt(-2.0 * log(uniform())) * sin(2.0 * M_PI * uniform());
  return mu + sigma * z;
}

// if two particles are too close, then throw exception
double eval(vector<double>& state){
  int n = state.size() / 2;
  double energy = 0.0;
  for(int i=0;i<n;++i){
    for(int j=i+1;j<n;++j){
      double dst = sqrt(
        pow(sin(state[2*i]) * cos(state[2*i+1]) - sin(state[2*j]) * cos(state[2*j+1]), 2)
        + pow(sin(state[2*i]) * sin(state[2*i+1]) - sin(state[2*j]) * sin(state[2*j+1]), 2)
        + pow(cos(state[2*i]) - cos(state[2*j]), 2));
      if(dst == 0){
        throw range_error("Diveded by zero.");
      }else{
        energy += 1.0 / dst;
      }
    }
  }
  return energy;
}

double temperature(double r){
  double startTemp = 10.0;
  double alpha = 0.90;
  return startTemp * pow(alpha, r);
}

double probability(double e1, double e2, double temperature){
  if(e1 >= e2){
    return 1.0;
  }else{
    return exp((e1 - e2) / temperature);
  }
}

vector<double> simulated_annealing(int n, int maxIter){
  vector<double> startState(2 * n);
  while(1){
    for(int i=0;i<n;++i){
      startState[2*i] = 2.0 * M_PI * uniform();
      startState[2*i+1] = M_PI * uniform();
    }
    try{
      eval(startState);
      break;
    }catch(range_error& exception){
      continue;
    }
  }
  vector<double> state = startState;
  vector<double> bestState = state;
  double e = eval(state);
  double bestE = e;
  for(int i=0;i<maxIter;++i){
    double temp = temperature((double)i / maxIter);
    
    // change 1 variable to follow normal distribution
    int index = (int)(uniform() * 2 * n);
    double previous = state[index];
    if(index % 2 == 1){
      state[index] = fmod(rand_normal(state[index], sqrt(temp / 2.0)), M_PI);
    }else{
      state[index] = fmod(rand_normal(state[index], sqrt(temp / 2.0)), 2.0 * M_PI);
    }
    
    try{
      double nextE = eval(state);
      if(nextE < bestE){
        bestState = state;
        bestE = nextE;
      }
      if(uniform() <= probability(e, nextE, temp)){
        e = nextE;
      }else{
        state[index] = previous;
      }
    }catch(range_error& exception){
      continue;
    }
  }
  return bestState;
}

int main(){
  ofstream ofs("thomson_step40000.txt");
  for(int n = 1; n <= 50; ++n){
    vector<double> result = simulated_annealing(n, 40000);
    ofs << n << " " << eval(result) << endl;
    cout << n << ": " << eval(result) << endl;
  }
  return 0;
}

結果

N=50までを反復回数20000, 40000のそれぞれで探索しました。最小値は冒頭のwikiに載っています。
真値との相対誤差をプロットしたものが下図
f:id:lan496:20150913192841p:plain
f:id:lan496:20150913192854p:plain

まとめ

お手軽に近似解を出すことができました。焼きなまし法の練習にもよい題材だったと思います。

それと、Thomson問題を焼きなまし法で解く論文はあるようなのですが、有料だったので読めてません。
Generalized simulated annealing algorithm and its application to the Thomson model