Hi guys,
I have done quite some digging in the docs and here in the forum, but I can't seem to find an elegant solution to this seemingly simple operation:
At the beginning of my program, I have a big m*m integer matrix M_dense, initialized entirely to 0 (m > 5.000).
In the end, this matrix is going to be really dense, there should be not a single 0.
This is how it is filled:
Repeat thousands of times:
1. Compute a _sparse_ m*n matrix M_sparse (dont ask how; n can range from 10 to 10000)
2. Multply this M_sparse by its transposed self => new m*m matrix
3. Add this new m*m matrix to the m*m matrix M_dense (above)

I want to do steps 2-3 in Eigen, but I can't seem to find an elegant (memory, cpu-friendly) way to do it.
As Matrix types, I use for M_dense:

Code: Select all

typedef Eigen::Matrix<int,Dynamic,Dynamic,RowMajor> DenseMat;
typedef Map<DenseMatTmp> DenseMap;

and for M_sparse:

Code: Select all

typedef Eigen::SparseMatrix<int,RowMajor> SparseMat;

(RowMajor, because the data comes RowMajor from a C program)

I initialize the sparse matrix via a list of triplets (see docs), and M_dense like this:

Code: Select all

DenseMap M_dense(pointer_to_int_array, m, m);

(where pointer_to_int_array is a pointer to an int array with m*m elements.)

Now, ideally, I would like to do steps 2 and 3 like this:

Code: Select all

M_dense += M_sparse * SparseMat(M_sparse.transpose());

But this line doesn't work and results in a huge error message (traceback) with these last errors from the compiler:

g++ -c eigen_wrapper.cpp -o eigen_wrapper.o -I.
...
./Eigen/src/SparseCore/SparseSparseProductWithPruning.h:69: error: class Eigen::Matrix<int, -0x00000000000000001, -0x00000000000000001, 1, -0x00000000000000001, -0x00000000000000001> has no member named startVec
./Eigen/src/SparseCore/SparseSparseProductWithPruning.h:71: error:class Eigen::Matrix<int, -0x00000000000000001, -0x00000000000000001, 1, -0x00000000000000001, -0x00000000000000001> has no member named insertBackByOuterInner
./Eigen/src/SparseCore/SparseSparseProductWithPruning.h:73: error: class Eigen::Matrix<int, -0x00000000000000001, -0x00000000000000001, 1, -0x00000000000000001, -0x00000000000000001> has no member named finalize

Then I introduce an intermediate result (what I actually wanted to avoid), but now the operator is undefined:

Code: Select all

SparseMat M_sparse_squared(m,m);
   M_sparse_squared= SparseMat(M_sparse * SparseMat(M_sparse.transpose()));
   M_dense = M_dense + M_sparse_squared;

error:

g++ -c eigen_wrapper.cpp -o eigen_wrapper.o -I.
eigen_wrapper.cpp: In function void selfMultiplyMatrix(int*, int, int*, int, int):
eigen_wrapper.cpp:51: error: no match for operator+ in M_dense + M_sparse_squared
/usr/include/c++/4.4/bits/stl_bvector.h:264: note: candidates are: std::_Bit_iterator std::operator+(ptrdiff_t, const std::_Bit_iterator&)
/usr/include/c++/4.4/bits/stl_bvector.h:350: note: std::_Bit_const_iterator std::operator+(ptrdiff_t, const std::_Bit_const_iterator&)
make: *** [compile_wrapper] Error 1

What am I doing wrong here? What can I do to avoid a) an intermediate sparse matrix that first stores the result of the multiplication and b) an extra loop over this sparse matrix to add the values to the dense matrix?

Thanks a lot for your help!

Best
new2you

P.S:
Here is the complete code (latter version):

Code: Select all

#include <Eigen/Sparse>
#include <iostream>
#include <vector>

using namespace Eigen;
using namespace std;

typedef Eigen::SparseMatrix<int,RowMajor> SparseMat;
typedef Eigen::Matrix<int,Dynamic,Dynamic,RowMajor> DenseMat;
typedef Map<DenseMat> DenseMap;
typedef Eigen::Triplet<int> T;

extern "C" void selfMultiplyMatrix(int *triplet_array, int num_triplets, int *M_dense_array, int m, int n)
{
   vector<T> triplet_list;
   triplet_list.reserve(num_triplets);
   int i;
   int array_size = num_triplets*3;
   for(i=0; i<array_size; i=i+3)
   {
      triplet_list.push_back(T(triplet_array[i],triplet_array[i+1],triplet_array[i+2]));
   }

   SparseMat M_sparse(m,n);
   SparseMat M_sparse_squared(m,m);
   DenseMap M_dense(M_dense_array, m, m);
   M_sparse.setFromTriplets(triplet_list.begin(), triplet_list.end());
   
   M_sparse_squared = SparseMat(M_sparse * SparseMat(M_sparse.transpose()));
   M_dense = M_dense + M_sparse_squared;

}


dense += sparse is not supported yet, you can workaround with dense = dense + sparse though this will be a n^2 op.

I have a similar question. I have sparse matrix and I need to multiply it by itself and assign to a dense matrix. If I try to directly assign the product to a dense matrix, it fails at compile time:

Code: Select all

Eigen::SparseMatrix<double,Eigen::ColMajor> A;
// ... fill out A ...
// This line fails!
Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic,Eigen::ColMajor> Q = A.transpose() * A;


If I assign to a temporary, however, the code compiles fine:

Code: Select all

// This code works
Eigen::SparseMatrix<double,Eigen::ColMajor> temp_matrix = A.transpose() * A;
Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic,Eigen::ColMajor> Q = temp_matrix;


Can I avoid this temporary?

Thank you!

This doesn't work for me with gcc-4.2.1 on Mac 10.7.5 with Eigen 3.1.2:

$ g++ -IEigen -c t2.cpp -o t2.o
t2.cpp: In function ‘void selfMultiplyMatrix(int*, int, int*, int, int)’:
t2.cpp:30: error: no match for ‘operator+’ in ‘M_dense + M_sparse_squared’

Any ideas?