连续特征张量

Question

我和本征中的张量一起工作，需要把两个张量连在一起，形成一个高阶张量。例如：

C_ijkl = A_ik B_jl

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

编辑:按照下面的答案，我将其实现为代码。使用g++编译器和-O3，并将显式循环与下面所示的收缩进行比较，我在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 empty_index_list = {}；Tensor A_ij(4，4)；Tensor B_kl(4，4)；Tensor C_ijkl = A_ij.contract(B_kl，empty_index_list)；