Very slow JuMP + Knitro for nonlinear problem

[ulneva]

I started using Julia with Knitro through Julia's JuMP for mathematical programming. I noticed that a small nonlinear problem with 33 variables and just over 300 non-zeros in Hessian takes quite a lot of time to generate and solve. While slow generation can be attributed to Julia's JuMP, Knitro itself takes as much as 15 minutes to solve the above problem, though it takes only 12 iterations? Would you be able to have a look? I can send you the code.

The analogous problem in AMPL+KNITRO takes a few seconds to generate and solve.

[joehuchette]

Please do share the model and any other code needed to fully run the example.

Is this a mixed-integer problem? If not, could you try it with Ipopt as well for comparison?

[mlubin]

JuMP is a bit more sensitive than AMPL to the formulation of the nonlinear expressions. It's likely that with a small change you'll see much improved performance. Can't say more until we see the code.

[ulneva]

I'll try with Ipopt, thanks joehuchette. Here is the code and data. As mlubin mentioned, I suspect the issue might be indeed with the nonlinear expressions that I embed one into another. I also tried using auxiliary variables as follows. Put an aux variable into the nonlinear expression and then define that variable in the constraint (but not defining a nonlinear expression and then using it in the constraint, would it make a difference?). In the case with aux variables Knitro was showing a very heavy model, as JuMp does not do the substitutions that AMPL usually does and the resulting model as shown by Knitro was having way too many variables and non-zeros in J and H, and estimation was unusually slow as well. What would be the efficient way to go around the issue? Thank you for your help!

TestOptProblem7.txt levels.txt a_choices.txt mean_l2.txt gstate.txt joining.txt

[mlubin]

Lots of things going on here. In order of importance:

• Don't use data frames to index variables inside of nonlinear expressions. For example, in:

@defNLExpr(lind_def[ct=CT,i=N],sum{
log(cprob_m_def[ct,mact[t,i],lev[t,i],mstate[t,i],round(mean_l[t,1])])+
log(cprob_a_def[ct,ach[t,i],lev[t,i],mstate[t,i],round(mean_l[t,1])]),t in 1:606})

mact, lev, mstate, and mean_l are all data frames which provide no type information on the elements. The Julia compiler has no idea that ach[1,1] is an integer (look at @code_warntype ach[1,1]), so the code generated will be amazingly slow. You can read in the data using DataFrames if you want, but move it over to a julia Matrix{Int} or Matrix{Float64} before using it inside of JuMP expressions.

• JuMP variables (and NL expressions) should be declared as indexed over explicit ranges if possible. Instead of:

CT=[1:NS]
M=[0:1]
@defVar(llh,coef_k[CT,M,M])

write

@defVar(llh,coef_k[1:NS,0:1,0:1])

JuMP can generate faster code if it knows that variables are indexed over linear ranges instead of arbitrary sets.

• Separate out linear expressions. It looks like util_m_def and util_a_def are linear in the variables. It's better to declare these as variables and add linear equality constraints @addConstraint to define them.

[ulneva]

thank you, mlubin, for the prompt response. So the cascade of nonlinear expressions embedded one into another is not really a problem for JuMP+Knitro?

[mlubin]

I don't think that's the bottleneck here, I've seen JuMP models with much deeper expression chains. What does help though is separating out linear subexpressions as I mentioned.

[ulneva]

I am not sure, if the following issue is related to the speed, but the code generates a warning "Base.Nothing is deprecated, use Void instead" and then proceeds to estimation

[mlubin]

I can't reproduce that warning locally, it shouldn't be generated by any of JuMP's code.

[ulneva]

Thank you so much for the hints with coding! Interestingly, in my case, switching from data frames to arrays did not help that much, however using explicit indexing for variables (I also did it for all expressions) helped tremendously - the solution time went down from 15 min to 2 min and 45 s. This timing is still noticeably slower than with AMPL+KNITRO, though. I then tried to implement the 3d point of the hints (defining util_m_def and util_a_def as variables and adding linear expression) but it slowed down the generation step and pushed non-zeros in H to around 50,000 instead of 300 something and solution became very slow again. Please, let me know if something else might help and many thanks again! (Just in case, here is the updated code) TestOptProblem7_update.txt

[mlubin]

I get a 25% speedup by giving proper types to pseg and ssize:

pseg=Float64[1/NS for ct=1:NS,i=N]#probability of belonging to a consumer segment (unobs heterogeneity)
ssize=Float64[1/NS for ct=1:NS] #segment sizes

[mlubin]

Next point is to not use dictionaries for conditional values:

