Adding infinity norm to VectorBase

[edrumwri]

I want to add an infinity norm to VectorBase:

  virtual T NormInf() const {
    using std::abs;
    T nrm(0);
    const int sz = size();
    for (int i=0; i< sz; ++i){
      T val = abs(GetAtIndex(i));
      if (val > nrm)
        nrm = val;
    }
    return nrm;
  }

The compiler chokes on this when building the basic_vector_test, complaining that '>' comparison can't be applied to FunctionalForm and 'abs' can't be applied to Polynomial (god forbid I try to use the max function).

[david-german-tri]

From Slack discussion: This can probably be resolved with SFINAE on the drake::is_numeric trait

[sherm1]

I don't think it's a good idea for us to use SFINAE in routine code. Code like Evan wrote above ought not to be controversial! We are writing mountains of numerical code that is going to be full of basic numeric operations -- shouldn't all legal types T allow us to write such code unimpeded? It would be fine if they want to throw a runtime error under some circumstances, but as it is the code won't compile if we instantiate it with the less-numerical types. That seems like an unreasonably large tax on Drake programmers.

/cc @RussTedrake, @bradking, @jwnimmer-tri (please comment)

[david-german-tri]

This is a proposed addition to the framework; it's by definition not routine, and the work that goes into it can't reasonably be characterized as a "tax on Drake programmers". If Drake users want to implement a NormInf on some specific vector outside of the framework, and they're comfortable assuming/requiring that type will never be instantiated on a Scalar that is not part of the real line, nothing is requiring them to use SFINAE.

[sherm1]

Evan could easily have decided to compute that function locally rather than add it as a VectorBase method (although NormInf does seem like a useful addition here). In that case it would have broken in his garden-variety numerical code. Won't these problems crop up repeatedly in even the most mundane numerical code, for example in some limited or one-sided force element where a max or abs function is needed? Our plan at the moment as I understand it is to pre-instantiate all the allowable T variants in the Drake libraries -- that would mean all that mundane stuff is going to get instantiated with FunctionalForm and Polynomial even if no one wants to use it that way.

[david-german-tri]

In that case it would have broken in his garden-variety numerical code. Won't these problems crop up repeatedly in even the most mundane numerical code, for example in some limited or one-sided force element where a max or abs function is needed?

Not if that code is only instantiated for real-valued scalars, as would be appropriate.

Our plan at the moment as I understand it is to pre-instantiate all the allowable T variants in the Drake libraries

That hasn't been decided. It's currently muddled and controversial. I have been meaning to try to clarify it, and I guess this thread is another good reason to start. I'll open a separate issue.

Even if we do have strict pre-instantiation, I do not believe that every Drake library has to support every Drake scalar type. (For instance, Simulator makes no sense on non-real-valued types, as far as I can tell!) However, deep-core libraries like VectorBase certainly should.

[jwnimmer-tri]

TL;DR: We should add min and max as required scalartype operations, assuming that Autodiff supports them. Then you can write NormInf on VectorBase in a straightforward way.

(1) I concur with what David has written so far -- the SystemsFramework library is central to systems design and analysis, so must support the maximal variety of scalartypes that Drake is involved with. That does not mean that the Simulator library (or other libraries) needs to also support the same wide variety of scalartypes.

(2) It is valid to ask whether our scalar types for SystemsFramework should be required to sanity-check at compile-time, or instead merely throw exceptions at runtime -- that is related to the recent discussions of whether our scalartypes should carry along some "nominal" double value like Autodiff does, so they can flow through arbitrary branching, loops, or bespoke other double-related operations. I do not think that NormInf is a good motivating example for that discussion.

(3) FunctionalForm is a valid concept for NormInf -- it should not throw an exception. The result should be like abs -- a vector of Constants stay constant, a single Arbitrary or Unknown trumps, but otherwise it's just Differentiable (which means mostly-differentiable, per the FunctionalForm docs). Making this happen would either require (3.1) an enable_if for FunctionalForm on NormInf with a hand-crafted implementation, or else (3.2) phrasing the implementation in terms of abs (already supported), max (should arguably be overloaded already, with similar semantics to abs), and perhaps a std::accumulate-like function that respects our scalartypes (or else just write a loop by hand in NormInf). I prefer (3.2), because I think accumulate-with-min/max will be a common pattern in different Systems, and across different scalartypes.

(4) It's been suggested that we replace Polynomial with SymbolicExpression in SystemsFramework (and systems code in general) as the symbolic-expressive scalartype. We can still have a System<SymbolicExpression> and ask it for its EvalOutput expressions, and then (at runtime) if that expression tree has polynomial form, we can realize the congruent Polynomiald expression and efficiently use that for mathematical programs and optimization problems. I think this is generally a good and widely-supported idea, but we haven't done it yet because SymbolicExpression has only just now appeared. While SymbolicExpression doesn't support everything needed for NormInf yet, I believe it would be in-scope to have it support max and thus be closed under accumulate-with-max-abs as is needed for NormInf.

