Hi all !

I have used BiCGSTAB solver for my sparse matrix algebra to solve computer vision problem, X = A \ B. matrix A size is 24000 x 24000 (it varies depends on number of 3D points). as you see, its size is too big. so, I had to use sparse matrix instead of dense matrix.

so far, so good. BiCGSTAB was a very reliable solution. however, only one problem is remained, speed. it consumes 400-500 ms for solving my matrices. however I need less than 200 ms. thus, I planned to test UmfPack because I believed it was well-known and fast, and it can be easily compiled on windows. (BTW, UmfPack is really faster than BiCGSTAB ?)

but, I encountered a runtime error. experts! please check the following code.

Code: Select all

typedef Eigen::SparseMatrix<double> SpMatd;

   SpMatd AA(4*NV,4*NV);      
   AA.reserve(Eigen::VectorXi::Constant(AA.outerSize(), A.nonZeros()/AA.outerSize()));
   AA = (A.transpose()*A);      // 4*NV x 4*NV (square matrix)
   AA.finalize();
   AA.makeCompressed();

   Eigen::VectorXd AB_col1(4*NV);
   Eigen::VectorXd AB_col2(4*NV);
   Eigen::VectorXd AB_col3(4*NV);
   AB_col1 = AB.col(0);
   AB_col2 = AB.col(1);
   AB_col3 = AB.col(2);

   Eigen::VectorXd Xcol1(4*NV);
   Eigen::VectorXd Xcol2(4*NV);
   Eigen::VectorXd Xcol3(4*NV);

   Eigen::UmfPackLU< Eigen::SparseMatrix<double> > solver;
   solver.compute(AA);            // <-- stopped here  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   Xcol1 = solver.solve(AB_col1);
   Xcol2 = solver.solve(AB_col2);
   Xcol3 = solver.solve(AB_col3);


in order to know the exact line that causes error, I changed the above command "solver.compute(AA)" like as follows based on UmfPackSupport.h.

Code: Select all

        double *null = (double *) NULL ;
   int i ;
   void *Symbolic, *Numeric ;
   (void) Eigen::umfpack_symbolic (4*NV, 4*NV, AA.outerIndexPtr(), AA.innerIndexPtr(), AA.valuePtr(), &Symbolic, null, null) ;
   (void) Eigen::umfpack_numeric (AA.outerIndexPtr(), AA.innerIndexPtr(), AA.valuePtr(), Symbolic, &Numeric, null, null) ;       // <-- stopped here !!!!


so, finally I was able to know which line causes the error. the criminal was umfpack_numeric function !
however I still don't know why. is there a mistake ? or is the matrix too big ?
I need your help to solve this problem.

Thanks in advance.
Best Regards,
Yeonchool.

p.s. FYI, I'm attaching simple test code. the following code was successfully performed.

Code: Select all

   int    n = 5 ;
   int    Ap [ ] = {0, 2, 5, 9, 10, 12} ;
   int    Ai [ ] = { 0,  1,  0,   2,  4,  1,  2,  3,   4,  2,  1,  4} ;
   double Ax [ ] = {2., 3., 3., -1., 4., 4., -3., 1., 2., 2., 6., 1.} ;
   double b [ ] = {8., 45., -3., 3., 19.} ;
   double x [5] ;

   double *null = (double *) NULL ;
   int i ;
   void *Symbolic, *Numeric ;
   (void) Eigen::umfpack_symbolic (n, n, Ap, Ai, Ax, &Symbolic, null, null) ;
   (void) Eigen::umfpack_numeric (Ap, Ai, Ax, Symbolic, &Numeric, null, null) ;
   umfpack_di_free_symbolic (&Symbolic) ;
   (void) Eigen::umfpack_solve (UMFPACK_A, Ap, Ai, Ax, x, b, Numeric, null, null) ;
   umfpack_di_free_numeric (&Numeric) ;
   for (i = 0 ; i < n ; i++) printf ("x [%d] = %g\n", i, x [i]) ;


Maybe your matrix is nearly singular? Have you tried the IncompleteLUT preconditioner for the BiCGSTAB? If not I recommend you to try it. I'm pretty sure you'll manage to fall below the 200ms. It has two threshold parameters to adjust it to your needs.

Thanks ggael for your guidance and quick reply !

Based on your comment, I set two parameters, setMaxIterations(1000) and setTolerance(1e-03). so now I got less than 100 ms performance.
however, it was only when I used DiagonalPreconditioner instead of IncompleteLUT.
On the contrary, I got more than 60 seconds when I used IncompleteLUT.

Regarding this problem, I have three questions..
First, could you guess what was the problem ?
Second, what's the difference between DiagonalPreconditioner and IncompleteLUT ?
Lastly, is direct solver generally faster than iterative solver ?

For other members, I'm attaching my code as follows.

Code: Select all

   Eigen::BiCGSTAB<Eigen::SparseMatrix<float>, Eigen::DiagonalPreconditioner<float> > solver;
   Eigen::BiCGSTAB<Eigen::SparseMatrix<float>, Eigen::IncompleteLUT<float> > solver;
   solver.setMaxIterations(1000);
   solver.setTolerance(1e-03);

   Xcol1 = solver.compute(AA).solve(AB_col1);


sorry for many questions.

thanks in advance and have a good day.
Yeonchool

Regarding question 1,

I found how to set the two parameters you said as follows.

solver.preconditioner().setDroptol(0.1);
solver.preconditioner().setFillfactor(5);

by the way, how to choose optimal parameters ? hmm... that would be fourth question..

I'm back..

I found a good information in terms of comparison between direct and iterative solvers.
http://www.mathworks.com/matlabcentral/newsreader/view_thread/299821

According to the comments on the site, direct solver is generally faster than iterative solvers.
however, direct solver has one big drawback. it is memory limitation.
so, I think if there are not many zeros and one doesn't use 64 bit machine, iterative solver is better than directive.

ggael, is this correct ? do you agree ?

regarding IncompleteLUT, it is usually better to set a large fill-factor (100), and play with the drop-tolerance. In most cases, for BiCGSTAB the IncompleteLUT outperform the DiagonalPreconditioner. You could also decompose the call to compute(A) into a call to analyzePattern(A) followed by factorize(A) and track the time of each step. For BiCGSTAB, the analyzePattern step consist in computing the preconditioner. It it takes too long, increase the droptol.