つまり、$\mathbb{F}$の全ての固有値が複素平面上で1を中心とする単位円に入っている必要があります。つまり、どう頑張っても任意行列$\mathbb{A}$のExponential $MatrixExp[\mathbb{A}]$に対してLogarithmの値を定めるには Taylor展開では役不足と言うことです。

ただし、$\mathbb{F}$のすべての固有値が右半面にあれば、適当なスケール変換を施した行列$\lambda\mathbb{F}$は収束半径に入るので次の式で級数展開による定義が可能になります。

$MatrixLog[\mathbb{F}] = MatrixLog[\lambda\mathbb{F}] - Log[\lambda]\mathbb{I}$

したがって、一般の行列$\mathbb{F}$に対するLogarithmの計算では、

式 $\mathbb{F} = MatrixExp[\mathbb{A}]$を解く
行列 $\mathbb{F}$を変形して級数展開の収束半径に入れる

という方針になります。

旧MatrixLog[]の実装

(1)の方針に相当し、コード中ではquadratic modeなどと書かれているあたりが該当します。

行列$\mathbb{F}_s$のLogarithmの近似値$\mathbb{A}_n$を次の処方で改良しています。

$\mathbb{X}_n = MatrixExp[-\mathbb{A}_n].\mathbb{F}_s$

$MatrixExp[\mathbb{A}_{n+1}] = MatrixExp[\mathbb{A}_n].X_n$

もちろん、このままでは$MatrixExp[\mathbb{A}_{n+1}]$の逆関数が解けないので、箱の鍵は箱の中状態です。ここで、近似が十分に良い($\mathbb{x}_n$が単位行列$\mathbb{I}$に近い)と仮定すると、$\mathbb{X}_n$のLogarithmを Taylor級数展開で与えることが可能になり、

$\mathbb{W}_n = MatrixLogTaylor[\mathbb{X}_n]$

$MatrixExp[\mathbb{A}_{n+1}] = MatrixExp[\mathbb{A}_n].MatrixExp[\mathbb{W}_n]$

を得ます。ここで、行列$\mathbb{X}_n$が単位行列に十分近ければ、行列$\mathbb{W}_n$は零行列に近いので、 Baker-Campbell-Hausdorff級数展開(Zassenhaus式)を用いて

$\mathbb{A}_{n+1} = \mathbb{A}_n + \mathbb{W}_n + \frac{1}{2}\{\mathbb{A}_n, \mathbb{W}_n\} + \frac{1}{12}(\{\mathbb{A}_n, \{\mathbb{A}_n, \mathbb{W}_n\}\} + \{\{\mathbb{A}_n, \mathbb{W}_n\}, \mathbb{W}_n\}) - \frac{1}{24}\{\mathbb{A}_n, \{\mathbb{W}_n, \{\mathbb{A}_n, \mathbb{W}_n\}\}\} + ...$

が得られます。

で、実際のデータ(乱数生成した行列のExponential)を流して $MatrixLogTaylor[\mathbb{X}_n]$の入力を見ていると 2次収束していた$\mathbb{X}_n$が突然発散に転じてあっという間に無限大になってしまうとい現象に遭遇することになります。 (まあ、それ以前にオリジナルのコードにはバグが有るのだが、それとは独立な現象)

結論から言うと、この手法では解は2次収束なのだが、数値的に不安定なので演算誤差や打ち切り誤差を種に指数関数的に発散が発生している。つまり、数列 $\mathbb{A}_1,\mathbb{A}_2,\mathbb{A}_3,...$に指数関数的増大列が含まれることが諸悪の根源なので、直せません。Orz (発散の直前で止めることで、それなりに有意な解を得られるが、精度がそこで制限されてしまう)

新MatrixLog[]の実装

(2)の方針に相当し、$\mathbb{F}$を変形して単位行列に近づけます。

アイデア自身は、MatrixExp[]の実装でノルムの大きな行列を計算する際に使っている次の恒等式そのものです。(ExponentialのTaylor級数展開は絶対収束だが、ノルムの小さな行列ほど収束が早い)

$MatrixExp[\mathbb{A}] = MatrixExp[\frac{\mathbb{A}}{n}]^n$

実際のコードでは、$n$に$2^k$を使っています。

で、これに$\mathbb{F}=MatrixExp[\mathbb{A}]$を代入してLogarithmをとると

$\mathbb{F} = n MatrixLog[\mathbb{F}^{\frac{1}{n}}]$

が得られます。ここで、例えば$n = 2$を選ぶと

$\mathbb{F} = 2 MatrixLog[MatrixSqrt[\mathbb{F}]]$

