Minimisation problem in thousands of dimensions

Question

I need to find the minimum of a function (a log-likelihood from a Potts model) in tens of thousands of dimensions. The function evaluation is quite fast, takes about $10^{-3} s$, and there is a gradient available (not implemented yet, though, but I expect similar timings). The hessian can be computed analytically, but it can be tedious given the number of dimensions.

Which algorithm will yield the best results in this case? All the documentation I have been able to find gets at most to 100 dimensions.

The process will be repeated several times for different configurations, so I will be able to benefit from a warm restart. In fact, as I am going to run them in parallel, each individual process doesn't have to be parallelisable on itself.

Answer 1

Matrix-free methods will work fine on problems with thousands of variables. Specifically, a properly implemented Newton-CG or trust-region code does not need to form a dense Hessian in order to compute a Newton iteration, but rather requires only the action of the Hessian on a vector. These methods rely on a truncated Krylov solver such as truncated-CG or truncated-MINRES in order to partially solve the linear systems. In my experience, a matrix-free, trust-region algorithm using truncated-CG works the best. It tends to work better than Newton-CG, which is a line-search method, because the trust-region helps determine how many Krylov iterations are required per iteration. In Newton-CG, you have to guess and set this number a priori. For a short description of Newton-CG, see Nocedal and Wright's book Numerical Optimization on page 169. That book also contains a description of a trust-region algorithm on page 69. Use this with truncated-CG, which is described on page 171. In addition, the big blue Trust Region Methods book has the most complete description of the trust-region algorithm. That book is a little overwhelming. Start with the basic trust-region algorithm on page 116. In step 2, use truncated-CG, which is described on page 202.

Now, these methods are not perfect. In general, you'll get better performance than a first-order method for a variety of reasons, but a high-order convergence rate is only guaranteed when close to the actual solution and only if we solve the Newton system accurately enough. If we're using a truncated Krylov solver such as truncated-CG, this is dictated by the clustering of the eigenvalues of the Hessian, which we likely can't compute. However, like I said above, in general it works well. If you want me to be frank, I've never seen a situation where a quasi-Newton method worked better than a good matrix-free optimization algorithm. I'm sure it happens, but not that often.

Finally, in terms of a warm-start, the only real issue with warm-starts would come into play if you had inequality constraints present in the problem. In unconstrained optimization, I'm not aware of any algorithm that has a problem with warm-starts. Though, there's a bit of a numerical conundrum if you actually start at the optimal solution. In any case, inequality constrained algorithms like primal-dual interior-point methods can have trouble with warm starts. However, other inequality constrained algorithms such as the Coleman-Li reflective Newton algorithm have zero problems with this. If you're close to the solution, just ignore the reflective bits of the code and it works great. As a note, good equality constrained algorithms also have no issues with warm-starts. These are things like sequential quadratic programming (SQP), or better, composite step SQP methods.

Finally, all of these algorithms are already implemented and freely available in Optizelle. It's BSD licensed and has hooks for C++, Python, and MATLAB. If you want a limited memory quasi-Newton method, it also has BFGS and SR1. Certainly, you can implement these algorithms yourself. And, to be honest, the codes for unconstrained optimization aren't that bad. However, Optizelle already does it and has been used on problems with upwards of half a billion variables.

Answer 2

In general, for something with maybe 100 variables, I'd probably try sequential quadratic programming. For very large problems, an interior point method would be better. Use an L-BFGS approximation; in practice, it's usually fine.

Most constrained nonlinear programming solvers will be based on sparse direct linear algebra; theory for preconditioning KKT systems is on the bleeding edge of research, so implementations for preconditioned iterative interior point method are on the horizon.

If your problem is box-constrained or unconstrained, and your problem is convex, CG is viable, too, and any quasi-Newton-like method (L-BFGS, etc.) or nonlinear iterative methods like NCG, NGMRES, etc. should work; use it with a globalization strategy (line search or trust region) for better results.

Answer 3

Is your Hessian sparse? If so, what I write below doesn't really apply.

The Hessian of a general function in tens of thousands of variables will have at least hundreds of millions of entries. That is, it'll take upwards of a gigabyte of RAM. Evaluating it will suck enough, and then you'll have to do linear algebra with the thing. No thanks!

You might try using nonlinear conjugate gradient or a limited-memory BFGS variant like L-BFGS here. Straight BFGS will also probably be too expensive, since it'll store an approximate Hessian inverse that's a gigabyte plus.

Answer 4

Depending on your problem details you may be able to use a Levenberg-Marquardt style algorithm. If its not amenable to transformation into this kind of problem then I'd stick with CG, possibly with some kind of preconditioner if you can find one. I've used both successfully in cases with 10's of thousands of dimensions.

Answer 5

I have good experiance with BOBYQA.

Although I don't know if it is the best algorithm, but it was sufficient for my needs (~250 parameter)

There are already some implementations (Fortran, C, python) with good descriptions

As the author states:

The user can substitute modern computation power for mathematical effort.

Answer 6

I suppose your problem (log-likelihood from a Potts model) is convex. For very large scale problems, using Hessian (second order) information may not be scalable and/or efficient. There is a growing body of work on first order methods for large scale optimization problems. They are designed to be scalable and efficient. Most first order methods are parallelizable. Refer to these surveys as a starting point:

Convex optimization for big data

Playing with duality: An Overview of Recent Primal-Dual Approaches for Solving Large-Scale Optimization Problems

Answer 7

Perhaps try some genetic algorithms if you are satisfied with a solution that is not perfect?