Damp initial iterations of Broyden's method used in ID parameter control

[iago-mendes]

Proposed changes

This PR modifies the initial data parameter control so that it doesn't go off-bounds in the q3 run of arXiv:1506.01689. We tried different modifications to Broyden's method, summarized in the plot below. The one that seems to work best for now is to damp the initial iterations of the Jacobian updates. We plan to make this damping more robust in the future.

Upgrade instructions

Code review checklist

• The code is documented and the documentation renders correctly. Run
  make doc to generate the documentation locally into BUILD_DIR/docs/html.
  Then open index.html.
• The code follows the stylistic and code quality guidelines listed in the
  code review guide.
• The PR lists upgrade instructions and is labeled bugfix or
  new feature if appropriate.

Further comments

[nilsdeppe]

Sorry, it's not clear which of the ones you mean. There are 4 curves that damp using the Jacobian. It also seems like delaying the P & C updates 2 iterations is even better, so I'm surprised by your recommendation. Naively I'd expect combining P&C delay, and then after 3 iterations damping with 1-e^{-k} to prove the best combination of the ideas presented.

[iago-mendes]

Apologies for the non-clarity. Out of the 4 curves that damp the Jacobian, this PR uses the one with -k/2 in the exponent (pink curve).

Indeed, the initial delay of Padm and CoM seems to perform better in this case. My understanding from conversations I've had with @nilsvu is that it would be more straightforward to modify the damping approach in the future to make the control scheme more robust, instead of changing it to work for each different configuration.

I am curious about how @nilsdeppe's suggestion performs. I'll experiment with it and share updates here as soon as possible. I'm also happy to explore other variations.

[teukolsky]

@iago-mendes I've had good success in other problems using Broyden's method by initializing the Jacobian using numerical finite differences instead of using the identity matrix. The idea is that using the identity matrix, the first N iterations are building up an approximation to the Jacobian (N = # of parameters to be found), so you need at least N iterations before you start converging. Instead, using finite differencing also takes N function evaluations to get the Jacobian, and it should be more accurate that the first N Broyden steps. I'm not sure this is worth exploring in your case, since the number of iterations is reasonably small, but I wanted to mention it for completeness.

[iago-mendes]

Thanks for the suggestion, @teukolsky! Here, we've carefully chosen our free data so that the Jacobian is approximately an identity matrix. In fact, in the Newtonian approximation, it is exactly an identity matrix. If interested, I explain how it works on this report (page 7). Hence, I believe it's not a problem to initialize the Jacobian in this way.

I've experimented with @nilsdeppe's suggestion, but 1-exp(-k) approaches 1 soon after the first couple of iterations, so the damping isn't doing much. While doing these tests again, I realized that what I had previously called "delay first 2 iterations" was actually iterations < 2, so it should be called "delay first iteration".

After talking to @nilsvu about the different options, we agree that none of them are as robust as we eventually want them to be, so it makes sense to go with the one that converges the fastest for now. Hence, I've updated this PR to the delay approach.

Our current focus is on controlling Eadm and Jadm, which will be relevant for hyperbolic encounters. Once that's done, we want to do a more complete study in parameter space. At that time, we'll want to find a more robust approach for this.

[teukolsky]

Ah! The fact that your initial Jacobian is a good approximation explains why you're in the convergent regime immediately. I didn't understand why there was good convergence before N iterations.

[nilsdeppe]

Sounds great! Thanks for the very thorough update and looking into that! I think the plan sounds good!

@nilsvu since this is python code and in things you wrote, would you mind reviewing the code itself?