1
\$\begingroup\$

I would like to improve the naive matrix transpose algorithm by using a cache friendly algorithm.

So far there are 4 variants.

Iterate over the source matrix:

  MatrixXd transposeSrc() const
  {
    MatrixXd At(m_cols, m_rows);

    for (int i = 0; i < m_rows; i++) {
      for (int j = 0; j < m_cols; j++) {
        At(j, i) = (*this)(i, j);
      }
    }
    return At;
  }

Iterate over the destination matrix:

  MatrixXd transposeDst() const
  {
    MatrixXd At(m_cols, m_rows);

    for (int i = 0; i < m_cols; i++) {
      for (int j = 0; j < m_rows; j++) {
        At(i, j) = (*this)(j, i);
      }
    }
    return At;
  }

A "tiling" version found on StackOverflow:

  MatrixXd transposeTilingSO(int tileSize = 16) const
  {
    MatrixXd out(m_cols, m_rows);

    for (int i = 0; i < m_rows; i += tileSize) {
      for (int j = 0; j < m_cols; ++j) {
        for (int b = 0; b < tileSize && i + b < m_rows; ++b) {
          out.m_data[j*m_rows + i + b] = (*this).m_data[(i + b)*m_cols + j];
        }
      }
    }

    return out;
  }

A block tiling version:

  MatrixXd transposeTiling(int tileSize = 16) const
  {
    MatrixXd At(m_cols, m_rows);

    for (int i = 0; i < m_rows;) {
      for (; i <= m_rows - tileSize; i += tileSize) {
        int j = 0;
        for (; j <= m_cols - tileSize; j += tileSize) {
          for (int k = i; k < i + tileSize; k++) {
            for (int l = j; l < j + tileSize; l++) {
              At(l, k) = (*this)(k, l);
            }
          }
        }

        for (int k = i; k < i + tileSize; k++) {
          for (int l = j; l < m_cols; l++) {
            At(l, k) = (*this)(k, l);
          }
        }
      }

      for (; i < m_rows; i++) {
        for (int j = 0; j < m_cols; j++) {
          At(j, i) = (*this)(i, j);
        }
      }
    }

    return At;
  }

Experimentally, I thought I have found how to choose which variant is preferable depending on the input matrix size:

  MatrixXd transposeOptim(int tileSize = 16) const
  {
    if (m_rows > 2 * m_cols && m_cols <= 64) {
      return transposeSrc();
    } else if (m_cols > 2 * m_rows && m_rows <= 64) {
      return transposeDst();
    } else if (m_rows % tileSize == 0) {
      return transposeTilingSO();
    } else {
      return transposeTiling();
    }
  }

And the essential information about the MatrixXd class (full source code):

class MatrixXd {
public:
  MatrixXd() :
    m_data(), m_rows(0), m_cols(0) {}

  MatrixXd(int row, int col) :
    m_data(row*col), m_rows(row), m_cols(col) {}

  double& operator() (int row, int col)
  {
    return m_data[row*m_cols + col];
  }

  double operator() (int row, int col) const
  {
    return m_data[row*m_cols + col];
  }

  std::vector<double> m_data;
  int m_rows;
  int m_cols;
};

On one computer, I have found that:

  • for Nx6 it was preferable to iterate over the source matrix
  • for 6xN the opposite
  • when the first dimension was divisible by 16 (the tile size), the SO tiling variant was better
  • otherwise the block tiling was better

It is interesting for me to optimize Nx6 size that corresponds to problem that involves N samples and 6 variables.

On another computer, the results are a little bit different that makes me think I have probably overengineered the problem (block tiling seems to be the best overall).

Moreover, on quick-bench, I don't have the same conclusion (block tiling performs the worst).

Any idea how to have a fast matrix transpose that performs well on any matrix size and specifically for matrices of Nx6 size?

How to choose the optimal tile size? I have used 16 experimentally.

I target classical architecture that has 32 kB of L1-cache, 8-way associative and 64 bytes of cache line size.

Full source code and log are here. I use Catch2 for micro-benchmarking.

\$\endgroup\$
  • \$\begingroup\$ The "best performance" is always bound to the specific platform you're running code on, so you will have to benchmark. That said, there are some easy general optimizations you can do if N is known at compile time. Is that the case? Or can N only be known at runtime? \$\endgroup\$ – hoffmale Sep 30 at 20:15
  • \$\begingroup\$ N is not known at compile time. An example is a computer vision application where you extract some features from a video stream. I am interested to know what optimization could be done with for instance N=1000 or a large number. \$\endgroup\$ – Catree Oct 1 at 11:43
  • \$\begingroup\$ The problem is that this code is very platform/compiler/hardware dependent. Normally, there is no need to optimize matrix transposition lest initially a big mistake was made - as this is neither a bottle neck nor is a frequently used operation. At most, I'd consider adding support for MatrixExpression a class or rather a concept that behaves matrix-like. So instead of transposing a matrix you make a MatrixExpression with reference to the original matrix, and make its operator () swap the indexes (i,j) when you access data. \$\endgroup\$ – ALX23z Oct 12 at 6:36
  • \$\begingroup\$ @ALX23z I am not convince about using a MatrixExpression and swaping the indexes. At least for my use case where I can potentially need to iterate over the matrix. MatrixExpression à la Eigen for lazy evaluation is for sure a good solution. But I am more interest in the theoretical optimisation for this specific operation, to understand more about memory cache topic. \$\endgroup\$ – Catree Oct 26 at 23:51
  • \$\begingroup\$ @Catree speed of memory manipulations depends heavily on the cache lines. If you somehow manage to transform in into exchange of cache lines or close to it then the code should be faster. So, if matrix rows memory are 64bit aligned then it should be most efficient to use tiled-like version with parameter 4 or 8. But at long as it is not aligned you'll surely get random results. \$\endgroup\$ – ALX23z Oct 27 at 4:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.