この記事は何か

これは僕の同人誌の内容をqiita向けに書き直したものです。
対象はある程度逆問題を知っている人向け。
内容的にはMNEの公式サイトを噛み砕いたやつ。僕にとっては難しかったからね。
https://mne.tools/stable/overview/implementation.html#the-minimum-norm-current-estimates

eLORETAという脳波研究用の計算方法を解説した所を一部改変の上、転載します。
脳の研究に応用されますが、それ以外の応用はよく知らないです。
あと、間違っていないという保証はないです。

TLDR;

連立方程式$Y=AX$が成立し、AとYが既知とする。
取り敢えず$A^\dagger$という逆問題が算出された場合、MAP推定となるが、
そもそも事前分布ないじゃん？だから、下記の式から再帰的に事前分布を更新して収束させる方法がある。
$\Sigma$を事前分布とした時に$\Sigma = \sqrt{AA^\dagger}$になると仮定して$\Sigma$を計算し、
その$\Sigma$から$A$を計算し、以下無限ループ

波に関する逆問題

Xが原因、Yが結果、Aが途中の計算式として、下記の式が成立するとします。
$$Y = AX$$
ここで、何が未知のパラメタであるかによって、使える手法が変わろうかと思います。

例えば、XとYが既知の場合、Aを推定するのはニューラルネットワークの「学習」でしょう。
カーネル法も使えますね。今回問題とするのは、YとAが既知の場合です。
さて、Aが既知であるというのには、例えば音波、地場、電流等、
一種の波の性質を持ったものが挙げられます。
これらは必ずどこかに伝わる時の伝わり方は物理の法則に基づくはずです。
即ち、Aを既知とすることが出来ます。

波に関する逆問題

$$Y=AX$$
この連立方程式を解きます。まぁ、普通の割り算というか、最小二乗法だとか、
最小ノルム解、ムーアペンローズの逆行列とか言われるやつです。

ここではそれらの内最小ノルム解に関する話になります。

横と縦を掛け算するので、行列の掛け算は連立方程式になるのは自明です。
これがどうしても連立方程式に見えない人は、行列代数の入門書でも読んで下さい。

解を$A^{\dagger}$として
$$X \Rightarrow A^{\dagger}Y$$
という感じにするのを目指します。

条件$Y=AX$のもとで$X$を小さくしたいのです。$X$が小さいってことは下記の式が小さいってことです。

$$||X||^2$$

ここで、$||X||$とはXの対角成分の事を指します。
Xが何らかのセンサーで捉えた波、例えばXの縦をチャンネル、横を時間軸とします。
${||X||}^2$がどうなるかと言うと、$||X^TX||$となります。何故これの対角行列かと言うと、
Xの1行目…例えばチャンネル1としましょう。これとXの2列目…これはチャンネル2です。
これらの掛け算はチャンネル間の相関を表しているのであって、振幅を表しません。
つまり、対角行列のみが大きさを表すのです。

まず、ラグランジュの未定乗数法を使います。$f(x,y)$が極値になるx,yは

$$L(x,y,\lambda)=f(x,y)-\lambda g(x,y)$$
の時に

$$\frac{\partial L}{\partial \lambda}= \frac{\partial L}{\partial y}= \frac{\partial L}{\partial x}=0$$
の解、または
$$\frac{\partial g}{\partial x}= \frac{\partial g}{\partial y}$$
の解。…という感じの公式ですね。
行列の微分をしないといけないので、行列の微分の仕方を確認しておきます。
$$\frac{\partial a^tx}{\partial x}=a$$
$$\frac{\partial x^ta}{\partial x}=a$$

今回はこれを微分します。
$$L=\frac{X^2}{2}-\lambda^T (AX-Y)$$

$$\frac{\partial L}{\partial X} = X - A^T\lambda$$
これが0になるので
$$X = A^T\lambda $$
$$Y = AX = A\lambda A^T = AA^T\lambda$$
$$X = A^T(AA^T)^{-1}Y$$

これで無事$X$を$A$と$Y$で表せました。
まぁ、ここまでは普通の最小二乗の最適解なのだけれど。

もっと良くする

もっと良くします。

$$Y = AX$$

という感じの割り算をときます。
が、今回はL2ノルムを使います。L2ノルムについてはググって下さい。
L2ノルムは要するに二乗をちっさくしようと言うことです。

$$argmin L_{x} = {||Y - AX||}^2 + \lambda {||WX||}^2$$

