【逆問題】脳磁逆問題をProny法をといてPythonで実装【理論編】
はじめに
初めましての人は初めまして。そうでない人はお久しぶりです、なるみやです。
自分は、東京大学計数学科というところに所属しているのですが、これはその計数工学科と物理工学科からなる応用物理学門(略称AP)で企画されたAPアドベントカレンダー2018/12/08分のものです。
Prony法とは、大雑把に言えば、逆問題の一種に対する解法です。得られている結果から要因を求めてやろうというやつの磁気バージョン。例えば、脳磁場のポテンシャルエネルギーが測定できてる時に、その脳磁の源となる磁気双極子の個数と、その場所、そして磁気モーメントを求めてやろうというときに用いられる一つの解法です。
これができたら、頭のどこに双極子の源(「ソース」と言います)があるかわかるため、人の頭の中の活動が、頭の中を侵さずに計測できてすごくハッピーというわけです！！！すごい！！！！
構成
この記事は以下のように進めて行きます。
- 準備編：磁場の「複素ポテンシャル」という概念の導入と、それを用いた双極子の表し方。
- 脳磁逆問題の設定
- 逆問題の解析(多重極子展開と双極子パラメタの導入)
- Prony法を用いた解法(変数変換による解法)
- 3の数値計算(N次方程式の解の求め方)
- 考察(笑)
今回の記事では、1-4まで書いて行きます。長くなりそうなので。5-6は、次の記事で書いて行きます。(記事が出来次第リンクを貼ります。)
そうなんです！！今回プログラムは一切出てきません！！！タイトルにPythonとかあるのに！！！！
つまるところ、この記事は理論編と実装編との２部構成となり、これは、その理論編です。タイトルにも書いてあるしね。そのまんま。
導入が長くなりましたね、さて初めて行きましょー。
1. 準備:磁場の「複素ポテンシャル」という概念の導入
磁場のポテンシャル
最初に断っておきますと、当方、電磁気が苦手なので、変なところがあったらすいません。反省します。
話は、まず、Maxwell方程式の第一式、ガウスの法則から出発します。
一様空間中で原点に点電荷$q$がある時は、ガウスの法則より、
$$
\boldsymbol{\nabla} \cdot \boldsymbol{E} = \frac{q\delta(\boldsymbol{r})}{\varepsilon_0}
$$
にしたがって電場が生じます。ただし、$\delta(*)$はディラックのデルタ関数です。この時、ポテンシャル$V$、つまり、$\boldsymbol{E} =- \boldsymbol{\nabla} V$を考えると、
$$
- \boldsymbol{\nabla}^2 V = \frac{q\delta(\boldsymbol{r})}{\varepsilon_0}
$$
となります。(これは、ポアソン方程式。)
このポアソン方程式の解は3次元だと
$$
V = \frac{q}{4\pi \varepsilon_0 r}
$$
2次元だと、
$$
V = - \frac{1}{2\pi \varepsilon_0} \log r
$$
と表されます。ここでは、この2次元平面でのポテンシャルについて話して行きます。
複素ポテンシャルについて
さて、先ほどのポテンシャル$V$を複素ポテンシャルに拡張して行きましょう。まず、点電荷の位置を$z=x+iy$とします。なお、ここから先、$ \varepsilon_0=1$と仮定して進めて行きます。いちいちこれを書いていくのも面倒だし。
このとき先ほど求めたポテンシャルは
$$
V = -\frac{1}{2\pi} \log |z|
$$
となります。このポテンシャルの実部として、複素数として拡張します。
この時に、虚部としてちょうどいいのってなんですかね。....そう、共役調和な関数(コーシーリーマンの関係式を満たす関数)だね！！そうすれば、拡張されたポテンシャルが正則なものになって都合がいい！！
というわけで、
$$
f(z) = u + iv
$$
で、係数はひとまず置いておくと$u(z) = \log |z|$とし、
$$
\frac{\partial u}{\partial x} = \frac{\partial v}{\partial x} , \,
\frac{\partial u}{\partial y} = - \frac{\partial u}{\partial x}
$$
となる$v(z)$を求めると、それは$v = \arg z$なので、求める複素ポテンシャルは
$$
f(z) = - \frac{q}{2\pi} ( \log |z| + \arg z ) = - \frac{q}{2\pi} \log z
$$
となります。
複素ポテンシャルの物理的な意味
複素ポテンシャルの物理的意味を考えると、
- $\log |z|$=一定...等電位線
- $\arg |z|$=一定...電気力線
と考えられます。すなわち、複素ポテンシャルは、等電位線と電気力線の２つを１で表してるものと言えます。
双極子の複素ポテンシャル
次に、原点に$+q$を, $z=-\epsilon \in \mathbb{C}$に$-q$の電荷をおくと、この時に生成される複素ポテンシャルは、重ね合わせの法則から、
$$ \begin{array}{cl} f(z) &= - \frac{q}{2\pi} \log z + \frac{q}{2\pi} \log (z+\epsilon) \\ &= - \frac{q}{2\pi} \log z + \frac{q}{2\pi} (\log z + \log (1 + \frac{\epsilon}{z}))\\ &= \frac{q}{2\pi} \log (1 + \frac{\epsilon}{z}) \\ &= \frac{q}{2\pi} \left( \frac{\epsilon}{z} - \left( \frac{\epsilon}{z} \right)^2 + ... \right) \end{array} $$
ただし、最後の等式は$\log (1+x)$のテイラー展開を用いました。これを$\epsilon \to 0$とする、すなわち、二つの単位電荷を無限に近づけると、
$$ f(z) \to \frac{m}{2\pi z} $$
となります。これが、双極子の複素ポテンシャルです。ただし、$ q \epsilon =m $とした。$m$は物理的にはモーメントを意味します。
同様にして、
逆向きの双極子を無限に近づけると四重極子$1/z^2$が生じ、
逆向きの四重極子を無限に近づけると八重双極子$1/z^4$が生じ...となっていきます。
一般に、2n重極の複素ポテンシャルは$\frac{1}{z^n}$と書けます。
原点が$z_0$の場所に移ったときは、$z \to z-z_0$とすればいいです。
2. 脳磁逆問題の設定
さて、準備が終わったところで問題設定です。これは、最初に言ったのと同じです。
逆問題:
観測場所$z$として、二次元領域(人の頭とかを考えてください)$\Omega$中の$z_1,...,z_N$の中にモーメント$m_1,...,m_N$の双極子がある。 このとき、$ \Omega $の境界線、$ \partial \Omega $ 上で、複素ポテンシャル $$ f(z) = \sum_{k=1}^{N} \frac{m_k}{2\pi (z-z_k)} $$ を観測した。このとき、$z_1,...,z_N, m_1,...,m_N,N$を推定せよ。
観測場所$z$として、二次元領域(人の頭とかを考えてください)$\Omega$中の$z_1,...,z_N$の中にモーメント$m_1,...,m_N$の双極子がある。 このとき、$ \Omega $の境界線、$ \partial \Omega $ 上で、複素ポテンシャル $$ f(z) = \sum_{k=1}^{N} \frac{m_k}{2\pi (z-z_k)} $$ を観測した。このとき、$z_1,...,z_N, m_1,...,m_N,N$を推定せよ。
というものです。
さて、この解き方の発想としては、いくつかあります。
今回はその中でもProny法をご紹介します。
まずは、この逆問題が解きやすくなるように少し工夫をして行きましょー。
3. 逆問題の解析(多重極子展開と双極子パラメタの導入)
本当に、$f(z)$だけで、$N$個のソースの位置とかモーメントとか求められるのお？？とお思いになるでしょう？
でもできるんです！(ジョン・カビラ風)
問題が難しいのなら、解きやすい方向に持っていってあげればいいんです！！
では、そのやり方を見て行きましょう。
ローラン展開と複素ポテンシャル
さて、この逆問題を解いていくために、まず、$f(z)$のローラン展開をします。なんでローラン展開を考えるかの理由は、すぐ後にわかります。
ローラン展開とは、簡単に言っちゃえば、テイラー展開の拡張版です。(ローラン展開についてわからない人は、この記事(EMANの物理学)をどうぞ。)
さて、$f(z)$を原点周りでローラン級数展開をすると、
$$
\begin{array}{cl}
f(z) &= \frac{1}{2\pi} \sum_{k=1}^{N} \frac{m_k}{z(1-z_k / z)} \\\
&= \frac{1}{2\pi} \sum_{k=1}^{N} \frac{m_k}{z} \left( 1 + \frac{z_k}{z} + (\frac{z_k}{z})^2 +\dots \right) \\\
&= \frac{1}{2\pi} \left( \frac{\sum_{k=1}^{N} m_k}{z} + \frac{\sum_{k=1}^{N} m_k z_k}{z^2} + \dots \right)
\end{array}
$$
2個目の等式には、$|z_k /z| <1$を使いました。$z$はセンサの観測位置で、$z_k$はソースの位置です。
さて、最後にできてた式を見ると、結局、右辺は、
となる。つまり、ローラン展開は、多重極展開であると言えます。
観測した$f(z)$から双極子パラメタを定める
ソースの位置$z_1,...,z_N$やモーメント$m_1,..,m_N$を直接求めるのは難しいので、$f(z)$を観測して双極子パラメタを決めて、そのパラメタから、求めて行きたいと思います。
まずは、$f(z)$を$\partial \Omega$上でそのまま積分してみることとします。今回、$\partial \Omega$を単位円と仮定します。
このとき、
$$
\int_{\partial \Omega}\frac{dz}{z} = \int_0^{2\pi} \frac{i \mathrm{e}^{i\theta}}{\mathrm{e}^{i\theta}}d\theta = \int_0^{2\pi} d\theta = 2\pi i \\\
\int_{\partial \Omega}\frac{dz}{z^2} = \int_0^{2\pi} \frac{i \mathrm{e}^{i\theta}}{\mathrm{e}^{i2\theta}}d\theta =\int_0^{2\pi} \frac{i }{\mathrm{e}^{i\theta}}d\theta = 0
$$
以下、同様にすると、
$$
\int_{\partial \Omega}\frac{dz}{z^2} = \int_{\partial \Omega}\frac{dz}{z^3} = \dots = 0
$$
となるので、結局、
$$
\begin{array}{cl}
\int_{ \partial \Omega} f(z) dz &= \frac{1}{2\pi} \int_{\partial \Omega} \left( \frac{\sum_{k=1}^{N} m_k}{z} + \frac{\sum_{k=1}^{N} m_k z_k}{z^2} + \dots \right) dz\\\
&= i \sum_{k=1}^{N} m_k
\end{array}
$$
となります。積分しただけで、モーメントの和の$i$倍が出ちゃいました！！！
それに、この計算過程をみると、要は、$1/z$の項が残り、それ以外は0になるっぽい！！！
というわけで、お次は、$f(z)$に$z$をかけたものを積分にしてみると...
$$
\begin{array}{cl}
\int_{ \partial \Omega} f(z) z dz
&= \frac{1}{2\pi} \int_{\partial \Omega} \left( \sum_{k=1}^{N} m_k z + \frac{\sum_{k=1}^{N} m_k z_k}{z} + \frac{\sum_{k=1}^{N} m_k z_k}{z^2} \dots \right) dz\\\
&= i \sum_{k=1}^{N} m_k z_k
\end{array}
$$
次は、モーメントにソースの位置をかけたものの和が出てきました！！！いい感じですね！！！
これを繰り返していくと、結局
$$
c_n := \frac{1}{i} \int_{ \partial \Omega} f(z) z^n dz = \sum_{k=1}^{N} m_k z_k^n
$$
が得られます。この$c_n$をパラメータとして、逆問題を解いていけそうですね。
3章のまとめ
$c_n$というパラメータのおかげで逆問題は以下のように書き換えられました！
逆問題の違う表現:
$$ c_n := \frac{1}{i} \int_{ \partial \Omega} f(z) z^n dz \hspace{0.5cm} (n=0,1,...,N) \tag{1} $$ が与えられたとき、 $$ \sum_{k=1}^{N} m_k z_k^n = c_n \tag{2} $$ を満たす、$z_k,m_k \, (k=0,1,..,N)$および$N$を求めろ。
$$ c_n := \frac{1}{i} \int_{ \partial \Omega} f(z) z^n dz \hspace{0.5cm} (n=0,1,...,N) \tag{1} $$ が与えられたとき、 $$ \sum_{k=1}^{N} m_k z_k^n = c_n \tag{2} $$ を満たす、$z_k,m_k \, (k=0,1,..,N)$および$N$を求めろ。
というわけで、次はこれを解くことを目指して行きます。
Prony法を用いた解法(変数変換による解法)
さて、では、書き換えられた問題を解いて行きましょー。解き方は色々あるみたいです。知ってるだけで4通りあります。(全て演習の授業で説明されましたものです。)
$N$個を厳密に推定する方法はありますが、計算コストが高いことなどがあり、$N=10$などと仮定して位置やモーメントを求めたのちに、モーメントが小さいものはカウントしないという手法を取るのが実用的みたいです。ここでも、その手法に乗っ取り、$N$は既知として進めて行きます。
まず単純に考えてみる
単純に書き換えられたものを考えてみると、それは以下の非線形連立方程式を解くこととなります。 例えば$N=2$のとき
$$
\left(
\begin{array}{ccc}
1 & 1 \\\
z_1 & z_2 \\\
z_1^2 & z_2^2 \\\
z_1^3 & z_2^3
\end{array}
\right) \left( \begin{array}{c}
m_1 \\\ m_2
\end{array} \right)
= \left( \begin{array}{c}
c_0 \\\ c_1 \\\ c_2 \\\ c_3
\end{array} \right)
$$
でも、この連立方程式をそのまま解くのは大変なので、ここは少し工夫を試みます。
視点を変えてみる(変数変換)
さて、ここからトリッキーな感じになってきます。
突然ですが、
$$g(x)= \prod_{i=0}^N (x- z_i)$$
とおくと、$g(x)=0$の根は$z_i$となるため、
$$
g(x)= x^N + a_{N-1} x^{N-1} + \dots + a_1 x + a_0
$$
としたときの係数$a_i$が定まれば、そのN次の方程式$g(x)=0$の根を求めることで$z_i$が求められることとなります。
これは、実質、$z_i$を$z_i$の基本対称式$a_i$に置き換えたことに相当します。(それは、解と係数の関係からわかることです。)
そして解法
次なるターゲットは$a_i$です。この$a_i$さえ求めることができれば、話はほぼ片付きます。Prony法も最後の段階です。
さて、$z_i$は$g(x)=0$の根なので、全ての$i=1,...,N$で、当然$g(z_i)=0$です。
よって、
$$
\sum_{i=1}^{N} g(z_i) m_i z_i^{k}= 0 \hspace{0.5cm} (k=0,1,...,N-1) \tag{*}
$$
です。左辺を考えると、
$$
\begin{array}{cl}
(\mbox{左辺}) &= \sum_{i=1}^{N} m_i ( z_i^{N+k} + a_{N-1} z_i^{N+k -1} + \dots + a_0 z_{i}^{k} ) \\\
&= \sum_{i=1}^{N} m_i z_i^{N+k} + a_{N -1 } \sum_{i=1}^{N} m_i z_i^{N+k -1 } + \dots + a_0 \sum_{i=1}^{N} m_i z_i^{k} \\\
&= c_{N+k} + a_{N-1}c_{N+k -1 } + \dots + a_0 c_{k}
\end{array}
$$
となります。よって式(*)は、
$$
c_{ N +k } + a_{ N-1} c_{ N+ k -1 } + \dots + a_0 c_k = 0 \hspace{0.5cm} (k=0,1,...,N-1)
$$
となります。これを$k$について連立させると...
$$
\left(
\begin{array}{ccccc}
c_0 & c_1 & c_2 & \dots & c_{N -1 } \\\
c_1 & c_2 & & & \vdots \\\
c_2 & & \ddots & & \vdots \\\
\vdots & & & \ddots & \vdots \\\
c_{N -1 } & \ldots & \ldots & \ldots & c_{2 N -2 }
\end{array}
\right) \left( \begin{array}{c}
a_0 \\\ a_1 \\\ \vdots \\\ a_{N -1 }
\end{array} \right)
= \left( \begin{array}{c}
c_N \\\ c_{N +1 } \\\ \vdots \\\ c_{2 N -1 }
\end{array} \right)
$$
という、線形連立方程式ができます。ということで、この連立方程式を解くことで$g(x)$の係数$a_i$が求まります。
ちなみに、ここで出てきた行列、
$$
\left(
\begin{array}{ccccc}
c_0 & c_1 & c_2 & \dots & c_{N -1 } \\\
c_1 & c_2 & & & \vdots \\\
c_2 & & \ddots & & \vdots \\\
\vdots & & & \ddots & \vdots \\\
c_{N -1 } & \ldots & \ldots & \ldots & c_{2 N -2 }
\end{array}
\right)
$$
はハンケル行列という名前がついており、$N \times N$の行列ですが、$2N -1$個のパラメータだけで行列が確定されます。
これによって、$g(x)$が求められたので、$g(x)=0$の方程式を解くことで$z_i$が求まります。
また、逆問題の元の式(2)より、
$$
\left(
\begin{array}{cccc}
1 & 1 & \dots & 1 \\\
z_1 & z_2 & & z_N\\\
\vdots & \vdots & & \ddots & \vdots \\\
z_1^{N -1 } & z_2^{N -1 } & \ldots & z_N^{N -1 }
\end{array}
\right) \left( \begin{array}{c}
m_1 \\\ m_2 \\\ \vdots \\\ m_{N}
\end{array} \right)
= \left( \begin{array}{c}
c_0 \\\ c_1 \\\ \vdots \\\ c_{N -1 }
\end{array} \right)
$$
の関係式がわかっていて、先ほど、$z_i$が求められたので、この線形連立方程式を解いてしまば$m_i$が求められるということになります。
以上がProny法の理論です。(多分)
最後に
というわけで、今回はここまでです。
お疲れ様でした！！！(僕も疲れた！)
この記事としては、珍しくデスマス調にした気がしますが、なれませんね。このブログ、あくまでも備忘録として使ってるので。
元々は、我らが学科のコース長で、逆問題の権威でいらっしゃる奈良先生が、授業でお話して宿題にもなった「Prony法を数値計算でとけ」というものに対して、「東大多浪」という冊子の執筆をしていたため、満足にできなかったことから、そのリベンジの意味で書いてます。あと、ネットで調べた時にあまりにも日本語の記事が少なかったため、残して置いたら、あとあとProny法を実装する人の力になれるかなと思ったのもあります。(ほんとか？)
なお、ここでは理論よりも、実装(数値計算)の法に重きを置きますので、正直、Prony法については大雑把に説明して終わります。どうしても、Prony法について気になる方は「数学セミナー」の5-7月の数理のクロスロードを読むことをお勧めします。
さて、今回の記事はどうだったでしょうか？僕ら計数工学科は普段こんな勉強しています！
...ごめんなさい、嘘つきました。全てが全てこんな内容やってるわけじゃないです。これは授業で取り扱った中でもトップレベルに面白いことでした。しかし、それでも、計数工学科の本質とする学問はこういうものにあるところにあると思います。
そもそも、計数工学科の「計数」とは、計測と数理を足し合わせてできた言葉です。言うなれば、計測数理工学科です。
だから、成り立ちは、計測対象に数理モデルを立ててエンジニアリングをしていく学問を扱うところ、それが計数工学科の期限です。
今では、それに加えてアクチュエータなど、計測だけでなく制御なども行なっていますが、それの本質も、働きかける対象を計測して数理モデルを立てることです。
僕も最近、全ての本質はセンシングにありと思うようになりました。というのも、僕の標榜とする、「人工現実感」は人に現実ではないものを現実と感じさせることがゴールだと思っています。ということは、人がどのようにして「現実」を感じ、その人が「現実感」を得る感覚をいじれるようになれば、 「人工現実感」は達成されると思うからです。すなわち、人のセンシング機構を知り、それを制御できれば現実は作れる！！！...はずです。ですので、この僕の考えと、計数工学科の理念はすごくマッチしていて、ここ最近は、特に、計数にきてよかったなあ...と思っています。いやぁ、過去の俺、グッジョブ！！！
というわけで、なんか、計数についての所感を語ってしまいました。隙あれば自分語り。恥ずかしいですね。
長さも1ブログ記事としては、かなり長くなってきたので、ここら辺で引き上げたいと思います！！
あくまでも、この記事、自分としては次の記事となる数値計算がメインとなるのですがそこまでちゃんと書けるかな？（）
それでは、ここで退散します！！ではでは、またどこかで！！！
