6
$\begingroup$

I would like to implement a Numerical Optimizer in my favorite language Go. It shall find solution(s) to this problem: maximize a function f(x) where f is non-linear and x is a real vector of dimension 10 or 20. f is real-valued.

What is the best method in terms of:

  • simplicity to implement
  • opportunities to paralellize

The algorithm is supposed to run on a x64 CPU with multiple cores. Maybe one needs to look at different cases, so it would be okay for me to implement 2 or 3 algorithms.

Interesting cases:

  • f is rational or f is a composition of rational functions, exp and maybe log
  • f has multiple maximums

This is not intended as an open question. So I am not looking for the answer that lists as many algorithms as possible. Or an endless discussion of what algorithm might be the nicest. Choosing an algorithm can be somehow subjective, so I want to emphasize on a decent solution, rather than a perfect solution.

If have heard that simulated annealing is considered state-of-the-art when it comes to multiple maximums and a fixed domain. Although this sounds to me like: if I do not fix the domain, although I could fix it a priori, simulated annealing might be inferior to some other algorithm.

Update: suboptimal solutions are fine, however I want to make sure the algorithm doesn't concentrate to easily on a single maximum.

About the function: I want to do maximum likelihood fits on functions that estimate time series. Something like y_{n+1} = a_1 * y_n + a_2 + y_{n-1} + ... + b_1 * epsilon_n + ... + epsilon_n in a simple case. (epsilon is i.i.d. noise, y some time series, a_i and b_i real parameters .) Currently I'm doing this by assuming epsilon to be normal distributed but I want to change to the Cauchy-distribution. Moreover I want to add some non-linear terms and extend this from y real value to y being a real vector. Unfortunately I don't know yet exactly what I want, so the optimizer should not be too scoped.

Improve this question
$\endgroup$
5
  • 1
    $\begingroup$ It depends on whether you want a certificate of optimality, or you're happy with converging to a suboptimal local minimum (or maximum, as the case may be). $\endgroup$
    – Geoff Oxberry
    Nov 14, 2012 at 4:56
  • $\begingroup$ Yes, the fact that $f$ has multiple maximums indicates it is non-convex. You are always better off in this case if you can say something additional about $f$ or how important it is to achieve the global maximum. $\endgroup$
    – Aron Ahmadia
    Nov 14, 2012 at 5:14
  • $\begingroup$ @GeoffOxberry: no, a certificate is not necessary. The problem domain consists of functions that might have many suboptimal local minimums. I have only experience with some basic optimizers from Maple and Matlab, when minimizing non-linear, multivariate function, I easily get stuck in uninteresting local minima. (Eventhough I know that there are even deepeer minima elsewhere.) $\endgroup$
    – Philip
    Nov 14, 2012 at 6:06
  • $\begingroup$ @Aron Ahmadia: I've added some more details to the question $\endgroup$
    – Philip
    Nov 14, 2012 at 6:21
  • $\begingroup$ See code.google.com/p/par-sa. The code is parallel (uses MPI) and small/simple (less than 200 lines). But SA requires many thousands of evaluations of the function so it must be cheap to evaluate. Computationally expensive functions can be approximated with Taylor series but then you're effectively doing 'local' optimization. $\endgroup$
    – stali
    Nov 14, 2012 at 11:33

2 Answers 2

Reset to default
6
$\begingroup$

