行列演算を普通の数式のままコーディングしたい

数値計算処理をするとき、やたらにforループのネストが増えたりする。いかにも計算している気になるかもしれないが、元の数式が複雑になってくると、それをコードに落とすのがどんどん大変になってくる。

そうなると「数式をそのままプログラム中に書きたい」となってくる。

完全にそのままは使える記号の種類とかで無理があるが、数式にできるだけ近い形でプログラムしたくなる。そうすると、プログラムも短くなり、バグの入る余地が減る。それ以上に、プログラムが圧倒的に見やすくなるのだ。

ほぼ数式どおりに書ける数値計算用のプログラミング言語がいつくかあるようだが、前回紹介したOctaveもその１つである。
人工知能で機械学習などをやろうとすると、大量のデータを行列にぶち込んで、行列演算をごちゃごちゃやることが非常に多くなる。
ということで、今回は、Octaveの行列関連の演算を紹介しようと思う。

行列は、[]の中に、行単位で数字を空白またはカンマ(,)で区切りながら並べ、行をセミコロン(;)で区切る。 転置行列は、肩にTの代わりにアポストロフィ(‘)を付けるだけでOK。

>> A = [1 2 1; 2 1 0; 1 3 2]
A =

1   2   1
2   1   0
1   3   2

>> B = A'
B =

1   2   1
2   1   3
1   0   2

以下にAの逆行列を３つの方法で計算してみた。行列の-1乗がちゃんと逆行列になってくれる。
eye(3)は、サイズ3の単位行列で、これを行列Aで割っても逆行列が求まる。
 

>> inv(A)
ans =

-2.00000   1.00000   1.00000
4.00000  -1.00000  -2.00000
-5.00000   1.00000   3.00000

>> A ^ -1
ans =

-2.00000   1.00000   1.00000
4.00000  -1.00000  -2.00000
-5.00000   1.00000   3.00000

>> eye(3) / A
ans =

-2.00000   1.00000   1.00000
4.00000  -1.00000  -2.00000
-5.00000   1.00000   3.00000

本当に逆行列になっているか、行列Aとの積を計算してみよう。
 

>> (A ^ -1) * A
ans =

1.00000   0.00000   0.00000
0.00000   1.00000   0.00000
0.00000   0.00000   1.00000

いきなり逆行列を求めてしまったが、ランク、行列式を求めてみよう。 さらに、固有ベクトル、固有値を求めてみる。
 

>> rank(A)
ans =  3
>> det(A)
ans = -1.0000
>> [V,L] =eig(A)
V =

0.50393   0.47733   0.10741
0.34334  -0.64817  -0.46439
0.79257   0.59332   0.87909

L =

Diagonal Matrix

3.93543         0         0
0  -0.47283         0
0         0   0.53740

こんな感じで、行列演算が簡単にできる。
ここでは、説明のために、１つ１つ示したが、これらの演算を連続してやっている数式は、Octaveなら１つの数式として書ける筈だ。

Octaveは、数値計算が便利にできるだけではなく、計算結果を３次元表示することもできる便利な道具である。なので、今後数回に渡って、色々紹介しようと思っている。