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.