在C++中最快的矩阵转置方法是什么？

Question

我有一个（相对较大的）矩阵需要转置。例如，假设我的矩阵是

a b c d e f
g h i j k l
m n o p q r 

我希望结果如下：

a g m
b h n
c I o
d j p
e k q
f l r

什么方法是最快的？

- mans

---

最快的方式不是旋转数组，而是在访问数组时仅仅交换索引顺序。 - High Performance Mark

无论多快，您都必须访问矩阵的所有元素。 - taocp

我猜这要看情况，如果你之后想要按行重复访问矩阵，有一个“转置”标志会对你造成很大影响。 - Matthieu M.

矩阵转置因其对内存缓存的影响而臭名昭著。如果您的数组足够大，以至于转置的性能非常重要，并且您无法通过提供具有交换索引的接口来避免转置，则最好的选择是使用现有的大型矩阵转置库例程。专家已经完成了这项工作，您应该使用它。 - Eric Postpischil

这个问题中有一些有用的信息（包括：使矩阵更大可以加快转置速度）。 - Eric Postpischil

Answer 1

这是一个好问题。实际上在内存中转置矩阵而不仅仅交换坐标有很多原因，例如在矩阵乘法和高斯模糊中。

首先让我列出我用于转置的函数之一（编辑：请参见我的答案末尾，我找到了一个更快的解决方案）。

void transpose(float *src, float *dst, const int N, const int M) {
    #pragma omp parallel for
    for(int n = 0; n<N*M; n++) {
        int i = n/N;
        int j = n%N;
        dst[n] = src[M*j + i];
    }
}

现在让我们看看转置为什么有用。考虑矩阵乘法 C = A*B。我们可以这样做。

for(int i=0; i<N; i++) {
    for(int j=0; j<K; j++) {
        float tmp = 0;
        for(int l=0; l<M; l++) {
            tmp += A[M*i+l]*B[K*l+j];
        }
        C[K*i + j] = tmp;
    }
}

然而，那种方式会有很多缓存未命中。一个更快的解决方案是先对B进行转置。

transpose(B);
for(int i=0; i<N; i++) {
    for(int j=0; j<K; j++) {
        float tmp = 0;
        for(int l=0; l<M; l++) {
            tmp += A[M*i+l]*B[K*j+l];
        }
        C[K*i + j] = tmp;
    }
}
transpose(B);

矩阵乘法的时间复杂度为O(n^3)，转置的时间复杂度为O(n^2)，因此对于大规模的n，进行转置应该对计算时间影响不大。在矩阵乘法中，循环分块比转置更有效，但这更加复杂。

我希望我知道一种更快的转置方法（编辑：我找到了一种更快的解决方案，请参见我的答案结尾）。当Haswell/AVX2在几周后推出时，它将具有gather函数。我不知道这在这种情况下是否有用，但我可以想象收集一列并写出一行。也许这会使转置变得不必要。

对于高斯模糊，您需要先进行水平模糊，然后进行垂直模糊。但是，垂直模糊存在缓存问题，因此您需要进行以下操作：

Smear image horizontally
transpose output 
Smear output horizontally
transpose output

这里有一篇由英特尔撰写的论文，解释了使用英特尔高级向量扩展实现IIR高斯模糊滤波器http://software.intel.com/en-us/articles/iir-gaussian-blur-filter-implementation-using-intel-advanced-vector-extensions。

最后，在矩阵乘法（以及高斯模糊）中，我实际上并不是精确地进行转置，而是按照某个向量大小的宽度进行转置（例如SSE / AVX的4或8）。 这是我使用的函数：

void reorder_matrix(const float* A, float* B, const int N, const int M, const int vec_size) {
    #pragma omp parallel for
    for(int n=0; n<M*N; n++) {
        int k = vec_size*(n/N/vec_size);
        int i = (n/vec_size)%N;
        int j = n%vec_size;
        B[n] = A[M*i + k + j];
    }
}

编辑：

我尝试了几个函数来寻找大矩阵的最快转置方法。最终，最快的结果是使用循环分块，block_size = 16（编辑：我找到了一种更快的解决方案，使用 SSE 和循环分块，请参见下文）。此代码适用于任何NxM矩阵（即矩阵不必为正方形）。