となります。ここで注目すべきは、平方根の主値をうまく選ぶと $MatrixSqrt[\mathbb{F}]$の固有値分布では$\mathbb{F}$の固有値分布と比べて各固有値の偏角$\theta$が半分になる点です。つまり、複素平面全域にあった固有値が右半面に集めることが出来ます。ようは、平方根を複数回作用させて固有値分布を右半面実軸近傍に集めれば、Taylor展開を使ってLogarithmを計算できるわけです。

ここで問題となるのは平方根の実装で、ExponentialやLogarithmを使わずに求めなければ、意味が有りません。これを愚直に次のような Newton法で実装してみました。

$\mathbb{A} = MatrixSqrt[\mathbb{B}]$

$\mathbb{A}_{n+1}=\frac{1}{2}(\mathbb{A_n}+\mathbb{A_n}^{-1}\mathbb{B})$

$\mathbb{A} = \lim_{n\rightarrow\infty}\mathbb{A}_n$

これによる前処理で、だいぶマシな状態になるのですが、この漸化式自体が数値的不安定性を持っているという罠が有りました。

つまり、精度を上げようと反復を増やすと発散してしまう問題は本質的には改善していないと言うことです。

で、何か手段はないかと Wikipedia(英語版)などをあさっていると... Denman-Beavers square root iterationsという 数値的に安定な行列の開平法を見つけました。

Denman-Beavers square root iterations

$\mathbb{M}_{n+1} = \frac{1}{2}\left(\mathbb{I} + \frac{\mathbb{M}_n + \mathbb{M}_n^{-1}}{2}\right), \mathbb{M}_0 = \mathbb{A}$

$\mathbb{Y}_{n+1} = \frac{1}{2}\mathbb{Y}_n(\mathbb{I}+\mathbb{M}_k^{-1}), \mathbb{Y}_0 = \mathbb{A}$

$\mathbb{Z}_{n+1} = \frac{1}{2}\mathbb{Z}_n(\mathbb{I}+\mathbb{M}_n^{-1}), \mathbb{Z}_0 = \mathbb{I}$

この漸化式は数値的に安定に2次収束し、その極限は

$\lim_{n\rightarrow\infty} \mathbb{M}_n = \mathbb{I}$

$\lim_{n\rightarrow\infty} \mathbb{Y}_n = \mathbb{A}^{\frac{1}{2}}$

$\lim_{n\rightarrow\infty} \mathbb{Z}_n = \mathbb{A}^{-\frac{1}{2}}$

に収束します。

これの周辺論文を調べていたら、応用数学分野の2001年ごろの論文でそのものズバリな下記の論文を見つけました。

APPROXIMATING THE LOGARITHM OF A MATRIX TO SPECIFIED ACCURACY,

S.H.Cheng, N.J.Higham, C.S.Kenney, and A.J.Laub,

SIAM J. Matrix Anal. Appl. Vol.22 No.4 pp.1112-1125

で、現在開発版でのMatrixLog[]はこの論文に基づいて実装しています。

速度は旧実装の半分ぐらい
旧実装のような数値不安定性がない
- テスト行列(Random[N, 6, 6]で生成)に対しては、(MatrixLog[MatrixExp[#]]-#&)のノルムを$10^{-10}$程度まで抑えた状態で安定に計算できる。

最後に行う単位行列に近い行列のLogarithmも論文と同様に Taylor級数展開ではなく、Pade近似による有理級数展開で実装しています。 (2/2次のPade近似よりも旧実装の MatrixLog[]を突っ込んだ方が精度が悪かった罠)

現在、内部実装のチューニングと精度評価部の調整中 (論文のアルゴリズム自身は、任意精度で計算できるが実装は倍精度浮動小数点のマシンイプシロンで制約されている)

Pade近似

関数$f(x)$を次の形の有理和で近似する。

$r_{mm}(x) = \sum_{j=1}^{m}\frac{\alpha_j^{(m)}x}{1+\beta_j^{(m)}x}$

$\log{(1-x)}$の近似式で低次のものは

$r_{11}(x) = \frac{-x}{1-\frac{1}{2}x}$

$r_{22}(x) = \frac{-\frac{1}{2}x}{1-\frac{1}{2}\left(1+\frac{1}{\sqrt{3}}\right)x} + \frac{-\frac{1}{2}x}{1-\frac{1}{2}\left(1-\frac{1}{\sqrt{3}}\right)x}$

$r_{33}(x) = \frac{-\frac{4}{9}x}{1-\frac{1}{2}x} + \frac{-\frac{5}{18}x}{1-\frac{1}{2}\left(1 - \sqrt{\frac{3}{5}}\right)x} + \frac{-\frac{5}{18}x}{1-\frac{1}{2}\left(1 + \sqrt{\frac{3}{5}}\right)x}$

[ツッコミを入れる]