Global sparsity detection with control flow

[gdalle]

We need to handle both branches of ifelse so that global sparsity patterns may be reused at different points x

Maybe with a "branch analysis" option

[gdalle]

Example where the sparsity pattern depends on whether x[1] <= 0:

https://github.com/JuliaSmoothOptimizers/OptimizationProblems.jl/blob/main/src/ADNLPProblems/AMPGO07.jl

The trouble is that we need to explore both branches to give a global pattern that holds for different values of x

[adrhill]

We can support ifelse, but I don't yet see how we could handle direct comparisons x <= 0 and hit all branches of the program using just operator overloading.

[adrhill]

The real issue is having to decide what comparisons like tracer > 42 return on non-dual tracers. To support ifelse, like I did in #94, the answer has to be a tracer. Otherwise we can't overload and dispatch on ifelse.

However, having comparisons return tracers could lead to all kinds of problems, including code invalidations in Julia Base. And without comparisons on non-dual tracers, control flow is basically useless.

Note that all of this is just to support ifelse. Regular if-else-blocks can't be supported via operator overloading.

[gdalle]

We can support ifelse, but I don't yet see how we could handle direct comparisons x <= 0 and hit all branches of the program using just operator overloading.

Indeed we can't, but supporting ifelse is enough in some cases for people to rewrite their code in a compatible way.

The real issue is having to decide what comparisons like tracer > 42 return on non-dual tracers. To support ifelse, like I did in #94, the answer has to be a tracer. Otherwise we can't overload and dispatch on ifelse.

I actually had a different overloading scheme in mind:

• <(x::Tracer, y::Tracer) = missing
• ifelse(::Missing, x::Tracer, y::Tracer) = union(x, y)

However, having comparisons return tracers could lead to all kinds of problems, including code invalidations in Julia Base.

That is indeed a valid worry, and it goes against the specification of comparison operators:

isless(x, y)


  Test whether x is less than y, according to a fixed total order (defined
  together with isequal). isless is not defined for pairs (x, y) of all types.
  However, if it is defined, it is expected to satisfy the following:

    •  If isless(x, y) is defined, then so is isless(y, x) and isequal(x,
       y), and exactly one of those three yields true.

    •  The relation defined by isless is transitive, i.e., isless(x, y)
       && isless(y, z) implies isless(x, z).

  Values that are normally unordered, such as NaN, are ordered after regular
  values. missing values are ordered last.

  This is the default comparison used by sort!.

  Implementation
  ≡≡≡≡≡≡≡≡≡≡≡≡≡≡

  Non-numeric types with a total order should implement this function. Numeric
  types only need to implement it if they have special values such as NaN.
  Types with a partial order should implement <. See the documentation on
  Alternate Orderings for how to define alternate ordering methods that can be
  used in sorting and related functions.

[gdalle]

https://discourse.julialang.org/t/how-bad-is-it-if-x-y-doesnt-return-a-bool/114611

[gdalle]

I think your way of returning a tracer would be better, because it also allows the quantity x < y to be used as a number, like they do in https://github.com/JuliaSmoothOptimizers/OptimizationProblems.jl/blob/main/src/ADNLPProblems/AMPGO07.jl

[adrhill]

https://discourse.julialang.org/t/how-bad-is-it-if-x-y-doesnt-return-a-bool/114611

I mirror a lot of the worries that have been expressed in here.
I really view this package as binary second-order ForwardDiff and not as Symbolics.jl.

I suggest we try out overloading < & friends in #94 and see whether this fixes issues in ADNLPProblems or whether they actually require abstract interpretation code like Symbolics.

[adrhill]

I suggest we try out overloading < & friends in #94 and see whether this fixes issues in ADNLPProblems

Done. AMPGO07.jl now passes, but it's a trivial scalar-in-scalar-out function. Would be interesting to see how many other problems in #83 this solves @gdalle.

[gdalle]

I mirror a lot of the worries that have been expressed in here.

I think the distinction between ForwardDiff and Symbolics is arbitrary. We can do just as much damage with operator overloading (and to some extent that's also how Symbolics works). The question is only how much we want to overload.

whether they actually require abstract interpretation code like Symbolics.

As you see in the Discourse thread, "abstract interpretation" is ill-defined. We are already doing abstract interpretation, by reinterpreting scalar operations to work on index sets. So it's not a good concept to use for clarification purposes.

[gdalle]

Would be interesting to see how many other problems in #83 this solves @gdalle.

Ideally we want all of the problems in the OptimizationProblems.jl test suite to be amenable to global sparsity detection, because the sparsity pattern must be reusable for various x.
That might require a bit of rewrite on their end, but nothing that goes out of the Base syntax. Typically, we don't want to introduce a new comparison operator that is specific to SCT, and that people would have to import for their code to work.

At the moment we don't see these problems because I fall back on local tracing every time. The real test would be to remove these fallbacks in #83.

[gdalle]

Upsides:

• we allow global sparsity detection with some control flow
  • they have to rewrite if blocks as ifelse

Downsides:

• invalidations
• slow instead of erroring (meh)
• less informative error (TypeError with failed Bool conversion instead of nice MethodError)
  • register an error hint

[gdalle]

The result of the comparison should be a tracer, at least for ConnectivityTracer, because:

• for connectivity, the influencing indices of ifelse(b, x, y) are union(b.inputs, x.inputs, y.inputs)
• for gradient, the gradient of ifelse(b, x, y) is union(x.gradient, y.gradient)

In a way, ifelse is just another ternary operator. We just don't have any other one