The numerical derivatives computed by the optimizer currently don't distinguish between the linear parameters (for which derivatives are trivial) and nonlinear parameters (for which they're hard), because we don't pass the information that distinguishes them to the object that computes the derivatives. If we move the computation of derivatives from the Optimizer class to the Objective class, we should be able to compute the derivatives much more efficiently.

While this doesn't matter much when fitting single component galaxy models (because there's only one linear parameter in that case), it should matter quite a bit when fitting high-order shapelets to PSF models.