If you have explicit gradients, BFGS with More-Thuente line searches is the method of choice. It is a serial algorithm, fairly easy to implement, and quite robust. (A Matlab implementation of a variant of BFGS/More-Thuente is described in http://prod.sandia.gov/techlib/access-control.cgi/2010/101422.pdf. The code is publicly available and can guide your reimplementation.)

Parallelization needs some modifications, that you can easily invent yourselves after having understood how the algorithm works. You need multiple start options (with far away starting points) to find multiple mimina with wide domain of attraction.

Branch and bound using convex relaxations would give you assurance to find the global minimum, but it is far from easy to make efficient.

If you don't have gradients, I propose that you reimplement my SNOBFIT package, being a very robust robust package for derivative-free optimization of expensive functions, with excellent parallelization capability. See http://www.mat.univie.ac.at/~neum/software/snobfit/

Improve this answer
$\endgroup$
4
  • $\begingroup$ Thanks a lot for your answer. Unfortunately I don't understand your point about convexity. Will it help me to prove convexity at some point? $\endgroup$
    – Philip
    Nov 14, 2012 at 10:45
  • $\begingroup$ @Philip: A convex relaxation is a convex function underestimating the original function. Its minimizer will be a lower bound for the minimum of the original function. Such information (if it is good enough) can be used to prove approximate global optimality at some poiint. $\endgroup$
    – Arnold Neumaier
    Nov 14, 2012 at 10:57
  • $\begingroup$ I see, thanks again, this is even more than I expected. $\endgroup$
    – Philip
    Nov 14, 2012 at 17:48
  • $\begingroup$ @Philip: The nontrivial part, however, is the construction of good convex relaxations. I would recommend not to implement a b&b approach without first having a good knowledge of global optimization methodology, as the simplest schemes tend to be quite inefficient. $\endgroup$
    – Arnold Neumaier
    Nov 16, 2012 at 9:22
3
$\begingroup$

Your problem is called estimation by maximum likelihood of an ARMA(p,q) model. You are going to extend it to a general VAR (vector auto-regression) case. This problem has been solved many times, there are existing Time Series Analysis software suites that do just that and what's more, output statistics for related tests (ML, Wald and likelihood ratio, for starters). For instance, there exist R packages for this; commercial statistics software has various implementations of ML routines. Quite unfortunately, to tap into the wealth of existing knowledge you have to name your problem in the most specific way.

That said, a few ideas on the implementation. First off, I have absolutely no idea on computational peculiarities and IEEE-754 conformance in the Go programming language. With (say) Fortran, the nitty-gritty details are well-known. There are some languages/math libraries that make my head ache on the computational front (Java, I'm looking at you!).

Second, parallelization is the last step you'd like to take from an already-working code. The natural approach to parallelize here would be running identical routines from widely spaced starting points to convergence.

Third, if you are keen on implementing an algorithm on your own and not getting a pre-packaged solution, would suggest experimenting with Powell's algorithm for differentiable functions (BOBYQA). You can port it from NLOPT.

Fourth, your idea on Cauchy errors sounds very interesting.

Some books and articles to keep you busy:

  • Hamilton, J.D. Time Series Analysis, 1994. Princeton.
  • Box, G.E.P., Jenkins, G.M. and Reinsel, G.C., 1994. Time Series Analysis: Forecasting and Control. Holden-Day.
  • Box, G.E.P. and Luceño, A., 1997. Statistical Control by Monitoring and Feedback Adjustment. Wiley.
  • E. J. Hannan and A. J. McDougall. Regression Procedures for ARMA Estimation. Journal of the American Statistical Association Vol. 83, No. 402, Jun., 1988. pp.490-498

Please note this sort of questions could be better answered at Stats SE (Cross-Validated).

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thanks a lot for your answer. In fact I am already using R for this, but I am kind of unsatisfied. The results are not as good as I expected and I miss options for customization. In particular VAR or even defining models with custom, non-standard terms seems overly complicated... Having ported already a part of my custom construction to C, I'm already faster than code in the TSA or tseries package... ;) So writing the optimizer in C or Go, I expect things to fly. (+ The built-in non-linear R optimizer refuses to work for me, it has only few options.) $\endgroup$
    – Philip
    Nov 14, 2012 at 10:42
  • 3
    $\begingroup$ To your last remark: Stats SE may have the better statistical expertise, but has probably not the best expertise on optimization methods. $\endgroup$
    – Arnold Neumaier
    Nov 14, 2012 at 16:51

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.