Provide support for the Kronecker product

[Klaus Iglberger]

Description

The Blaze library should provide support for the computation of the Kronecker product. More specifically, it should be possible to ...

• ... compute the Kronecker product between two dense matrices;
• ... compute the Kronecker product between a dense and a sparse matrix;
• ... compute the Kronecker product between a sparse and a dense matrix;
• ... compute the Kronecker product between two sparse matrices.

Conceptually, the computation of the Kronecker product should work as follows:

Conceptual Example

// Creating the matrix
//    ( 1 2 3 )
//    ( 4 1 2 )
blaze::StaticMatrix<int,2UL,3UL> A{ { 1, 2, 3 }
                                  , { 4, 1, 2 } };

// Creating the matrix
//    (  1 -2 )
//    ( -3  1 )
blaze::StaticMatrix<int,2UL,2UL> B{ {  1, -2 }
                                  , { -3,  1 } };

// Evaluating the Kronecker product
blaze::StaticMatrix<int,4UL,6UL> C;
C = kronecker( A, B );  // Results in (   1 -2  2 -4  3 -6 )
                        //            (  -3  1 -6  2 -9  3 )
                        //            (   4 -8  1 -2  2 -4 )
                        //            ( -12  4 -3  1 -6  2 )

Tasks

• implement the Kronecker product between two dense matrices
• implement the Kronecker product between a dense and a sparse matrix
• implement the Kronecker product between a sparse and a dense matrix
• implement the Kronecker product between two sparse matrices
• provide a full documentation of the Kronecker product
• ensure compatibility with all existing matrix classes
• ensure compatibility with all existing matrix expressions
• guarantee maximum performance for the operations
• add test cases for the entire Kronecker product functionality

[Klaus Iglberger]

The following function implements the Kronecker product for both dense and sparse matrices:

template <typename MT1  // Type of the left-hand side dense matrix operand
         ,bool SO1      // Storage order of the left-hand side dense matrix operand
         ,typename MT2  // Type of the right-hand side dense matrix operand
         ,bool SO2      // Storage order of the right-hand side dense matrix operand
         ,typename MT3  // Type of the target dense matrix
         ,bool SO3>     // Storage order of the target dense matrix
void kronecker( const blaze::Matrix<MT1,SO1>& left
              , const blaze::Matrix<MT2,SO2>& right
              , blaze::Matrix<MT3,SO3>& out)
{
  using namespace blaze;

  using Left  = If_t< IsMatMatMultExpr<MT1>::value, const typename MT1::ResultType, const MT1& >;
  using Right = If_t< IsExpression<MT2>::value, const typename MT2::ResultType, const MT2& >;

  Left  A( ~left  );
  Right B( ~right );

  for (std::size_t i = 0; i < rows(A); ++i) {
    for (std::size_t j = 0; j < columns(A); ++j) {
      const std::size_t row_lower = i * rows(B);
      const std::size_t col_lower = j * columns(B);
      auto block = submatrix(out, row_lower, col_lower,
                             rows(B), columns(B));
      (~block) = A(i, j) * B;
    }
  }
}

Please note that this function provides the required functionality, but is not optimized for speed! Also, the signature is different from the desired signature of a kronecker() function.

Please feel free to use this function until the functionality has been integrated into Blaze and is officially supported.

[Klaus Iglberger]

Description

The Blaze library should provide support for the computation of the Kronecker product. More specifically, it should be possible to ...

• ... compute the Kronecker product between two dense matrices;
• ... compute the Kronecker product between a dense and a sparse matrix;
• ... compute the Kronecker product between a sparse and a dense matrix;
• ... compute the Kronecker product between two sparse matrices.

Conceptually, the computation of the Kronecker product should work as follows:

Conceptual Example

// Creating the matrix
//    ( 1 2 3 )
//    ( 4 1 2 )
blaze::StaticMatrix<int,2UL,3UL> A{ { 1, 2, 3 }
                                  , { 4, 1, 2 } };

// Creating the matrix
//    (  1 -2 )
//    ( -3  1 )
blaze::StaticMatrix<int,2UL,2UL> B{ {  1, -2 }
                                  , { -3,  1 } };

// Evaluating the Kronecker product
blaze::StaticMatrix<int,4UL,6UL> C;
C = kronecker( A, B );  // Results in (   1 -2  2 -4  3 -6 )
                        //            (  -3  1 -6  2 -9  3 )
                        //            (   4 -8  1 -2  2 -4 )
                        //            ( -12  4 -3  1 -6  2 )

Tasks

• implement the Kronecker product between two dense matrices
• implement the Kronecker product between a dense and a sparse matrix
• implement the Kronecker product between a sparse and a dense matrix
• implement the Kronecker product between two sparse matrices
• provide a full documentation of the Kronecker product
• ensure compatibility with all existing matrix classes
• ensure compatibility with all existing matrix expressions
• guarantee maximum performance for the operations
• add test cases for the entire Kronecker product functionality

[Klaus Iglberger]

Summary

The feature has been implemented, tested, optimized, and documented as required. It is immediately available via cloning the Blaze repository and will be officially released in Blaze 3.6.

Kronecker Product

Vector Kronecker Product

The Kronecker product of two vectors with the same transpose flag can be computed via the kron() function:

using blaze::DynamicVector;
using blaze::CompressedVector;

DynamicVector<double>   v1( 28UL );
CompressedVector<float> v2( 17UL );

// ... Initialization of the vectors

CompressedVector<double> v3 = kron( v1, v2 );

Both dense and sparse vectors can be used for a Kronecker product. It is possible to multiply two vectors with different element type, as long as the element types themselves can be multiplied.

Matrix Kronecker Product

The Kronecker product of two matrices can be computed via the kron() function:

using blaze::DynamicMatrix;
using blaze::CompressedMatrix;

DynamicMatrix<double>   M1( 28UL, 35UL );
CompressedMatrix<float> M2( 17UL, 11UL );

// ... Initialization of the matrices

CompressedMatrix<double> M3 = kron( M1, M2 );

Both dense and sparse matrices can be used for a Kronecker product. The storage order of the two matrices poses no restrictions on the operation, all variations are possible. It is also possible to multiply two matrices with different element type, as long as the element types themselves can be multiplied.