ジョイジョイジョイ

ジョイジョイジョイジョイジョイ

Parikhの定理

この前 Parikh の定理(パリークの定理)を知ってびっくりしたので紹介します。

形式言語界隈では常識らしいです。

三行で説明して

文脈自由言語と正規言語は、単語を記号の度数で同一視(つまり記号の順番を無視する)と、同じクラスになるというものです。

具体的には準線形というクラスになります。

びっくりしませんか?僕はびっくりしました。

続きを読む

櫟井唯ちゃんの画像を無限に生成する話りぴーと

joisino.hatenablog.com ▲昔の記事

前回の試みから二年以上経ちましたがまだゆゆ式二期は発表されません。(*1)

二期が発表されないためこの記事のタイトルものんのんびよりさんから拝借することになりました。

やはり今話題のディープラーニングでなんとかするしかなさそうです。

三行で説明して

Progressive Growing GAN (Tero Karras et. al. 2017) を chainer で実装して櫟井唯さん(ゆゆ式)の画像を生成しました。

唯

こんな感じのができました。(下の方にスクロールするともっとたくさんあります。)

github.com

レポジトリです。

Progressive Growing GAN とは

浅いネットワークで小さい画像を生成・識別することからはじめ、段階的にネットワークと画像を大きくしていく手法です。太古に流行った AE の事前学習のような感じです。

こうすることで本来不安定な GAN の学習が安定します。

GAN の枠組みは基本的に WGAN-GP を使っています。

この論文では他にも GAN を安定させたり、多様性をもたせるためのいくつかの技術を提唱しているので紹介します。

バッチごとの統計

GAN の generator は教師データ 1 枚を覚えてそれだけ生成する過学習に陥りがちです。そこで、この論文では、 discriminator の中間層に、バッチ全体の統計量を用いることを提唱しています。

具体的には、

  1. その層の入力 (batch x channel x height x width 次元) について、各位置・チャネルごと (channel x height x width 箇所) に、バッチを通して(各 batch 個ある値についての)標準偏差を求める。
  2. それらの値 (channel x height x width 個ある) の平均 m を求める。
  3. m を height x width 個並べたものを(それぞれのバッチについて)入力に append して、出力とする(batch x (channel + 1) x height x width 次元)

というものです。こうすることで、データ 1 枚だけ覚えるのでは、分散が小さすぎて discriminator に見破られてしまい、結果として generator は多様な画像を生成するようになります。

論文では、この層を最終層付近に一箇所挿入しています。

論文によると、さらにリッチな統計量も試したが、このシンプルな方法が一番良かったとのことです。

論文では 1. で標準偏差を用いていましたが、後述の実験では、安定しなかったので分散を用いました。(chainer の F.sqrt は現状二回微分に対応していなくて、自分で書いたので何か間違えていたのかも)

初期化と動的補正

畳み込みの重みの初期化に He normal (平均 0 標準偏差 sqrt(2/N) の正規分布で初期化)を使うことが多いですが、この論文では標準正規分布で初期化し、畳み込むときに sqrt(2/N) 倍することを提唱しています。

こうすることで、様々な畳み込み層の学習が同じペースで進み安定するそうです。

ピクセルごとの正規化

generator の畳み込み層の出力の際に、毎回ピクセルごとに正規化をかけます。

つまり、出力は batch x channel x height x width 次元 あるわけですが、 batch x height x width 個の場所(各場所には channel 個値がある)について、標準偏差で割ってから出力します。

こうすることで値が発散するのを防ぐ効果があるそうです。

交互に学習

WGAN-GP では discriminator(critic) 5 回学習するごとに generator 1 回学習などとすることが一般的ですが、この論文では交互に学習 (n_critic = 1) させたそうです。

特に理由は書いてませんでした。

(これでもうまくいくから)簡単のため?

discriminator の出力にペナルティ

この論文では WGAN-GP の loss に加えて、 discriminator の出力値の二乗に比例する項を設けています。

これは、この論文では discriminator の最終層の活性化関数に恒等関数を使っていて、discriminator が強い時に出力がどんどん大きくなるのを防ぐためだと思います。

