Iterative Solvers? - Githubissues

Dear Alec,

I spent a while going down the rabbithole of trying to evaluate iterative solvers as well. Having some more careful, geometry processing oriented, benchmarks of iterative methods could be useful too. gmres is popular in graphics in practice as far as I understand, but there aren't great references for how it should be used as far as I know and if you start reading the numerical analysis literature it's quite easy to get lost by all the different tricks they try.

Anyway in my experience from personally trying a few solvers in a much less principled manner than you do here, 🐐 is: https://github.com/thanospol/MARIE/blob/master/src_numeric/src_iterative/pgmresDR.m

There's also a hard earned trick that I figured out that I'll document here. When using this sort of solver for the problem I've been working on, it's not clear how to set the "relative residual" term. The trick that I found that works pretty well is to do an initial approximate solve to determine the appropriate stopping accuracy for the given problem (which in my understanding depends on the conditioning of the input matrix), and then doing a second "warm start" solve, now with a well informed guess for the stopping criterion.

In particular this is useful for alternating minimization style algorithms where maybe at each step you don't care about necessarily having a super high accuracy result, but a priori it's not clear what the correct accuracy is during your solve.

My solver code to document solver parameters, YMMV:

      inner_steps = 25;
      deflate_basis = 2;

      % run once to make some progress and estimate tolerance
      [x_after_step1,flag,relres,iter,resvec] = pgmresDR(lagrangian,b, inner_steps, 1e-8, 0, deflate_basis, pre_fac',[],pre_fac,[], warmstart);
      init_gmres_flag_relres_iter_residue = [ flag, relres, iter, resvec(end) ] 

      solver_conf.gmres_cur_tol = relres;

      res_after_step1 = norm(b - lagrangian*x_after_step1);

      [fac,flag,relres,iter,resvec] = pgmresDR(lagrangian,b, 100, relres * .65, 4, 20, pre_fac',[],pre_fac,[], x_after_step1);

Other anecdata:

In my experience gurobi is also pretty great. Not free software but worth evaluating.

Eigen::SpQR, which binds Cholmod is much better than Eigen SparseQR, which is unreliable and not worth using.

Iterative solvers can produce meaningfully different results on high order problems from exact solvers and should be used with caution. Also much more sensitive to choice of preconditioner, which is it's own piece of dark magic where it is quite difficult to understand best practices.

Best,

Josh Vekhter

Interesting! I'll be a bit surprised if gmres is competitive with CG or cholesky for a lot of the sparse symmetric problems in geometry processing when high accuracy is required. You make a nice point about alternating minimizations. This is likely true for other descents like newton etc. where the system is changing until convergence.

Haha, I suppose my question then is: Suppose you've got a really big sparse matrix, something very high res, it seems like gmres and variants are the best bets?

Cholesky is the best, always totally agree! \ is very reliable, until the matrix is ~ 1^10x1^10 block sparse matrix when it starts running out of memory on most any machine. How should one then solve a system? You prefer CG?

In my experiments, per justin solomon, robert bridson and a maybe a few other peoples's suggestions, tried regularizing with block laplacian preconditioner with gmres. Specifically in the code above:

pre_fac = ichol( block_tet_mesh_laplacian )

A*pre_fac ~ I
pre_fac' * A~ I

This is then used as a symmetric preconditioner. At some point the linear algebra starts to take some thinking, but for big systems where you have resolution it seems to work pretty well actually. It's not actually faster then cholesky too, slower if you aren't fairly liberal with the parameters, (i.e. bit of early stopping needed to get convenient performance, though with the two step scheme the convergence isn't horrible... depending on what metric you are specifically looking at ;). )

P.S. More concretely:

In the problem this code is from, we have like 22 variables per element in a tet mesh. We alternate:

1 step of some convex program in these variables --- that we want to solve using any matrix factorization method that scales to sufficiently a large sparse matrix. 1 step local projection, 2nd order accurate ( i.e. newtonstep with hessian, implemented by Paul Zhang 🍿 ).

The second step is imminently scalable, and could even in principle be implemented on a GPU. For the former, is there a better option than a good version of gmres? The one I used has some mix of deflating and random restarts??? On top of that I did this thing where I take like 50 steps and see to what level it converges, then take that initialization, and run it again with an increased accuracy. This is necessary because it influences restarts.....

To get more symmetry we work in reduced coordinates. the 22 variables actually encode a 1 + 9 + 81 = 91 dimensional tensor.

On some level it's a question of what you are interested in, physical correctness, or plausible and scalable performance. I should investigate further sometime I suppose! I've learned the hard way that one must always use cholesky whenever possible though!!! That's the good stuff :).

On Mon, Jun 6, 2022 at 11:49 AM Alec Jacobson @.***> wrote:

Interesting! I'll be a bit surprised if gmres is competitive with CG or cholesky for a lot of the sparse symmetric problems in geometry processing when high accuracy is required. You make a nice point about alternating minimizations. This is likely true for other descents like newton etc. where the system is changing until convergence.

— Reply to this email directly, view it on GitHub https://github.com/alecjacobson/sparse-solver-benchmark/issues/3#issuecomment-1147602070, or unsubscribe https://github.com/notifications/unsubscribe-auth/AARSGZ5ECLNXCWKCCELL25LVNYMZJANCNFSM5X52GT6Q . You are receiving this because you authored the thread.Message ID: @.***>

alecjacobson / sparse-solver-benchmark

Iterative Solvers? #3