Thanks for the feedback, Thomas! Many great points, and I will definitely look into reaching out to the GPyTorch community.
Regarding integrating the various bits of the optimizer into optimizer.step, we have chosen not to do this in order to ensure flexibility in the sample selection when defining the updates and gradient differences in L-BFGS. This (I think) is the easiest way to make it compatible with both multi-batch and full-overlap L-BFGS. We also suspect that as stochastic quasi-Newton methods improve and other potential mechanisms are introduced, one may need to make use of the two-loop recursion outside of simply defining the search direction, which is why we had chosen to make this separate for now. (One example of this is in what our group calls “progressive batching” or “adaptive sampling” tests where adaptive tests are used for automatically increasing the batch size in the algorithm.)