连接自己的张量
Concatenating Eigen Tensors
我正在使用 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";
}
我认为这个问题与:How do I do outer product of tensors in Eigen? 相同,我也回答了以下问题:
You have to contract over no indices.
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);
我正在使用 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";
}
我认为这个问题与:How do I do outer product of tensors in Eigen? 相同,我也回答了以下问题:
You have to contract over no indices.
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);