0

I am working in a code that makes use of large sparse matrices. So far I have been using the Eigen3 SparseMatrix class to store linear operators that degrade an image. In the operator, each row corresponds to one ouput pixel in the degraded image and each column describes the contribution from other pixels in the field-of-view.

I need to calculate the product Op.transpose()*Op, but it is taking more than 10 minutes to calculate that. Given that the operator uses real numbers, the result should be symmetric and I could save 1/2 of the time by only computing the upper or lower triangle of the matrix.

I thought I could do:

    Eigen::SparseMatrix<T,RowMajor,ind_t> b(npix,npix);
    b.triangularView<Eigen::Lower>() = Op.transpose() * Op;

I am just getting compilation errors related to the operator= not being overloaded, but several posts and the original documentation suggest that this should be possible. I am not sure I understand what am I doing wrong. Any tips to speed up even more that calculation?

Here is an example code that illustrates the problem with a small matrix.

#include <iostream>
#include <Eigen/SparseCore>
#include <Eigen/Core>
#include <vector>

using namespace Eigen;

int main(){
  using T = double;

  // --- Fill a test sparse matrix from triplets --- //
  
  SparseMatrix<T,RowMajor,long> a(4,4), c(4,4);
  std::vector<Triplet<T,long>> p(16);
  
  for(int ii=0; ii<4; ii++)
    for(int jj=0; jj<4;++jj)
      p[ii*4+jj] = Triplet<T,long>(ii,jj,double(ii*4+jj));


  a.setFromTriplets(p.begin(), p.end());


  
  // --- try calculating only 1/2 of the matrix-matrix product --- //
  
  c.triangularView<Upper>()  = a.transpose() * a;


  
  // --- print results --- //

  std::cout<<c<<std::endl;

}

Edit: Thanks to one of the comments, I could make it work with:

c.selfadjointView<Lower>().rankUpdate(a.transpose());
JdlCR
  • 155
  • 7
  • 2
    Did you try the `rankUpdate()` method, as suggested in this answer: https://stackoverflow.com/a/60270593/6870253? – chtz Jun 20 '23 at 11:29
  • 1
    Also, how large and how dense are your matrices? Spending 10minutes for a matrix product is quite a lot. Maybe you don't need this explicitly, e.g., just apply `A` and `A.transpose()` to a vector after each other. Or if you want to solve a least-squares system, try an iterative solver. – chtz Jun 20 '23 at 11:33
  • The sparsity of the matrix depends quite a bit on the image degradation process. In the example that I am using now, the operator matrix has dimensions 11664x46656 with 64256256 non-zero elements (~12%). Another example has dimensions 46656x46656 with 140991876 non-zero elements (~6.5%). – JdlCR Jun 20 '23 at 13:00
  • I suspect the matrix is still too dense. If you have a uniform distribution of nonzeros with 6% probability over `n=46656`, the chance of a scalar in the output being zero is `(1-0.06**2)**n`, effectively zero. Meaning the result is fully dense. You probably have to dispatch between dense and sparse routines and use parallelization. (Eigen can only parallelize [dense rank updates](https://eigen.tuxfamily.org/dox/classEigen_1_1SparseSelfAdjointView.html#a843bdc3844b4b9f5bc2daa84520d7bb7) through [external BLAS libraries](https://eigen.tuxfamily.org/dox/TopicUsingBlasLapack.html) like OpenBLAS) – Homer512 Jun 20 '23 at 15:20

1 Answers1

1

Thanks to one of the comments, the code compiles now with the correct result. For completeness:

#include <iostream>
#include <Eigen/SparseCore>
#include <Eigen/Core>
#include <vector>

using namespace Eigen;

int main(){
  using T = double;

  // --- Fill a test matrix from triplets --- //
  
  SparseMatrix<T,RowMajor,long> a(4,4), c(4,4);
  std::vector<Triplet<T,long>> p(16);
  
  for(int ii=0; ii<4; ii++)
    for(int jj=0; jj<4;++jj)
      p[ii*4+jj] = Triplet<T,long>(ii,jj,double(ii*4+jj));


  a.setFromTriplets(p.begin(), p.end());


  
  // --- try calculating only 1/2 of the matrix-matrix product --- //
  
  
  //c.triangularView<Lower>()  = a.transpose() * a;
  c.selfadjointView<Lower>().rankUpdate(a.transpose());

  
  // --- print results --- //

  std::cout<<c<<std::endl;
  std::cout<<a.transpose()*a<<std::endl;
   
}
JdlCR
  • 155
  • 7