@defNLExpr(util_a_def[ct=1:NS,a=0:nA,l=2:nA,ms=0:nM,ml=4:by:20],sum{ coef_a0[ct]+coef_a1[ct,a]+(l-ml)*(sum{ coef_ce[ct,1] ; ml-l > 0}+sum{ coef_ce[ct,0] ; ml-l <= 0}) + (a-l)*( sum{ coef_g[ct,1]; a-l > 0} + sum{ coef_g[ct,0]; a-l <= 0}); a > 0})

We're working on some internal restructuring which will make these transformations less necessary in the future. For the time being, if the expressions use slow data structures then evaluation will be slow.

[ulneva]

Huge thanks, Miles!

[ulneva]

Unfortunately, I have to come back to discussing the speed issue. I added one more nonlinear expression while developing the model and the model generation and evaluations became unbearably slow again (it takes 15-20 min to generate the model and then 10-15 minutes per iteration). The new line is as follows (conditional definition using the sum{} approach). Before that, the model was running much faster than when I opened the issue, but still rather slow: out of the total 222.153 seconds spent on estimation, 221.805 spent in evaluations (in Julia/JuMP, I guess). But again, after the introduction of the new NLExpr it became much worse as described above. Would you have any hints on how the performance could be improved and on where a possible bottleneck is? I attach the updated code, the data files are as before, they are attached above.

Total program time (secs) = 222.153 ( 222.083 CPU time) Time spent in evaluations (secs) = 221.805

@defNLExpr(cprob_a_obs[p=1:4,ct=1:NS,a=0:maxa[p],l=2:maxa[p],ms=0:nM,ml=4:maxa[p]],sum{obsp_a[a]_cprob_a[p,ct,a,l,ms,ml];a>0}+sum{cprob_a[p,ct,0,l,ms,ml] + sum{(1-obsp_a[b])_cprob_a[p,ct,b,l,ms,ml],b in 1:maxa[p]};a==0}) TestOptProblem9.txt

[mlubin]

Try using IpoptNLSolver() from the AmplNLWriter package. If that works faster then you can use the same file-based interface to call KNITRO.

[ulneva]

As we chatted on the gmail groups, using AmplNLWriter is even slower, unfortunately, so I am back to using JuMP+Knitro and looking for way to speed up the code. Are there any approaches I can try to make it faster? It looks like it spends too much time in derivatives evaluations. Thank you so much!

[mlubin]

From a glance it looks like the problem is convex, is that correct?

[ulneva]

I believe it is convex at this point in the code development process, but as i introduce a dynamic programming part and will start optimizing with respect to value functions (treating them as variables to optimize over), it won't be any more

[mlubin]

Okay, I did a bit of digging into this. Even though there are only 58 variables, the nonlinear expressions involve summations with tens of millions of terms. JuMP doesn't currently save the intermediate values of expressions, so any time an expression is used its value is recomputed. (This should change with the NLP rewrite.) If you look at the definition of cprob_a, this means that the value of v_a_def gets computed a lot of times.

This might call for some clever reformulating, e.g., define a term for

sum{exp(v_a_def[p,ct,a_prime,l,ms,ml]),a_prime=0:maxa[p]}

and then pull it up a few levels out of the sum using the property exp(x+y) = exp(x)*exp(y). Then this term will be computed only once (per evaluation) instead of maxa[p] times. That should easily give a 10x speedup.

[ulneva]

Miles, thank you very much for the hint that JuMP does not remember the results for the intermediate expressions, I rewrote the code taking it into account, now runs much faster!

[mlubin]

Confirmed fixed by #638, because we will store the values of intermediate expressions. Before (30 iterations of ipopt):

Total CPU secs in IPOPT (w/o function evaluations)   =      2.220
Total CPU secs in NLP function evaluations           =    146.544

After:

Total CPU secs in IPOPT (w/o function evaluations)   =      0.676
Total CPU secs in NLP function evaluations           =      2.040

[ulneva]

Miles, thanks a lot for letting me know. Is the fix already available through JuMP package update? The fix will work with KNITRO as well, right?

Regards, Yulia

On Tue, Dec 29, 2015 at 6:27 PM, Miles Lubin notifications@github.com wrote:

Confirmed fixed by #638 https://github.com/JuliaOpt/JuMP.jl/pull/638, because we will store the values of intermediate expressions. Before (30 iterations of ipopt):

Total CPU secs in IPOPT (w/o function evaluations) = 2.220 Total CPU secs in NLP function evaluations = 146.544

After:

Total CPU secs in IPOPT (w/o function evaluations) = 0.676 Total CPU secs in NLP function evaluations = 2.040

— Reply to this email directly or view it on GitHub https://github.com/JuliaOpt/JuMP.jl/issues/614#issuecomment-167904414.

[mlubin]

No, it's not released nor ready to test yet, you can follow julia-opt for announcements. Yes, the behavior will be the same for KNITRO.