快速的方式來遍歷不同形狀的多維數組？

向量操作可以通過在多個索引上廣播每個索引操作來簡潔和快速地編寫。例如，將矢量複製到更大的矢量中（例如，在FFT卷積之前進行零填充時）。但是，如果數組具有不同的形狀，則張量（多維數組）具有不同的連續整數索引。例如，2x4矩陣中的元素（1,3）處於平整整數索引1 * 3 + 4 = 7，但3x5矩陣中的相同元素具有索引1 * 5 + 3 = 8（請參見下面的示例）。快速的方式來遍歷不同形狀的多維數組？

所以複製矩陣到一個更大的矩陣是更棘手。如果你知道在編譯時的形狀，你可以只寫嵌套的循環：

typedef unsigned long*__restrict const tup_t; 
typedef const unsigned long*__restrict const const_tup_t; 
void nested_for_loops(const_tup_t shape) { 
    // Writing x.dimension separate nested for loops: 
    for (unsigned long i0=0; i0<shape[0]; ++i0) 
    for (unsigned long i1=0; i1<shape[1]; ++i1) 
     for (unsigned long i2=0; i2<shape[2]; ++i2) 
     // ... 
     { 
      // Inside innermost loop: 
      unsigned long x_index = ((i0*shape[1] + i1)*shape[2] + i2)*shape[3] /* + ... using each loop variable once */ ; 
     // Perform operations on x.flat[x_index] for some tensor x 
     // in global scope: func() 
     } 
}

在此之前的方案：

x_index = T ·s的·s的 ·秒·s的 ...小號_d-1 +噸·s的·s的·s的 ...小號_d-1 +噸·s的·s的 ...小號_d-1 + ... +噸_{D- 2}·s _d-1 + t _d-1。

但是，如果您不知道編譯時的維數（因爲您需要知道for循環的數量），這是不可能的。解決這個問題的一個方法是使用元組索引，在每次迭代中增加列，並在碰到張量邊界（形狀值）後執行進位操作。與形狀（2,2,2）張量的一個例子可能看起來像這樣：

但是，代碼涉及if語句，該語句轉換爲代碼中的分支以執行進位操作。

另一種方法是簡單地用形狀X的張量的平坦整數索引重新映射到在具有不同形狀的張量扁平索引Y（這可以用模運算來完成）：

inline unsigned long reindex(unsigned long index, const_tup_t shape, 
          const_tup_t new_shape, unsigned int dimension) { 
    unsigned long new_index = 0; 
    unsigned long new_axis_product_from_right = 1; 
    for (int i=dimension−1; index>0 && i>=0; −−i) { 
    unsigned long next_axis = shape[i]; 
    unsigned long new_next_axis = new_shape[i]; 

    unsigned long next_value = index % next_axis; 

    new_index += next_value * new_axis_product_from_right; 
    index /= next_axis; 

    new_axis_product_from_right *= new_next_axis; 
    } 
    return new_index; 
}

這消除如果語句，但它確實有模和除法操作，這不會像加法或乘法那樣快。當張量具有所有軸爲2的冪的形狀時，可以通過位旋轉來加速，用&和/用>>代替％操作。

現在的問題是，這些方法中的哪一種在實踐中更快？當然，有多維數組的庫（例如boost），但它們似乎要求在編譯時已知數組的維數，並且當張量具有不同的形狀時，像scala或go這樣的某些映射函數是非常棘手的。

來源

2017-04-11 Flow

玩了一段時間後導致了另一種方法，使我們可以用模板元編程結合C++ 11的可變參數模板和lambda功能展開的for循環所需的數字：

template <unsigned int DIMENSION> 
inline unsigned long tuple_to_index_fixed_dimension(const_tup_t tup, const_tup_t shape) { 
    unsigned long res = 0; unsigned int k; 
    for (k=0; k<DIMENSION−1; ++k) { 
    res += tup[k]; 
    res *= shape[k+1]; 
    } 
    res += tup[k]; 
    return res; 
} 

template <unsigned int DIMENSION, unsigned int CURRENT> 
class ForEachFixedDimensionHelper { 
public: 
    template <typename FUNCTION, typename ...TENSORS> 
    inline static void apply(tup_t counter, const_tup_t shape, FUNCTION function, TENSORS & ...args) { 
    for (counter[CURRENT]=0; counter[CURRENT]<shape[CURRENT]; ++counter[CURRENT]) 
     ForEachFixedDimensionHelper<DIMENSION−1, CURRENT+1>::template apply<FUNCTION, TENSORS...>(counter, shape, function, args...); 
    } 
}; 