実験

実際に python 3.5.2 + chainer 3.0.0 で実装して実験しました。(レポジトリは上の方に貼ってあります)

データは唯の画像 4077 枚です。(同じ・ほぼ同じ画像も多く含んでいます)

4x4 の画像からはじめて 256x256 の画像を出力するようにしました(論文では 1024x1024 までやっています)

学習にかかった時間は p3.2xlarge (GPU は Tesla V100) を使っておおよそ丸 2 日です。

生成は以下に載せるような gif アニメを一本作るくらいなら c5.large (GPU なし)で 1 分くらいでできます。

結果をどうぞ。

ギャラリー

割といい感じに生成できてるのから、雑に二枚の画像補完しただけやんけみたいなのまで様々です。

唯

唯

唯

唯

唯

唯

唯

唯

唯

唯

唯

唯

唯

唯

考察

体感的に他の GAN の手法に比べてかなり安定しているように思います。

論文のデモビデオほどはうまくいかなかったですが、思ったよりいい感じに生成できました。

論文のデモビデオでは、見たら分かるとおり、目や鼻や口の位置などをかなり丁寧に正規化しています(論文の appendix に正規化する手法も書いてあります)

今回はそこまで丁寧にデータを洗ったわけではないので、この辺をもっと丁寧にやるともっといい感じになりそうです。

データももっと必要そう。

そもそも唯の画像がそれほど大量には無いので、一般のアニメキャラクターの生成モデル作って唯でファインチューニングしたりすると良いかもしれません。

もっといい感じにしてまたどこかで発表できたらいいな。

おまけ

f:id:joisino:20171106173720p:plain

▲データに少しゆずこが混ざってしまったせいでできた合体画像

最後に

何か間違いや質問などがあればお願いします。

最後まで読んでいただきありがとうございました。

*1 OVA は発表されました!めでたい!!

参考文献

[1] Progressive Growing of GANs for Improved Quality, Stability, and Variation | Research

[2] GitHub - dhgrs/chainer-WGAN-GP: A Chainer implementation of WGAN-GP.

[3] GitHub - joisino/chainer-PGGAN: Progressive Growing of GANs implemented with chainer

Convolutional LSTM

大学の実験で必要になって実装したのでメモしておきます。

Convolutional LSTM の説明

名前で完全にネタバレしてる感が否めないですが、Convolutional LSTM とは、LSTM の結合を全結合から畳み込みに変更したものです。

例えば画像を RNN に食わすときに、位置情報が失われないので便利です。

動画の次フレームの予測や天気予報などに使えます。

具体的には、以下のような構造になっています。

f:id:joisino:20171027192558j:plain

x は要素ごとの掛け算、 * は畳み込みを表します。

通常の LSTM との差分を赤で書きました。といっても、一部の掛け算が畳み込みになっているだけですが。

peephole の部分だけ要素ごとの掛け算になっていることに注意してください。

実装

Convolutinoal LSTM と、それをビルディングブロックとして使って画像予測のネットワークを実装しました。

github.com

レポジトリです。

python3.5.2 + chainer 3.0.0 で開発しました。

画像予測のネットワークは以下のような構造をしています。([1] より抜粋)

f:id:joisino:20171027151517p:plain

左側のネットワーク(Convolutional LSTM を 2 層並べたもの)に、入力となる画像(例えば 1 フレーム目, 2 フレーム目, 3 フレーム目)を流し込んでいって、データの表現を得ます。

次に、左側の LSTM たちの内部状態を右側の LSTM たちにそれぞれコピーします。

最後に、右側のネットワークにゼロを入力して、出てきた値を concat して、それらを 1x1 で畳み込んで出力を得ることを繰り返します(例えばこれが 4 フレーム目の予測画像, 5 フレーム目の予測画像, 6 フレーム目の予測画像になります)

右側のネットワークの二回目以降の入力は、前回の出力を入れるという実装もありそうです。論文を読んだ限り特に書いてなさそうだったので、今回はゼロを入力するようにしました。

実験

Moving MNIST [2] で実験しました。

