Optim.jl : ParticleSwarm Multithreading at the function level

Specific Domains Optimization (Mathematical)
1

Hi,

I’m using the PSO algorithm in Optim.jl to minimise a certain loss function, which is a positive multinomial of very high degree (over a constraint domain, a product of several simplexes), and the optimisation is done in BigFloat precision.

The loss function itself consists of recursive computations that are not suited to parralelisation, so i thought I’ll parallelise at the Swarm level: after all, the Particle Swarm algorithm computes in parallel objectives values for many points, hence parralelisation can be done there.

I noticed that they were several issues on github talking about that :

I’ve been using Julia for only one week now so I’m not really sure of what’s the good way of doing what I want to do. My current plan is to update the first PR code to build my own version of Optim.jl and make julia work with it. The PR is old (~1year), and i guess that the stable version of Optim.jl has moved since.

Last but not least, I’d like the possibility to tweak myself some of the details of the ParticleSwarm implementation I’m using, e.g including a Scout phase as in Hasan Koyuncu & Rahime Ceylan’s ScPSO. How hard would it be to export through Optim.jl layers a tweakable Boolean to enable / disable this modification to the classical PSO ?

TL;DR:

  • Do you know about other implementation of a PSO that I missed (and that are pure-julia for BigFloat support) ?
  • How hard would it be to bring the old PR back to life ?
  • What documentation should i read before diving into making my own PR / testing the package / documenting? What are julia’s standard tools for testing / documentations / etc (again, only a few days since i started learning julia).

Thanks for any remarks you might have :slight_smile:

1 Like
2

Edit: sorry, I see that the PR mentioned above, and farther down in the thread suggests exactly what I wrote.

The basic idea of doing particle swarm in parallel seems perfectly feasible. A first stab at it could probably be made using Threads.@threads before the loop where the cost function is evaluated at each particle. In the current code in Optim.jl master, there is

function compute_cost!(f,
                       n_particles::Int,
                       X::Matrix,
                       score::Vector)

    for i in 1:n_particles
        score[i] = value(f, X[:, i])
    end
    nothing
end

I don’t see why you couldn’t just using threads for that for loop.

A second comment is that you might switch to BigFloat only in a final optimization, once the neighborhood of the global optimizer is located (supposing that can be successfully done with lower precision).

1 Like
3

Please open an issue for the scout phase, but I already did the rest for you :wink:

julia> using NLSolvers

julia> using ThreadsX

julia> using Base.Threads

julia> function himmelblau!(x)
           println(Threads.threadid())
           fx = (x[1]^2 + x[2] - 11)^2 + (x[1] + x[2]^2 - 7)^2
           return fx
       end
himmelblau! (generic function with 1 method)

julia> function himmelblau_batched!(F, X)
           ThreadsX.map!(himmelblau!, F, X)
           return F
       end
himmelblau_batched! (generic function with 1 method)

julia> x0 = big.([3.0,1.0])
2-element Array{BigFloat,1}:
 3.0
 1.0

julia> obj = ScalarObjective(himmelblau!, nothing, nothing, nothing, nothing, nothing,  himmelblau_batched!, nothing)
ScalarObjective{typeof(himmelblau!),Nothing,Nothing,Nothing,Nothing,Nothing,typeof(himmelblau_batched!),Nothing}(himmelblau!, nothing, nothing, nothing, nothing, nothing, himmelblau_batched!, nothing)

julia> prob = OptimizationProblem(obj, (x0.-1, x0.+1))
OptimizationProblem{ScalarObjective{typeof(himmelblau!),Nothing,Nothing,Nothing,Nothing,Nothing,typeof(himmelblau_batched!),Nothing},Tuple{Array{BigFloat,1},Array{BigFloat,1}},NLSolvers.Euclidean{Tuple{0}},Nothing,NLSolvers.InPlace,Nothing}(ScalarObjective{typeof(himmelblau!),Nothing,Nothing,Nothing,Nothing,Nothing,typeof(himmelblau_batched!),Nothing}(himmelblau!, nothing, nothing, nothing, nothing, nothing, himmelblau_batched!, nothing), (BigFloat[2.0, 0.0], BigFloat[4.0, 2.0]), NLSolvers.Euclidean{Tuple{0}}(), nothing, NLSolvers.InPlace(), nothing)

