连接特征张量

Question

我正在使用 Eigen 中的张量，需要将两个张量连接在一起以形成一个更高阶的张量。举个例子：

C_ijkl = A_ik B_jl

我知道我可以使用一些 for 循环来做到这一点，但有更好的方法吗？

编辑： 按照下面的答案，我将其实现到代码中。使用带有 -O3 的 g++ 编译器并将显式循环与下面显示的收缩进行比较，我得到了 10,000 次迭代的以下结果：

Average time (ns) for explicit loop: 196.782
Average time (ns) for contraction  : 790.439

因此，在我的测试设置中，天真的显式循环似乎比收缩更快（尽管更混乱）。

代码：

#include <iostream>
#include <Eigen/Core>
#include <Eigen/Dense>
#include <unsupported/Eigen/CXX11/Tensor>

#include<chrono>
typedef std::chrono::high_resolution_clock Clock;

using namespace Eigen;

int main() {

//Initialize
Tensor<double, 2> A_ij(4, 4);
Tensor<double, 2> B_kl(4, 4);
Tensor<double, 4> C_ijkl(4, 4, 4, 4);

A_ij.setValues({{ 1, 2, 3, 4},
                { 5, 6, 7, 8},
                { 9,10,11,12},
                {13,14,15,16}});

B_kl.setValues({{17,18,19,20},
                {21,22,23,24},
                {25,26,27,28},
                {29,30,31,32}});

Tensor<double, 4> Ctrue_ijkl(4, 4, 4, 4);

int nrepeat = 10000;
double avg_time_explicit;
double avg_time_contract;

auto t1 = Clock::now();
auto t2 = Clock::now();

array<Eigen::IndexPair<long>,0> empty_index_list = {};

for (int r=0; r<nrepeat; r++){
    t1 = Clock::now();
    for (int i=0; i<4; i++){
        for (int j=0; j<4; j++){
            for (int k=0; k<4; k++){
                for (int l=0; l<4; l++){
                    Ctrue_ijkl(i,j,k,l) = A_ij(i,j)*B_kl(k,l);
                }
            }
        }
    }
    t2 = Clock::now();
    avg_time_explicit += std::chrono::duration_cast<std::chrono::nanoseconds>(t2-t1).count()/static_cast<double>(nrepeat);
}



for (int r=0; r<nrepeat; r++){
    t1 = Clock::now();
    C_ijkl = A_ij.contract(B_kl, empty_index_list);
    t2 = Clock::now();
    avg_time_contract += std::chrono::duration_cast<std::chrono::nanoseconds>(t2-t1).count()/static_cast<double>(nrepeat);
}

//Check to make sure the two methods produce the same result
Map<VectorXd> mt(C_ijkl.data(), C_ijkl.size());
Map<VectorXd> mr(Ctrue_ijkl.data(), Ctrue_ijkl.size());

std::cout << "Ctrue_ijkl == C_ijkl: " << mt.isApprox(mr) << "\n";

std::cout << "Average time (ns) for explicit loop: " << avg_time_explicit << "\n";
std::cout << "Average time (ns) for contraction:   " << avg_time_contract << "\n";

}

Answer 1

我相信这个问题与：How do I do outer product of tensors in Eigen? 相同，我也回答了以下问题：

你必须在没有索引的情况下收缩。

Eigen::array<Eigen::IndexPair<long>,0> empty_index_list = {};
Tensor<double, 2> A_ij(4, 4);
Tensor<double, 2> B_kl(4, 4);
Tensor<double, 4> C_ijkl = A_ij.contract(B_kl, empty_index_list);