Я работаю над библиотекой C ++, которая опирается на тензорные сокращения. Я не буду публиковать полное заявление здесь, но я разобрал его до следующего.
Мы определяем игрушечный тензор ранга 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)
, С моим текущим пониманием модуля Tensor в Eigen, мы запишем сжатую пару как
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
Я тестировал те же игрушечные тензоры в numy, используя 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
которые совпадают с делом A(ib)*T(ijkl)=N(bjkl)
в Эйгене. Что заставляет Эйгена не давать одинаковый результат в обоих случаях?
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
Других решений пока нет …