RangeとExpression Templateと行列ライブラリ

D言語といえば、コンパイル時処理がものすごい。
様々な処理をコンパイル時に終わらせることができ、逆にコンパイル時に動かないものは実装が悪いというD言語erまでいる。
[要出典]

今日は、そんなコンパイル時処理、というよりもtemplateを活用して私が作成したGigueという行列ライブラリを解説していきたい。

Gigueとは

Gigueは、D言語の行列ライブラリである。
特徴として、uBLASやEigenといったC++の行列ライブラリでも採用されている式テンプレート(Expression Template)というもので実装されている。
また、インターフェイスはそれらのライブラリよりも圧倒的に簡素であり、テンプレートがある程度扱えるD言語ユーザーであれば簡単に改造が可能である。

import std.stdio;
import std.random;

auto m33 = matrix!(3, 3, (i, j) => uniform(0, i + j + 1));           // 3 行 3 列の行列
auto v3 = (ones!int * 4.5).matrix!(3, 1);

v3 = m33 * v3;                              // 行列とベクトルの積

// v3 = v3.transpose * m33;                 // コンパイル時エラー
                                            // 左辺(3, 1), 右辺(1, 3)

v3 = (v3.transpose * m33).transpose;        // OK

writeln(v3);

ETはRangeから学べ！！

Rangeは素晴らしい。
一時オブジェクトを生成せず、最適化次第では最高のパフォーマンスを発揮する。
例として、std.algorithm.mapを単純化したRangeを以下に示す。
この記事はRangeを解説するためのものではないので、Input rangeにのみ対応している。

auto map(alias f, R)(R r)
if(isInputRange!R)
{
    static struct Map
    {
        auto ref front() @property { return f(_r.front); }
        bool empty() @property { return _r.empty; }
        void popFront() { _r.popFront(); }

      private:
        R _r;
    }

    return Map(r);
}

この実装がまさにExpression Templateである。
Rangeの特徴である「遅延評価」はExpression Templateの特性である。

Rangeには、大きく2つに分けるとInput rangeとOutput rangeがある。
さらにInput rangeには、Forward range, Bidirectional range, Random access range, Infinite rangeといったような区別がある。
これらは、インターフェイスによってのみ定義され、あるinterfaceやclassを継承しているわけではない。
Input rangeになるためには、rangeの先頭要素を返すr.front, rangeを進めるr.popFront(), rangeが空かどうか調べるr.emptyの3つのインターフェイスが最低でも必要である。

Gigueでは同様の仕組みを行列とベクトルに対して適用した。
行列の場合に重要となることは、「行列の大きさ」と「行列の要素へのアクセス」である。
つまり、行列の大きさを取得できるだけのなにかと、行列のそれぞれの要素へアクセスできるだけの何かが必要である。

古くから様々な行列ライブラリでは、rowsやcolsといった名前で前者が実装され、後者は演算子オーバーロードなどを用いて解決されてきた。
残念ながら、私には古くからのしきたりを破壊できるような力はない。
Gigueでもm.rowsやm.colsで行列の大きさにアクセスでき、m[i, j]によって行列の要素へアクセス可能とした。

/**
行列型かどうかを判定する
*/
template isMatrix(T)
{
    enum bool isMatrix = is(typeof((const T m){
            size_t rsize = m.rows;
            size_t csize = m.cols;

            auto e = m[0, 0];
        }));
}

行列とベクトル

数学では、1階のテンソルをベクトルといい、2階のテンソルを行列というようだ。
私は、数学についてはあまりわからないので、テンソルとは「スカラー, ベクトル, 行列を一般化した何か」という、多分間違った認識をしている。
なので、テンソルについては考えず、単純に行列とベクトルという2つについて考える。

ベクトルと行列を別々に考え、それらについて個々に実装することは、大変な労力の無駄である。
つまり、行ベクトルもしくは列ベクトルと考えた方が、実装は簡素になり、発生するバグも少なく、無駄な努力がない。

しかし、ベクトルはベクトルで扱うことも重要である。
ベクトルの要素にアクセスするためにv[i, 0]と書かなければいけないのは不便である。

Gigueでは、Input rangeに対してForward rangeなどがあるように、行列に対していくつかのインターフェイスを追加したものをベクトルとした。
まず、ベクトルはv[i]によって要素にアクセス出来る必要があり、v.lengthによってベクトルの長さを取得可能としなければならない。

/**
ベクトル型かどうか判定する
*/
template isVector(V)
{
    enum isVector = isMatrix!V && !isInferableMatrix!V
                 && is(typeof((const V v){
                        static if(hasStaticRows!V)
                        {
                            static if(hasStaticColumns!V)
                                static assert(V.rows == 1 || V.cols == 1);
                            else
                                static assert(V.rows == 1);
                        }
                        else
                            static assert(V.cols == 1);

                        size_t size = v.length;

                        auto a = v[0, 0];
                        auto b = v[0];
                        static assert(is(typeof(a) == typeof(b)));
                    }));
}

このインターフェイスは、Random access-ibleでhas-lengthなrangeと親和性が高いように思えるかもしれない。
しかし、ベクトルは行列として扱える必要があり、v[i, j]のコンパイルが通る必要がある。
つまり、Gigueではrangeをベクトルとして扱えない。
しかし、制約としてisRandomAccessRange!R && hasLength!Rを課すことによってrangeをベクトルに変換する関数を提供している。
注意として、以下の実装は、実際よりも少しだけ簡略化されている。

