以前の記事がSVD分解を用い、また点群も平均移動させている方法を使いました。
今回は絶対値の方の損失関数で、また定数項用の次元を追加する考え方でフィッティングしてみます。
理論は回帰の誤差、２乗するか絶対値をとるか 補足のおまけをご覧ください。
絶対値誤差最小化による平面回帰
まず必要なモジュールのインポート。
import random
import math
import numpy as np
from scipy.optimize import minimize
答えの平面パラメータ v0 = [a, b, c], d0 を三角関数を使って決めます。 phi0 と theta0 で平面の３次元中での角度を決めます。
phi0 = random.random()*np.pi
theta0 = random.random()*np.pi
v0 = np.array([math.sin(phi0)*math.cos(theta0), math.sin(phi0)*math.sin(theta0), math.cos(phi0)])
d0 = random.random()*10.0
このパラメータで平面上の点群を作ります。とりあえず２０００点あたり。
def create_sample(nv, d, N):
""" nv: np.array of (a, b, c), |nv| = 1 """
""" nv*r + d = 0 """
dxyz = np.random.random(3*N).reshape(3, N)*0.005
x = np.random.random(N)*2 - 1.0 + dxyz[0, :]
y = np.random.random(N)*2 - 1.0 + dxyz[1, :]
z = - (x*nv[0] + y*nv[1] + d)/nv[2] + dxyz[2, :]
return (x, y, z)
xs, ys, zs = create_sample(v0, d0, 2000)
点群から平面をフィットさせる関数は以下のようになります。
def find_plane(xs, ys, zs):
r = np.c_[xs, ys, zs, np.ones_like(xs)].T
p0 = np.ones(4)
def loss(p):
return np.sum(np.abs(np.dot(p, r)))
res = minimize(loss, p0, method="Nelder-Mead")
print(res)
nvd = res["x"]
nv = nvd[:3] # パラメータ４つのうち平面の向きに関する３つだけをとりだす
n = np.linalg.norm(nv) # 平面の向きのベクトルの大きさを１にしたい
print(nvd/n)
find_plane(xs, ys, zs)
最小化のアルゴリズムとしてNelder-Meadだとうまく行きました。実はなぜかよくわかっていません。やっぱ絶対値関数の導関数の性質のせいですかね？
最後の nvd/n で、平面の法線ベクトル $a, b, c$ の長さを１にしつつ、変位 $d$ を正しいスケールにすることができます。
真値は以下のようなものの場合、
0.49957542861731424 x + 0.007143327762107544 y + -0.866240938763754 z + 9.262962179997787 = 0
minimize関数の出力は print(res) で出しています。
final_simplex: (array([[ 1.21480601e-01, 1.76214791e-03, -2.10653182e-01,
2.25312000e+00],
[ 1.21477207e-01, 1.75002916e-03, -2.10653604e-01,
2.25312674e+00],
[ 1.21482217e-01, 1.76227612e-03, -2.10653652e-01,
2.25311870e+00],
[ 1.21489275e-01, 1.75442357e-03, -2.10653048e-01,
2.25311743e+00],
[ 1.21475212e-01, 1.74415270e-03, -2.10653912e-01,
2.25312772e+00]]), array([0.51681127, 0.516823 , 0.51683176, 0.51684715, 0.51684861]))
fun: 0.5168112711864259
message: 'Optimization terminated successfully.'
nfev: 298
nit: 170
status: 0
success: True
x: array([ 1.21480601e-01, 1.76214791e-03, -2.10653182e-01, 2.25312000e+00])
得られた答えは print(nvd/n)で出力しています。
[ 4.99554811e-01 7.24633774e-03 -8.66251973e-01 9.26532238e+00]
絶対値誤差でもフィットできていることが確認できました。
SVD分解
せっかくなので最小二乗法、もしくはその同値なものでの性能と比較してみます。
SVD分解による、平面の平行移動ではなく１の次元を追加する場合でのフィッティングは
def find_plane_SVD(xs, ys, zs):
r = np.c_[xs, ys, zs, np.ones_like(xs)]
u, s, v = np.linalg.svd(r)
nvd = v[-1, :]
nv = nvd[:3]
n = np.linalg.norm(nv)
return nvd/n
となります。平面の向きに関するところだけ大きさを１にしてやると $d$ の成分も調度良くスケールされるのは先程の時と同じですね。
そしてテストとしてサンプル生成時に
def create_sample(nv, d, N):
""" nv: np.array of (a, b, c), |nv| = 1 """
""" nv*r + d = 0 """
dxyz = np.random.random(3*N).reshape(3, N)*0.005
x = np.random.random(N)*2 - 1.0 + dxyz[0, :]
y = np.random.random(N)*2 - 1.0 + dxyz[1, :]
z = - (x*nv[0] + y*nv[1] + d)/nv[2] + dxyz[2, :]
# 外れ値 追加ここから
Nnoise = 20 # 外れ値の個数
noise_width = 10.0 # 外れ値が外れる距離、この範囲で一様分布
x = np.append(x, (np.random.rand(Nnoise) - 0.5)*noise_width)
y = np.append(y, (np.random.rand(Nnoise) - 0.5)*noise_width)
z = np.append(z, (np.random.rand(Nnoise) - 0.5)*noise_width)
# 追加ここまで
return (x, y, z)
と、一様分布した外れ値を追加してやります。大体のサンプルが0.005くらいのズレに対して、外れ値は最大で10近くずれているものになります。
これを用いると、SVD分解が誤差が大きくなっていくのに比べて絶対値誤差最小化だとあまり影響を受けないことが確かめられるかと思います。
例えば正解が
-0.0898563576016248 x + 0.6389776072060572 y + 0.7639590646675976 z + 8.472614388573925 = 0
という平面だった場合、絶対値誤差の場合
-0.08994579 0.63910422 0.76384262 8.46941575
が得られましたが、SVDの場合
-0.14249166 0.52318171 0.84022439 9.33522365
となりました。絶対値誤差は$a, b, c$ の誤差が大体 0.0001、$d$だと 0.003の違いですが、SVDでは$a, b, c, d$ の誤差が大体 0.05〜0.1 で、絶対値誤差最小化のほうがはるかに正確な値を得ることができました。
なお、外れ値がない場合はどちらも正確な値を出してくれます。
正解
0.033990801030594014 x + 0.09240041507798535 y + 0.9951415923067002 z + 0.9028615255295003 = 0
絶対値誤差
0.03386859 0.09227052 0.99515781 0.90037193
SVD
0.0339317 0.09233158 0.99515 0.90034869