Skip to content

Instantly share code, notes, and snippets.

@lelemmen
Created November 30, 2017 09:49
Show Gist options
  • Save lelemmen/d01a32387363376d00191d10208a5440 to your computer and use it in GitHub Desktop.
Save lelemmen/d01a32387363376d00191d10208a5440 to your computer and use it in GitHub Desktop.
Undocumented behavior of Eigen::Tensor contractions
#include <iostream>
#include <unsupported/Eigen/CXX11/Tensor>
// Inspired by my question and its accepted answer on StackOverflow (https://stackoverflow.com/questions/47556726), I've decided to create this in a little gist.
// This Gist is about some undocumented (https://bitbucket.org/eigen/eigen/src/default/unsupported/Eigen/CXX11/src/Tensor/README.md?fileviewer=file-view-default) behavior of Eigen::Tensor contractions.
/** Print the contents of a rank-four tensor in a fashionable way
*/
void print(const Eigen::Tensor<double, 4>& T) {
auto dim = T.dimension(0);
for (size_t i = 0; i < dim; i++) {
for (size_t j = 0; j < dim; j++) {
for (size_t k = 0; k < dim; k++) {
for (size_t l = 0; l < dim; l++) {
std::cout << i << ' ' << j << ' ' << k << ' ' << l << " " << T(i, j, k, l) << std::endl;
}
}
}
}
}
int main () {
// Let's create a toy rank-4 tensor (which is equal to np.arange(15).reshape(2, 2, 2, 2))
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;
}
}
}
}
// and a toy rank-2 tensor (which is equal to np.arange(1, 5).reshape(2, 2))
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;
}
}
// The goal of this Gist is to research the behavior of tensor contractions. Let's therefore create a contraction pair:
Eigen::array<Eigen::IndexPair<int>, 1> contraction_pair = {Eigen::IndexPair<int>(0, 0)};
// This contraction pair allows us to do a contraction of the form
// T(ijkl) A(ib)
Eigen::Tensor<double, 4> M = T.contract(A, contraction_pair);
std::cout << "M is" << std::endl;
print(M);
std::cout << std::endl;
// M is
// 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
// Since we have two equal indices in this contraction pair, we could as well have done the contraction of the form
// A(ib) T(ijkl)
Eigen::Tensor<double, 4> N = A.contract(T, contraction_pair);
std::cout << "N is" << std::endl;
print(N);
std::cout << std::endl;
// N is
// 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
// As we can see, the results are not equal. The accepted answer (https://stackoverflow.com/a/47558349/7930415) provides insight in how the contractions are done in Eigen:
// The surviving axes are kept in order (from left to right)
// According to this rule, we have
// T(ijkl) A(ib) = M(jklb)
// A(ib) T(ijkl) = N(bjkl)
// We can check the validity of this rule by permuting (shuffling) M's axes ...
Eigen::Tensor<double, 4> M_shuffled = M.shuffle(Eigen::array<int, 4> {3, 0, 1, 2}); // this means that M_shuffled's first axis becomes N's fourth axis
// ... and checking if the result is true for every value in N
std::cout << (M_shuffled == N).all() << std::endl;
// 1
};
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment