Created
November 30, 2017 09:49
-
-
Save lelemmen/d01a32387363376d00191d10208a5440 to your computer and use it in GitHub Desktop.
Undocumented behavior of Eigen::Tensor contractions
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#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