計算資源の都合上論文の実装より小さめのネットワークです。

具体的には 2 -> 5x5 Conv -> 128 -> 5x5 Conv -> 64 -> 5x5 Conv -> 64 -> 1x1 Conv -> 2 (右側の forecasting network) という感じです。

以下は生成された画像たちです。

f:id:joisino:20171027145711p:plain

f:id:joisino:20171027143117p:plain

f:id:joisino:20171027143126p:plain

f:id:joisino:20171027143133p:plain

f:id:joisino:20171027143145p:plain

まあいいんじゃないでしょうか。

f:id:joisino:20171027145210p:plain

Loss です。

もう少しちゃんとチューニングする必要がありそうですかね。

だいたい 20 epoch = 8750 iteration = 140000 images の訓練を p2.xlarge (GPU は K80) 上で行って 22 時間くらいでした。

最後に

何か間違いや質問などがあればお願いします。

最後まで読んでいただきありがとうございました。

参考文献

[1] Xingjian Shi, et. al. Convolutional LSTM network: A machine learning approach for precipitation nowcasting. In NIPS, 2015

[2] Moving MNIST

[3] わかるLSTM ~ 最近の動向と共に - Qiita

神絵師がtwitterに上げた神絵を収集する

色んな人がやってそうだけど自分用に作ってみました。

リポジトリです。

github.com

使い方

  • リポジトリをクローンして初期設定します
  • 神絵師のリストを twitter で作って登録します
  • 動かします
  • slack に神絵が流れてきます

動作例

f:id:joisino:20171016205716p:plain

詳しく

動作の流れは以下の通りです。

  1. twitter API を使って、リストから画像を含むツイートを抽出して保存する
  2. 保存した画像から二次元イラストを抽出する
  3. 抽出した画像を slack incoming-webhook をつかってチャンネルに流す

2 が難しいところです。

神絵師が画像を投稿したからといってそれが必ずしも神絵であるとは限りません。

例えば今回の場合、神絵師が食べた焼肉の写真は slack に流れてきてほしくありません。

▲流れてきてほしくないツイートのイメージ(このアカウントは神絵師ではありません)

写真とイラストを区別するのは意外と難しかったです。(何か良いアルゴリズムがあったら教えてください)

ディープラーニングでなんとかしたりできるんでしょうか。

今回は OpenCV の Haar Cascade で nagadomi さんのアニメ顔検出器を使わせてもらいました。

github.com

これでアニメ顔が 1 箇所以上で検出されると、イラストという判定にしました。

体感 accuracy は 7 割くらいですが、 false negative はほとんど無く、slack に流れてくるのは神絵ばかりです。

今回の場合は、網羅的に収集したいというよりは、よさそうな画像をなんとなく眺めたい用途なので、これでよしとします。

今後の課題

ソシャゲガチャ問題

神絵師がソシャゲガチャの結果画像を上げることはよくあります。

今回の場合、ソシャゲガチャの結果画像は slack に流れてきてほしくないですが、アニメ顔が含まれるため防げません。

10 連ガチャの一覧結果はまだ例外でなんとかなりそうですが、特定の一枚の結果を上げている場合、その絵を描いたのも(別の)神絵師なので、弾くのはなかなかむずかしそうです。

最後に

何か間違いや質問などがあればお願いします。

最後まで読んでいただきありがとうございました。

ArtClass(IOI2013)をディープラーニングで解く

はじめに

IOI 2013 オーストラリア大会に Art Class という問題があります。

この問題は、画像データが与えられるのでその画像が

様式1(新造形主義の現代芸術) f:id:joisino:20171003140837p:plain

様式2(印象派の風景画) f:id:joisino:20171003140845p:plain

様式3(表現派のアクション・ペインティング) f:id:joisino:20171003140855p:plain

様式4(カラーフィールド・ペインティング) f:id:joisino:20171003140907p:plain

のいずれであるかを判定する問題です。

正答率が 0.9 以上になると満点が得られます。

IOI にしては珍しい機械学習的な問題であることと、ジャッジが壊れて結果が返ってこなくなったことなどで有名なので、知っている人も多いかもしれません。

