Feature request: block matrix view

[Mikhail Katliar]

In some cases it is desirable to have a block view of a numeric matrix, and at the same time to be able to use it as a normal numeric matrix. With the approach proposed here, inconsistent block sizes are possible. For example, the following code will run:

    using namespace blaze;

    using DynM = DynamicMatrix<double>;
    DynamicMatrix<DynM> A{ { DynM(1, 2), DynM(3, 4) },
                           { DynM(5, 6), DynM(7, 8) } };

A can be seen as a matrix of matrices, but not as a matrix of numbers. The main problem with this is that the assignment to a numeric matrix will not work:

    DynamicMatrix<int> A1;
    A1 = A; // compiler error

I propose introducing a BlockMatrixView<> class. The BlockMatrixView<> ctor would take 3 arguments: a reference to a matrix m, and two integer arrays si and sj defining the row and column sizes of the blocks. The operator()(size_t i, size_t j) would return a Submatrix representing the block of m with si[i] rows and sj[j] columns.

Another possible implementation would be a function returning a matrix of submatrices:

template <typename MT, typename SI, typename SJ>
DynamicMatrix<Submatrix<MT>> blockView(Matrix<MT>& m, SI const& row_size, SJ const& col_size);

Related issues: #127

[Klaus Iglberger]

Hi Mikhail!

Thanks a lot for the proposal. You correctly summarize the shortcomings of the current block matrix approach: It is possible to set up blocks with inconsistent dimensions (as each block is an independent matrix object) and it is not possible to assign a matrix of matrices to a matrix of numbers (even if the setup of the block matrix would fit).

Both of your proposed solutions sound interesting, but are unfortunately also not without disadvantages: It is not possible to exploit knowledge about the matrix structure (i.e. mix sparse and dense matrix types) and the block matrix is not a native type, but only a view (a proxy object). However, we understand that depending on the context of the application your solution can be better suited. Therefore a choice between both approaches would be desirable. Therefore we are thinking about how to properly integrate this into Blaze.

Until we integrate a solution, the following blockMatrixToNumberMatrix () conversion function might be of use. It properly performs the conversion from a matrix of blocks to a matrix of numbers under the assumption that the block dimensions are consistent (i.e. it doesn't check the consistency of the block dimensions of the given block matrix):

template< typename MT1  // Type of the left-hand side matrix
        , bool SO1      // Storage order of the left-hand side matrix
        , typename MT2  // Type of the right-hand side matrix
        , bool SO2 >    // Storage order of the right-hand side matrix
void blockMatrixToNumberMatrix( const Matrix<MT1,SO1>& blocks, Matrix<MT2,SO2>& numbers )
{
   const size_t M( (~blocks).rows() );
   const size_t N( (~blocks).columns() );

   if( M == 0UL || N == 0UL ) {
      clear( ~numbers );
      return;
   }

   size_t totalRows( 0UL );
   for( size_t i=0UL; i<M; ++i ) {
      totalRows += (~blocks)(i,0UL).rows();
   }

   size_t totalColumns( 0UL );
   for( size_t j=0UL; j<N; ++j ) {
      totalColumns += (~blocks)(0UL,j).columns();
   }

   if( totalRows != numbers.rows() || totalColumns != numbers.columns() ) {
      BLAZE_THROW_INVALID_ARGUMENT( "Invalid conversion between block matrix and number matrix" );
   }

   size_t rowOffset( 0UL );

   for( size_t i=0UL; i<M; ++i )
   {
      size_t columnOffset( 0UL );

      for( size_t j=0UL; j<N; ++j )
      {
         auto& block( (~blocks)(i,j) );
         submatrix( ~numbers, rowOffset, columnOffset, block.rows(), block.columns() ) = block;
         columnOffset += block.columns();
      }

      rowOffset += (~blocks)(i,0UL).rows();
   }
}

Thanks again for the proposal,

Best regards,

Klaus!

[Mikhail Katliar]

Thank you for the prompt reply, Klaus!

In my case, I want to avoid copying elements from a block matrix to a normal matrix. One reason is that my software can be run on embedded platforms, where memory is often a scarce resource; the other reason is avoiding unnecessary copying. I still want to have a block-view of the matrix, because it simplifies the algorithm a lot. My current solution is implementing my own BlockMatrixView class.

Best regards,

Mikhail

[Emil Fresk]

This would be an amazing feature indeed! I run into this problem a lot when doing Kalman Filter implementations, especially on embedded systems where all performance matters and most of my matrices have something like this form:

Where being able to define a matrix based on the block structure would do magic to run time performance (and I manually write out the block operations now).

Is there anything in the TODO list for Blaze that would add this?

BR Emil

[Klaus Iglberger]

Hi Emil!

Since Mikhail has made a convincing argument and since we also see the advantages of such a view we are indeed planning to implement it. We cannot give you a release date, though. The only thing I can tell you is that we are planning to implement it for Blaze 3.4 or 3.5.

Till the feature is available you might consider to use a block-structured matrix. Since you have a lot of zeros in the matrix, I believe that it will work pretty well.

Best regards,

Klaus!

[Emil Fresk]

Nice, thanks for taking it under consideration! I have been poking Gael about it (Eigen) for a while but there was very little interest, which gives me a lot of hope for Blaze. :)

When it comes to the block-structured matrix, sorry for my ignorance but, does it optimize away matrices only containing zeros? I could not find a reference for this in the documentation.

BR Emil

[Klaus Iglberger]

Hi Emil!

If you consider the second example in the wiki, then all matrices that are not explicitly added to the sparse matrix are not stored and implicitly considered to be zero matrices. Only those matrices that you explicitly add are stored. Thus the matrices that only contain zeros are effectively "optimized away". I hope this helps,

Best regards,

Klaus!

[Emil Fresk]

Hi Klaus!

Thank you, however sadly it does not apply for me. The CompressedMatrix does dynamic allocation which is sadly forbidden in embedded systems, this is why I am looking for a way to do this at compile time. :)

I have envisioned something like this (functionality wise, not syntax wise):

using mySparseMatrix = StaticCompressedMatrix<
    BlockSize<3, 3>, // block size (rows/cols)
    NumBlocks<5, 5>, // number of blocks (rows/cols) - now in total a sparse 15x15 matrix
    BlockPos<0, 0>,  // list of existing blocks
    BlockPos<1, 1>, 
    BlockPos<2, 2>, 
    BlockPos<3, 4>,
    BlockPos<4, 4>,
    BlockPos<0, 2>,   
  >;

// Access example
mySparseMatrix m{};
m.atPos<0, 0>() = { { 1, 2, 3 }, { 4, 5, 6 }, { 7, 8, 9 } };

If something like this existed, code generators such as CVXGen would not need to exist as we could do the same in the type system of C++.

Thank you for the help though! For "normal" systems the CompressedMatrix can clearly be used!

BR Emil

Edit: Moreover, if something like this existed, optimized kalman filters could easily be generated by scripts when only needing to know the block structure of the system. I see a lot of possibilities with a feature like this. :)

[Klaus Iglberger]

Issue #186 was marked as a duplicate of this issue.