template <unsigned int CURRENT> 
class ForEachFixedDimensionHelper<1u, CURRENT> { 
public: 
    template <typename FUNCTION, typename ...TENSORS> 
    inline static void apply(tup_t counter, const_tup_t shape, FUNCTION function, TENSORS & ...args) { 
    for (counter[CURRENT]=0; counter[CURRENT]<shape[CURRENT]; ++counter[CURRENT]) 
     function(args[tuple_to_index_fixed_dimension<CURRENT+1>(counter, args.data_shape())]...); /* tensor.data_shape() is an accessor for returning the shape member. */ 
    } 
}; 

template <unsigned char DIMENSION> 
class ForEachFixedDimension { 
public: 
    template <typename FUNCTION, typename ...TENSORS> 
    inline static void apply(const_tup_t shape, FUNCTION function, TENSORS & ...args) { 
    unsigned long counter[DIMENSION]; 
    memset(counter, 0, DIMENSION*sizeof(unsigned long)); 
    ForEachFixedDimensionHelper<DIMENSION,0>::template apply<FUNCTION, TENSORS...>(counter, shape, function, args...); 
    } 
};

還要注意元組值和形狀可以安全地聲明爲__restrict，這意味着它們指向不同的內存位置，因爲它們將專門構建用於迭代，然後解除分配。當另一個指針被解除引用和更改時（「指針別名」問題），由這些指針索引的值不需要從內存中重新讀取。當調用ForEachFixedDimension :: template apply時，可以在編譯時根據張量參數的內容推斷出typename FUNCTION（可能是一個lambda函數）和模板參數包typename ... TENSORS（variadic support）參數類型的功能。的展開for循環

所需的數量可以在運行時擡頭：

typedef unsigned int TEMPLATE_SEARCH_INT_TYPE; 
template <TEMPLATE_SEARCH_INT_TYPE MINIMUM, TEMPLATE_SEARCH_INT_TYPE MAXIMUM, template <TEMPLATE_SEARCH_INT_TYPE> class WORKER> 
class LinearTemplateSearch { 
public: 
    template <typename...ARG_TYPES> 
    inline static void apply(TEMPLATE_SEARCH_INT_TYPE v, ARG_TYPES && ... args) { 
    if (v == MINIMUM) 
     WORKER<MINIMUM>::apply(std::forward<ARG_TYPES>(args)...); 
    else 
     LinearTemplateSearch<MINIMUM+1, MAXIMUM, WORKER>::apply(v, std::forward<ARG_TYPES>(args)...); 
    } 
}; 

template <TEMPLATE_SEARCH_INT_TYPE MAXIMUM, template <TEMPLATE_SEARCH_INT_TYPE> class WORKER > 
class LinearTemplateSearch<MAXIMUM, MAXIMUM, WORKER> { 
public: 
    template <typename...ARG_TYPES> 
    inline static void apply(TEMPLATE_SEARCH_INT_TYPE v, ARG_TYPES && ... args) { 
    assert(v == MAXIMUM); 
    WORKER<MAXIMUM>::apply(std::forward<ARG_TYPES>(args)...); 
    } 
};

注意，在這裏，儘管模板遞歸時，尺寸需要不直到運行時是已知的。這基本上是通過使用模板作爲即時（JIT）編譯形式，爲所有感興趣的維度預先計算策略，然後在運行時查找正確的策略來實現的。

所以這些方法用Benchmark進行了測試。在基準1中，數據被從形狀的張量（2 ，2 ，）複製到形狀的張量（2 ，2 ，2 ）。在基準2，和形狀的兩個張量之間的內積（2 ，2 ，2 ）（2 ，2 ，2 ）被計算（僅來訪由兩者共享的元組索引）。將模板遞歸的實現與其他替代方法進行比較：元組迭代;元組迭代，其中維在編譯時已知;整數重新索引;整數重新索引，其中軸限制爲2的冪; numpy的; C風格for循環（硬編碼）;矢量化的Fortran代碼; Go中的循環。

事實證明模板遞歸比元組索引，並且該方法通過升壓使用更快：

灰色數字表示平均運行時和誤差棒的分和最大。這裏是本方法的那些實施了用於基準1爲每個方法：

// Tuple iteration (DIMENSION must be compile−time constant): vector<unsigned long> t(DIMENSION); 
t.fill(0); 
unsigned long k; 
for (k=0; k<x.flat.size(); advance_tuple_fixed_dimension<DIMENSION>(&t[0], &x.data_shape()[0]), ++k) 
    x[k] = y[tuple_to_index_fixed_dimension<DIMENSION>(&t[0], &y.data_shape()[0])]; 

// boost: 
x[ boost::indices[range(0, x.shape[0])][range(0,x.shape[1])][range(0,x.shape[2])] ] = y[ boost::indices[range(0,x.shape[0])][range(0,x.shape[1])][range(0,x.shape[2])] ]; 

