Take the 2-minute tour ×
Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It's 100% free, no registration required.

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.

share|improve this question
    
Have you looked at Monte Carlo methods? Those are typically used in very high dimensional problems. –  H H Aug 27 at 23:40
    
@HH I have never heard of MC used in minimisation. Yes, you can explore the space, but in very high dimensions you get to sample only a small section of the search space. –  Davidmh Aug 28 at 0:04

6 Answers 6

up vote 6 down vote accepted

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.

share|improve this answer
    
This answer deserves several careful readings. Matrix free methods seem the best way forward, and probably I will be able to increase the maximum size of my problem. Also, Optizelle looks a great asset. –  Davidmh Aug 28 at 12:25

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.

share|improve this answer

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.

share|improve this answer
    
The hessian is as dense as it gets, good point about the RAM. Do you know how L-BFGS or CG scale with number of dimensions? –  Davidmh Aug 28 at 0:02
    
@Davidmh: Linearly. But they're sensitive to the shape of the objective functions. It is at this point that I point out that both, especially CG, are really easy to code up; try CG, see if it works for your problem, and report back. –  tmyklebu Aug 28 at 0:04
    
On second thought, even if memory was not a problem, unless there is some structure you can make use of, just computing it would take hours per iteration! –  Davidmh Sep 1 at 9:52

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.

share|improve this answer

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.

share|improve this answer

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

share|improve this answer
1  
The region around the minimum is quite flat, so I don't think it would get good results. It may be good for a first iteration, though, but I think a warm restart would give better starting points. –  Davidmh Aug 28 at 13:45

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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