はじめに

いままでpythonに触る機会がなかったものの、機械学習方面に強くテンソルの計算なんかが得意らしいということぐらいは知っている。
しかもアインシュタインの縮約記法に則った計算が出来ると聞いた。

アインシュタインの縮約記法とは、例えばn次元の行列$A$とベクトル$x$の掛け算$y = Ax$を考えた時に、
$$
y_{i} = \sum^n_j{A_{ij} x_{j}} = A_{ij} x_{j}
$$
というように同じ添え字($j$)が掛け算の中に現れたらn次元和算しようという規則だ。
これを用いることでとても式を簡潔に記載することができる。

これをプログラム中で書けたら、とても便利だ。

pythonのnumpyというライブラリにあるeinsumという機能を使うと

import numpy as np
A = np.arange(9).reshape(3,3)
x = np.array([1,2,3])
y = np.einsum("ij,j -> i",A,x)

print(y)
# [ 8 26 44]

と書けるらしい。

y = np.einsum("ij,j -> i",A,x)

この部分がそうである。

D言語でも是非欲しい機能である。
これをそのままD言語で再現してもいいのだが（すでに実装している方はいる）、~~見た目がちょっとアレである~~ せっかくならもっと数式っぽく自然に書きたい。

というわけでggedライブラリ内に実装した。この記事はその紹介だ。

使い方

行列とベクトルの積

pythonの例をと同じものを書いてみよう。

import std.stdio : writeln;
import ggeD;
void main()
{
	auto A = tensor([0.,1.,2,3,4,5,6,7,8],3,3);
	auto x = tensor([1.,2,3],3);

	auto y = Einsum | A.ji*x.j;
	foreach(ijk;y)
	{
		writeln(ijk ," | ", y[ijk]);
	}
}
// 0 | 8
// 1 | 26
// 2 | 44

となる。

auto y = Einsum | A.ji*x.j;

は
$y_{i} = A_{ij} x_{j}$

と対応しており似たような書き方になっている。が、添え字がi,jと逆だ。これは仕様である。

(どうしようか悩み中)

$A_{xy}$と書いたときに添え字のイメージと横縦が合わないのに納得いかないし、そもそも$A_{ijk}$と増えていったら、関係ないし...でも世にある行列の式と齟齬があるのは問題があるし...
万が一、自分のライブラリを使う人が居て今の書き方だと不都合だなんてリクエストが来たら考える。

行列のトレースや内積

auto trA = Einsum | A.ii;    // 12
auto dot = Einsum | x.i*x.i; // 14

足し算や関数適用も

Einsum | A.ij*x.i + x.j;    // [9 28 48]
Einsum | BroadCast!sin(x.i); // [0.841471 0.909297 0.14112]

クロネッカーのデルタやエディントンのイプシロン

クロネッカーのデルタやエディントンのイプシロンはこんな感じで書ける。

auto delta = tensor!double(3,3);
auto epcilon = tensor!double(3,3,3);
foreach(i,j;delta)
{
    delta[i,j] = i == j ? 1. : 0.;
}
foreach(i,j,k;epcilon)
{
    epcilon[i,j,k] = sgn((1.*j-i)*(1.*k-j)*(1.*k-i));
}
writeln(Einsum | delta.ij*A.ij);
writeln(Einsum | epcilon.ijk*A[0,0..$].i*A[1,0..$].j*A[2,0..$].k );

まとめ

アインシュタインの縮約記法は便利。