optimizationjuliahessian

How to use Optim.jl optimize function with a summation


I'm currently working on a research project where I need to find the correct value for $\beta$ that minimizes the function $\sum_i (p_i-\avg(n_c)_i)^2$. Here, $\avg(n_c)_i$ is an array that is the function I'm trying to match $p_i$ to. $p_i$ depends on $\beta$ but also has some other dependencies, like $\chi^2=\sum_i p_i=N$ where $N$ is the number of particles in the system.

So essentially, I'm trying to run my code to take a starting $\beta$ value, solve for $p_i$ so it matches the summation condition, and check if $\chi^2$ is minimized. From there, I want to update $\beta$ to a new value until I get within a reasonable cutoff value.

This is where I started trying to implement the optimize function in optim.jl. I've run through a few examples and tests, but I can't seem to get it to work. The main issue is that all the functions that calculate $p_i$ rely on other variables. Specifically, my function is called pcalc(β,μ,E,cutoff). As you can see, I can't really cast this in a format similar to the examples in the Optim.jl documentation. I'd like to write it as pcalc(P,β) like is in the documentation, but it's really not feasible.

Currently, I have a sort of hessian type solver, but it takes too long to run, which is why I was hoping to implement Optim.jl.

Any advice or help will be appreciated, otherwise, I'm probably going to have to write my own optimizer function.


Solution

  • I have to admit I struggle a little bit with following your problem, but the gist of it seems to be that you have a function pcalc(β,μ,E,cutoff) which you want to optimize with respect to β, and are struggling to write this in a way that allows you to pass it to optimize which requires a single input.

    This is dealt with in the manual in the section tips and tricks: essentially you define an anonymous function or closure to optimize, i.e. you do

    optimize(x -> pcalc(x, μ, E, cutoff), ...)
    

    I don't understand how the summation part comes into this, but if you're saying there's an additional constraint on the solution to pcalc which says your solution for β should add to some value, you might want to look into JuMP which has special syntax for constraints, see e.g. here for an example. Alternatively you can hack this in Optim by doing something like

    optimize(x -> pcalc(x, μ, E, cutoff) + 1e6 * (sum(x) - 1.0)^2, ...)
    

    i.e. just adding a penalty term to your objective function.