julia> solve(prob, x0, ParticleSwarm(), OptimizationOptions(maxiter=5))
6
2
5
1
3
1
1
2
6
4
1
1
6
2
4
3
5
1
5
4
6
3
4
1
5
4
6
2
3
1
5
3
4
6
3
1
Results of minimization

* Algorithm:
  Adaptive Particle Swarm

* Candidate solution:
  Final objective value:    1.14e-01

  Initial objective value:  1.00e+01

* Convergence measures

* Work counters
  Seconds run:   2.33e-01
  Iterations:    6

edit: by the way, just to show that it is indeed bigfloats


julia> solve(prob, x0, ParticleSwarm(), OptimizationOptions(maxiter=5)).info.X
5-element Array{Array{BigFloat,1},1}:
 [1.472835050193244429866528557757068401216554971469976874333357342070522307611754, 1.884851874326613485201770475071946518492050863729267730128441497983851351087877]
 [3.044800398162936776498703643433726991716450681598752750575706905503580284267865, 2.027991525583957853734345605408913381706412690691050521940793203231028729464288]
 [2.443207626133586676667050210022595623254037417956167541930171908405951768308145, 1.826333595882882102297845394332928078730587674111063119891896293103183632536135]
 [3.425374922497072526350988729574577628392837441317390961634126980687446459936601, 2.160115732289009819057009308557057457401693643132817943282186120744482283748051]
 [3.075501200630894081231377926805973854639330995303417744810659757900593837237966, 2.097209845779889338937521859608940285351207056338943956833012290980295280787417]

julia> solve(prob, x0, ParticleSwarm(), OptimizationOptions(maxiter=5)).info.minimum
0.01387799604048359887762015986372579531022827515101120461744260610507974173300061
1 Like
4

@lrnv I created that PR. As I mentioned in the comments there, I thought it would be best to just get something into Optim that made PSO actually feasible and on par with solutions in other languages. Seems I was not convincing enough.

I think I most recently used it perhaps 6 months ago. I believe it should probably work fine if you just checkout that branch if that is the only thing you are trying to get from Optim. The package did do some breaking releases though, so I can’t guarantee anything.

It does seem like Patrick gave you a solution above that doesn’t use Optim, perhaps give that one a shot.

5

Sorry if you felt like your contributed was not valued. I just think that once you go parallel you might as well decide how you’re doing it yourself. In the example I posted above I think it’s very easy to do especially with ThreadsX. I just need to adjust the ScalarObjective to be constructed from keyword args so you don’t need all the nothings.

NLSolvers will be the backend code for Optim sometimes in the beginning of next year when I am again free to do the final push.

6

No I don’t feel slighted, but like I said, a PSO implementation that isn’t parallelized in any fashion just isn’t competitive. I provided the bare minimum to get one that was.

I’m not sure that things like ThreadsX even existed when I made that PR, so it is quite possible that there are much better ways to do it now than there were a year ago.

7

I will ASAP. Thanks for sharing this piece of code ! Indeed, NlSolvers is not yet documented and i understood that you were in the process of replacing optim with it. Very nice work btw.

8

That is correctly understood. I cannot promise stability of the API either, but I decided to publish it anyway, and I will follow semver in terms of breaking changes of course. I just wanted to point you to an implementation you could actually use.

9

And you did, thanks. Maybe one last question : running your code, I have no output(appart for the threadid). I tried adding show_trace=true in the Optimizationoption() call but no success, I think I missed the right parameter. Could you give me a hint ?

10

Yeah, that’s one of the places where NLSolvers.jl is still rough around the edges… I have not yet implemented progress printing. You could try keeping a variable of lowest value found so far and print it out on an objective call.