C++
メタプログラミング
設計
テンソル

概要

少しだけ時間が取れたためC++でテンソル型を作成することで思いついたことをメモしていきたいと思います。

共変性や反変性の扱い

ベクトルは反変ベクトルであれば、ベクトル空間$V$によるテンソル空間$T^{1}(V)$の元として定義できます。同様にして行列は$T^{1}_{1}(V)$の元として定義できたかと思います。
そして、別のベクトル空間$W$とその元$w^{j}\in W$となる線型変換は$v^{i}\in T^{1}(V)$と$\alpha^{j}_{k}\in T^{1}_{1}(V)$および縮約を示す写像$C^{i}_{k}$によって
$$
w^{j}=C^{i}_{k}(\alpha^{j}_{k}v^{i})
$$
とあらわされたかと思います。また、ほとんどの場合は$C^{i}_{k}$を省略して
$$
w^{j}=\alpha^{j}_{i}v^{i}
$$
とあらわすことが多いと思います(個人的にはテンソル積があたかも可換にあらわされることがあるため嫌いです、ただ縮約としてだけなら便利)。
さて、ここまでテンソルの基底を意識なんてしていませんでしたが$V$の双対空間$V^{\ast}$と$V$によって$T^{1}_{1}(V)\cong V\otimes V^{\ast}$となることを考えれば双対基底と反変基底を考慮しなければならないことがわかります。
ここで、プログラムによる実装を考えると、基本的に正規標準基底以外は使わないと考えられます(使いたいなら基底変換すればいいだけだし実質いらないはず)。このとき、反変ベクトルと共変ベクトルは一致するためテンソルの反変と共変は区別しなくてもいいとすることができます。
というわけでそれに合わせた仕様を考えていきます。

次元と添え字の定義

テンソル型では反変と共変を意識しなくてもいいため、意識すればいいのはテンソル空間の構成因子となるベクトル空間の次元と添え字の記号だけでいいはず。この構築には前のやつも含めて考えます。
というわけで基本形は

namespace iml {

    //順序対の定義(何故かエイリアスで"using set = index_tuple<Indices...>"としたときエラー)
    template <size_t...>
    struct set;

    //テンソル型(Dims:次元,Suffixes:添え字)
    template <class, class...>
    class tensor;
    template <class T, size_t... Dims, size_t... Suffixes>
    class tensor<T, set<Dims, Suffixes>...> {
        //多次元配列
        typename multi_array<T, Dims...>::type x;

    public:

    };
}

といった感じの実装になるでしょう。
これに対して、メタを与えます。まずは配列の次元を取得するメタです。

namespace iml {

    //配列の次元の取得
    template <size_t, class>
    struct _Dimension;
    template <size_t Dim, size_t N>
    struct _Dimension<Dim, index_tuple<N>> {
        //次元0の配列はいけない
        static_assert(Dim*N > 0, "0 parameter should not exist.");
        static constexpr size_t value = Dim * N;
    };
    template <size_t Dim, size_t First, size_t... Indices>
    struct _Dimension<Dim, index_tuple<First, Indices...>> : _Dimension<Dim*First, index_tuple<Indices...>> {};
    template <size_t First, size_t... Indices>
    struct dimension : _Dimension<First, index_tuple<Indices...>> {};
}

そこまで難しいものでもなく、ただ単にそれぞれの配列の次元の積をとっているだけです。
次に考えるべきことは添え字です。テンソル内で同じ添え字は存在しないべきであるため、そのときにはstatic_assertによるコンパイルエラーを出力させるのではなく根本的に定義不可能という風にしたいと思います。というわけで、まずは同じ記号が存在するかの探索するメタを定義します。

namespace iml {

    //同じ記号が存在するか
    template <bool f, size_t... Indices>
    struct _Is_exists_same_signs {
        static constexpr bool value = f;
    };
    //要素が2個以上存在するかつ同じ要素が見つかっていないとき
    template <size_t First, size_t Second, size_t... Indices>
    struct _Is_exists_same_signs<false, First, Second, Indices...>
        : _Is_exists_same_signs<(First == Second) || _Is_exists_same_signs<false, Second, Indices...>::value, First, Indices...> {};
    template <size_t... Indices>
    struct is_exists_same_signs : _Is_exists_same_signs<false, Indices...> {};
}

これを用いることでクラスの構成をできなくしたいのですが、部分特殊化されたテンプレートに対してデフォルトのテンプレートを定義することができないので、その場しのぎ的なものでとりあえずコンパイルエラーを引き起こすことにします。

namespace iml {

    //条件式が真のときのみ型が有効
    template<bool, class T = void>
    struct enable_if {};
    template<class T>
    struct enable_if<true, T> {
        using type = T;
    };

    //条件式に応じてクラスの構成でコンパイルエラーを引き起こす補助
    template <bool f, class = typename enable_if<f>::type>
    class class_if {};

