在 Eigen::Tensor 收缩中交换张量会导致不同的结果?
What causes different results by interchanging the tensors in an Eigen::Tensor contraction?
我正在开发一个依赖于张量收缩的 C++ 库。我不会在此处 post 完整的应用程序,但我已将其提炼为以下内容。
我们定义了一个 toy rank-4 张量,它只是 (0, 1, ..., 15) 重塑:
Eigen::Tensor<double, 4> T (2, 2, 2, 2);
for (size_t i = 0; i < 2; i++) {
for (size_t j = 0; j < 2; j++) {
for (size_t k = 0; k < 2; k++) {
for (size_t l = 0; l < 2; l++) {
T(i, j, k, l) = l + 2 * k + 4 * j + 8 * i;
}
}
}
}
和一个要收缩的 2 阶张量,它只不过是 (1, 2, 3, 4) 重塑:
Eigen::Tensor<double, 2> A (2, 2);
for (size_t i = 0; i < 2; i++) {
for (size_t j = 0; j < 2; j++) {
A(i, j) = 1 + j + 2 * i;
}
}
要在 Eigen 中收缩两个张量,我们必须指定一个收缩对。我们的目标是收缩张量的前两个指标,如 T(ijkl)*A(ib)=M(bjkl)
。以我目前对Eigen中Tensor模块的理解,我们将收缩对写成
Eigen::array<Eigen::IndexPair<int>, 1> contraction_pair = {Eigen::IndexPair<int>(0, 0)};
不过,我认为应该可以使用完全相同的收缩对来进行收缩A(ib)*T(ijkl)=N(bjkl)
。不幸的是,事实并非如此,M
的元素是
0 0 0 0 24
0 0 0 1 32
0 0 1 0 28
0 0 1 1 38
0 1 0 0 32
0 1 0 1 44
0 1 1 0 36
0 1 1 1 50
1 0 0 0 40
1 0 0 1 56
1 0 1 0 44
1 0 1 1 62
1 1 0 0 48
1 1 0 1 68
1 1 1 0 52
1 1 1 1 74
而N
中的这些是
0 0 0 0 24
0 0 0 1 28
0 0 1 0 32
0 0 1 1 36
0 1 0 0 40
0 1 0 1 44
0 1 1 0 48
0 1 1 1 52
1 0 0 0 32
1 0 0 1 38
1 0 1 0 44
1 0 1 1 50
1 1 0 0 56
1 1 0 1 62
1 1 1 0 68
1 1 1 1 74
我已经在 numpy 中使用 einsum 测试了相同的玩具张量:
T = np.arange(16).reshape(2, 2, 2, 2)
A = np.arange(1, 5).reshape(2, 2)
contraction1 = np.einsum('ijkl,ia->ajkl', integrals, C)
contraction2 = np.einsum('ia,ijkl->ajkl', C, integrals)
并且contraction1
和contraction2
都是
0 0 0 0 24
0 0 0 1 28
0 0 1 0 32
0 0 1 1 36
0 1 0 0 40
0 1 0 1 44
0 1 1 0 48
0 1 1 1 52
1 0 0 0 32
1 0 0 1 38
1 0 1 0 44
1 0 1 1 50
1 1 0 0 56
1 1 0 1 62
1 1 1 0 68
1 1 1 1 74
这与 Eigen 中的情况 A(ib)*T(ijkl)=N(bjkl)
一致。是什么导致 Eigen 在这两种情况下不给出相同的结果?
Eigen
界面似乎只采用收缩轴进行规范。因此它必须自行决定如何安排非收缩轴。显而易见的方法是保持原始轴的顺序,首先是第一个参数,然后是第二个参数。
为了确认这一点,我们可以
- 使用
np.einsum
指定不同的输出布局并与 Eigen
: 的输出进行比较
.
import numpy as np
T = np.arange(16).reshape(2, 2, 2, 2)
A = np.arange(1, 5).reshape(2, 2)
print(np.einsum('ijkl,ia->ajkl', T, A))
# [[[[24 28]
# [32 36]]
# [[40 44]
# [48 52]]]
# [[[32 38]
# [44 50]]
# [[56 62]
# [68 74]]]]
print(np.einsum('ijkl,ia->jkla', T, A))
# [[[[24 32]
# [28 38]]
# [[32 44]
# [36 50]]]
# [[[40 56]
# [44 62]]
# [[48 68]
# [52 74]]]]
- 或直接在
Eigen
中随机播放(感谢@lelemmen):
.
M.shuffle(Eigen::array<int, 4> {3, 0, 1, 2})
(M_shuffled == N).all()
# 1
我正在开发一个依赖于张量收缩的 C++ 库。我不会在此处 post 完整的应用程序,但我已将其提炼为以下内容。
我们定义了一个 toy rank-4 张量,它只是 (0, 1, ..., 15) 重塑:
Eigen::Tensor<double, 4> T (2, 2, 2, 2);
for (size_t i = 0; i < 2; i++) {
for (size_t j = 0; j < 2; j++) {
for (size_t k = 0; k < 2; k++) {
for (size_t l = 0; l < 2; l++) {
T(i, j, k, l) = l + 2 * k + 4 * j + 8 * i;
}
}
}
}
和一个要收缩的 2 阶张量,它只不过是 (1, 2, 3, 4) 重塑:
Eigen::Tensor<double, 2> A (2, 2);
for (size_t i = 0; i < 2; i++) {
for (size_t j = 0; j < 2; j++) {
A(i, j) = 1 + j + 2 * i;
}
}
要在 Eigen 中收缩两个张量,我们必须指定一个收缩对。我们的目标是收缩张量的前两个指标,如 T(ijkl)*A(ib)=M(bjkl)
。以我目前对Eigen中Tensor模块的理解,我们将收缩对写成
Eigen::array<Eigen::IndexPair<int>, 1> contraction_pair = {Eigen::IndexPair<int>(0, 0)};
不过,我认为应该可以使用完全相同的收缩对来进行收缩A(ib)*T(ijkl)=N(bjkl)
。不幸的是,事实并非如此,M
的元素是
0 0 0 0 24
0 0 0 1 32
0 0 1 0 28
0 0 1 1 38
0 1 0 0 32
0 1 0 1 44
0 1 1 0 36
0 1 1 1 50
1 0 0 0 40
1 0 0 1 56
1 0 1 0 44
1 0 1 1 62
1 1 0 0 48
1 1 0 1 68
1 1 1 0 52
1 1 1 1 74
而N
中的这些是
0 0 0 0 24
0 0 0 1 28
0 0 1 0 32
0 0 1 1 36
0 1 0 0 40
0 1 0 1 44
0 1 1 0 48
0 1 1 1 52
1 0 0 0 32
1 0 0 1 38
1 0 1 0 44
1 0 1 1 50
1 1 0 0 56
1 1 0 1 62
1 1 1 0 68
1 1 1 1 74
我已经在 numpy 中使用 einsum 测试了相同的玩具张量:
T = np.arange(16).reshape(2, 2, 2, 2)
A = np.arange(1, 5).reshape(2, 2)
contraction1 = np.einsum('ijkl,ia->ajkl', integrals, C)
contraction2 = np.einsum('ia,ijkl->ajkl', C, integrals)
并且contraction1
和contraction2
都是
0 0 0 0 24
0 0 0 1 28
0 0 1 0 32
0 0 1 1 36
0 1 0 0 40
0 1 0 1 44
0 1 1 0 48
0 1 1 1 52
1 0 0 0 32
1 0 0 1 38
1 0 1 0 44
1 0 1 1 50
1 1 0 0 56
1 1 0 1 62
1 1 1 0 68
1 1 1 1 74
这与 Eigen 中的情况 A(ib)*T(ijkl)=N(bjkl)
一致。是什么导致 Eigen 在这两种情况下不给出相同的结果?
Eigen
界面似乎只采用收缩轴进行规范。因此它必须自行决定如何安排非收缩轴。显而易见的方法是保持原始轴的顺序,首先是第一个参数,然后是第二个参数。
为了确认这一点,我们可以
- 使用
np.einsum
指定不同的输出布局并与Eigen
: 的输出进行比较
.
import numpy as np
T = np.arange(16).reshape(2, 2, 2, 2)
A = np.arange(1, 5).reshape(2, 2)
print(np.einsum('ijkl,ia->ajkl', T, A))
# [[[[24 28]
# [32 36]]
# [[40 44]
# [48 52]]]
# [[[32 38]
# [44 50]]
# [[56 62]
# [68 74]]]]
print(np.einsum('ijkl,ia->jkla', T, A))
# [[[[24 32]
# [28 38]]
# [[32 44]
# [36 50]]]
# [[[40 56]
# [44 62]]
# [[48 68]
# [52 74]]]]
- 或直接在
Eigen
中随机播放(感谢@lelemmen):
.
M.shuffle(Eigen::array<int, 4> {3, 0, 1, 2})
(M_shuffled == N).all()
# 1