ガウス・ニュートン法の実装について

私は今，ガウス・ニュートン法の実装を行っており，測定データ(x,y)のデータ群から楕円形の式のパラメータを推定しようとしています．
方法として2乗の誤差関数を用いようと思っていますが，序盤で躓いてしまって，先に進めない現状です．

楕円形の基本形は
ax^2+bxy+cy^2+dx+ey+f=0 ・・・(1)
なのですが，ここから式(1)を式(2)の形に変形したいです．
Y=f(X,a,b,c,d,e,f) ・・・(2)

説明が足りない部分もあると思うので．その際は遠慮せずに聞いてきてほしいです．
よろしくお願いします

質問者からの補足コメント

• 

  最小二乗法を用いることでa~fまでのパラメータを算出まではできてはいるのです．しかし，データによって双曲線の式( a>0 && c<0など )になってしまう場合もあるためにどうすればいいかわからなくなり，この方法でやってみようかと考えた感じです．ちなみにデータは非線形のものです．

    補足日時：2019/04/10 10:12

No.4 ベストアンサー 回答者： stomachman 回答日時： 2019/04/11 22:48

No.3への追加です。

> データによって双曲線の式( a>0 && c<0など )になってしまう

　かなりノイズが大きくて、しかも楕円の少数の部分だけに偏ったデータなんでしょうかね。ならば、x, yそれぞれについて最大値と最小値が±1になるようにデータを一次変換してからやってみるのも良いかと思います。
　それでもダメなら、「楕円形とデータとのズレの尺度（評価関数）をきちんと定義した上で、その尺度を最小化する」という、最もマトモな（したがって非線形の）問題にガッツリ取り組むしかないでしょう。その時に適当なのがストレートなガウスニュートン法なのか、もうちょっと安定している（Marquadt法など）最急降下法とのミックスなのか、その辺はやってみないとわからんところがあります。この場合、双曲線を排除するという条件も加えられる。それには、尺度にペナルティ関数を含める（「a>0 && c<0など」になったら大罰金を加える）か、あるいは、実数の範囲では楕円にしかなりようがない表式を工夫して、そのパラメータを探索するかです。

No.3 回答者： stomachman 回答日時： 2019/04/10 12:34

No.2へのコメントについてです。

> 縦と横を正規化する（=原点中心に平行移動し，縦と横を最大値で割ることで半径1の円上のデータに変換）

回転も必要です。なので、平行移動した後で2×2の行列を掛ける、という演算で表せます。

No.2 回答者： stomachman 回答日時： 2019/04/10 08:25

データを (x[i], y[i])として

　　r[i] = p[1]x[i]^2 + p[2]x[i]y[i] + p[3]y[i]^2 + p[4]x[i] + p[5]y[i] + p[6]
を考えると、「データが何らかの楕円でよく近似できる」のであれば、r[i]は絶対値が小さい値になるでしょう。そして右辺は パラメータ p[1]〜p[6]の一次式になっています。つまり行列で書けば
　　r = J p
ただし
　　J[i,1] = x[i]^2, J[i,2] = x[i]y[i], J[i,3] = y[i]^2, J[i,4] = x[i], J[i,5] = y[i], J[i,6] = 1
です。Jの転置をJ’として
　　J’ r = J’ J p
を解けばpが決まる。これは線形問題ですから一発で解けて、ガウスニュートン法の出番はない。

　さて、この計算は 「Σ(r[i]^2) を最小化する最小二乗法」になっていて、「データ(x[i],y[i])と楕円との『近さ』」そのものの最小化になっているとは限りません。てか、「データ(x[i],y[i])と楕円との『近さ』」が具体的にきちんと定義されていないからこそ、そこに乗じて、「Σ(r[i]^2) を最小化するんでもいいでしょ、これなら簡単に計算できるしィ〜」とやっているわけです。

　で、実際、たとえばデータ(x[i],y[i])を平行移動した(x[i]-A,y[i]-B)を使うと、算出される楕円は別のものになり、平行移動を元に戻しただけでは、前の問題の結果と一致しません。が、その不一致は、もし「データが何らかの楕円でよく近似できる」のであればごくわずかです。
　じゃあどうするか。
　以上のやり方で楕円がひとつ決まれば、その楕円を単位円に写す座標変換T[1]が決まります。そこで、データを全部T[1]で変換する。そして上記の「Σ(r[i]^2) を最小化する最小二乗法」を使います。（このとき、変換されたデータはほぼ単位円上に乗っているわけですから、今度のr[i]は「(x[i],y[i])から原点までの距離と、単位円の半径との差」を表しています。）こうして変換T[2]を決定してやる。これを（実用上は1〜2度ぐらい）繰り返し適用すれば、T[k]はほとんど恒等変換になるでしょう。そしたら収束とする。
　何をやったことになっているかというと、「単位円を伸ばしたり拡大したり回転したり平行移動して楕円に写す座標変換(T*[1]T*[2}…T*[k])（ただしT*はTの逆変換）で、単位円で近似される点たちを変換したものが、与えられたデータ(x[i],y[i])になっている」ということによってデータ(x[i],y[i])と楕円との『近さ』」を定義したわけです。

　さて、このやり方で行くのなら、変換T[1]を出す前に簡単な前処理をしておくべきです。すなわち、x[i], y[i] それぞれの平均が0になるように平行移動する変換T[0]を適用しておきます。こうすると、変換T[1]以降に含まれる「平行移動の成分」が小さくなり、収束が速くなります。

この回答へのお礼

丁寧な回答ありがとうございます．
知識が足らなくて少し疑問があるのですが，「その楕円を単位円に写す座標変換T[1]」についてですが，縦と横を正規化する（=原点中心に平行移動し，縦と横を最大値で割ることで半径1の円上のデータに変換）という認識でいいですか？

お礼日時：2019/04/10 10:35

No.1 回答者： Tacosan 回答日時： 2019/04/10 00:21

要は

2次方程式を解きたい
ってこと?