85

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

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

最快的方法是什么？

4

11 回答 11

141

这是一个很好的问题。您希望在内存中实际转置矩阵而不只是交换坐标的原因有很多，例如在矩阵乘法和高斯拖尾中。

首先让我列出我用于转置的一个函数（编辑：请参阅我的答案的结尾，在那里我找到了一个更快的解决方案）

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 几周后问世时，它将具有收集功能。我不知道这在这种情况下是否有帮助，但我可以想象收集一列并写出一行。也许它会使转置变得不必要。

对于高斯涂抹，您所做的是水平涂抹，然后垂直涂抹。但是垂直涂抹有缓存问题，所以你要做的是

Smear image horizontally
transpose output 
Smear output horizontally
transpose output

这是英特尔的一篇论文，解释了 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);
                }
            }
        }
    }
}
于 2013-05-24T20:33:02.150 回答
41

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

于 2013-05-24T14:28:24.697 回答
5

关于用 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);

一个有趣的观察是，两个 shuffle 可以像这样转换为一个 shuffle 和两个混合（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 条指令。优点是它降低了端口压力，这在某些情况下可能有好处。目前，所有的 shuffle 和 unpacks 只能去一个特定的端口，而 blends 可以去两个不同的端口中的任何一个。

我尝试使用像 MSVC 这样的 8 次随机播放并将其转换为 4 次随机播放 + 8 次混合，但它不起作用。我仍然不得不使用 4 个解包。

我对 8x8 浮点转置使用了相同的技术（请参阅该答案的末尾）。 https://stackoverflow.com/a/25627536/2542702。在那个答案中，我仍然不得不使用 8 次解包，但我设法将 8 次随机播放转换为 4 次随机播放和 8 次混合。

对于 32 位整数，没有什么shufps比这更好的了（除了 128 位的 AVX512 洗牌），所以它只能通过我认为不能转换为混合（有效）的解包来实现。使用 AVX512 时，除了 4 个整数的 128 位通道而不是 32 位浮点数之外，AVX512 的vshufi32x4行为有效，因此在某些情况下shufps可能使用相同的技术。vshufi32x4对于 Knights Landing，shuffle 的速度（吞吐量）是混合的四倍。

于 2016-12-28T16:00:29.260 回答
4

如果事先知道数组的大小，那么我们可以使用联合来帮助我们。像这样-

#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;
}
于 2019-08-01T03:04:15.580 回答
1

将每一行视为一列，将每一列视为一行..使用 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;
}
于 2013-05-25T05:28:56.073 回答
1

没有任何开销的转置（类不完整）：

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);
     }
}

可以这样使用：

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

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

于 2014-03-17T20:55:39.320 回答
0 
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];
        }
    }
}
于 2013-05-24T14:29:52.487 回答
0

现代线性代数库包括最常见操作的优化版本。其中许多包括动态 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);
}
于 2019-09-21T10:28:06.800 回答
0

英特尔 mkl 建议就地和就地转置/复制矩阵。这是文档的链接。我会建议尝试以更快的 10 就地实现不适当的实现，并且在最新版本的 mkl 的文档中包含一些错误。

于 2019-10-18T15:39:57.580 回答
-1

我认为最快速的方法不应该高于 O(n^2)，这样你也可以只使用 O(1) 空间：
这样做的方法是成对交换，因为当你转置矩阵时，你做的是： M[i][j]=M[j][i] ，所以将 M[i][j] 存储在 temp 中，然后 M[i][j]=M[j][i] 和最后一步：M[j][i]=temp。这可以通过一次完成，所以它应该需要 O(n^2)

于 2013-05-29T10:02:51.737 回答
-6

我的答案是 3x3 矩阵的转置

 #include<iostream.h>

#include<math.h>


main()
{
int a[3][3];
int b[3];
cout<<"You must give us an array 3x3 and then we will give you Transposed it "<<endl;
for(int i=0;i<3;i++)
{
    for(int j=0;j<3;j++)
{
cout<<"Enter a["<<i<<"]["<<j<<"]: ";

cin>>a[i][j];

}

}
cout<<"Matrix you entered is :"<<endl;

 for (int e = 0 ; e < 3 ; e++ )

{
    for ( int f = 0 ; f < 3 ; f++ )

        cout << a[e][f] << "\t";


    cout << endl;

    }

 cout<<"\nTransposed of matrix you entered is :"<<endl;
 for (int c = 0 ; c < 3 ; c++ )
{
    for ( int d = 0 ; d < 3 ; d++ )
        cout << a[d][c] << "\t";

    cout << endl;
    }

return 0;
}
于 2013-12-25T06:50:03.707 回答