Orz日記 by Akio Morita

_ [SAD][数学]MatrixLog[]再実装

C.K.Allenの実装した MatrixLog[]なのだが、精度が悪い＆収束性がよろしく ないので、直せないかゴソゴソ調査していたのだが...

結論、アルゴリズム的に数値不安定性が有って無理！

と言うわけで、まったく別のアルゴリズムで再実装しました。

MatrixLog[]の基礎

実数の場合と同様に Logarithmを Exponentialの逆関数で定義する。

行列$\mathbb{F}$の Logarithmである行列 $\mathbb{A}$は、 $MatrixExp[\mathbb{A}] = \mathbb{F}$を満足する行列。

MatrixExp[]の基礎

Exponentialは実数の場合の定義をそのまま拡張して定義する。 例えば、積展開の極限

$MatrixExp[\mathbb{A}] = \lim_{n\rightarrow\infty}\left(1 + \frac{\mathbb{A}}{n}\right)^n$

もしくは、Taylor級数展開

$MatrixExp[\mathbb{A}] = \sum_{n=0}^{\infty}\frac{\mathbb{A}^n}{n!}$

実際の計算では、任意の精度で計算を打ち切るために Taylor級数展開が 使われます。重要なのは、このTaylor級数展開においては 行列のノルム $||\mathbb{A}||$が有界ならば級数が絶対収束し、 級数展開の収束半径は無限大であることです。

複素数まで拡大して、一般化すると有界な複素正方行列$\mathbb{A}$の Exponentialは上記の Taylor展開で一意に定義できます。

MatrixLog[]の級数展開

一般に逆関数を解くのは骨なのでExponentialと同様に Logarithmを級数展開してみると、 Taylor級数展開では

$MatrixLog[\mathbb{I} - \mathbb{G}] = -\sum_{n=1}^{\infty}\frac{\mathbb{G}^n}{n}, \mathbb{G} = \mathbb{I} - \mathbb{F}$

となります。この例では収束半径は、$||\mathbb{G}|| < 1$です。

つまり、$\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の計算では、

1. 式 $\mathbb{F} = MatrixExp[\mathbb{A}]$を解く
2. 行列 $\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}$