python でテンソル代数．モードn積の実装について．

テンソル代数を用いた研究をしています．この度，テンソルのタッカー分解のpython実装が求められています．このタッカー分解の実装にあたって，テンソルのモードn積を実行する必要があります．テンソル代数では頻繁に登場する「テンソルと行列」の演算ですが，NumPyを調べたところテンソルのモード積については，自前で用意しないといけないようです．

モードn積の定義

モードn積は，テンソル X と行列 M の間に成立する二項演算です．
テンソルと行列の次元を次のようにします．
X ∈ R^{l_1×l_2×…×l_N}
M ∈ R^{J×l_n}

ここで，Rは実数．l_1,l_2,…,l_N と J はゼロを含まない自然数です．nをN以下の自然数とするとき，テンソルXと行列Mの間のモードn積 X ⊗_n M は次のように定義されます．

( X ⊗_n M ) ∈ R^{l_1×l_2×…×l_(n-1)×J×l_(n+1)×…×l_N}
( X ⊗_n M ){i_1,…,i{n-1},j,i_{n+1},…,i_N} = Σ_{i_n}^{l_n} X_{i_1,…,i_{n-1},i_n,i_{n+1},…,i_N} m_{j,i_n}

これをてっとり早くpythonで実装するにはどうしたらよいでしょうか．どうも大量のfor文ループになるようで困惑しています．

参照

モードn積の定義やRでの実装は以下を参考にしています．
https://www.alexejgossmann.com/tensor_decomposition_tucker/

np.tensordot で対応できると良いのですが...

実装

たとえば，n=3 だあれば，こんな感じかと思います．
絶対にもっと良い実装方法があるとは思うのですが...

Python
1l1 = 3
2l2 = 3
3l3 = 3
4l4 = 3
5J = 2
6X = np.random.uniform(0,1,(l1,l2,l3,l4))
7M = np.random.uniform(0,1,(J,l3))
8
9def mode_3_prodcut(X,M,mode=3):
10    A = np.zeros((l1,l2,J,l4))
11    for i1 in range(l1):
12        for i2 in range(l2):
13            for j in range(J):
14                for i4 in range(l4):
15                    term = 0
16                    for i3 in range(l3):
17                        term += X[i1,i2,i3,i4] * M[j,i3]
18                        
19                    A[i1,i2,j,i4] = term
20                    
21    return A

行動規範の内容に同意します

回答1件

ベストアンサー

モードn積については知らなかったのですが、定義に従うと以下のコードになると思います。

python
1def mode_n_product(x, m, mode):
2    x = np.asarray(x)
3    m = np.asarray(m)
4    if mode <= 0 or mode % 1 != 0:
5        raise ValueError('`mode` must be a positive interger')
6    if x.ndim < mode:
7        raise ValueError('Invalid shape of X for mode = {}: {}'.format(mode, x.shape))
8    if m.ndim != 2:
9        raise ValueError('Invalid shape of M: {}'.format(m.shape))
10    return np.swapaxes(np.swapaxes(x, mode - 1, -1).dot(m.T), mode - 1, -1)

補足:

np.swapaxes(x, m, n) :
テンソルxのm個目の添字とn個目の添字を交換する関数です。
ただし添字の番号は0始まりで、-1は末尾の添字を表します。
したがって、np.swapaxes(x, mode - 1, -1)で i_n と i_N が交換されて
X'{i_1, i_2, …, i{n-1}, i_N, i_{n+1} …, i_{N-1}, i_n}
:= X_{i_1, i_2, …, i_{n-1}, i_n, i_{n+1} …, i_{N-1}, i_N}
が得られます。
m.T:
Mの転置です。M'{i_n, j} := M{j, i_n} が得られます。
c = a.dot(b):
b が行列であれば以下の関係が成り立ちます。
C_{…, i, j} = Σ_k A_{…, i, k} B_{k, j}
したがって np.swapaxes(x, mode - 1, -1).dot(m.T) と書くと
A'{i_1, …, i{n-1}, i_N, i_{n+1}, …, i_{N-1}, j}
:= Σ_{i_n}^{l_n} X'{i_1, i_2, …, i{n-1}, i_N, i_{n+1} …, i_{N-1}, i_n} M'{i_n, j}
= Σ{i_n}^{l_n} X_{i_1, i_2, …, i_{n-1}, i_n, i_{n+1} …, i_{N-1}, i_N} M_{j, i_n}
が得られます。
最後に、再び np.swapaxes(x, mode - 1, -1) を適用することで
( X ⊗_n M ){i_1, …, i{n-1}, j, i_{n+1}, …, i_N}
= A'{i_1, …, i{n-1}, i_N, i_{n+1}, …, i_{N-1}, j}
を得ます。