inline void transpose_scalar_block(float *A, float *B, const int lda, const int ldb, const int block_size) {
    #pragma omp parallel for
    for(int i=0; i<block_size; i++) {
        for(int j=0; j<block_size; j++) {
            B[j*ldb + i] = A[i*lda +j];
        }
    }
}

inline void transpose_block(float *A, float *B, const int n, const int m, const int lda, const int ldb, const int block_size) {
    #pragma omp parallel for
    for(int i=0; i<n; i+=block_size) {
        for(int j=0; j<m; j+=block_size) {
            transpose_scalar_block(&A[i*lda +j], &B[j*ldb + i], lda, ldb, block_size);
        }
    }
}

lda和ldb的值是矩阵的宽度。这些值需要是块大小的倍数。为了找到这些值并为例如3000x1001的矩阵分配内存，我会做如下操作：

#define ROUND_UP(x, s) (((x)+((s)-1)) & -(s))
const int n = 3000;
const int m = 1001;
int lda = ROUND_UP(m, 16);
int ldb = ROUND_UP(n, 16);

float *A = (float*)_mm_malloc(sizeof(float)*lda*ldb, 64);
float *B = (float*)_mm_malloc(sizeof(float)*lda*ldb, 64);

对于3000x1001，此函数返回ldb = 3008和lda = 1008 编辑： 我使用SSE内置函数找到了一个更快的解决方案：

inline void transpose4x4_SSE(float *A, float *B, const int lda, const int ldb) {
    __m128 row1 = _mm_load_ps(&A[0*lda]);
    __m128 row2 = _mm_load_ps(&A[1*lda]);
    __m128 row3 = _mm_load_ps(&A[2*lda]);
    __m128 row4 = _mm_load_ps(&A[3*lda]);
     _MM_TRANSPOSE4_PS(row1, row2, row3, row4);
     _mm_store_ps(&B[0*ldb], row1);
     _mm_store_ps(&B[1*ldb], row2);
     _mm_store_ps(&B[2*ldb], row3);
     _mm_store_ps(&B[3*ldb], row4);
}

inline void transpose_block_SSE4x4(float *A, float *B, const int n, const int m, const int lda, const int ldb ,const int block_size) {
    #pragma omp parallel for
    for(int i=0; i<n; i+=block_size) {
        for(int j=0; j<m; j+=block_size) {
            int max_i2 = i+block_size < n ? i + block_size : n;
            int max_j2 = j+block_size < m ? j + block_size : m;
            for(int i2=i; i2<max_i2; i2+=4) {
                for(int j2=j; j2<max_j2; j2+=4) {
                    transpose4x4_SSE(&A[i2*lda +j2], &B[j2*ldb + i2], lda, ldb);
                }
            }
        }
    }
}

- user2088790

---

不错的尝试，但我不确定“矩阵乘法是O(n^3)”这个说法是否正确，我认为它应该是O(n^2)。 - ulyssis2

@ulyssis2说：这是O(n^3)，除非你使用Strassen矩阵乘法（O(n^2.8074)）。 user2088790说：做得非常好，我会把它保存在我的个人收藏中。 :) - saurabheights

如果有人想知道谁写了这个答案，那就是我。我曾经离开过Stack Overflow（SO），但后来克服了困难又回来了。 - Z boson

朴素矩阵乘法的时间复杂度肯定是O(n^3)，据我所知，计算核实现了朴素算法（我认为这是因为Strassen算法需要执行更多的操作（加法），如果你可以进行快速的乘法，那么这种情况就不好，但我可能是错的）。目前有一个开放问题，即矩阵乘法是否可以达到O(n^2)。 - étale-cohomology

通常来说，依赖于线性代数库来完成工作是更好的选择。现代化的库例如Intel MKL、OpenBLAS等提供了动态CPU分派功能，可以为您的硬件选择最佳实现（例如，可能有比SSE更宽的向量寄存器：AVX AVX2、AVX512...），因此您不需要编写非可移植程序来获得快速程序。 - Jorge Bellon

请注意，如果行数和列数不是4的倍数，则最后一个SSE代码片段无法正常工作。它将保留边框单元格不变。 - Sopel

Answer 2

这将取决于您的应用程序，但通常转置矩阵的最快方式是在查找时反转坐标，然后您不必实际移动任何数据。

- Shafik Yaghmour

---