    //テンソル型(Dims:次元,Suffixes:添え字)
    template <class, class...>
    class tensor;
    template <class T, size_t... Dims, size_t... Suffixes>
    class tensor<T, set<Dims, Suffixes>...> : class_if<!is_exists_same_signs<Suffixes...>::value> {
        //多次元配列
        typename multi_array<T, Dims...>::type x;

    public:

    };
}

多分機能を追加するうちにclass_ifを別の機能と統合するとは思いますが、現状はこんなもんでいいでしょう。

テンソル積に関するメタ

テンソル積や縮約をするには戻り値型が必要になりますが、それを決定するメタを定義する必要が(多分)あります。例えば、添え字に対して縮約を実行するならば

namespace iml {

    //テンソルの縮約(違う添え字かつSuffixesに存在しなければならない)
    template <size_t M, size_t N, class T, size_t... Dims, size_t... Suffixes
        , class = typename enable_if<(M != N) && is_exists_same_signs<M, Suffixes...>::value && is_exists_same_signs<N, Suffixes...>::value>::type>
    auto contraction(const tensor<T, set<Dims, Suffixes>...>& t) {

    }
}

といった感じになります。外部でメタ処理をしてautoは自動決定するかもとかと考えたりもしましたが、それをまだ思いついていないのでとりあえず戻り値型を取得できるメタを定義しておきましょうみたいな感じです。多分そのうち役に立つときが来ると思うというか思いたいですし。
色々と実装は面倒そうですが、よく考えると同じ添え字のペアが存在したときそれを除去するような操作をするだけということがわかります。例えば、テンソルの添え字の組$(i,j,k)$と$(l,m,j)$と与えられれば、一旦$(i,j,k,l,m,j)$といった組を生成して、空の組$()$に順に探索対象の添え字が元データに対して2つ以上同じものが存在しないならば順に入れていく感じになります。
$$
\begin{align}
(i,j,k,l,m,j) \{\} \ &\Rightarrow\ () \\
(i,j,k,l,m,j) \{i\} \ &\Rightarrow\ (i) \\
(i,j,k,l,m,j) \{j\} \ &\Rightarrow\ (i) \\
(i,j,k,l,m,j) \{k\} \ &\Rightarrow\ (i,k) \\
(i,j,k,l,m,j) \{l\} \ &\Rightarrow\ (i,k,l) \\
(i,j,k,l,m,j) \{m\} \ &\Rightarrow\ (i,k,l,m) \\
(i,j,k,l,m,j) \{j\} \ &\Rightarrow\ (i,k,l,m) \\
\end{align}
$$
ちなみに$\{\}$の中身は現在操作対象の添え字です。そのため、まずは添え字探索のためのメタを与えます。

namespace iml {

    //Nと同じ記号のカウント
    template <size_t, size_t, class>
    struct _Find_same_signs;
    template <size_t Cnt, size_t N>
    struct _Find_same_signs<Cnt, N, set<>> {
        static constexpr size_t value = Cnt;
    };
    template <size_t Cnt, size_t N, size_t First, size_t... Indices>
    struct _Find_same_signs<Cnt, N, set<First, Indices...>> : _Find_same_signs<Cnt + (N == First), N, set<Indices...>> {};
    template <size_t, class>
    struct find_same_signs;
    template <size_t N, size_t... Indices>
    struct find_same_signs<N, set<Indices...>> : _Find_same_signs<0, N, set<Indices...>> {};
}

これを用いることで2つのテンソルのテンソル積から縮約を考慮した新しい型を取得するメタを作成することができます。

namespace iml {

    //任意の位置の整数パラメータの取得
    template <size_t N, size_t First, size_t... Args>
    struct at_sign : at_sign<N - 1, Args...> {
        //テンプレート引数の範囲内に存在しない
        static_assert(sizeof...(Args) + 1 > N, "N is greater than the number of template arguments.");
    };
    template <size_t First, size_t... Args>
    struct at_sign<0, First, Args...> {
        static constexpr size_t value = First;
    };


