Hello everybody,

I started to use Eigen library a couple of weeks ago, so I'm pretty much a newbie here.
My project relies on the efficient computation of sparse matrices operations, therefore I'm very focused into efficency.

For the project I'm developing, I need to compute the kronecker product (kp) between a small (3x3) sparse and a dense (10000x10000) matrix, and store its result in a sparse matrix(30000x30000). I tried the built-in kronecker product which, If I understood correctly, to perform such computation translates the dense matrix into a sparse one, and a sparse-sparse kronecker product is performed through kroneckerProduct_sparse.

I don't know whether I messed up something, but the performance I had were very disappointing (it took around 70 seconds to compute the result), therefore I tried to implement my own function, which rationale is as follows. Given resMat = KP(sMat, dMat), with sMat of size (is,js) and dMat of size (id,jd), I progressively compute the kp between each column of sMat and dMat, and then store the results in resMat inserting them column by column.

Here's the function code:

Code: Select all

   void Miscellanea::kroneckerProd(const SparseMatrixXd& m1, const MatrixXd& m2, SparseMatrixXd& res)
   {
      res.resize(m1.rows()*m2.rows(), m1.cols()*m2.cols());
      double tolerance = std::numeric_limits<SparseMatrixXd::Scalar>::epsilon();

      for (int jM1=0; jM1<m1.outerSize(); ++jM1)
      {
         int numElements = 0;
         int colIdx = jM1*m2.cols();
         std::vector<int> rowIdx;
         std::vector<MatrixXd> KpCol;
         for (SparseMatrixXd::InnerIterator it1(m1,jM1); it1; ++it1)
         {
            if(it1.value()>tolerance)
            {
               KpCol.push_back(it1.value()*m2);
               rowIdx.push_back(it1.row()*m2.rows());
               numElements++;
            }
         }

         for(int jM2=0; jM2<m2.cols(); ++jM2)
         {
            res.startVec(colIdx+jM2);
            if(numElements>0)
               for(int k=0; k<numElements; ++k)
                  for(int iM2=0; iM2<m2.rows(); ++iM2)
                     if(KpCol[k](iM2,jM2) > tolerance)
                        res.insertBack(rowIdx[k]+iM2,colIdx+jM2) = KpCol[k](iM2,jM2);
         }
      }
      res.finalize();
   }

Now the computation required less than 0.2 seconds to compute the same results. As I said, I may have forgot something. However, I decided to post the code so that (if it makes sense) it may be of some use to others as well. Cheers,

Francesco

Thanks for sharing. I think the KP in *unsupported* does not explicitly handle sparse*dense and so threat it as a sparse*sparse...

Hi Ggael,

yes, i can imagine that there is a lot of stuff "under construction", particularly for the sparse matrices which are quite new. That's why I wanted to share my code so that, if it can be of any help for further developing the library, someone may incorporate it.

I only tested it for small sparse matrices, meaning a small number of temporary matrices that need to be created and stored in the vector, therefore I don't know how performance degrades for bigger mats, however... hope it helps! Cheers!