Hi,

I recently started using Eigen (3.1.0-alpha1) more seriously and am fairly impressed. I'm trying to use it to replace CHOLMOD in my code and am running into a few issues. I'm solving a nonlinear least squares problem with the Gauss-Newton algorithm, where the inner loop involves solving a sparse linear least squares problem in the form of Jx = r, where J is tall (m >> n) and sparse (most rows have only 6 non-zero entries). I'm solving it using the method of normal equations: Jt J x = Jt r. Eigen nicely provides a routine for computing the product of two sparse matrices, which would normally be a giant pain. In my case, A = Jt J is guaranteed to be symmetric (and positive definite). I'd like to feed the result of A = Jt J to PARDISO from Intel MKL, but PARDISO requires that A be:

- stored in compressed sparse row
- (since it's symmetric) only the upper triangle be stored
- (ideally) 1-based indexing (otherwise it makes a copy)

Is there a way to get the CSR form of the upper triangle of a sparse-sparse product in a not-incredibly-expensive fashion? Better yet, is there a way to tell Eigen that the result of the multiplication will be symmetric to make it faster?

In a related question, since this is a Gauss-Newton solver, I run many iterations of the sparse solver, where the sparsity structure stays the same. What is the recommended way of clearing the matrix and avoiding reallocating memory? I'm currently creating a new matrix every time and doing insert(). Should I use insert() on the first iteration, and use coeffRef() on later iterations? Finally, is there a way to reset (i.e., clear internal memory) a sparse matrix besides constructing a totally new one?

Thanks,

Jiawen

In short:

SparseMatrix<double,RowMajor> copy(n,n);
copy = A.triangularView<Upper>();

or directly:

SparseMatrix<double,RowMajor> A(n,n);
A.selfadjointView<Upper>().rankUpdate(J);

This last operation still does the full product but it will be optimized in the future.

If you are sure the pattern stay exactly the same, then coeffRef() is not bad. If you "insert" the element in a coherent order you might even use an iterator:

for (int k=0; k<mat.outerSize(); ++k)
for (SparseMatrixType::InnerIterator it(mat,k); it; ++it)
{
it.valueRef() = ...;
}

Hi, thanks for the fast reply. There are a few issues. I wrote the following tiny piece of code:

int main( int argc, char* argv[] )
{
int m = 4;
int n = 3;
Eigen::SparseMatrix< float, Eigen::RowMajor > j( m, n );
j.insert( 0, 0 ) = 1;
j.insert( 0, 2 ) = 1;
j.insert( 1, 2 ) = 1;
j.insert( 2, 0 ) = 1;
j.insert( 2, 2 ) = 1;
j.insert( 3, 0 ) = 1;
j.insert( 3, 1 ) = 1;

j.makeCompressed();

auto jt = j.transpose();

Eigen::SparseMatrix< float, Eigen::RowMajor > jtj( n, n );

#if 0
// does not work!
auto triangularView = ( jt * j ).triangularView< Eigen::Upper >();
jtj = triangularView;
#else
auto selfadjointView = ( jt * j ).selfadjointView< Eigen::Upper >();
jtj = selfadjointView;
#endif

std::cout << "j = " << std::endl << j << std::endl;

std::cout << "jtj = " << std::endl << jtj << std::endl;

return 0;
}

The triangular view trick doesn't seem to work, it complains that:


Error 3 error C2039: 'InnerIterator' : is not a member of 'Eigen::SparseSparseProduct<LhsNested,RhsNested>' c:\work\libs\cpp\include\eigen\src\sparsecore\sparsetriangularview.h 79 1 TestEigen
Error 5 error C2039: 'InnerIterator' : is not a member of 'Eigen::SparseSparseProduct<LhsNested,RhsNested>' c:\work\libs\cpp\include\eigen\src\sparsecore\sparsetriangularview.h 81 1 TestEigen
Error 4 error C2499: 'Eigen::SparseTriangularView<MatrixType,Mode>::InnerIterator' : a class cannot be its own base class c:\work\libs\cpp\include\eigen\src\sparsecore\sparsetriangularview.h 80 1 TestEigen

But the self adjoint view does work, so I went with that.

The other problem is that the resulting matrix jtj is a very strange: when I print it, it has no non-zero entries! Prunning the product doesn't help:

auto selfadjointView = ( jt * j ).pruned().selfadjointView< Eigen::Upper >();
jtj = selfadjointView;

However, I can prune the result:

auto selfadjointView = ( jt * j ).selfadjointView< Eigen::Upper >();
jtj = selfadjointView;
jtj.prune( 0.f );

I've only experimented with relatively small matrices (~3k x ~3k, 6 entries per row) in debug mode, but I suspect that this sparse-sparse product followed by copying is is very slow. Is this what you meant by "the last operation does the full product"? I can imagine that being the case for rankUpdate(), but am I correct in assuming that the same thing happens when I assign a selfadjointView to a SparseMatrix?

Thanks again for all the help,

Jiawen

sparse-sparse product is not a lazy operation, it evaluates into a temporary. So:

jtj = ( jt * j ).selfadjointView< Eigen::Upper >();

is equivalent to:

tmp = ( jt * j );
jtj = tmp.selfadjointView< Eigen::Upper >();

and the last copy has actually no effect since it means "take the upper triangular part, form a virtual selfadjoint matrix from it, and copy the full selfadjoint matrix into jtj". What you want is:

jtj = tmp.triangularView< Eigen::Upper >();


Currently, jtj.selfadjointView<Upper>().rankUpdate(J), does exactly that, i.e., evaluate into a temporary, and copy the upper triangular part only. In the future, it will call a special product that computes only the respective triangular part.


Regarding pruning, rankUpdate() does not allow that currently. However, the following should really works:

tmp = ( jt * j ).pruned();
jtj = tmp.triangularView< Eigen::Upper >();

or in one line:

jtj = ( jt * j ).pruned().eval().triangularView< Eigen::Upper >();


Note the .eval() to workaround the issue you hit. It will be fixed soon.

the evaluation is fixed, and I confirm the following does what it is expected to do:

A = (M*M.transpose()).pruned().triangularView<Upper>();

i.e., it removes explicit zeros and copy only the upper half into A.

That all makes sense. It still doesn't work though (with 3.1.0-alpha1, should I try using the development branch instead?)

// I need Jt * J, not J * Jt
Eigen::SparseMatrix< float, Eigen::RowMajor > jtj3( n, n );
jtj3 = ( j.transpose() * j ).triangularView< Eigen::Upper >();

Gives me the same error as before:

Error 3 error C2039: 'InnerIterator' : is not a member of 'Eigen::SparseSparseProduct<LhsNested,RhsNested>' c:\work\libs\cpp\include\eigen\src\sparsecore\sparsetriangularview.h 79 1 TestEigen
Error 5 error C2039: 'InnerIterator' : is not a member of 'Eigen::SparseSparseProduct<LhsNested,RhsNested>' c:\work\libs\cpp\include\eigen\src\sparsecore\sparsetriangularview.h 81 1 TestEigen
Error 4 error C2499: 'Eigen::SparseTriangularView<MatrixType,Mode>::InnerIterator' : a class cannot be its own base class c:\work\libs\cpp\include\eigen\src\sparsecore\sparsetriangularview.h 80 1 TestEigen

If I stick in an eval:

jtj3 = ( j.transpose() * j ).eval().triangularView< Eigen::Upper >();

I get an assert:

Assertion failed: (m_outerIndex[outer+1]-m_outerIndex[outer]==0 || m_data.index(m_data.size()-1)<inner) && "Invalid ordered insertion (invalid inner index)", file c:\work\libs\cpp\include\eigen\src\sparsecore\sparsematrix.h, line 385

Adding a pruned() doesn't help, same error:

jtj3 = ( j.transpose() * j ).pruned().eval().triangularView< Eigen::Upper >();

Finally, the rankUpdate trick almost works, but gives me the wrong answer:

jtj2.selfadjointView< Eigen::Upper >().rankUpdate( j.transpose() );

=============== printout ===============
j =
Nonzero entries:
(1,0) (1,2) (1,2) (1,0) (1,2) (1,0) (1,1)

Column pointers:
0 2 3 5 $

1 0 1
0 0 1
1 0 1
1 1 0

Right answer:

jtj =
Nonzero entries:
(3,0) (1,1) (2,2) (1,0) (1,1) (2,0) (3,2)

Column pointers:
0 3 5 $

3 1 2
1 1 0
2 0 3

Wrong answer:

jtj2 =
Nonzero entries:
(3,0) (1,1) (2,2) (1,1) (3,2)

Column pointers:
0 3 4 $

3 1 2
0 1 0
0 0 3

I guess it ignored the lower triangle somewhere...

Thanks for the detailed replies!

--Jiawen

you need the "devel branch" to get the fix:
http://eigen.tuxfamily.org/index.php?ti ... e#Download

Then I thought that you wanted a sparse matrix with only the upper half stored? All my suggestion were to achieve that precise goal!