(Edited to fix: s/Abstract/Symbolic/g typo in the first draft.)

[david-german-tri]

(3) FunctionalForm is a valid concept for NormInf -- it should not throw an exception.

This is a great point, and makes the rest of the discussion moot IMO.

..., but otherwise it's just Differentiable (which means mostly-differentiable, per the FunctionalForm docs).

For a while, I struggled with whether the nondifferentiable points in NormInf/max are really measure zero for all vectors of Differentiable, but I'm pretty sure I believe it now. Cool!

[bradking]

FunctionalForm is a valid concept for NormInf -- it should not throw an exception. ...max (should arguably be overloaded already, with similar semantics to abs)

This is a great point, and makes the rest of the discussion moot IMO.

So FunctionalForm should gain min and max operations?

[david-german-tri]

For the record, here's the argument I've used to convince myself. I'm not a mathematician and this is probably garbage.

As Jeremy pointed out above, if the max of Differentiable functions is Differentiable, then so is the NormInf, so let's just consider max.

Let {f_1(x)...f_n(x)} be a set of almost-everywhere-differentiable functions from R^N to R^1. The max of these functions is only nondifferentiable at the points where one functions themselves is nondifferentiable (which are measure zero by the premises), and at some subset of the points where one function becomes greater than another, which are {Q : (f_i(q) == f_j(q)) && (f_i'(q) != f_j'(q))} for all (i, j) in n X n. By the definition of derivatives, q + ε is not in Q for any |ε| > 0, because f_1(q + ε) != f_2(q + ε). Therefore, Q is measure zero, and the union of measure zero sets is still measure zero, so the set of points where max is nondifferentiable is a subset of a measure zero set, so max is almost-everywhere-differentiable, QED.

[david-german-tri]

So FunctionalForm should gain min and max operations?

Yep, I think so. @edrumwri, would you agree that solves the original problem?

[SeanCurtis-TRI]

Generally, I don't believe the max of a vector is generally differentiable. Even if each of the individual components varies smoothly, when one value pops up higher than the previous maximum value, you get a discontinuity in the rate of change of the max; in other words, at that moment, the derivative coming from negative infinity has a different value than the derivative coming from positive infinity. So, the max of a smooth vector is only C0 continuous. Am I missing something?

[david-german-tri]

Am I missing something?

Only that Differentiable, for FunctionalForm purposes, is defined as "almost-everywhere-differentiable". Those points where one value pops up higher than the other are measure zero. (I think. See above.)

[SeanCurtis-TRI]

Again, my point is not about the data type that goes into max, but the max function itself. Even if the inputs were infinitely smooth, max itself is only C0 continuous. As long as that fact is accounted for, we're good. But if we expect it to be differentiable we'll be disappointed.

[david-german-tri]

It seems like we're talking past each other, but I'm not sure what the disconnect is. I agree, max is not everywhere-differentiable on smooth inputs, but it is almost-everywhere-differentiable, which is what matters here.

[SeanCurtis-TRI]

I was only responding to the quote:

As Jeremy pointed out above, if the max of Differentiable functions is Differentiable, then so is the NormInf, so let's just consider max.

I interpreted "Differentiable" as differentiable and was pointing out that the max function is not going to be differentiable; so the dependent clause is false, rendering the conclusion moot. If "Differentiable" means piecewise differentiable, then I'm good.

[edrumwri]

Aside from @SeanCurtis-TRI's point (mathematically, I don't expect the max operator to be differentiable either, but I don't feel that I'm yet sufficiently aware of the programming implications), what about comparison operators and the abs operation? I'm concerned that the infinity norm is just the tip of the iceberg for vector arithmetic operations (i.e., that we'll have to keep revisiting this discussion).

[david-german-tri]

You don't need comparison for NormInf if you have min and max, right?

As for Polynomial, @jwnimmer-tri is suggesting that we can drop support for it in drakeSystemsFramework, which sounds reasonable to me. Feel free to delete the one test case in drake/systems/framework/test/basic_vector_test.cc that exercises it.

[edrumwri]

I'd need abs. I could make do with comparison and negation. Still concerned that we'll be revisiting similar issues repeatedly.

[jwnimmer-tri]

We already have abs throughout (just not Polynomial, which is why we are dropping it).

[edrumwri]

I see. Yes, this current problem will be resolved if I have a max operator.

[bradking]

Please see PR #3851 for an implementation of max and min for FunctionalForm.

[edrumwri]

Ok, am now getting a compile problem due to missing operator< (!!!) in Expression for symbolic::Expression.