如果矩阵很小或者只需要读取一次，那么这样做是很好的。然而，如果转置后的矩阵很大并且需要被多次重复使用，你可能仍然需要保存一个快速转置版本以获得更好的内存访问模式。（顺便加1分） - Agentlien

@Agentlien：为什么A[j][i]会比A[i][j]慢？ - beaker

如果你有一个大矩阵，不同的行或列可能占据不同的缓存行或页面。在这种情况下，你需要以一种访问相邻元素的方式迭代元素。否则，每个元素访问都可能成为缓存未命中，这会完全破坏性能。 - Agentlien

@beaker: 这与 CPU 层面的缓存有关（假设矩阵是一大块连续内存），缓存行实际上是矩阵的行，预取器可以预取下几行。如果你改变访问模式，CPU 缓存/预取器仍然按行工作，而你却按列访问，性能会急剧下降。 - Matthieu M.

@taocp 基本上，您需要某种标志来指示它已经转置，然后请求 (i,j) 将被映射到 (j,i)。 - Shafik Yaghmour

Answer 3

关于使用x86硬件转置4x4方阵浮点数的一些细节（我稍后会讨论32位整数），这对于转置更大的方阵如8x8或16x16非常有帮助。

_MM_TRANSPOSE4_PS(r0, r1, r2, r3)在不同编译器中的实现方式不同。GCC和ICC（我还没有检查Clang）使用unpcklps，unpckhps，unpcklpd，unpckhpd，而MSVC仅使用shufps。我们可以将这两种方法组合在一起，像这样。

t0 = _mm_unpacklo_ps(r0, r1);
t1 = _mm_unpackhi_ps(r0, r1);
t2 = _mm_unpacklo_ps(r2, r3);
t3 = _mm_unpackhi_ps(r2, r3);

r0 = _mm_shuffle_ps(t0,t2, 0x44);
r1 = _mm_shuffle_ps(t0,t2, 0xEE);
r2 = _mm_shuffle_ps(t1,t3, 0x44);
r3 = _mm_shuffle_ps(t1,t3, 0xEE);

有一个有趣的观察是，两次洗牌可以转换为一次洗牌和两次混合（SSE4.1），就像这样。

t0 = _mm_unpacklo_ps(r0, r1);
t1 = _mm_unpackhi_ps(r0, r1);
t2 = _mm_unpacklo_ps(r2, r3);
t3 = _mm_unpackhi_ps(r2, r3);

v  = _mm_shuffle_ps(t0,t2, 0x4E);
r0 = _mm_blend_ps(t0,v, 0xC);
r1 = _mm_blend_ps(t2,v, 0x3);
v  = _mm_shuffle_ps(t1,t3, 0x4E);
r2 = _mm_blend_ps(t1,v, 0xC);
r3 = _mm_blend_ps(t3,v, 0x3);

这实际上将4个洗牌转换为2个洗牌和4个混合。这比GCC、ICC和MSVC的实现多使用了2个指令。优点是减少了端口压力，在某些情况下可能有好处。 目前所有的洗牌和解包都只能发送到一个特定的端口，而混合可以发送到两个不同的端口之一。

我尝试使用8个洗牌，就像MSVC那样，并将其转换为4个洗牌+8个混合，但没有成功。我仍然需要使用4个解包。

我在一个8x8浮点数转置中使用了相同的技术（请参阅答案末尾）。 https://dev59.com/cV8e5IYBdhLWcg3wja_5#25627536。在那个答案中，我仍然需要使用8个解包，但我设法将8个洗牌转换为4个洗牌和8个混合。

对于32位整数，除了AVX512的128位洗牌之外，没有像shufps这样的操作码，所以只能用解包来实现，但我认为无法高效地将其转换为混合。使用AVX512 vshufi32x4可以有效地执行与shufps类似的操作，但它是针对4个整数的128位通道而不是32位浮点数，因此在某些情况下可能可以使用这种技术。在骑士着陆芯片上，洗牌的吞吐量比混合慢四倍。

- Z boson

---

你可以在整数数据上使用 shufps。如果你要进行大量的洗牌操作，最好使用 FP 领域中的 shufps + blendps 进行全部操作，特别是当你没有同样高效的 AVX2 vpblendd 可用时。此外，在英特尔 SnB 系列硬件上，使用 shufps 在整数指令（如 paddd）之间没有额外的旁路延迟。（根据 Agner Fog 的 SnB 测试，使用 blendps 与 paddd 混合会有旁路延迟。） - Peter Cordes

