Python でテンソル代数．行列化．

以前の質問に関連した内容です．

今度は，テンソルXの第n軸に沿って，Xを行列Y_nに直す，「テンソルのモードn展開」の実装に手間取っています．Xとnを与えられたときにY_nを返す関数をPythonでNumPyで設計したいです．

モードn展開の定義

テンソルX∈R^{I_1×I_2×…×I_N}をモードnで展開することで，行列Y_n∈R^{I_n×(Π_{m≠n}^{N} I_m)}にします．

テンソルXの(i_1,i_2,…,i_N)要素は行列Y_(n)の(i_n,j)に移ります．ここで，jは
j = 1 + Σ_{k=1,k≠n}^N (i_k-1)J_k
J_k = Π_{m=1,m≠n}^{k-1} I_m
を満たす整数です．

参考

上記の定義は http://www.kolda.net/publication/TensorReview.pdf の p460を参照しています．

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

回答1件

ベストアンサー

ご確認ください。

python
1def mode_n_expansion(x, n):
2    x = np.asarray(x)
3    if n <= 0 or n % 1 != 0 or x.ndim < n:
4        raise ValueError((x.shape, n))
5    return np.moveaxis(x, int(n) - 1, -1).T.reshape(x.shape[int(n) - 1], -1)

テストコード:

python
1def test_mode_n_expand():
2    X = np.random.random(np.random.randint(1, 5, np.random.randint(3, 5)))
3    for n in range(X.ndim):
4        import itertools as it
5        for i in it.product(*map(range, X.shape)):
6            Y = mode_n_expansion(X, n=n + 1)
7            j = int(sum(
8                i_k * np.product([I_m for m, I_m in enumerate(X.shape[:k]) if m != n])
9                for k, i_k in enumerate(i)
10                if k != n
11            ))
12            assert X[i] == Y[i[n], j]