問題文やデータは、 http://www.ioinformatics.org/locations/ioi13/contest/ から手に入ります。

普通の解法

例えば 3x365x65 の大きさの窓を作って分散を計算して、それらを使って手で決定木を作るなどすると解けます。

C++ での実装例

gist.github.com

正答率は 99/102 = 0.97... で満点が獲得できました。

ディープラーニングで解く

手で特徴量を作るのは面倒なのでディープラーニングでなんとかして欲しいですよね。

今回は VGG の pretrained model (with ImageNet) を使ってこの問題を解いてみました。(念の為に言っておくと、 IOI の本番では chainer も pretrained model も GPU も使うことはできません)

最終層だけ 4096 x 4 の FC 層に変えて学習しなおしました。

学習データが 4 * 9 = 36 枚しかないのでかなり厳しかったですが、augmentation(crop) してなんとか正答率 92/102 = 0.90... で満点を獲得できました。

適当に作った決定木より精度低い。悲しい。

ImageNet で前学習したモデルを使ったのが筋が悪かったかもしれません。絵画的なデータで前学習したモデルがあれば、もっと楽に満点を達成できたと思います。

テストデータに対しても少しチューニングしたので、汎化性能はさらにもう少し低いと思います。体感 accuracy は 0.8 ~ 0.85 くらいです。IOI 本番でもテストデータに対する正答率は見える(はずだった)のでよしとしましょう。

python + chainer での実装例

環境は python 3.5.2 + chainer 2.1.0 です。

gist.github.com

最後に

正直今回の場合、人力特徴量抽出とハイパーパラメータチューニングの面倒さは同じくらいでした。(本末転倒)

何か間違いや質問などがあればお願いします。

最後まで読んでいただきありがとうございました。

参考文献

画像や問題文は全て http://www.ioinformatics.org/locations/ioi13/contest/ からの引用です。

符号付き整数を二冪で割る

符号付き整数を 2 冪で割るときに右算術シフトを使ってしまって痛い目にあったことのある皆さんこんにちは。僕はあなたの仲間です。

コンピュータで割り算を行うのは一般に重いので、できるだけシフト等の軽い操作で済ませたいですよね。

最近のコンパイラは賢いので適当に書けば最適化してくれますが、何かしらの事情で自分で最適化しなければならいないときのためにやり方を紹介しておきます。

C99 で 32 ビット符号付き整数 x1<<n で割ることを考えます。

符号なし整数なら二冪の割り算 x/(1<<n) と右シフト x>>n は一致するのですが、符号あり整数については、割り算は 0 方向への切り捨てなのに対して、右算術シフトは切り下げとなります。

例えば (-3) / 2 の結果が -1 なのに大して (-3) >> 1 の結果は -2 です。

これは困りました。

簡単な対策として、x が負のときシフトを行う前に (1<<n)-1 を足すという方法があります。

こうすると、負の数のとき切り上げになり、割り算と同じになります。

ですが、if は一般に重いので、できるだけ使わずに済ませたいです。

そこで、全体として以下の手順を踏みます

  1. x31 だけ右算術シフトしたものを y とする。
  2. y32-n だけ右論理シフトしたものを z とする。
  3. xz を足す。
  4. xn 右論理シフトする。

1, 2 が重要です。これによって zx が正のとき 0 に、負のとき (1<<n)-1 になっています。

たとえば n=6 のとき

x = 0xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
y = 00000000000000000000000000000000
z = 00000000000000000000000000000000

x = 1xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
y = 11111111111111111111111111111111
z = 00000000000000000000000001111111

となっています。

正の数のときは最上位ビットが 0 なので、算術シフトによって 0 が引き出され全て 0 で埋まります。論理シフトをしても結果は変わりません。

対して負の数のときは最上位ビットが 1 なので、算術シフトによって 1 が引き出され全て 1 で埋まり、続いて論理シフトによって上位ビット 32-n ビットが全て 0 で埋まります。 1 は全部で n 個並んでいるのでこの値は (1<<n)-1 となります。

例えば以下のように実装できます。

gist.github.com

