1

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.

TBA
  • 1,921
  • 4
  • 13
  • 26
Casey Sobecks
  • 43
  • 2

2 Answers2

1

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.

Nils Gudat
  • 13,222
  • 3
  • 39
  • 60
  • The key here is, as Nils said, you need to define a function (he used an anonymous one which is convenient). The key addition is that the values of \mu and E and such will be taken from the environment of the call to optimize. This presumes that these are constants rather than, say, functions of \beta. – Ted Dunning Nov 06 '21 at 21:12
  • Okay this clears up a lot, I was looking at the Optim documentation, but I seem to have passed over the part where you can pass other variables to the function so thank you. As for the summation, I was able to get it working, I just needed to cast it into the form that optimize is looking for. So essentially, it ends up being a for loop with the optimization function inside. Again, thank you for your help – Casey Sobecks Nov 07 '21 at 20:19
  • I was able to get it working now. Thank you, your comment helped out a lot. – Casey Sobecks Nov 21 '21 at 06:13
1

This is also a great application for LsqFit.jl, which solves exactly for sums of squares where the elements may be interrelated like in your case. I've found it useful in "generalized method of moments" problems where I also want to impose bounds for beta. (The alternative is Optim.jl with FminBox but I have found that to be very slow.)

https://github.com/JuliaNLSolvers/LsqFit.jl

Similar to Optim.jl you want to define a function that takes both beta and your data, so that MYFUNCTION(ALLDATA, beta) returns a vector where the i'the element is p_i-\avg(n_c)_i. Then call

fit = curve_fit(MYFUNCTION, ALLDATA, zeros(N), beta0)

Are there other parameters you need to estimate, that p_i depends on? Then add them alongside beta in the second argument (in a vector).

GEK
  • 61
  • 5