! Fortran 95 
x = y(1:2**5,1:2**9,1:2**9) 

// Hard−coded for loops in C: unsigned long k; 
for (k=0; k<x.data_shape()[0]; ++k) { 
    for (unsigned long j=0; j<x.data_shape()[1]; ++j) { 
    unsigned long x_bias = (k*x.data_shape()[1] + j)*x.data_shape()[2]; 
    unsigned long y_bias = (k*y.data_shape()[1] + j)*y.data_shape()[2]; 
    for (unsigned long i=0; i<x.data_shape()[2]; ++i) 
     x[x_bias + i] = y[y_bias + i]; 
    } 
} 

// Integer reindexing: 
unsigned long k; 
for (k=0; k<x.flat.size(); ++k) 
    x[k] = y[reindex(k, &x.data_shape()[0], &y.data_shape()[0], DIMENSION)]; 

// Integer reindexing (axes are powers of 2): 
unsigned long k; 
for (k=0; k<x.flat.size(); ++k) 
    x[k] = y[reindex_powers_of_2(k, &x_log_shape[0], &y_log_shape[0], DIMENSION)]; 

// Tuple iteration (DIMENSION unknown at compile time): 
vector<unsigned long> t(DIMENSION); 
t.fill(0); 
unsigned long k; 
    for (k=0; k<x.flat_size(); advance_tuple(&t[0], &x.data_shape()[0], DIMENSION), ++k) 
    x[k] = y[t]; 

# numpy (python): 
x_sh = x.shape. 
x = np.array(y[:x_sh[0], :x_sh[1], :x_sh[2]]) 

// Go: 
for i:=0; i<1<<9; i++ { 
    for j:=0; j<1<<9; j++{ 
    for k:=0; k<1<<5; k++{ 
     x[i][j][k] = y[i][j][k] 
    } 
    } 
} 

// TRIOT (DIMENSION unknown at compile time): 
apply_tensors([](double & xV, double yV) { 
    xV = yV; 
}, 
x.data_shape(), 
x, y);

令人驚訝地，整數重新索引（即使軸是2的冪）基本上不是使一個元組計數器慢。與模板遞歸版本有時更快（包括比boost更快30％，即使boost :: multi_array必須知道編譯時的維數）。

這裏是你將如何使用這個嵌套的循環招用模板遞歸一個又一個例子：

double dot_product(const Tensor & x<double>, const Tensor<double> & y) { // This function written for homogeneous types, but not unnecessary 
    double tot = 0.0; 
    for_each_tensors([&tot](double xV, double yV) { 
    tot += xV * yV; 
    }, 
    x.data_shape(), /* Iterate over valid tuples for x.data_shape(); as written, this line assumes x has smaller shape*/ 
    x, y); 
    return tot; 
}

，並通過元組迭代多維卷積的實現，該版本與模板遞歸和numpy的進行了比較通過卷積兩個矩陣，每個矩陣具有形狀（2 ，2 ）。

Tensor<double> triot_naive_convolve(const Tensor<double> & lhs, const Tensor<double> & rhs) { 
    assert(lhs.dimension() == rhs.dimension()); 

    Tensor<double> result(lhs.data_shape() + rhs.data_shape() − 1ul); 
    result.flat().fill(0.0); 
    Vector<unsigned long> counter_result(result.dimension()); 

    enumerate_for_each_tensors([&counter_result, &result, &rhs](const_tup_t counter_lhs, const unsigned int dim_lhs, double lhs_val) { 
    enumerate_for_each_tensors([&counter_result, &result, &rhs, &counter_lhs, &lhs_val](const_tup_t counter_rhs, const unsigned int dim_rhs, double rhs_val) { 
     for (unsigned int i=0; i<dim_rhs; ++i) 
     counter_result[i] = counter_lhs[i] + counter_rhs[i]; 
     unsigned long result_flat = tuple_to_index(counter_result, result.data_shape(), dim_rhs); 
     result.flat()[result_flat] += lhs_val * rhs_val; 
    }, 
    rhs.data_shape(), 
    rhs); 
    }, 
    lhs.data_shape(), lhs); 
    return result; 
}

這些基準是定時的2.0 GHz的英特爾Core i7的芯片上的優化（-std = C++ 11 -Ofast -march =天然 - mtune中=天然-fomit幀指針）。所有Fortran實現都以相反順序使用軸，並以緩存優化方式訪問數據，因爲Fortran使用列主數組格式。可在此small journal article中找到詳細信息和源代碼（一個簡單的多維數組庫，其中在編譯時不需要知道維數）。

來源

2017-04-11 21:57:54 Flow

快速的方式來遍歷不同形狀的多維數組？

回答

相關問題