ここで、WXとは何かという疑問が生じます。
XはまぁXなんですが、WXは「これがちっさくなって欲しいな？」というものです。
WがIでも計算は出来るのですが、WがIの場合どうなるかと言うと
「なるべく全てが、平等に正則化されるように」です。
Wは重みとでも言えるもので、Xの振幅の逆数とでも言えるようなものです。
つまり、Xが小さくなってほしいと言うより、Xの大きい部分は大きく、
小さい部分は小さく予測したいのです。$W^TW$の逆数は分散になるような感じです。

ここで、$||X||$はXのトレースです。
$$ \hat L= tr((Y^T - X^TA^T)(Y - AX)) + \lambda tr(X^TW^TWX)$$
$$=tr(Y^TY - X^TA^TY - Y^TAX + X^TA^TAX) + \lambda tr(X^TW^TWX)$$
ここで、微分します。
$$\frac{\partial \hat L}{\partial X} = - 2A^TY + 2 A^TAX + 2 \lambda W^TWX= 0$$

$$-A^TY + A^TAX + \lambda W^TWX = 0$$
$$X = (A^TA + \lambda W^TW)^{-1}A^TY$$
さて、$W^TW$は分散の逆数なのですが、
$(W^TW)^{-1}$が長いので$\Sigma$って書き換えさせて下さい。
$$X = (A^TA + \lambda \Sigma ^{-1})^{-1}A^TY$$
さて、ここで以下の謎の式に持っていきます。
$$(A^TA + \lambda \Sigma)A^T = \Sigma A^T(A\Sigma A^T + \lambda I)^{-1}$$
上の式をじっと見て下さい。成立することが解ると思います。
故に、こうなります。

$$X = \Sigma A^T(A\Sigma A^T + \lambda I)^{-1}Y$$
だけど、まだ$\lambda I$が怪しいですね？ここで、はじめの式を見て下さい。
$$argmin L_{x} = {||Y - AX||}^2 + \lambda {||WX||}^2$$
これ、あからさまに第一項がYにまつわるもの、第二項がXにまつわるものですね？
ということは、XにWがあるならYにもそれに相当するものがあって良いはずです。
こんどはYで捉えた結果の分散共分散行列をCとしてあげます。
…ノイズと言いましたが「背景」と言ったほうがわかりやすいかも知れません。
数学得意な人は「重み付き最小ノルム解」と言います。
数学得意な人は一々強そうな名前言いますよね…。式はこんな感じ。
$$argmin L_{x} = C{||Y - AX||}^2 + \lambda {||WX||}^2$$

実際計算するとこうなるようです。
$$X = \Sigma A^T(A\Sigma A^T + C)^{-1}Y$$
で、この値は正則化されていないわけですから、
グニョングニョンの値になってしまいます。なので、結局正則化のパラメータ入れるのです。

$$X = \Sigma A^T(A\Sigma A^T + \lambda^2 C)^{-1}Y$$
ここで$\Sigma$の正体が何なのかと言うと、ベイズで言う所の事前分布です。
読者の皆さんはつかれましたか？僕も疲れました。

だけど、まだ$\Sigma$が未知です。まだ続きます(´・ω・｀)

白色化

Yの各成分、それぞれ独立でいて欲しいですよね？
こういう風に、それぞれの相関を打ち消す操作を白色化と言います。
まず、想定されるノイズに関して、分散共分散行列を作ります。

毎回観測される波は同じようなものであるという想定のもとで推定をしたいので、
この分散共分散行列を式に入れ込みます。

$$Y = AC^{1/2}C^{-1/2}X$$
$$\hat A^{\dagger}Y = C^{-1/2}X$$
$$\hat A = A^{\dagger}C^{1/2}$$
の元で計算できる$A^\dagger$を$\hat {A^\dagger}$とすると
$$ \hat{A^{\dagger}}C^{1/2}Y = X$$
こんなかんじ。
$$ \hat{A^{\dagger}}C^{1/2} $$
$$ = \Sigma A^TC^{1/2}(AC^{1/2}\Sigma C^{1/2}A^T + C)^{-1}C^{1/2}$$
$$ = \Sigma A^T(A\Sigma A^T + \lambda ^2 I)^{-1}Y$$
これが上記の最小ノルム解の白色化バージョンです。
せっかくここまでCを入れてたのに脱力感半端ないっすね。
まぁ、あとでしつこく出てきます。

これが、例えば脳の研究ではMNEとして実用化されています。

MAP推定

今回はラグランジュの未定乗数法で二乗したやり方を書きましたが、
ここまででも書きましたが、ベイズ統計学の視点からの見方もあります。
繰り返しになりますが、$\Sigma$は事前分布ですです。
しかし、ここまで書いて疲れました。これはここに書くのが超絶面倒いので書きません。

