# Trust-region Newton: implementation issue with Conjugate Gradient calculations

UPDATE: The problem turned out to be the step (refer penultimate paragraph below) where I was factoring out a small value from the vectors of the numerator and denominator and then computed dot products/norm-square. Computing these terms in straightforward manner resulted in CG steps converging with almost same number of steps as Eigen Solver.

I am implementing Steihaug method in C++, for large scale unconstrained convex optimization problems. The current instance, involves around 33000 variables; ideally, the problem sizes will be half-a-million to a million variables. The Hessian matrix is symmetric and quite sparse (~30% non-zero entries) and has a block structure about the main diagonal. The Hessian is PSD (never strictly PD) and has a huge (~$10^8$) condition number. I am using a diagonal preconditioner. The scientific problem is about applying log-sum-exp based smoothing to LP relaxation of combinatorial optimization problems.

As the method approaches the optimum it takes increasingly more internal CG iterations per outer iteration: ~7700 CG iterations in the last outer iteration which reaches the global optimum.

When I plug in the data of the last outer iteration to Eigen::BiCGSTAB solver, with diagonal preconditioner, it converges in 100 iterations to the desired optimum.

I would appreciate some pointers on making the implementation numerically robust.

edit: I suppose an important issue is that all non-zero numbers in the gradient and hessian are quite small in magnitude. Since, we are nearing the global optimum.

For the calculation of $\alpha$ (steepest descent step) and $\beta$ (coefficient in direction update), I am factoring out the smallest non-zero magnitude (and cancel it between the numerator and denominator) and computing the dot products. This is not helping. (UPDATE: this turned out to be an unnecessary step that caused the error in the implementation.)

Also, I am computing the residual (as $Ax - b$) every fifty iterations to compensate for drift.

33000 decision variables is still within the realm of direct solvers, so for at least that problem instance, you could try a sparse LU or sparse Cholesky factorization.

For larger problem instances, where direct solvers are not an option, you could try an incomplete Cholesky conditioner. Separate papers by Gondzio and by Waechter, Curtis, et al use left-ILU-preconditioned GMRES to solve the KKT system in their interior point methods. I would think that you could try something similar with CG.

• Ideally, the problem sizes will be half-a-million to a million variables. I will try incomplete Cholesky preconditioner. Thanks. Will the intermediate solutions from GMRES behave similar to CG for a trust region method? So that I can truncate the iterations based on a trust region condition?
– Hari
Jun 19 '15 at 7:21
• When do small numerical values create problems? can you point to some references?
– Hari
Jun 19 '15 at 7:22
• Any gut feelings about how Eigen managed to make it work with a diagonal preconditioner?
– Hari
Jun 19 '15 at 7:27
• Can you check the residual? Maybe Eigen just failed after 100 iterations. Jun 19 '15 at 12:44
• Thanks for your inputs. Please, refer the updated statement of the question.
– Hari
Jun 19 '15 at 16:21