[C++] [Python] change constraint satisfaction behavior for nonlinear optimizers

[suntzu86]

Note: the current satisfaction behavior, while hacky, does work and will converge on the correct solution.

Currently when an update would move us out of the domain, I do one (or more) of:

1. halve the step size till we're inside the domain again
2. project the point back to the nearest point on the domain
3. find the intersect of the update vector with the plane; reduce the update so that it remains inside

After speaking w/Prof. Frazier, I learned #2 is the only legitimate solution for domains with planar boundaries and problems that are locally convex. It's provably convergent and is equivalent to optimizing for having 0 directional derivative in all non-normal directions; i.e., it'll keep going until you're on the constraint plane and the only derivative component is in the direction of the normal (i.e., the direction that you cannot travel in). (This is definitely true for linear constraints).

That analogy doesn't hold up at corners/edges (i.e., violating multiple constraints at once) but projecting back to the nearest point on the constraint space should still be correct (again at least for linear constraints).

So let's standardize the behavior of my code. In particular, since I do the right thing and some other things, there are probably bits of code that are never exercised in practice.

For problems with nonlinear boundaries or 'very' nonlinear behaviors, other heuristics/techniques may be needed.

[suntzu86]

It occurs to me that this also implies a much better testing technique. Currently I do hacky things to make it so that the true optima lies inside the domain in most cases. Instead, I don't have to be as worried about it: I can instead test that in the event that we hit constraints, the direction of the gradient is parallel to the normal vector of the active constraint. (Again, we might need special handling for corners and edges.)

I think it'll go something like this:

Optimize.
Compute gradient at optima, call it GRAD.
IF GRAD is nonzero: (things did not converge; let's investigate the possibility that we've hit a constraint)
  Get list of bounding planes that we are "close to" (need to define close)
  Compute dot product of GRAD with outward normals. Keep only bounding planes where this dot prod is positive
  IF only one bounding plane remains: (only 1 constraint, verify we've gone as far as possible in this direction)
    Assert dot product of outward normal with GRAD is "nearly" 1.0
  ELSE: (multiple bounding constraints, try edge/corner logic)
    Verify that we are in a corner and that GRAD points "out" of the corner (not 100% sure how to do this, ugh, geometry.)
ELSE:
  everything converged, yay

Edges will be like corners except 1 or more directions of travel will have to have 0 gradient component. I'm not sure how detect this robustly (maybe it won't be worth the trouble).