３次元の近似直線

質問

こんにちは。２次元で実験データなどの点列から近似直線を求めるのは、最小二乗法の基本問題ですが、３次元の点群から直線の方程式(x-x0)/a=(y-y0)/b=(z-z0)/cを求めるにはどんなアルゴリズムを使いますか？スマートな方法があれば教えていただけたら幸いです。よろしくお願いします。

回答No.2

stomachman · ベストアンサー

３次元空間の曲面ではなく、直線に乗ると仰るのだから、
(1) x, y, zのどれかを与えて、残りの２つを推定する問題。
(2) <x[i],y[i],z[i]>と直線との距離d[i]の二乗和が最小になる直線を求める問題。
と分類すべきでしょう。

(1)の場合は、たとえばzを与えてx,yを求めたいのであれば、
・zからxを求める問題。
・zからyを求める問題。
の二つを別々に解けばおしまいです。
それぞれの解は（x=Az+B, yは任意）という平面と、（y=Cz+D, xは任意）という平面を定めますから、この二つの平面の交線が、求める直線ということですね。

(2)の場合はやっかいです。
[1]ちょっと手抜きしながらも、まともにやってみましょう。
(i) 直線をどう表すか。
ご質問の式を見ると、この直線はx軸、y軸、z軸のどれとも平行でも垂直でもないことが仮定されています。
ですから、zをパラメータとして
x=az+c
y=bz+d
と書いても良いでしょう。a,b,c,dが決められれば良い訳です。
(ii) 点<p,q,r>と直線との最短距離を求める。
直線上の任意の点<az+c,bz+d, z>と点<p,q,r>の距離をdとすると
d^2 = (az+c-p)^2+(bz+d-q)^2+(z-r)^2
= (az)^2+2az(c-p)+(c-p)^2+(bz)^2+2bz(d-q)+(d-q)^2+z^2-2rz+r^2
です。これが最小になるzを求めると、
0=∂(d^2)/∂z = 2(az+c-p)a+2(bz+d-q)b+2(z-r)
ゆえに
z=(ap+bq+r-ac-bd)/(a^2+b^2+1)
であって、このときの最短距離h(p,q,r)は
h(p,q,r)^2 = (ap+bq+r-ac-bd)^2/(a^2+b^2+1)+2(ap+bq+r-ac-bd)(ac-ap+bd-bq-r)/(a^2+b^2+1)+(c-p)^2+(d-q)^2+r^2
わあ、とんでもないですね。
(iii) じゃあ、直線を求めるには？
S=Σ(h(x[i],y[i],z[i]))^2　　（i=1,2,...,N)
を最小化するには
∂S/∂a = 0
∂S/∂b = 0
∂S/∂c = 0
∂S/∂d = 0
を解く必要があります。言い換えれば
∂(h(x[i],y[i],z[i]))/∂a
∂(h(x[i],y[i],z[i]))/∂b
∂(h(x[i],y[i],z[i]))/∂c
∂(h(x[i],y[i],z[i]))/∂d
を求めておいて
Σh(x[i],y[i],z[i]) (∂(h(x[i],y[i],z[i]))/∂a)=0
Σh(x[i],y[i],z[i]) (∂(h(x[i],y[i],z[i]))/∂b)=0
Σh(x[i],y[i],z[i]) (∂(h(x[i],y[i],z[i]))/∂c)=0
Σh(x[i],y[i],z[i]) (∂(h(x[i],y[i],z[i]))/∂d)=0
という連立方程式を解くことになります。
これがa,b,c,dについて非線形である（一次式でない）ことは言うまでもありません。一筋縄では行かず、反復計算で徐々に収束させていくしかありません。

[2]手抜き
もうすこし手抜きの方法を考えてみましょう。
この座標系を回転・平行移動した座標系をX-Y-Zとします。そして、求めたい直線がZ軸と一致するようにしたとします。回転と平行移動は行列を使って
Ｘ　＝　Ｒ　ｘ　＋　ｐ
Ｙ　　　　　ｙ　　　ｑ
Ｚ　　　　　ｚ　　　ｒ
と表せます。Ｒは３×３の行列で Ｒの転置をＲ'とすると ＲＲ' = Ｒ' Ｒ = 単位行列
となる行列です。各点<x[i],y[i],z[i]>をこの変換で<Ｘ[i],Ｙ[i],Ｚ[i]>に写したとすると、
直線、すなわちＺ軸との最短距離はＸ[i]^2 + Ｙ[i]^2ですから、他のどんな回転・平行移動の仕方に比べても
U=Σ(Ｘ[i]^2 + Ｙ[i]^2)　　（i=1,2,....,N)
が最小になっている筈で、しかも
S=U
です。
　さて、UはＺ[i]の値とは無関係ですからＺ[i]を求める必要はない。さらに座標系をＺ軸の周りで回転してもUは変化しません。従って、
Ｘ　＝　Ｒ　ｘ　＋　ｐ
Ｙ　　　　　ｙ　　　ｑ
　　　　　　ｚ
Ｒ ＝Ｐ（α）Ｑ（β）
Ｐ（α）＝ｃｏｓα　　０　　－ｓｉｎα
　　　　　　　０　　　１　　　　０
Ｑ（β）＝　１　　０　　　　　０
　　　　　　０　ｃｏｓβ　－ｓｉｎβ
　　　　　　０　ｓｉｎβ　　ｃｏｓβ
とすれば良いのです。展開すれば
Ｘ[i] = x[i]ｃｏｓα－y[i]ｓｉｎαｓｉｎβ－z[i]ｓｉｎαｃｏｓβ＋ｐ
Ｙ[i] = y[i]ｃｏｓβ－z[i]ｓｉｎβ＋ｑ
ですね。
ここでα、β、ｐ、ｑを決めたい訳です。