    //同じ記号が存在するとき除去して新しいテンソル型を得る
    //縮約するか,次の探索対象が存在するか(Temp_S...の第1が現在の探索対象),その他
    template <bool, bool, class, class, class, class>
    struct _Create_tensor_type;
    //Res:結果,Temp:一時データ,Org:オリジナル
    //縮約対象
    template <class T, class... Res_set, size_t Temp_D_F, size_t Temp_S_F, class Org>
    struct _Create_tensor_type<true, true, T
        , _Arg_types<Res_set...>
        , _Arg_types<set<Temp_D_F, Temp_S_F>>
        , Org> {
        using type = tensor<T, Res_set...>;
    };
    //縮約対象でない
    template <class T, class... Res_set, size_t Temp_D_F, size_t Temp_S_F, class Org>
    struct _Create_tensor_type<false, true, T
        , _Arg_types<Res_set...>
        , _Arg_types<set<Temp_D_F, Temp_S_F>>
        , Org> {
        using type = tensor<T, Res_set..., set<Temp_D_F, Temp_S_F>>;
    };
    //縮約対象
    template <class T, class... Res_set, size_t Temp_D_F, size_t Temp_S_F, size_t... Temp_D, size_t... Temp_S, class Org>
    struct _Create_tensor_type<true, false, T
        , _Arg_types<Res_set...>
        , _Arg_types<set<Temp_D_F, Temp_S_F>, set<Temp_D, Temp_S>...>
        , Org>
        :_Create_tensor_type<find_same_signs<at_sign<0, Temp_S...>::value, Org>::value == 2, sizeof...(Temp_S) < 2, T
        , _Arg_types<Res_set...>
        , _Arg_types<set<Temp_D, Temp_S>...>
        , Org> {};
    //縮約対象でない
    template <class T, class... Res_set, size_t Temp_D_F, size_t Temp_S_F, size_t... Temp_D, size_t... Temp_S, class Org>
    struct _Create_tensor_type<false, false, T
        , _Arg_types<Res_set...>
        , _Arg_types<set<Temp_D_F, Temp_S_F>, set<Temp_D, Temp_S>...>
        , Org>
        :_Create_tensor_type<find_same_signs<at_sign<0, Temp_S...>::value, Org>::value == 2, sizeof...(Temp_S) < 2, T
        , _Arg_types<Res_set..., set<Temp_D_F, Temp_S_F>>
        , _Arg_types<set<Temp_D, Temp_S>...>
        , Org> {};
    //0階のテンソルを含む場合と分けて定義
    template <class, class>
    struct create_tensor_type;
    template <class T, size_t... Dims1, size_t... Suffixes1, size_t... Dims2, size_t... Suffixes2>
    struct create_tensor_type<tensor<T, set<Dims1, Suffixes1>...>, tensor<T, set<Dims2, Suffixes2>...>>
        : _Create_tensor_type<(find_same_signs<at_sign<0, Suffixes1..., Suffixes2...>::value, set<Suffixes1..., Suffixes2...>>::value == 2)
        , (sizeof...(Suffixes1) + sizeof...(Suffixes2) < 2), T
        , _Arg_types<>
        , _Arg_types<set<Dims1, Suffixes1>..., set<Dims2, Suffixes2>...>
        , set<Suffixes1..., Suffixes2...>> {};
    template <class T>
    struct create_tensor_type<tensor<T>, tensor<T>> {
        using type = tensor<T>;
    };


    //テンソルの縮約(違う添え字かつSuffixesに存在しなければならない)
    template <size_t M, size_t N, class T, size_t... Dims, size_t... Suffixes
        , class = typename enable_if<(M != N) && is_exists_same_signs<M, Suffixes...>::value && is_exists_same_signs<N, Suffixes...>::value>::type>
    typename create_tensor_type<tensor<T, set<Dims, Suffixes>...>, tensor<T, set<1, M>, set<1, N>>>::type
        contraction(const tensor<T, set<Dims, Suffixes>...>& t) {
        //中身はまだない
        return typename create_tensor_type<tensor<T, set<Dims, Suffixes>...>, tensor<T, set<1, M>, set<1, N>>>::type();
    }
}

テンソルのパラメータは多かったりするためそれ相応に長くなっていますが、実際は結構単純だったりします。それでも十分に複雑なのは否めませんが。多分もう少し単純化できる気がするけど思いつかない。
あと、プログラムを組んでいる最中にコンパイラのバグ(?)に遭遇してそれを回避するために色々と書き換えたりと思ったより大変でした。
とりあえずざっとイメージ的なコードです。添え字は整数パラメータですが多分charを使う方が多いと思う。

int main() {

    iml::tensor<float, iml::set<1, 'i'>, iml::set<3, 'j'>, iml::set<1, 'k'>, iml::set<3, 'l'>> temp;
    auto temp2 = iml::contraction<'j', 'k'>(temp);

    return 0;
}

ちなみに環境はVC++2017です。まぁ趣味程度なら比較的使いやすくて便利です。C++11の頃から改善していないコンパイラのバグがいくつかあったりしますが回避は可能なのでいいでしょう(実はVC++2015で回避不可なバグがあってVC++2017に移行したりとか)。

終わりに

実際のコードはもっと自作のライブラリ依存部が多いため色々とまとめる分には問題ない程度に省略していますが、まぁそれはそれでいいでしょうみたいな。
演算とか中身の実装は何もしていませんがそれはまた追々考えていこうかと思います。今回はあくまでもイメージ程度です。
テンソル積を縮約をしながら演算をする効率的なアルゴリズムは思いついていないので、今後はそれをのんびりと考えていきたいと思います。多分バタフライ演算あたりが参考になるかな。
縮約なら縮約する添え字と次元に応じたイテレータ範囲を取得するメタを作れば結構簡単に作ることはできる気はします(メモリが連続しているのでランダムアクセスができて楽という点で)。というか縮約ができればテンソル積に縮約を適用するだけでいいのですが速度がゴミになること間違いなしなので前述のように考えていきたいと思います。

追記

2018/10/12
誤字の修正