@PeterCordes，我需要再次审查域更改。是否有某个表格（也许在SO上的答案）可以总结Core2-Skylake的域更改惩罚？无论如何，我已经更加深入地思考了这个问题。现在我明白为什么wim和你一直提到我的16x16转置答案中的vinsertf64x4而不是vinserti64x4。如果我读取然后写入矩阵，那么使用浮点域或整数域确实没有关系，因为转置只是移动数据。 - Z boson

Agner的表格列出了Core2和Nehalem（以及我认为是AMD）每个指令的域，但不包括SnB系列。 Agner的微架构指南只有一个段落，说SnB的速度降至1c，并且通常为0，带有一些示例。 Intel的优化手册有一个表格，但我没有尝试理解它，所以我不记得它有多少细节。 我记得某些情况下并不明显给定指令属于哪个类别。 - Peter Cordes

即使您不仅仅是写回内存，整个转置也只需要1个额外的时钟。每个操作数的额外延迟可以并行发生（或交错方式），因为转置的使用者开始读取由混洗或混合写入的寄存器。乱序执行允许前几个FMA或其他操作开始，而最后几个混洗正在完成，但没有dypass延迟链，只有最多一个额外的延迟。 - Peter Cordes

不错的回答！英特尔64-ia-32架构优化手册，表2-3列出了Skylake的旁路延迟，也许这对您有兴趣。Haswell的表2-8看起来完全不同。 - wim

我认为在Skylake上，vinsertf64x4和vinserti64x4是可以互换的。我没有理由提到其中一个或另一个。我只是在考虑64x4位数据。 - wim

Answer 4

如果数组的大小是已知的，那么我们可以使用union来帮助我们。像这样-

#include <bits/stdc++.h>
using namespace std;

union ua{
    int arr[2][3];
    int brr[3][2];
};

int main() {
    union ua uav;
    int karr[2][3] = {{1,2,3},{4,5,6}};
    memcpy(uav.arr,karr,sizeof(karr));
    for (int i=0;i<3;i++)
    {
        for (int j=0;j<2;j++)
            cout<<uav.brr[i][j]<<" ";
        cout<<'\n';
    }

    return 0;
}

- Sandeep K V

---

我对C/C++还很陌生，但这看起来很厉害。 因为union使用共享内存位置来存储其成员，所以可以以不同的方式读取该内存。因此，您可以获得一个转置矩阵，而无需进行新的数组分配。我是正确的吗？ - Doğuş

我认为这不正确。这只是以不同的行大小打印元素。转置需要交换行和列。@Doğuş所指的可以通过主帖子评论中描述的方式实现，“只需在访问数组时交换索引顺序即可”。 - jezza

Answer 5

将每一行视为一列，每一列视为一行..使用j,i而不是i,j。

示例：http://ideone.com/lvsxKZ

#include <iostream> 
using namespace std;

int main ()
{
    char A [3][3] =
    {
        { 'a', 'b', 'c' },
        { 'd', 'e', 'f' },
        { 'g', 'h', 'i' }
    };

    cout << "A = " << endl << endl;

    // print matrix A
    for (int i=0; i<3; i++)
    {
        for (int j=0; j<3; j++) cout << A[i][j];
        cout << endl;
    }

    cout << endl << "A transpose = " << endl << endl;

    // print A transpose
    for (int i=0; i<3; i++)
    {
        for (int j=0; j<3; j++) cout << A[j][i];
        cout << endl;
    }

    return 0;
}

- Khaled.K

Answer 6

无额外开销的转置（类未完成）：

class Matrix{
   double *data; //suppose this will point to data
   double _get1(int i, int j){return data[i*M+j];} //used to access normally
   double _get2(int i, int j){return data[j*N+i];} //used when transposed

   public:
   int M, N; //dimensions
   double (*get_p)(int, int); //functor to access elements  
   Matrix(int _M,int _N):M(_M), N(_N){
     //allocate data
     get_p=&Matrix::_get1; // initialised with normal access 
     }