C99 では論理シフトになるか算術シフトになるかは処理系依存であることに注意してください。

上記のコードは gcc 5.4.0 で動作を確認しています。

おまけ

ちなみに最近のコンパイラでも、ここで紹介した感じになります。

以下は clang 3.8.0x / 64MIPSコンパイルした結果(の一部)です。

sra $2, $1, 31
srl $2, $2, 26
addu $1, $1, $2
sra $1, $1, 6

最後に

何か間違いや質問などがあればお願いします。

最後まで読んでいただきありがとうございました。

参考文献

[1] Henry S. Warren. (2012). Hacker’s Delight (2nd Edition)

シュトラッセンのアルゴリズムとその導出の仕方

 O(n^{3}) よりも小さい時間計算量  O(n^{\log_{2} 7}) \approx O(N^{2.807}) で行列積を計算するシュトラッセンアルゴリズムとその導出方法を紹介します。

シュトラッセンアルゴリズム

 n \times n の行列積  A B = C を考えます。

まず  A, B, C

$$ \begin{pmatrix} A_{11} & A_{12} \\ A_{21} & A_{22} \end{pmatrix} \begin{pmatrix} B_{11} & B_{12} \\ B_{21} & B_{22} \end{pmatrix} = \begin{pmatrix} C_{11} & C_{12} \\ C_{21} & C_{22} \end{pmatrix} $$

というように  \frac{n}{2} \times \frac{n}{2} 行列に分解します。

 n が奇数のときは余った行、列を適当にどちらかに分けてあげてください。

まず  7 つの行列積

$$ \begin{align} P_{1} &= ( A_{11} + A_{22} ) ( B_{11} + B_{22} ) \\ P_{2} &= ( A_{21} + A_{22} ) B_{11} \\ P_{3} &= A_{11} ( B_{12} - B_{22} ) \\ P_{4} &= A_{22} ( B_{21} - B_{11} ) \\ P_{5} &= ( A_{11} + A_{12} ) B_{22} \\ P_{6} &= ( A_{21} - A_{11} ) ( B_{11} + B_{12} ) \\ P_{7} &= ( A_{12} - A_{22} ) ( B_{21} + B_{22} ) \end{align} $$

を計算します。

すると、 C の各ブロックは

$$ \begin{align} C_{11} &= P_{1} + P_{4} - P_{5} + P_{7} \\ C_{12} &= P_{3} + P_{5} \\ C_{21} &= P_{2} + P_{4} \\ C_{22} &= P_{1} + P_{3} - P_{2} + P_{6} \end{align} $$

となります。

この証明は両辺を頑張って展開して一致することを確かめればできます。

ここで、各  P_{i} を計算するときに、再帰的にこのアルゴリズムを適用することを考えます。

各ステップの計算は再帰部分以外は行列の加減算なので、 O(n^{2}) でできます。

再帰は大きさが半分の問題を  7 回呼び出しています。

よって、 Master theorem より、全体の計算量は  O(n^{\log_{2} 7}) \approx O(N^{2.807}) です。

導出方法

 P_{i} C_{ij} の式が与えられれば証明は頑張って計算するだけですが、 P_{i} C_{ij} の式がどうやって出てきたのか分かりません。

ここでは簡単な(とはいっても骨の折れる)導出方法を一つ紹介します。

流れは以下のとおりです。

  1.  2 \times 2 の実行列の掛け算で、要素同士の掛け算を  7 回することを目指す
  2. 平面上の  120 度回転を表す行列  D をとる
  3. 非零ベクトル  u をとる
  4.  u と直交して  u^{\bot} D u = 1 を満たす横ベクトル  u^{\bot} をとる
  5.  X = u u^{\bot} とする
  6.  \{ D, X, D^{-1}XD, DXD^{-1} \} \{ D^{-1}, X, D^{-1}XD, DXD^{-1} \} 2 \times 2 の行列の基底になっている
  7.  A, B をそれぞれ上の基底で表現してかけると  7 回ですむ
  8.  D は実は固有多項式が同じなら同じ議論が成り立つので、簡単な行列をとってみる
  9.  u も簡単なベクトルをとってみる
  10. すると実は 7 はシュトラッセンアルゴリズムと同じだった

