3
$\begingroup$

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.

Improve this question
$\endgroup$

1 Answer 1

Reset to default
1
$\begingroup$

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.

Improve this answer
$\endgroup$
6
  • $\begingroup$ 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? $\endgroup$
    – Hari
    Jun 19, 2015 at 7:21
  • $\begingroup$ When do small numerical values create problems? can you point to some references? $\endgroup$
    – Hari
    Jun 19, 2015 at 7:22
  • $\begingroup$ Any gut feelings about how Eigen managed to make it work with a diagonal preconditioner? $\endgroup$
    – Hari
    Jun 19, 2015 at 7:27
  • 1
    $\begingroup$ Can you check the residual? Maybe Eigen just failed after 100 iterations. $\endgroup$
    – Nick Alger
    Jun 19, 2015 at 12:44
  • $\begingroup$ Thanks for your inputs. Please, refer the updated statement of the question. $\endgroup$
    – Hari
    Jun 19, 2015 at 16:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.