Support Nonlinear Expressions

[pulsipher]

Well after months of work, nonlinear has finally arrived! This implements our own nonlinear expression abstraction NLPExprs which use a left-child right-sibling expression tree that can store constants, variables, affine expressions, and/or quadratic expressions as leaves. Moreover, we support NLPExpr as a 1st class citizen in InfiniteOpt having nearly all the functionality of affine/quadratic expressions! Closes #16.

Early adopters can try this out via:

julia> ]

pkg> add InfiniteOpt#nlp

The preview of the relevant new nonlinear documentation is available here: https://pulsipher.github.io/InfiniteOpt.jl/previews/PR161/guide/expression/#nlp_guide

Some of the many features are highlighted below.

Operator Definition
Nonlinear expressions can be defined outside of macros like their quad/affine counterparts:

julia> using InfiniteOpt; m = InfiniteModel(); @variable(m, x);

julia> x^2.3 + sin(x)
x^2.3 + sin(x)

Macro Definition
By making the necessary extensions to MutableArithmetics, we can define NLPExprs inside of @expression, @objective, and @constraint (i.e., no need for the special NL macros):

julia> @expression(m, x + x^3 * sin(x)^x)
x + x^3 * sin(x)^x

Linear Algebra Support
We can also do linear algebra operations:

julia> @variable(m, X[1:3, 1:2]); @variable(m, v[1:3]); @variable(m, w[1:2]);

julia> @expression(m, v' * X * w)
0 + (X[1,1]*v[1] + X[2,1]*v[2] + X[3,1]*v[3]) * w[1] + (X[1,2]*v[1] + X[2,2]*v[2] + X[3,2]*v[3]) * w[2]

Function Tracing
Because we support operator definition, we can trace functions in like manner to Symbolics.jl:

julia> my_func(y) = 2^y + exp(y^-2.3);

julia> my_func(x)
2^x + exp(x^-2.3)

Deletion
Variables can be deleted from NLPExprs and constraints with NLPExprs can also be deleted:

julia> @constraint(m, con, sin(x) == 0)
con : sin(x) - 0 == 0.0

julia> delete(m, con)

User Defined Functions
For user-defined functions that are not compatible with tracing, we allow them to be registered like in JuMP and Symbolics.jl.

julia> using SpecialFunctions;

julia> my_func(a) = sin(a) * eta(a);

julia> @register(m, my_func(a)); # the gradient/hessian functions can also be given as additional arguments

julia> my_func(x)
my_func(x)

Note that it must be a macro since we have to overload the function behind the scenes to create NLPExprs.

[codecov]

Codecov Report

Merging #161 (f42e9a7) into master (b6be9e3) will increase coverage by 0.03%.
The diff coverage is 100.00%.

@@            Coverage Diff             @@
##           master     #161      +/-   ##
==========================================
+ Coverage   99.78%   99.81%   +0.03%     
==========================================
  Files          32       33       +1     
  Lines        6522     7149     +627     
==========================================
+ Hits         6508     7136     +628     
+ Misses         14       13       -1     

Impacted Files	Coverage Δ
src/InfiniteOpt.jl	100.00% <ø> (ø)
src/macros.jl	100.00% <ø> (ø)
src/Collections/VectorTuple.jl	100.00% <100.00%> (ø)
src/MeasureToolbox/expectations.jl	94.20% <100.00%> (ø)
src/MeasureToolbox/integrals.jl	100.00% <100.00%> (ø)
src/MeasureToolbox/support_sums.jl	100.00% <100.00%> (ø)
src/TranscriptionOpt/model.jl	100.00% <100.00%> (ø)
src/TranscriptionOpt/transcribe.jl	99.73% <100.00%> (+0.01%)	⬆️
src/array_parameters.jl	99.79% <100.00%> (ø)
src/datatypes.jl	100.00% <100.00%> (ø)
... and 13 more

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update b6be9e3...f42e9a7. Read the comment docs.

[sshin23]

@frapac might be of interest for experimenting next-generation JuMP nonlinear

[odow]

julia> @expression(m, x + x^3 * sin(x)^x)
x + x^3 * sin(x)^x

I'm a little in awe that you made this work. Have you done any benchmarking to see how this scales?

E.g., how does it handle

@expression(m, sum(x^i for i in 1:100_000))

cc @blegat @mlubin

[pulsipher]

I'm a little in awe that you made this work. Have you done any benchmarking to see how this scales?

E.g., how does it handle

@expression(m, sum(x^i for i in 1:100_000))

I benchmarked a lot of different configurations as this was developed and this approach was the most performant. Admittedly, it is not always as fast as the @NL macros since building the tree via operators means we have to start from the bottom up and not the top down. Also, there are certain cases where this is faster outside @expression (likely because of a combination of my limited understanding of MutableArithmetics and it being intended for quadratic-affine expressions).

It does handle large sums reasonably well:

julia> @time @expression(m, sum(x^i for i in 1:1000)); 
0.001774 seconds (25.96 k allocations: 1.129 MiB)

julia> @time @expression(m, sum(x^i for i in 1:10000)); 
0.013832 seconds (268.96 k allocations: 11.429 MiB)

julia> @time @expression(m, sum(x^i for i in 1:100000)); # we start to hit gc at about this point
0.167060 seconds (2.70 M allocations: 114.426 MiB, 18.01% gc time)

julia> f(n) = sum(x^i for i in 1:n);

julia> @time f(1000);
0.001436 seconds (25.96 k allocations: 1.129 MiB)

julia> @time f(10000);
0.014032 seconds (268.96 k allocations: 11.429 MiB)

julia> @time f(100000);
0.176825 seconds (2.70 M allocations: 114.426 MiB, 22.59% gc time)