長いですね。

一つずつ見ていきましょう。

回転行列  D

 120 度 ( 2/3 \pi) の回転行列

$$ D = \begin{pmatrix} -\frac{1}{2} & -\frac{\sqrt{3}}{2} \\ \frac{\sqrt{3}}{2} & -\frac{1}{2} \end{pmatrix} $$

を取ります。

 tr(D) = -1,  det(D) = 1 なので、固有多項式 \lambda^{2} + \lambda + 1 です。

ケーリーハミルトンの定理より  D^{2} + D + E = 0 \cdots (1)が成り立ちます。

両辺に  D をかけると  D^{3} + D^{2} + D = 0 \cdots (2) が成り立ちます。

(1) と (2) の左辺どうしは等しいので  D^{3} = E \cdots (3) となります。これは、 D 120 度の回転行列であることを考えると当たり前のことです。

(3) より  D^{-1} = D^{2} \cdots (4) です。

よって、 (1) は  D^{-1} + D + E = 0 \cdots (5) と書き換えられます。

(5) より  tr(D^{-1}) = -tr(D) - tr(E) = 1 - 2 = -1 \cdots (6) となります。

非零ベクトル  u

非零ベクトル  u をとります。

 u と直交する非零の横ベクトル  u' をとります。

 D の固有多項式の実根は存在しないので、 u D固有ベクトルではありません。

よって、 Du u と別方向を向いているはずなので、  u'Du \neq 0 です。

 u^{\bot} = \frac{u'}{u'Du} とします。

 u^{\bot}u と直交し  u^{\bot} D u = 1 となります。

行列  X

 X = u u^{\bot} をとります。ここで、

 X^{2} = (u u^{\bot}) (u u^{\bot}) = u (u^{\bot} u) u^{\bot} = 0 \cdots(7)

 XDX = (u u^{\bot}) D (u u^{\bot}) = u (u^{\bot} D u) u^{\bot} = u u^{\bot} = X \cdots(8)

 XD^{-1}X = (u u^{\bot}) D^{-1} (u u^{\bot}) = u (u^{\bot} D^{-1} u) u^{\bot} = -u u^{\bot} = -X \cdots(9)

が成り立ちます。

一次独立性

前々節でも言ったように、 D の固有多項式の実根は存在しないので、 u D固有ベクトルではありません。

よって、 u Du は一次独立です。

同様に  u^{\bot} u^{\bot} D は一次独立です。

よって、これらの積  u \cdot u^{\bot}, u \cdot u^{\bot} D, Du \cdot u^{\bot}, Du \cdot u^{\bot} D 2 \times 2 行列の基底になります。

これらに両側から  D を掛けます。

 D^{-1} = D^{2} \cdots (4) を考えると、  DXD, DXD^{-1}, D^{-1}XD, D^{-1}XD^{-1} となります。

 D は正則なので、これらも  2 \times 2 行列の基底になります。

これらを足し合わせると  (D + D^{-1}) X (D + D^{-1}) = (-E) X (-E) = X となります。ひとつ目の等号は  D^{-1} + D + E = 0 \cdots (5) から従います。

よって、基底のうち一つを  X で置き換えても基底になります。

これより、  X, DXD^{-1}, D^{-1}XD は一次独立です。

基底表現

 X, DXD^{-1}, D^{-1}XD のトレースはそれぞれ、

 tr(X) = tr(u^{\bot}u) = tr(uu^{\bot}) = 0

 tr(DXD^{-1}) = tr(XD^{-1}D) = tr(X) = 0

 tr(D^{-1}XD) = tr(XDD^{-1}) = tr(X) = 0

となります。

よって、 \{ X, DXD^{-1}, D^{-1}XD \} から生成される  2 \times 2 行列の部分空間の要素のトレースは  0 です。

 tr(D) = -1, tr(D^{-1}) = -1 より、これらはそれぞれ  X, DXD^{-1}, D^{-1}XD と一次独立です。