正規化

さて…式をよく眺めてみましょう。これは
$$Y=AX$$
の変形であり、シンプルな掛け算なのは言うまでもありません。
そこで、式を次のように書き直してみます。
$$X=A^\dagger Y$$

ここで、ふとした疑問が出てきます。
$A^\dagger$のノルムが1じゃない場合に奇妙なことが起こります。
$A^\dagger$が1じゃない場合で、何も対象がない時にセンサーで計測してみたと仮定してみて下さい。

センサーが捉えるノイズと何もない空間の分散は同じ大きさのはずですが、
$A^\dagger$が1じゃなかったら空室の中のノイズが
大きくなったり小さくなったりしておかしいですね？？？
何もない所に何かが生じてしまう。

ということで、計算結果を補正してあげる必要があります。
通常、補正の方法は分散の1/2乗で割り算してあげる事によって為します。
分散って言うと、振幅の二乗ですから、つまり振幅で割ってる的なことです。
それぞれのソースを、そのソースの振幅で割り算してあげるのです。

振幅で割るってどうするのかって？対角成分で割り算してあげるのです。
答えを分散で割るとき、その方法が有名所が2つあります。
一つの方法は下記です。

$$Y = AX$$
$$A^{\dagger}Y = X$$
この$A^\dagger$を補正してあげるのですが、ソースの振幅を出さねばなりません。
センサーの振幅は分散の1/2乗なので$C^{1/2}$とします。
そこで、ソースの対角成分を考えます。Yの誤差をyとしておきましょう。
CはYの分散、即ちyの二乗に比例します。そして、縦横ともにYのチャンネル数に等しいです。
$||A^\dagger y||$の大きさ的な物で割り算したいんですよね？
なら、チャンネル数じゃなくてyのソース数分の大きさの行列にするべきです。
$A^\dagger$は縦がチャンネル数、横がソース数の行列になるはずです。
その、縦横を考慮して、ソースを割り算できる行列を作るなら、
$||A^\dagger y (yA^\dagger)^T||$ と言う感じの形になるはずです。

つまり$\sqrt{||A^\dagger C(A^\dagger)^T||^2} $となります。
これで割るわけです。だから、こうなります。

$$X =||A^\dagger C(A^\dagger)^T||^{-1/2} A^\dagger Y $$
これは脳の研究ではdSPMと呼ばれる方法です。

もう一つの正規化

ソースの分散を別の方法で計算してみます。
$$\sqrt{||X^TX||} \Rightarrow \sqrt{||Y^TA\dagger X||} = \sqrt{||Y^TA^\dagger AY||}$$
つまり、推定された世界ではXの分散は$||A^\dagger A||$倍になっているということです。
じゃ、それで割ればいいということになりますね。

$$X = ||A^\dagger A||^{-1/2}A^\dagger Y $$
この方法は脳の研究ではsLORETAです。

本題・事前確率を事後確率で更新する

もうほんとに疲れてきましたが、ここが本題です。

上記の方法は$A^\dagger A$で割り算していました。
だけど、本当に$A^\dagger A$でいいん？という問題があります。そもそも、$\Sigma$自体が怪しさ満点です。
$\Sigma$はベイズで言いますと事前分布です。そして、Xの実際の分散が事後分布にあたります。
これををどうにかしてみます。
…事前分布、という言い方をするから分かりにくいですね。
要するに、$\Sigma$は重み行列です。そして、これから重み行列を補正します。
ところで、sLORETAを見て下さい。割り算してるってことは、つまり重み付けをしているんです。
そう、実は$\Sigma$と$\sqrt{||A^\dagger A||}^{-1}$は似たようなものなのです。
似たようなものなら同じにしちゃえばいいじゃん？ってなります。
つまり、下記の等式が成立して欲しいのです。

$$\sqrt{A^T(A\Sigma A^T + \lambda ^2 C)^{-1}A}^{ｰ1} = \Sigma$$
この$\Sigma$を収束させるために、繰り返し上記の式を再帰的に計算してきます。

まず、$\Sigma$を適当に決めます。一寸ズレてます。
次に、 $\sqrt{A^T(A\Sigma A^T + \lambda ^2 C)^{-1}A}^{-1}$を計算します。
2を新しい$\Sigma$として、以下無限ループ。

こうすると、ベイズの事後分布でベイズの事前分布を更新していることが解ると思います。
脳の研究ではeLORETAという方法です。
という感じの理解で多分合ってると思うけど、間違ってたらごめん(´・ω・｀)