c++ - Armadillo(C++)中的快速数组置换(广义张量转置)

我有一个项目，涉及3D数组(arma::Cube<cx_double>)上的许多置换。特别地，所需的排列是切片之间的列互换。在Matlab中，这可以通过permute(cube,[1,3,2])有效地计算，而在Python中，可以通过numpy.transpose(cube,axis=[0,2,1])有效地计算。

不幸的是， Armadillo 本身没有permute函数。我尝试了不同的方法，但是与Matlab相比，所有方法都相当慢。 我想知道在Armadillo 中置换(相当大)多维数据集的更快方法是什么。使用gprof对代码进行性能分析，大部分时间都花在了我下面尝试过的permute函数中，而在Matlab中，对于同一个移植项目，大部分时间都花在了SVD或QR矩阵分解上(重塑和置换非常快)在matlab中)。

我想了解在Armadillo中进行排列的最快方法，以及为什么某些方法比其他方法更好。

选项1:原始排列(最快的选项)(是否有更快的方法？)

输入多维数据集到输出多维数据集的逐元素分配。

template <typename T>
static Cube<T> permute (Cube<T>& cube){

uword D1=cube.n_rows;
uword D2=cube.n_cols;
uword D3=cube.n_slices;

Cube<T> output(D1,D3,D2);

for (uword s = 0; s < D3; ++s){
for (uword c = 0; c < D2; ++c){
for (uword r = 0; r < D1; ++r){
    output.at(r, s, c) = cube.at(r, c, s);
    // output[ D1*D3*c + D1*s+ r ] = cube[ D1*D2*s + D1*c + r ];

}
}
}

return output;
}

选项2:填充切片(非常慢)

通过非连续的subcube View 填充输出多维数据集的切片。

template <typename T>
static Cube<T> permute (Cube<T>& cube_in){

    uword D1 = cube_in.n_rows;
    uword D2 = cube_in.n_cols;
    uword D3 = cube_in.n_slices;

    Cube<T> output;
    output.zeros(D1, D3, D2);

    for (uword c=0; c<D2; ++c) {
        output.slice(c) = cube_in.subcube( span(0,D1-1),span(c),span(0,D3-1) );
    }

    return output;
}

选项3:转置图层(比Raw排列慢，但具有可比性)

我们可以遍历输入多维数据集的各层(固定行)并进行转置。

template <typename T>
static Cube<T> permute (Cube<T>& cube_in){
    // in a cube, permute {1,3,2} (permute slices by columns)

    uword D1 = cube_in.n_rows;
    uword D2 = cube_in.n_cols;
    uword D3 = cube_in.n_slices;

    if(D3 > D2){
        cube_in.resize(D1,D3,D3);
    } else if (D2 > D3) {
        cube_in.resize(D1,D2,D2);
    }

    for (uword r=0; r<D1; ++r) {
        static cmat layer = cmat(cube_in.rows(r,r));
        inplace_strans(layer);
        cube_in.rows(r,r)=layer;
    }

    cube_in.resize(D1,D3,D2);

    return cube_in;
}

选项4:查询表
通过读取 vector 中的indes来获得非连续访问。

template <typename T>
arma::Cube<T> permuteCS (arma::Cube<T> cube_in){
    // in a cube, permute {1,3,2} (permute slices by columns)
    uword D1 = cube_in.n_rows;
    uword D2 = cube_in.n_cols;
    uword D3 = cube_in.n_slices;

    cx_vec onedcube = cube_in.elem(gen_trans_idx(cube_in));

    return arma::Cube<T>(onedcube.memptr(), D1, D3, D2, true ) ;
}

其中gen_trans_idx是用于生成置换后的多维数据集的indes的函数:

template <typename T>
uvec gen_trans_idx(Cube<T>& cube){

    uword D1 = cube.n_rows;
    uword D2 = cube.n_cols;
    uword D3 = cube.n_slices;

    uvec perm132(D1*D2*D3);
    uword ii = 0;
    for (int c = 0; c < D2; ++c){
    for (int s = 0; s < D3; ++s){
    for (int r = 0; r < D1; ++r){
        perm132.at(ii) = sub2ind(size(cube), r, c, s);
        ii=ii+1;
    }}}

    return perm132;
}

理想情况下，如果预先确定多维数据集的尺寸，则可以预先计算这些查找表。

选项5 (原位转置)非常慢，有效存储

// Option: In-place transpose
template <typename T>
arma::Cube<T> permuteCS (arma::Cube<T> cube_in, uvec permlist ){

    T* Qpoint = cube_in.memptr(); // pointer to first element of cube_in

    uvec updateidx = find(permlist - arma::linspace<uvec>(0,cube_in.n_elem-1,cube_in.n_elem)); // index of elements that change position in memory
    uvec skiplist(updateidx.n_elem,fill::zeros);

    uword rr = 0; // aux index for updatelix
    for(uword jj=0;jj<updateidx.n_elem;++jj){

        if(any(updateidx[jj] == skiplist)){ // if element jj has already been updated
            // do nothing
        } else {

            uword scope = updateidx[jj];
            T target = *(Qpoint+permlist[scope]);        // store the value of the target element

            while(any(scope==skiplist)-1){              // while wareyou has not been updated

                T local  = *(Qpoint+scope);                  // store local value
                *(Qpoint+scope) = target;
                skiplist[rr]=scope;
                ++rr;
                uvec wareyou = find(permlist==scope);        // find where the local value will appear
                scope = wareyou[0];
                target = local;
            }

        }
    }
   cube_in.reshape(cube_in.n_rows,cube_in.n_slices,cube_in.n_cols);

    return cub

e_in;
}

最佳答案

此代码是我对memcpy hack的评论的补充。另外，不要忘记尝试添加const reference以防止复制对象。

template <typename T>
static Cube<T> permute(const Cube<T> &cube){
    const uword D1 = cube.n_rows;
    const uword D2 = cube.n_cols;
    const uword D3 = cube.n_slices;
    const uword D1_mul_D3 = D1 * D3;
    const Cube<T> output(D1, D3, D2);

    const T * from = cube.memptr();
    T *to = output.memptr();

    for (uword s = 0; s < D3; ++s){
        T *to_tmp = to + D1 * s;
        for (uword c = 0; c < D2; ++c){
            memcpy(to_tmp, from, D1 * sizeof(*from));
            from += D1;
            to_tmp += D1_mul_D3;
        }
    }

    return output;
}

关于c++ - Armadillo(C++)中的快速数组置换(广义张量转置)，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/51842301/