I believe that I asked about this a few months ago but I cannot seem to find the posting now.

In my application I optimize a function of a vector. The elements of that vector determine the values in the "non-zero" positions in a square sparse matrix, Lambdat, which is multiplied by another sparse matrix, Ut. The function value requires the sparse Cholesky factor of

Code: Select all

Lambdat * Ut * Ut.adjoint() * Lambdat.adjoint() + Identity


where "Identity" is a sparse identity matrix of the same size as Lambdat.

Each evaluation of the function to be optimized requires an update of the numerical values of Lambdat * Ut (I use CHOLMOD functions for the Cholesky factorization and they can factor the expression above from Lambdat * Ut). Because I will do repeated evaluations I do the symbolic decomposition (analyzePattern) once only and perform just the numeric decomposition (factorize) in the evaluation.

This all works fine until I need to evaluate on the boundary of the parameter space, which can result in values of zero in some of the "non-zero" positions in the product. When this occurs these zeros are pruned from the product. I assume the pruning takes place in the calculation of the product, probably in the templated function Eigen::internal::sparse_product_impl. In fact, if I am not mistaken it occurs in the piece of code

Code: Select all

    res.startVec(j);
    for (typename AmbiVector<Scalar,Index>::Iterator it(tempVector); it; ++it)
      res.insertBackByOuterInner(j,it.index()) = it.value();


because the increment operator on an iterator skips positions whose values are zero.

What is the best workaround? For the time being I fall back on cholmod_ssmult when I detect that the number of nonzeros in the Eigen sparse-sparse product has shrunk and then convert the cholmod_sparse result to an Eigen sparse matrix. Given the structure of the code in SparseSparseProduct.h it does not seem easy to suppress the pruning without redefining the behavior of the Ambivector::Interator but that is clearly undesirable because the current behavior makes sense for most applications.

Another possibility I have considered is caching a copy of the matrix used in the analyzePattern method call then, when pruning has been detected, zeroing the values in the cached matrix and copying the non-zeros from the product to there.

Suggestions welcome. Also I am happy to elaborate if what I have described is too obscure.

there is an "experimental" second implementation of sparse-sparse products which does not prune zeros. You can try it like this:

res._experimentalNewProduct(A,B); // res = A*B

You have to uncomment the lines 86 to 107 in Eigen/src/Core/SparseSparseProduct.h to get it working.

It is currently slower than the default implementation based on the AmbiVector because this algorithm does not naturally generated sorted non zeros. The aforementioned lines does the sort. However, in some applications (e.g., for LLT solving) we don't really need sorted entries, and I'm currently implementing this possibility in SparseMatrix...

ggael wrote:there is an "experimental" second implementation of sparse-sparse products which does not prune zeros. You can try it like this:

res._experimentalNewProduct(A,B); // res = A*B

You have to uncomment the lines 86 to 107 in Eigen/src/Core/SparseSparseProduct.h to get it working.

It is currently slower than the default implementation based on the AmbiVector because this algorithm does not naturally generated sorted non zeros. The aforementioned lines does the sort. However, in some applications (e.g., for LLT solving) we don't really need sorted entries, and I'm currently implementing this possibility in SparseMatrix...

Thank you. I had seen that code and was wondering why it was commented out.

I eventually decided on another approach which is to cache the original product used for the symbolic decomposition and, in subsequent evaluations, create the product directly in that storage. Like most other sparse matrix operations some of the effort in the formation of the product (the use of the intermediate storage in the Ambivector, for example) is determining where the non-zeros in the result occur. If you already know that you can fall back on the numeric operations only.

The code is

Code: Select all

   fill(d_LamtUt._valuePtr(), d_LamtUt._valuePtr() + d_LamtUt.nonZeros(), Scalar());
   for (Index j = 0; j < d_Ut.cols(); ++j) {
       for (SpMatrixd::InnerIterator rhsIt(d_Ut, j); rhsIt; ++rhsIt) {
      Scalar     y = rhsIt.value();
      Index      k = rhsIt.index();
      for (SpMatrixd::InnerIterator lhsIt(d_Lambdat, k); lhsIt; ++lhsIt) {
          SpMatrixd::InnerIterator prdIt(d_LamtUt, j);
          Scalar                       x = lhsIt.value();
          Index                        i = lhsIt.index();
          while (prdIt.index() != i && prdIt) ++prdIt;
          if (!prdIt) throw runtime_error("logic error in updateLamtUt");
          prdIt.valueRef() += x * y;
      }
       }
   }



dmbates wrote:

ggael wrote:there is an "experimental" second implementation of sparse-sparse products which does not prune zeros. You can try it like this:

res._experimentalNewProduct(A,B); // res = A*B

You have to uncomment the lines 86 to 107 in Eigen/src/Core/SparseSparseProduct.h to get it working.

It is currently slower than the default implementation based on the AmbiVector because this algorithm does not naturally generated sorted non zeros. The aforementioned lines does the sort. However, in some applications (e.g., for LLT solving) we don't really need sorted entries, and I'm currently implementing this possibility in SparseMatrix...

Thank you. I had seen that code and was wondering why it was commented out.

I eventually decided on another approach which is to cache the original product used for the symbolic decomposition and, in subsequent evaluations, create the product directly in that storage. Like most other sparse matrix operations some of the effort in the formation of the product (the use of the intermediate storage in the Ambivector, for example) is determining where the non-zeros in the result occur. If you already know that you can fall back on the numeric operations only.

The code is

Code: Select all

   fill(d_LamtUt._valuePtr(), d_LamtUt._valuePtr() + d_LamtUt.nonZeros(), Scalar());
   for (Index j = 0; j < d_Ut.cols(); ++j) {
       for (SpMatrixd::InnerIterator rhsIt(d_Ut, j); rhsIt; ++rhsIt) {
      Scalar     y = rhsIt.value();
      Index      k = rhsIt.index();
      for (SpMatrixd::InnerIterator lhsIt(d_Lambdat, k); lhsIt; ++lhsIt) {
          SpMatrixd::InnerIterator prdIt(d_LamtUt, j);
          Scalar                       x = lhsIt.value();
          Index                        i = lhsIt.index();
          while (prdIt.index() != i && prdIt) ++prdIt;
          if (!prdIt) throw runtime_error("logic error in updateLamtUt");
          prdIt.valueRef() += x * y;
      }
       }
   }



Actually, I think I can tighten that code a bit by moving the declaration of prdIt outside the loop on lhsIt