limited support for embedded nonlinear expressions

[mlubin]

Needs docs and a bunch of testing.

The current limitation is that nonlinear subexpressions cannot appear inside of a sum{} or prod{} term, because we would need to add in a lot of new infrastructure to handle this case and ensure that subexpressions are homogeneous.

CC @tkelman

[tkelman]

Interesting. So could these expression objects be made to propagate symbolically through user-defined nonlinear functions? Like the Unknown() type in Sims.jl?

[mlubin]

We're not planning on implementing operator overloading for these in the short term, that would be a very big change. ReverseDiffSparse is too dependent on homogeneity of expressions (since a new derivative evaluation function has to be compiled for each unique expression graph) that it would be horribly inefficient to run with 1000 constraints each with its own expression graph derived from operator overloading. Once we have a few MPB "solvers" which take in ExprGraph input, it would probably be worthwhile to experiment with this idea outside of JuMP.

[tkelman]

I do expect to re-apply the same expression as a template across multiple index sets of variables (propagating dynamics over a time horizon, in the case of optimal control). But it's a lot to ask to expect the user to understand how to do that, and writing things this way is useless outside of JuMP (people will want to use the same nonlinear functions for forward simulation, etc). If a user provides a function as input, I want to apply that function to generate the constraint expression that should be applied multiple times over different index sets. Should be doable outside of JuMP, but people really want to be able to do @addNLConstraint(m, f(x) == 0) where f is defined somewhere else.

[mlubin]

I definitely recognize the utility of @addNLConstraint(m, f(x) == 0), but from my understanding of the literature and from discussions with AD people, applying reverse-mode AD transparently to a somewhat general user-defined function all while exploiting sparsity of second-order derivatives is beyond what exists in the state of the art and would almost certainly earn someone a PhD thesis.
My impression though is that it wouldn't be all that difficult in Julia, but it seems more appropriate to experiment with it outside of JuMP first, and we've intentionally designed MPB so that there's not a big startup cost for writing your own modeling interface.

[tkelman]

I guess therein lies the novelty in projects like casadi and acado, of being able to apply specific structural information about a class of problem to make it actually tractable. Or punt on second-order derivatives and prioritize convenience over performance the way yalmip does.

With appropriate overloads on an expression placeholder type, how would @addNLConstraint(m, f(x) == 0) where f(x) returns an expression for expression inputs be any different than manually writing out the operations performed in f(x)? Assuming model generation time/memory allocation doesn't really matter and is an offline cost for some application areas.

[IainNZ]

I tried to do

@defNLExpr(grav[j=0:n], g_0*(h_0/h[j])^2)

and got The value of j has changed unexpectedly

[mlubin]

@IainNZ, that happens when you reuse the index j in both the definition of the expression and in the constraint.

[mlubin]

I'll fix that by keeping each expression in its own namespace.

[mlubin]

@IainNZ, try again with ReverseDiffSparse master.
We could definitely use a fuzzer for this.

With appropriate overloads on an expression placeholder type, how would @addNLConstraint(m, f(x) == 0) where f(x) returns an expression for expression inputs be any different than manually writing out the operations performed in f(x)?

@tkelman, the only difference is that it requires a complete rewrite of JuMP's internal infrastructure. We wouldn't use ReverseDiffSparse in that case but maybe NL or osil.

[tkelman]

Makes sense that JuMP is designed for macros rather than operator overloading. I could probably write a simple scalar overload-based incremental osil modeling layer (I did the same in Matlab a while ago), but I don't trust OS' AD at this point. I also wouldn't be constructing Expr objects in that case, I'd probably just build up the XML tree in-memory as I go.

Why wouldn't ReverseDiffSparse still work? Or do you just expect it to scale more poorly than ASL/CppAD for generic, large expressions?

[mlubin]

@JackDunnNZ is looking at writing an NL writer based on expression graph output, so it would pay to go through MPB instead of tying yourself to OS.

My sense is that ReverseDiffSparse would not even work for these problems. If your expression graph has a flattened out sum with 10,000 terms, ReverseDiffSparse would try to generate and compile a function with 10,000+ lines. I did a few experiments with this before, and it crashes and burns. You would have to rewrite the AD to interpret the expression graph instead of compiling functions, which is what most AD tools do anyway.

[tkelman]

Any updates there @JackDunnNZ? If I get too frustrated/impatient with OS I might have to dig out the "Writing NL files" technical report again and actually try to understand it this time. ampl/mp#30 also may become relevant.

If your expression graph has a flattened out sum with 10,000 terms, ReverseDiffSparse would try to generate and compile a function with 10,000+ lines. I did a few experiments with this before, and it crashes and burns. You would have to rewrite the AD to interpret the expression graph instead of compiling functions, which is what most AD tools do anyway.

Ah, thanks. ReverseDiffSparse is a lot less like conventional AD tools than I thought it was. I do generally expect my problems to be a lot sparser than that, so maybe it wouldn't be completely impossible? I guess ReverseDiffSource would be the more traditional-style choice if staying in pure Julia, but has more limitations and less solver upsides than NL/osil would.

[JackDunnNZ]

I've been busy this week, but finally got some time today to work on the NL writer (thanks for the reminder!). It's a bit rough in parts, but is working and solving correctly for the few problems I've tested including a proper MINLP. It still needs a lot more tidy up and thorough testing, but glad it's working for now :)

[tkelman]

@JackDunnNZ cool, that's looking pretty good so far. We'd need to be sure we build the various solvers with ASL support, so the command-line executables get installed. And it'll need a longer more specific name once the package is ready to be registered.

[mlubin]

Nice @JackDunnNZ!

I've added documentation for this with some general hints for nonlinear modeling syntax. I think we can just merge this since it shouldn't break anything that already works.

[mlubin]

Also unfortunately the model I wanted to translate needs expressions inside of sum{}.