Simplify NonlinearExpr with + or * and 1 or 0

[odow]

One possible fix for #3498

@mitchphillipson could you try out this branch? (] add JuMP#od/nlp-simplify).

It'd be useful if you had the test cases using sparse matrices.

[codecov]

Codecov Report

Patch coverage: 95.00% and project coverage change: -0.02% ⚠️

Comparison is base (8bcea81) 98.13% compared to head (aa0d10a) 98.12%.
Report is 1 commits behind head on master.

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #3502      +/-   ##
==========================================
- Coverage   98.13%   98.12%   -0.02%     
==========================================
  Files          37       37              
  Lines        5536     5555      +19     
==========================================
+ Hits         5433     5451      +18     
- Misses        103      104       +1     

Files Changed	Coverage Δ
src/nlp_expr.jl	99.31% <95.00%> (-0.21%)	⬇️

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[mitchphillipson]

I'm still testing this, but have observed the following:

julia> using JuMP

julia> M = Model();

julia> A = [1, 0, 3]

julia> @variable(M, y[1:3])

julia> sum(y[i]^10*A[i] for i=1:3) #not as expected
(((y[1] ^ 10) * 1.0) + ((y[2] ^ 10) * 0.0)) + ((y[3] ^ 10) * 3.0)

julia> @expression(M, tmp, sum(y[i]^10*A[i] for i=1:3)) #As expected
(y[1] ^ 10) + ((y[3] ^ 10) * 3.0)

julia> T = sum(y[i]^10*A[i] for i=1:3)
(((y[1] ^ 10) * 1.0) + ((y[2] ^ 10) * 0.0)) + ((y[3] ^ 10) * 3.0)

julia> @expression(M, tmpT, T)
((y[1] ^ 10) * 1.0) + ((y[3] ^ 10) * 3.0) + ((y[2] ^ 10) * 0.0)

The simplification only takes hold inside an @expression. And as a second check, I constructed the expression first then placed it inside an @expression, no simplification.

I'm going to test this more tomorrow morning.

[odow]

The simplification only takes hold inside an @expression

Correct

[odow]

I don't really know how I feel about this. Is there any evidence that keeping the +0 and *0 is a problem? Do we just want to remove them because it looks nice?

The overhead of walking every expression just to see if there is a +0 might be larger than continuing with one. One problem here is that it makes * unstable, because it might return 0.0 or it might return a nonlinear expression.

[mitchphillipson]

I'm tending towards agreeing that it's mainly aesthetic and probably adds unnecessary overhead to model creation. I'm fine if this gets shelved. I may build a more specific solution for my use case.

…

On Mon, Sep 11, 2023, 8:56 PM Oscar Dowson ***@***.***> wrote: I don't really know how I feel about this. Is there any evidence that keeping the +0 and *0 is a problem? Do we just want to remove them because it looks nice? The overhead of walking every expression just to see if there is a +0 might be larger than continuing with one. One problem here is that it makes * unstable, because it might return 0.0 or it might return a nonlinear expression. — Reply to this email directly, view it on GitHub <#3502 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AEVAPXKDDMCG7RCHTTOJCN3XZ66NZANCNFSM6AAAAAA4UBXKQU> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[odow]

I do think this is a nice to have, and I can see people requesting it in future.

Maybe I should try what you suggested, and add drop_zeros!, which wouldn't be type stable and would be opt-in by users.

[mitchphillipson]

Having an optional function would be nice. 90% of the time that operation isn't necessary, but it's nice to have for that 10%. The type instability is intriguing to me. I'm guessing this comes from something like $0*x^5+x^2$ becoming quadratic after `drop_zeros!` is applied. Is there an issue forcing it to remain nonlinear?

…

On Mon, Sep 11, 2023, 9:23 PM Oscar Dowson ***@***.***> wrote: I do think this is a nice to have, and I can see people requesting it in future. Maybe I should try what you suggested, and add drop_zeros!, which wouldn't be type stable and would be opt-in by users. — Reply to this email directly, view it on GitHub <#3502 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AEVAPXPZYMNY2DZ6QIU5LKTXZ7BQZANCNFSM6AAAAAA4UBXKQU> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[odow]

The type instability is intriguing to me.

It would come from a definition like this:

Base.:*(x::Real, y::NonlinearExpr) = iszero(x) ? x : NonlinearExpr(:*, Any[x, y])

Now in a = 2.0; y = a * sin(x) Julia will infer y as Union{Float64,NonlinearExpr}.

[odow]

So one option would be:

_is_zero(x::Real) = iszero(x)
_is_zero(::Any) = false

_is_one(x::Real) = isone(x)
_is_one(::Any) = false

simplify(x) = x

simplify(x::Union{GenericAffExpr,GenericQuadExpr}) = drop_zeros!(x)

function simplify(x::GenericNonlinearExpr{V}) where {V}
    if x.head == :* && any(_is_zero, x.args)
        return zero(value_type(V))
    end
    for i in 1:length(x.args)
        x.args[i] = simplify(x.args[i])
    end
    if x.head == :*
        filter!(!_is_one, x.args)
    elseif x.head == :+
        filter!(!_is_zero, x.args)
    end
    if x.head == :+ || x.head == :* && length(x.args) == 1
        return x.args[1]
    end
    return x
end

So it's a fairly simple feature addition. But I think it could wait for a bit more usage. It feels like we might want a much more robust expression rewriting system that deals with more than just + 0, * 0, and * 1.

[mitchphillipson]

simplify resonates with me, but I would want it to be more robust. My background in Mathematics so my mind always goes to CAS, and I ask "What would Maple/Mathematica do with simplify?" The answer is that it would do much more. However, the "much more" is complicated and requires extensive manipulation of the expression tree and various heuristics of what is "simple".

The Symbolics.jl (or SymbolicUtils.jl) package has code that does this, simplify.jl. But then they have a helper, simplify_rules.jl, that has a long list of identities and special algebra rules. This is a lot to add right now, especially for cosmetics.

[odow]

Yeah the expectations that simplify might be CAS-like is what is stopping me from adding it, and why I called the other function flatten! instead of simplify.

I think we can wait for the current stuff to be out in the world for a while before we decide to add anything.

[odow]

I'm going to close this PR for now, but I'll leave the issue open.