Welcome! It’s more helpful if you can provide a fully reproducible example that someone can copy-paste to see what the problem is. But if I’ve understood your question correctly, I think you’re looking for something like this (I haven’t run, so there might be typos etc):
M = 20
N = 25
X = rand(N, 100)
locations = rand(1:M, N)
location_lookup = [Int[] for j in 1:M]
for (i, j) in enumerate(locations)
push!(locations_lookup[j], i)
end
model = Model()
@variable(model, θ[i=1:100])
@variable(model, y[1:N])
@constraint(model, y .== X * θ)
@NLexpression(model, denom[j=1:M], 1 + sum(exp(y[i]) for i in location_lookup[j]))
Your main problem is that location is a vector with 25,000,000 elements! So performing that look-up is expensive, and you are doing this 20,000,000 times! In situations like this, it can be helpful to step back and ask if there are ways to do less work, e.g., by doing the lookup once and caching the answer, instead of trying to use threading or similar.