   double get(int i, int j){
     //there should be a way to directly use get_p to call. but i think even this
     //doesnt incur overhead because it is inline and the compiler should be intelligent
     //enough to remove the extra call
     return (this->*get_p)(i,j);
    }
   void transpose(){ //twice transpose gives the original
     if(get_p==&Matrix::get1) get_p=&Matrix::_get2;
     else get_p==&Matrix::_get1; 
     swap(M,N);
     }
}

can be used like this:

Matrix M(100,200);
double x=M.get(17,45);
M.transpose();
x=M.get(17,45); // = original M(45,17)

当然，我在这里没有去关注内存管理，这是至关重要但不同的话题。

- Reza Baram

---

每次访问元素都需要遵循函数指针的额外开销。 - user877329

Answer 7

template <class T>
void transpose( const std::vector< std::vector<T> > & a,
std::vector< std::vector<T> > & b,
int width, int height)
{
    for (int i = 0; i < width; i++)
    {
        for (int j = 0; j < height; j++)
        {
            b[j][i] = a[i][j];
        }
    }
} 

- Rachel Gallen

---

我认为如果你交换两个循环，速度会更快，因为写入时缓存未命中的惩罚比读取时要小。 - phoeagon

这仅适用于方阵。矩形矩阵是一个完全不同的问题！ - NealB

问题要求寻找最快的方法，但这只是其中一种方式。你怎么能认为它快呢，更不用说最快了？对于大型矩阵来说，这会使缓存失效，导致性能极差。 - Eric Postpischil

@NealB：你是怎么想到那个的？ - Eric Postpischil

@EricPostpischil OP询问的是一个相对较大的矩阵，因此我认为他们想要“原地”操作以避免分配双倍的内存。当这样做时，源矩阵和目标矩阵的基地址相同。通过翻转行和列索引来转置只适用于方阵。有一些方法可以使矩形矩阵正确转置，但它们会更加复杂。 - NealB

Answer 8

Intel MKL建议使用就地和非就地转置/复制矩阵。这里是文档链接。我建议尝试使用非就地实现，因为它比就地实现更快，并且最新版本的MKL文档中包含一些错误。

- Gennady.F

Answer 9

现代线性代数库包括最常见操作的优化版本。其中许多包括动态CPU分发，可以在程序执行时选择最佳实现（不会影响可移植性）。

这通常比通过向量扩展内部函数手动优化您的函数更好。后者会将您的实现与特定的硬件供应商和型号绑定：如果您决定换成不同的供应商（例如Power、ARM）或更新的向量扩展（例如AVX512），则需要重新实现以获得最大效益。

例如，MKL转置包括BLAS扩展函数imatcopy。你也可以在其他实现中找到它，比如OpenBLAS：

#include <mkl.h>

void transpose( float* a, int n, int m ) {
    const char row_major = 'R';
    const char transpose = 'T';
    const float alpha = 1.0f;
    mkl_simatcopy (row_major, transpose, n, m, alpha, a, n, n);
}

对于一个 C++ 项目，你可以使用 Armadillo C++ 库：

#include <armadillo>

void transpose( arma::mat &matrix ) {
    arma::inplace_trans(matrix);
}

- Jorge Bellon

Answer 10

最快的转置是那个将在下一次操作中保留在缓存中的转置。

例如，不要一次性全部转置。只转置一个子矩阵。然后在需要转置数据的下一个算法的一部分中使用它。然后转置下一个子矩阵。然后计算。然后再转置另一个子矩阵。重复此过程，直到整个矩阵被转置。这样，数据会保持在缓存中。

如果你一次性完全转置一个128MB的矩阵，在一个2MB缓存的CPU上，那么在操作结束时，只有矩阵的最新位于缓存中。然后你最好从最新位开始乘以矩阵，以使用那2MB的热数据。

但是当你将工作分成较小的块，比如用子矩阵进行乘法运算时，你可以简单地进行懒惰转置，就像这样：

multiply:
  for all sub_matrices in mat1 row
  for all sub_matrices in mat2 column
    select sub_matrix1
    select sub_matrix2
    if sub_mat2 is not transposed
        transpose sub_mat2
    multiply sub_mat1 and sub_mat2 <---- data in cache!
    accumulate result

优势：

• 使用L1/L2缓存带宽进行下一步操作
• 转置延迟被隐藏在下一步操作的计算后面
• 适用于小缓存，最低可达64kB，取决于块大小

- huseyin tugrul buyukisik