/**
レンジからベクトルを作る
*/
auto toVector(Major mjr = Major.row, R)(R range)
if(isRandomAccessRange!R && hasLength!R)
{
    alias E = Unqual!(std.range.ElementType!R);

    static struct ToVector
    {
        enum rows = 1;

        auto cols() const @property
        {
            return _range.length;
        }


        auto ref opIndex(size_t i, size_t j) inout
        in{
            assert(i == 0);
            assert(j < this.cols);
        }
        body{
            return this._range[j];
        }


        mixin(defaultExprOps!(false));

      private:
        R _range;
    }

    return ToVector(range);
}

Gigueの素晴らしい点は、たとえベクトルであろうが行列としての最低限のインターフェイスであるrows, cols, opIndex(size_t, size_t)さえ実装しておけば、mixin(defaultExprOps!(false));によってベクトル用インターフェイスや演算子オーバーロードが自動的に生成されることである。
もちろん、mixin(defaultExprOps!(false));を捨てて、自ら全ての演算子と残りのベクトルとしてのlengthを提供してもよいが、それは必要に迫られた場合にのみに限定すべきである。

行列とレンジ

Gigueには、rangeとの親和性はほとんど無いことは先に述べた。
この理由は、Rangeとベクトルでは根本的に意味が異なり、またRangeをベクトルとして扱ったとしてもRandom access rangeしか上手くラップできないからである。

Rangeを行列として扱うことは、ほぼ不可能であろう。
行列はN×M行列というように、行の大きさと列の大きさが決まっている。
RangeのRange、つまり要素がRangeであるRangeは確かに行列のようであるが、要素となっているRange全てが同一の長さを持つかどうかは保証できない。

以上に加え、行列に対する通常の基本操作は、そこまで複雑ではない。
std.algorithm.mapのようなものは必要かもしれないが、std.algorithm.filterのようなものは不必要である。

Gigueの使い方―宣言―

今までの解説によって、GigueはExpression Templateによって実装された行列ライブラリであることは理解できたと思う。
Gigueで3x3行列を扱い場合には、次のように宣言する。

alias M33f = SMatrix!(float, 3, 3);

M33f m;

SMatrixは、固定長な行列をスタック上に構築する。
つまり、SMatrixはただ単に固定長配列をラップしているだけだ。

スタックに確保されるのが嫌、もしくはスタックでは容量が足りない場合には、newで確保するか、もしくはmatrix関数を使用することでDMatrixを作れば良い。

alias M33f = SMatrix!(float, 3, 3);

// ヒープ上に作成
auto m = new M33f;

// ヒープ上に作成
auto m = matrix!(float, 3, 3)();

今までの3つの例に出現した宣言方法では、行列の大きさは静的に確定しており、実行時に変更できない。
しかし、次のように行列を作成することで、動的な大きさを持つ行列を作成可能である。

// 行の数が動的
auto mdr = matrix!(float, dynamic, 3)(3);
assert(mdr.rows == 3);
assert(mdr.cols == 3);

mdr.rows = 12;
assert(mdr.rows == 12);

同様に列についても、もちろん行と列両方について動的な大きさを持つ行列を作ることが出来る。

ベクトルについては、行列の行か列の大きさを1とすれば良い。

Gigueの使い方―和, 差, 積―

演算子オーバーロードによって、思うように動作する。
しかし、Gigueは行列の大きさにうるさく、静的な大きさを持つ行列同士の演算についてその大きさに対して演算が適応できないのであればコンパイル時にエラーを出力する。

import std.stdio;
import std.random;

auto m33 = matrix!(3, 3, (i, j) => uniform(0, i + j + 1));           // 3 行 3 列の行列
auto v3 = (ones!int * 4.5).matrix!(3, 1);

v3 = m33 * v3;                              // 行列とベクトルの積

// v3 = v3.transpose * m33;                 // コンパイル時エラー
                                            // 左辺(3, 1), 右辺(1, 3)

v3 = (v3.transpose * m33).transpose;        // OK

Gigueの使い方―連立方程式を解く―

GigueにはLU分解法による連立方程式を解く関数が用意されている。
Gigueに実装されているLU分解法は、PA = LUと分解する方法である。

real[3][3] mStack = [[2, 4, 2],
                     [4, 10, 3],
                     [3, 7, 1]];

auto m = mStack.matrix();

SMatrix!(real, 3, 3) org = m;           // コピーを作る
auto plu = m.pluDecomposeInPlace();     // InPlaceで分解

auto v = matrix!(3, 1)(cast(real[])[8, 17, 11]);

plu.solveInPlace(v);                    // v について解く
foreach(i; 0 .. 3)
    assert(approxEqual(v[i], 1));       // チェック

Gigueに足りないもの

現在のGigueでは、上手く三角行列や帯行列, 疎行列といった特殊な行列を表現できない。
というのも、行列aとbの積を演算子オーバーロードによって行う場合に、a * bという式の演算量を、これらの行列を演算に使用したとしても減らすことができない。
理由は、単純にGigueは特殊な行列を理解しないからである。

私は、この機能のエレガントな実装をまだできないでいる。
思いつく方法としては、行列の非ゼロ要素の密度を行列型に提供させるように強制するか、UDAを用いる方法である。

さいごに

実は明日は電力系統工学のテストである。
私はまだ用語や電力潮流計算, 単位法について理解していない。
それにくわえて締め切りの迫った卒業研究である。
非常につらい。
追い込まれた心理状態で書いたこの記事は、支離滅裂で非常にひどい内容であろう。
そして私はこの記事を添削せずにQiitaに投稿しているだろう。
そのために、この記事はさらに粗悪な文章の固まりとなっているに違いない。