In file included from /usr/include/c++/5/memory:62:0, from /home/drum/drake/drake/../drake/systems/framework/primitives/gain.h:3, from /home/drum/drake/drake/../drake/systems/framework/primitives/gain-inl.h:8, from /home/drum/drake/drake/systems/framework/primitives/test/gain_scalartype_test.cc:1: /usr/include/c++/5/bits/stl_algobase.h: In instantiation of ‘constexpr const _Tp& std::max(const _Tp&, const _Tp&) [with _Tp = drake::symbolic::Expression]’: /home/drum/drake/drake/../drake/systems/framework/vector_base.h:139:16: required from ‘T drake::systems::VectorBase::NormInf() const [with T = drake::symbolic::Expression]’ /home/drum/drake/drake/systems/framework/primitives/test/gain_scalartype_test.cc:140:1: required from here /usr/include/c++/5/bits/stl_algobase.h:224:15: error: no match for ‘operator<’ (operand types are ‘const drake::symbolic::Expression’ and ‘const drake::symbolic::Expression’) if (a < b) ^ In file included from /usr/include/c++/5/bits/stl_algobase.h:64:0, from /usr/include/c++/5/memory:62, from /home/drum/drake/drake/../drake/systems/framework/primitives/gain.h:3, from /home/drum/drake/drake/../drake/systems/framework/primitives/gain-inl.h:8, from /home/drum/drake/drake/systems/framework/primitives/test/gain_scalartype_test.cc:1: /usr/include/c++/5/bits/stl_pair.h:220:5: note: candidate: template<class _T1, class _T2> constexpr bool std::operator<(const std::pair<_T1, _T2>&, const std::pair<_T1, _T2>&) operator<(const pair<_T1, _T2>& x, const pair<_T1, _T2>& y) ^ /usr/include/c++/5/bits/stl_pair.h:220:5: note: template argument deduction/substitution failed: In file included from /usr/include/c++/5/memory:62:0, from /home/drum/drake/drake/../drake/systems/framework/primitives/gain.h:3, from /home/drum/drake/drake/../drake/systems/framework/primitives/gain-inl.h:8, from /home/drum/drake/drake/systems/framework/primitives/test/gain_scalartype_test.cc:1: /usr/include/c++/5/bits/stl_algobase.h:224:15: note: ‘const drake::symbolic::Expression’ is not derived from ‘const std::pair<_T1, _T2>’ if (a < b)

[jwnimmer-tri]

The the new code calling max (using ADL, correct) or std::max (only defined for std types, so wrong)? I believe the NormInf body should look like using std::max; for ... { foo = max(foo, bar); } .... Happy to help with gists / slack if I can.

[edrumwri]

Yes, I should have updated the fragment at the top to match my current code. Apologies!

  virtual T NormInf() const {
    using std::abs;
    using std::max;
    T nrm(0);
    const int sz = size();
    for (int i=0; i< sz; ++i) {
      T val = abs(GetAtIndex(i));
      nrm = max(nrm, val);
    }

    return nrm;
  }

[jwnimmer-tri]

Aha. I think that code is correct, but we haven't added max to SymbolicExpression yet. I lost that tidbit in the upthread discussion. Sorry!

Here is a patch that gets you compiling:

--- a/drake/common/symbolic_expression.h
+++ b/drake/common/symbolic_expression.h
@@ -211,6 +211,7 @@ class DRAKE_EXPORT Expression {
   friend DRAKE_EXPORT Expression sinh(const Expression& e);
   friend DRAKE_EXPORT Expression cosh(const Expression& e);
   friend DRAKE_EXPORT Expression tanh(const Expression& e);
+  friend DRAKE_EXPORT Expression max(const Expression& e1, const Expression& e2) { return e1; }

   friend DRAKE_EXPORT std::ostream& operator<<(std::ostream& os,
                                                const Expression& e);

@soonho-tri Would you like to prioritize a max / min feature patch to SymbolicExpression? Otherwise I could get to it quickly too.

[soonho-tri]

Ok, am now getting a compile problem due to missing operator< (!!!) in Expression for symbolic::Expression.

There is operator<(const symbolic::Expression& e1, const Symbolic::Expression& e2), but its return type is not bool but symbolic::Formula.

Both of symbolic::Expression and FunctionalForm represent functions which can be evaluated later. Let's say that you have e1 = x + y and e2 = x - y, without having assignments on x and y, it's meaningless to ask a Boolean value for e1 < e2. In case of FunctionalForm, we deleted relational operations (<, <=, ==, !=, >, >=) over FunctionalForms to prevent this problem. In case of symbolic::Expression, we introduced symbolic::Formula to represent the results of relational operations over symbolic expressions and to evaluate the results afterwards with assignments on variables (it's called symbolic::Environment, a mapping from variable to value(double for now)).

[soonho-tri]

@soonho-tri Would you like to prioritize a max / min feature patch to SymbolicExpression? Otherwise I could get to it quickly too.

Sure thing. I'll write one soon.

[edrumwri]

Thanks guys- you rock!

[soonho-tri]

@edrumwri, https://github.com/RobotLocomotion/drake/pull/3878 is merged. Please rebase your commits and compile again.

[edrumwri]

Pull request #4006 shows that this has been addressed. Thanks!