Matrix multiplication optimization

Question

I am performing a series of matrix multiplications with fairly large matrices. To run through all of these operations takes a long time, and I need my program to do this in a large loop. I was wondering if anyone has any ideas to speed this up? I just started using Eigen, so I have very limited knowledge.

I was using ROOT-cern's built in TMatrix class, but the speed for performing the matrix operations is very poor. I set up some diagonal matrices using Eigen with the hope that it handled the multiplication operation in a more optimal way. It may, but I cannot really see the performance difference.

// setup matrices
int size = 8000;

Eigen::MatrixXf a(size*2,size);

// fill matrix a....

Eigen::MatrixXf r(2*size,2*size); // diagonal matrix of row sums of a

// fill matrix r

Eigen::MatrixXf c(size,size); // diagonal matrix of col sums of a

// fill matrix c

// transpose a in place
a.transposeInPlace();

Eigen::MatrixXf c_dia;
c_dia = c.diagonal().asDiagonal();

Eigen::MatrixXf r_dia;
r_dia = r.diagonal().asDiagonal();

// calc car
Eigen::MatrixXf car;
car = c_dia*a*r_dia;

Is something like this faster? Eigen::MatrixXf car = ((a.transpose().array().rowwise() * a.colwise().sum()).colwise() * a.rowwise().sum()).matrix(). — jdehesa
You honestly expect multiplying an 8000x8000 matrix by a 16000x8000 to be fast? Unless you know something about the structure of the matrix, and that something can be exploited, your chances of getting decent speed will be remote. — Peter
I did not expect it to be fast. I was hoping that I was overlooking something incredibly easy and obvious. It looks like I was. Thank you for the comment! — j_natt

chtz chtz · Accepted Answer · 2019-05-15T13:31:33

You are doing way too much work here. If you have diagonal matrices, only store the diagonal (and directly use that for products). Once you store a diagonal matrix in a square matrix, the information of the structure is lost to Eigen.

Also, you don't need to store the transposed variant of a, just use a.transpose() inside a product (that is only a minor issue here ...)

// setup matrices
int size = 8000;

Eigen::MatrixXf a(size*2,size);

// fill matrix a....
a.setRandom();

Eigen::VectorXf r = a.rowwise().sum(); // diagonal matrix of row sums of a
Eigen::VectorXf c = a.colwise().sum(); // diagonal matrix of col sums of a

Eigen::MatrixXf car = c.asDiagonal() * a.transpose() * r.asDiagonal();

Finally, of course make sure to compile with optimization enabled, and enable vectorization if available (with gcc or clang compile with -O2 -march=native).

Matrix multiplication optimization

1 Answers