　始めに（１）の問題を解けば、α、β、ｐ、ｑの大体の値を求めることができます。これを使ってU(α,β,ｐ,ｑ)を計算します。
　それから、U(α,β,ｐ,ｑ)が小さくなるようにα、β、ｐ、ｑをちょっとずつ改良して行けば良いでしょう。これには微小な角度Δα、Δβを使って、
Ｐ（Δα）＝ｃｏｓΔα　　０　　－ｓｉｎΔα
　　　　　　　　０　　　　１　　　　０
　　　　　　ｓｉｎΔα　　０　　　ｃｏｓΔα
Ｑ（Δβ）＝　１　　０　　　　　　０
　　　　　　　０　ｃｏｓΔβ　－ｓｉｎΔβ
　　　　　　　０　ｓｉｎΔβ　　ｃｏｓΔβ
を作り、Ｐ（α）、Ｑ（α）にそれぞれ掛け算すれば良い。
Ｐ（α＋Δα）＝Ｐ（Δα）Ｐ（α）
Ｑ（β＋Δβ）＝Ｑ（Δβ）Ｑ（β）
だからです。さらにここで、Δα、Δβは微小だから、
ｃｏｓΔα≒１、ｃｏｓΔβ≒１、ｓｉｎΔα≒Δα、ｓｉｎΔβ≒Δβ
(Δα)^2≒０、(Δβ)^2≒０、ΔαΔβ≒０
という近似をしても構わないでしょう。
この近似を利用すると計算は一層簡単になり、Uを最小にするようにΔα、Δβ、ｐ、ｑを求める問題は線形最小二乗法（一次式の最小二乗法）になってしまい、簡単に解けます。
それを解いてから、真面目にＰ（α）、Ｑ（α）を計算しなおし、また線形最小二乗法を解く。これを収束するまで繰り返せば良いのです。

なお、stomachmanは計算間違いの常習犯ですから、チェックは慎重に。

お礼

ありがとうございました。非常に理解しやすく、活用することができました。また、内容も大変興味深いので、仕事の合間を見つけて非線形最小二乗法も勉強しようと思います。今後とも、よろしくお願いいたします。よいお年をお迎えください。

回答No.1

stomachman

求めた直線の式をどう使う積もりなのかによって、やり方は違います。

用途に適した「近似」がどういうものであるか、をまず検討しなくてはなりません。さらに、たとえば「近似」を最小二乗の意味であると決めたとしても、まだ話は決まりません。

2次元の場合でも、話は単純ではないんです。
(1)「データ<x[i],y[i]> (i=1,2,....,N)から
f(x) = a x + b
というモデルについて、
Σ(f(x[i])-y[i])^2
が最小になるようにa,bを決める。」
というのが典型的な最小二乗法の問題で、こうして求めた
f(x) = a x + b
は、測定値xからyの値を推定するのに用いる式です。

しかしこれを
(y-b)/a
と変形して測定値yからxを推定するのに使ったら誤りなんです。それがやりたければ
(2)「データ<x[i],y[i]> (i=1,2,....,N)から
g(y) = A y + B
というモデルについて、
Σ(g(x[i])-y[i])^2
が最小になるようにA,Bを決める。」
という問題を解かなくてはなりません。そうやって求めたA,Bは
A y + B ≠(y-b)/a
なんです。

(3) ２次元平面に点<x[i],y[i]> (i=1,2,....,N)が打ってある。これらのなるべく近くを通る直線を求めよ。
という問題ですと、これまた別の話になります。(1),(2)の場合にはx軸、y軸を何倍に引き延ばそうとa,b,A,Bの値に変化はありませんが、(3)ではそうは行かない。
x=αt+x0
y=βt+y0
という直線と、各点との距離を小さくしたい。ですから、
点<x[i],y[i]>とこの直線の距離をd[i]とすれば
d[i]^2 = (αt[i]+x0-x[i])^2+(βt[i]+y0-y[i])^2
ここにt[i] = (α(x[i]-x0)+β(y[i]-y0))/ (α^2+β^2)
です。よって、
Σ(d[i]^2)
を最小にするα,β,x0,y0（たとえばα≧0, (α^2+β^2)=1, x0 y0 = 0という条件を追加すると一意的に決まります)を求める問題です。

かくて、ご質問は３次元の場合ですけど、どう使うのかを明示して戴かないと答は出ません。

お礼

丁寧なご説明ありがとうございます。（１）と（２）、（３）混同して考えてました。確かにそれぞれの解は異なりますよね。そうすると、３次元の場合では、ｘ、ｙからｚを推定する場合と、同様にｘ、ｙを推定する場合、各点からなるべく近くを通る直線を求める場合に分類されるわけですね。今回の私の求めたいのは、最後の問題（３）です。よろしくお願いします。もしご面倒でなければ、最初の問題の場合も教えていただけると非常に嬉しいです。すみませんが、お世話になります。