# Minimisation problem in thousands of dimensions

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.

• Have you looked at Monte Carlo methods? Those are typically used in very high dimensional problems. Aug 27, 2014 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. Aug 28, 2014 at 0:04

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.

• 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. Aug 28, 2014 at 12:25
• @wyer33 What are your thoughts if the Hessian is highly ill-conditioned or rank deficient, while it is still Positive Semidefinite?
– Hari
Aug 29, 2015 at 8:56
• @haripkannan Truncated Newton methods can still work well with ill-conditioned or rank-deficient Hessians. Really, what matters most in this context is the clustering of the eigenvalues of the Hessian. If they're well clustered, the methods still work well. If they are not well clustered, we work hard to generate a preconditioner that clusters them. Note, a Hessian can completely singular and the Krylov method only take two iterations. Clustering is the key. There's more discussion about this here: optimojoe.com/post/matrix-free-codes-affects-optizelle Aug 30, 2015 at 14:49
• @haripkannan There's a zoo of different Krylov methods and which one works best for a particular problem can be a little random. Mostly, your choice will boil down to what information you have available and what kind of solution you need. Honestly, the best reference I've seen for this is a flow chart from David Fong's Ph.D thesis on page 116, figure 7.1, (web.stanford.edu/group/SOL/dissertations/david-fong-thesis-online.pdf). Most of the time, I use MINRES or GMRES, but the diagram in David's thesis is the best reference. Aug 31, 2015 at 21:33
• @haripkannan As far as preconditioning, that's a black art and depends on the nature of the system. For example, for systems that arise from PDEs, there are physics specific preconditioners. For generic systems, we're stuck mostly with sparse direction factorizations. On this point, if the system is singular, the LU factorization still exists whereas things like a Choleksi do not, so I prefer that. Depending, either L or U will be singular, so we need a custom backsolver to skip 0s on the diagonal. Sometimes, incomplete factorizations are OK. Also, check out Saad's papers on multilevel ILU. Aug 31, 2015 at 21:42

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.

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.

• 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? Aug 28, 2014 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. Aug 28, 2014 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! Sep 1, 2014 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.

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.

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

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

• 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. Aug 28, 2014 at 13:45