分子構造の重ねあわせ (行列、最小二乗法など？)

質問

すいません、化学っぽい内容の質問なのですが、
アルゴリズムは数学になると思うのでこちらで質問させていただきます。
また、説明不足を補うために一部プログラムっぽい表記をさせていただきますのでご了承ください。

現在、分子構造を比較するためのプログラムを作成しています。
そのプログラムに構造の「重ね合わせ」機能を追加したいのですが
そのアルゴリズムに最小二乗法などが必要らしく、悩んでいます。
(詳細は長くなるので下に書きます)

分かる方がいましたらご教示お願いします。
よろしくお願いします。

■重ね合わせについて
分子構造のデータは、各原子がそれぞれX, Y, Z座標の3つの値を持っています。
例えば、200原子からなる分子のデータ構造は以下のようになります。

Atom { double x; double y; double z; }
Atom atom = new Atom[200];

重ね合わせる2つの構造は原子数は同じとします。なので、分子構造2つをa1, a2とし、
対応する点の誤差をとっていくと、誤差の平均を求める計算は以下のようになります。

Atom a1 = new Atom[200];
Atom a2 = new Atom[200];
double d = 0;
for (i = 0; i < a1.length; i++)
d += sqrt(pow((a1[i].x-a2[i].x),2) + pow((a1[i].y-a2[i].y),2) + pow((a1[i].z-a2[i].z),2));
d = d / a1.length;

このdの値が最小値をとるようにa2の原子の座標を変更することを「重ね合わせ」と定義します。

■重ね合わせのアルゴリズムについて
基準となる構造に対し、もう一つの構造を並進や回転によって合わせる感じになると思います。
今回質問する箇所はこのアルゴリズムになります。
この際に行列計算や最小二乗法を使うことになるかと思われます。

■質問者の考え
途中まで、しかも間違っている可能性がありますが私の考えたアルゴリズムを追記しておきます。
(i) a1, a2の重心を求める
double x1 = y1 = z1 = x2 = y2 = z2 = 0;
for (i = 0; i < a1.length; i++) {
x1 += a1[i].x; y1 += a1[i].y; z1 += a1[i].z;
x2 += a2[i].x; y2 += a2[i].y; z2 += a2[i].z;
}
x1 = x1 / a1.length; y1 = y1 / a1.length; z1 = z1 / a1.length;
x2 = x2 / a1.length; y2 = y2 / a1.length; z2 = z2 / a1.length;

a1の重心(x1, y1, z1)とa2の重心(x2, y2, z2)が求まります。

(ii) 重心の誤差を求め、誤差をa2の各原子の座標に反映させる
xd = (x2 - x1); yd = (y2 - y1); zd = (z2 - z1);
for (i = 0; i < a1.length; i++) {
a2[i].x -= xd;
a2[i].y -= yd;
a2[i].z -= zd;
}

この方法で誤差はかなり小さくなります。しかし、最小値にはなってないと思われます。
なので、この方法が正しいのかよく分かりません。

回答No.1

echoes_x86 · ベストアンサー

こんにちは．
ソースコードは良く見てませんが…
要するに合同であることが期待されるn点のデータ集合XとX'に関して，
Σ ||Xi' - (R*Xi+T) ||^2 -> min
なる回転と並進ベクタを計算したいと考えます．
ここで，XiとXi'はi番目の座標を意味する3次元ベクタであり，
総和演算は範囲(i=1,...,n)とします．

質問者様の方法では上記の変換(R,T)に関して，
並進量Tまでしか解決できません．
ここからさらに回転量Rを解決する必要があります．
基本的な戦略は相関行列を特異値分解してRを決定します．
Rの推定には質問者様の方法で重心を合わせたXとX'に関して
M = Σ Xi'*Xi^T
なる3*3相関行列を計算し，M -> U*D*V^T と特異値分解します(^Tは転置記号)．
最小二乗の意味で最適なRは
R = U*V^T
で計算できます．
ただし，回転の向きは適宜調整してください．

データ集合XとX'が相似の意味で同じである場合は
重心を会わせた後にスケール調整を行ってください．

お礼

回答ありがとうございます。
やはり行列計算が必要になるのですね。
ちょっと難しいですが調べてみます。
大変参考になりました。