Iterative Solvers?

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

Start by reading the referenced pgmresDR.m solver and the repository's existing benchmark entry points, though no in-repository files or tests are named. The issue discusses GMRES, tolerance estimation, warm starts, and preconditioners; completion would require a decided scope for iterative-solver benchmarks and documented comparative results.