Hi,

I need to compute kronecker product matrix. so, I coded like as follows.
however, its speed is quite slow. it consumes around 25000 ms (size of matrix M is 69756 x 6319 (sparse) and the value of nonzeros() is 69756, G: 4x4 sparse matrix).
could you let me know why ? is there a problem or a mistake ?
when I tested this in Matlab, it took around 800 ms for 48760 non-zeros. it's weird.

Code: Select all

typedef Eigen::SparseMatrix<float> SpMat;
   SpMat M(triplet_vec.size(), NV);      // 69756 x 6319
   M.setFromTriplets(triplet_vec.begin(), triplet_vec.end());
   M.finalize();
   M.makeCompressed();

        SpMat G(4,4);   // identity matrix ( G(i,i)=1 )

   SpMat MG(M.rows()*G.rows(), M.cols()*G.cols()) ;
   Eigen::kroneckerProduct(M, G, MG);


if you have any idea to speed-up the code, any comments are welcome.

Thanks in advance..
Cheers!
Yeonchool.

nobody knows ?

What if you change line 56 of unsupported/Eigen/src/KroneckerProduct/KroneckerTensorProduct.h to:

AB.reserve(VectorXi::Constant(AB.outerSize(), A.nonZeros()*B.nonZeros()/AB.outerSize()));

Thanks ggael,

I already tested "reserve" function. however its improvement is not sufficient to me.
it still consumes 8 sec for the same routine even though it is 3 times faster than previous method, from 25 sec.

so, I changed my plan to compute Kronecker product. I used Eigen's Triplet for multiplication.
after all, its performance is amazing ! it took only 20 ms. it is 400 times faster.

the code is as follows.

Code: Select all

   std::vector<Eigen::Triplet<float>> triplet_vec;   // sparse triplet vectors
   std::vector<Eigen::Triplet<float>> MGv;
        Eigen::MatrixXf G(4,4);
        G.setIdentity();    // identity matrix

        for (size_t i=0; i<triplet_vec.size(); ++i)
   {
      size_t ri = triplet_vec[i].row()*4;
      size_t ci = triplet_vec[i].col()*4;
      MGv.push_back(T(ri, ci, triplet_vec[i].value()*G(0,0) ));
      MGv.push_back(T(ri+1, ci+1, triplet_vec[i].value()*G(1,1) ));
      MGv.push_back(T(ri+2, ci+2, triplet_vec[i].value()*G(2,2) ));
      MGv.push_back(T(ri+3, ci+3, triplet_vec[i].value()*G(3,3) ));
   }

   SpMat MG(row, col);
   MG.setFromTriplets(MGv.begin(), MGv.end()); 


I think Eigen sparse still has more possibilities to improve and optimize performance.

anyway thanks a lot ggael ! I thought that nobody cares my post.

always appreciate your help.
Yeonchool.

Please excuse. I do not know any of your original problem *) but nevertheless or because of this I'm perplexed by some things:
a) A "4x4 sparse Matrix" : Is'nt that a contratiction in it self?
b) "G.setIdentity();" => "
float value = triplet_vec[i].value();
MGv.push_back(T(ri, ci, value ));
MGv.push_back(T(ri+1, ci+1, value ));
MGv.push_back(T(ri+2, ci+2, value ));
MGv.push_back(T(ri+3, ci+3, value ));
" should suffice?!
c)

actually, G is not a real identity matrix. it's diagonal but, last element is not 1. i.e., G is actually diag(1,1,1,r).
so, b) is not sufficient for my problem.
however when it is only limited to my problem, you are right. of course, last element can be changed to like MGv.push_back(T(ri+3, ci+3, value*r ));
and, what if G is not a diagonal matrix ? of course, the for-loop should be changed again.

FMD wrote:Please excuse. I do not know any of your original problem *) but nevertheless or because of this I'm perplexed by some things:
a) A "4x4 sparse Matrix" : Is'nt that a contratiction in it self?
b) "G.setIdentity();" => "
float value = triplet_vec[i].value();
MGv.push_back(T(ri, ci, value ));
MGv.push_back(T(ri+1, ci+1, value ));
MGv.push_back(T(ri+2, ci+2, value ));
MGv.push_back(T(ri+3, ci+3, value ));
" should suffice?!
c)

I see, you have a special case with only one element per row.

ggael wrote:What if you change line 56 of unsupported/Eigen/src/KroneckerProduct/KroneckerTensorProduct.h to:

AB.reserve(VectorXi::Constant(AB.outerSize(), A.nonZeros()*B.nonZeros()/AB.outerSize()));

Hi,

I am also experiencing that the Kronecker product for sparse matrices is (very) slow (I am using Eigen 3.2.0). It seems to be increasingly slow for large matrices. I suspect that it is entirely due to the poor estimation of the number of non-zeros. ggael suggests to change the reserve-statement to:

Code: Select all

AB.reserve(VectorXi::Constant(AB.outerSize(), A.nonZeros()*B.nonZeros()/AB.outerSize()));

This will distribute the number of non-zeros uniformly and it will almost always reserve too little for some inner vectors. It does not fix the problem.

I suggest to compute the exact number non-zeros by iterating over the A and B matrix. This requires minimal overhead and will always ensure that enough memory is allocated. I have tested this on a small set of matrices with varying sizes and it performs very well. To serve as an example I have included a simple, but working "pseudo code" below. I am sure the code can be improved.

Regards,
Bjarne

Code: Select all

// Calculate exact number of non-zeros for each inner vector
Eigen::VectorXi nnzA = Eigen::VectorXi::Zero(A.outerSize());
Eigen::VectorXi nnzB = Eigen::VectorXi::Zero(B.outerSize());
Eigen::VectorXi nnzAB = Eigen::VectorXi::Zero(AB.outerSize());

for (int jA = 0; jA < A.outerSize(); ++jA)
{
  int nnz = 0;
  for (Eigen::SparseMatrix<double>::InnerIterator itA(A,jA); itA; ++itA) nnz++;
    nnzA(jA) = nnz;
}

for (int jB = 0; jB < B.outerSize(); ++jB)
{
  int nnz = 0;
  for (Eigen::SparseMatrix<double>::InnerIterator itB(B,jB); itB; ++itB) nnz++;
    nnzB(jB) = nnz;
}

int innz = 0;
for (int i = 0; i < nnzA.rows(); ++i)
{
  for (int j = 0; j < nnzB.rows(); ++j)
  {
    nnzAB(innz) = nnzA(i)*nnzB(j);
    innz++;
  }
}

AB.reserve(nnzAB);


Done:

https://bitbucket.org/eigen/eigen/commi ... 883529a4e8

Hi,

I used your updated version (2014-02-14) of KroneckerProductSparse() with DynamicSparseMatrices.
The compiler remarks the way how reserve() is used with error: no matching function for call to reserve(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1> >::MapType)
Details see KroneckerTensorProduct.h:177:5
dst.reserve(VectorXi::Map(nnzAB.data(), nnzAB.size()));
instead of
void reserve(Index reserveSize = 1000)

Is this reproducible from your side?

--
Martin