*注意)この「べじぱみゅの学習メモ」のカテゴリー記事は、ワタクシ自身がこれまでに勉強したいろいろな項目について、テキストにあんまり書いてない内容などを勝手に妄想したメモです。
ワタクシ自身の備忘録のために書いており「初学者にわかりやすく説明する」というものではございません。
導入なしに唐突に話が始まり、おそらく意味不明な文章かもしれません。
しかし、せっかく考えたことなので、記事の内容がもし誰か1人でもお役に立てれば幸いです。
なんで「2乗」ばっかり?
ワタクシ自身の備忘録のために書いており「初学者にわかりやすく説明する」というものではございません。
導入なしに唐突に話が始まり、おそらく意味不明な文章かもしれません。
しかし、せっかく考えたことなので、記事の内容がもし誰か1人でもお役に立てれば幸いです。
なんで「2乗」ばっかり?
以前の記事で、スペクトル解析手法として交互最小二乗法を紹介しました。まさしくベースは「最小二乗法」なのですが、こんな疑問を持ったことはないでしょうか?
なぜ「最小1乗法」や「最小3乗法」は流行ってないのか?
結論から言うと「2乗以外は難しいから」なのですが、今回はノリで、あの「交互最小二乗法」を「交互最小x乗法」に拡張してみましょう。
交互最小二乗法はその名のとおり、もとのデータとフィッティングパラメータとの偏差の2乗和を最小化する方法なのですが、ここではこれ

偏差のx乗和、を考えてみます。なおxは0以上の実数としますが、xが整数でない場合
「a^x」はaが正のときしか定義されません。よって実際には「偏差の絶対値のx乗和」とする必要があります。
絶対値にすると微分した際に場合分けが必要となり、この後の表記がちょっとだけ複雑になって面倒なので、まずは絶対値でないx乗和を考え、最後に絶対値の効果を考慮します。
(2乗和のときはこの差は問題になりません)
2乗和がx乗和になっても、やることは一緒です。パラメータでの偏微分をゼロとします。

これをマジメに計算すると

こうなります。文字だらけで厄介なのでそれぞれFcik, Fsjkと置きました。このFcik(全部でI×K個)及びFsjk(全部でJ×K個)をすべてゼロにすることを目指します。
なお、行列Cと行列Sを「同時に」いじることも可能ですが、式が面倒になるのと計算がちょっと不安定になりがちなので、もとの「交互最小二乗法」にならいCとSを交互に更新する(更新している間、相手は「定数」とみなす)というスキームを採用します。
ちなみにx=1のときは特異点で、Fcikはsjkのみ、Fsjkはcikのみに依存するようになります。そのため「交互最小x乗法」のスキームが使えません。それは気が向いたら別の機会にお話するとして、ここではx=1は除外します。
上記は、x=2のとき以外は非線形方程式になります。非線形連立方程式を解く方法はいろいろありますが、ここではメジャーな「ニュートン法」を用います。ニュートン法自体の解説は別の場所を見ていただくとして、ニュートン法では以下の、ゼロにすべき関数を全パラメータで偏微分したものを用います。
なおありがたいことに、この2式はいずれもi=jでのみ値を持ち、それ以外ではゼロとなります。この「2階偏微分」を用いて、以下の「ヤコビ行列」を準備します。

これを用いて、ニュートン法のループは

このように表せます。右辺はすべてn回目のループにおける値で、左辺のベクトルを更新していき、収束するまで繰り返します。
これでもいいのですが、先ほど申し上げたようにヤコビ行列の成分のほとんどはゼロになるため、実際の計算はもっと簡単で

このように、iごと、jごとの小さなヤコビ行列を用いて

このような計算を順番にしていけばよいのです。

このように、iごと、jごとの小さなヤコビ行列を用いて

このような計算を順番にしていけばよいのです。
式として打ってしまうと形状は一緒に見えますが、もとのヤコビ行列はIK行IK列、JK行JK列の大きな行列であり大がかりです。こっちのスモールヤコビ行列はK行K列なので計算の負担はだいぶ小さくなります。
これで計算は可能です!なお説明の都合上Cの更新ループとSの更新ループを並列に書いてますが、実際は「交互最小」なので並列処理ではなく
これで計算は可能です!なお説明の都合上Cの更新ループとSの更新ループを並列に書いてますが、実際は「交互最小」なので並列処理ではなく
Cを更新→更新したCを用いてSを更新→更新したSを用いてCを更新→・・・
とやります。
はっぴぃ理系らいふ、いぇい
ヽ(・ε・)人(・ε・)ノ キミモナカマニナロウゼ
コメント