Eigen matrix with block diagonal structure

Question

This is a query regarding the c++ Eigen library.

In my calculation, I frequently come across matrices with a block diagonal form. A simple example of the structure is below.

Matrix6f M = Matrix6f::Zero(); //< M has a special structure
M.block<3, 3,>(0, 0) = Matrix3f::Random();
M.block<3, 3,>(3, 3) = Matrix3f::Random();

Matrix6f out = M * Matrix6f::Random(); //< Want to optimize this operation

Now, lets say I have a dense 6x6 matrix I and need to evaluate the product M*I, is there any way I can avoid having to perform redundant multiplications/additions with all the zeros that are there in M. Short of writing out the multiplication explicitly that is.

Eigen has this neat feature where if I do an auto u = v + w, I just get an operator without an actual evaluation. This also appears to work with auto z = Matrix6f::Zero() and if I ever did a Matrix6f out = z + Matrix6f::Random(), I could imagine writing an operator on the special Zero matrix so that I effectively just return the elements of the random matrix without attempting to add a zero to it.

Is something like this possible for the specific example I have list above?

If M is sufficiently large, you could gain by using a SparseMatrix instead of a normal one. — RHertel
How large is your matrix in practice (approximately). Are the sub-blocks all of the same size/is their size known at compile time? — chtz

asherikov asherikov · Accepted Answer · 2019-12-26T17:30:10

I did some experiments with block-wise sparse matrices implemented on top of dense Eigen matrices, the code is incomplete and not maintaned, but you can take a look at blockmatrix* headers here -> https://github.com/asherikov/eigenut/tree/master/include/eigenut. This allows you to do stuff like this:

eigenut::DiagonalBlockMatrix<3,3> diagonal;                
diagonal.setZero(2);                                      
diagonal(0) = Eigen::Matrix3d::Random();                  
diagonal(1) = Eigen::Matrix3d::Random();                

Eigen::MatrixXd matrix;                                        
matrix.setZero(6,6);                                           
matrix.diagonal().head(3).setRandom();                             
matrix.diagonal().tail(3).setConstant(1);              

Eigen::MatrixXd result1;                                                  
diagonal.multiplyLeft(result1, matrix);
std::cout << result1 << std::endl;                                 

eigenut::GenericBlockMatrix<3,3> generic;                              
generic.setZero(2,2);                                     
generic(1,0).setIdentity();                               
generic(0,1) = Eigen::Matrix3d::Random();

eigenut::GenericBlockMatrix<3,3> result2 = generic*diagonal;
std::cout << diagonal.getRaw() << std::endl;
std::cout << generic.getRaw() << std::endl;
std::cout << result2.getRaw() << std::endl;

Performance gain was substantial in my application.

Eigen matrix with block diagonal structure

2 Answers