よって、 \{ D, X, D^{-1}XD, DXD^{-1} \} \{ D^{-1}, X, D^{-1}XD, DXD^{-1} \} はそれぞれ  2 \times 2 行列の基底になっています。

$$ A = a_{1} D + a_{2} X + a_{3} D^{-1}XD + a_{4} DXD^{-1} $$

$$ B = b_{1} D^{-1} + b_{2} X + b_{3} D^{-1}XD + b_{4} DXD^{-1} $$

と表します。このとき基底同士の掛け算は

f:id:joisino:20170903160923p:plain

となります。ここでは (7), (8), (9) の結果をたくさん使っています。(表は [2] より引用)

同じ結果が同列、同行に現れているので、これらはまとめられます。

つまり、

$$ \begin{align} AB &=& a_{1} &\times& b_{1} &\times& id \\ &+& a_{2} &\times& ( b_{1} + b_{4}) &\times& XD^{-1} \\ &+& a_{3} &\times& ( b_{1} + b_{2}) &\times& D^{-1}X \\ &+& a_{4} &\times& ( b_{1} + b_{3}) &\times& DXD \\ &+& ( a_{1} - a_{4}) &\times& b_{2} &\times& DX \\ &+& ( a_{1} - a_{2}) &\times& b_{3} &\times& XD \\ &+& ( a_{1} - a_{3}) &\times& b_{4} &\times& D^{-1}XD^{-1} \cdots (10) \end{align} $$

となります。

もう少し簡単に

ここまで  D の性質としては  tr(D) = -1,  det(D) = 1 であることしか使っていません。よって、たとえば

$$ D = \begin{pmatrix} 0 & -1 \\ 1 & -1 \end{pmatrix} $$

ととっても構いません。このとき

$$ D^{-1} = \begin{pmatrix} -1 & 1 \\ -1 & 0 \end{pmatrix} $$

です。

また、 u にも任意性がありましたが、例えば  u = ( 1, 0 )^{{\rm T}} ととります。

このとき

$$ X = \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix} $$

となります。

このようにとることで上記の議論を  2 \times 2 行列の要素がブロック行列でも適用できることになります。( 1単位行列で、 0 は零行列です)

また、計算すると

$$ D^{-1}XD = \begin{pmatrix} -1 & 1 \\ -1 & 1 \end{pmatrix} $$

$$ DXD^{-1} = \begin{pmatrix} 0 & 0 \\ -1 & 0 \end{pmatrix} $$

となります。 A, B の要素ごとに考えると、

 X, DXD^{-1}, D^{-1}XD はトレースが  0 なので、  a_{1} = - tr(A) = - a_{11} - a_{22} となります。

また、 X, D^{-1}XD (2,2) 要素が  0 なので、 a_{3} = a_{22} + a_{1} = - a_{11} となります。

残りは、それぞれ異なる場所に  1, -1 が立っているだけなので、それぞれ

 a_{2} = a_{12} + a_{1} - a_{3} = a_{12} - a_{22}

 a_{4} = -( a_{21} - a_{1} + a_{3}) = - a_{21} - a_{22}

となります。

同様に

 b_{1} = -b_{11} - b_{22}

 b_{2} = b_{11} + b_{12}

 b_{3} = b_{22}

 b_{4} = b_{11} - b_{21}

となります。

また、瑣末な計算なのでここでは省きますが、 id, XD^{-1}, D^{-1}X, DXD, DX, XD, D^{-1}XD^{-1} もそれぞれ  1, 0, -1 からなる簡単な  2 \times 2 行列になっています。

最後に (10) で表したようにそれぞれ係数を掛けて足し合わせると、シュトラッセンアルゴリズムが出てくるというわけです。

最後に

何か間違いや質問などがあればお願いします。

もっと簡単な導出があれば教えていただけると嬉しいです。

最後まで読んでいただきありがとうございました。

参考文献

[1] Volker Strassen. (1969) Gaussian elimination is not optimal.

[2] Christian Ikenmeyer. Vladimir Lysikov. (2017) Strassen’s 2x2 matrix multiplication algorithm: A conceptual perspective.

[3] Master theorem - Wikipedia