Repository: astropy/astropy

Issue:
Implement deriv and fit_deriv for compound models
Compound models currently do not implement analytic derivatives (with respect to their parameters) through a `fit_deriv` method.  This hasn't been a showstopper, since all of the optimizers that use the derivative can use an estimate if an analytic expression is not provided.

However, this is much less straightforward than mere evaluation of compound models.  It will have compute partial derivatives with respect to each parameter in the compound model, and will have to take into account which operators each model is combined with, and implement the appropriate product rule, quotient rule, chain rule, etc.  There's nothing especially complicated about the math here, but getting the implementation right may be a bit tedious.  The most annoying part about this is that we'll have to implement derivatives with respect to the inputs for each model (`fit_deriv` implements partial derivatives w.r.t. the fittable parameters).  Perhaps there should be a separate issue for this.

To provide some hints for anyone who might like to tackle this, model _evaluation_ is currently implemented by passing a dictionary mapping operators (given as strings) to functions that implement that operator between two models:

https://github.com/astropy/astropy/blob/cee42f3923a4fa9ea4d20c206baa9a9789e9defd/astropy/modeling/core.py#L1379

For example, the `'+'` operator is simply implemented as "A(x) + B(x)" where A and B are the `evaluate` methods of two models (there is a little more complexity involved to support things like multiple inputs/outputs.  You can see how this is implemented in `_make_arithmetic_operator` above.  

Similarly the `'|'` operator is implemented to evaluate "B(A(x))" (again modulo some slight complications for generalization, and parameter handling--let me know if these need better explanation).

So for evaluating `fit_deriv`, if you have a model like "A + B", the `'+'` operation would have to be implemented like `A.fit_deriv(x) + B.fit_deriv(x)`.  That's easy enough.  But for implementing `'|'` we'll need a derivative method.  So "A | B" could be implemented like "B.deriv(A(x)) \* A.fit_deriv(x)".  Though things get more complicated for models of multiple arguments, etc.  Not also that in this example, if "B" is itself a compound model it will have to use a `B.deriv` for the entire compound model; i.e. this would have to work recursively.


Discussion:
github-actions[bot]: Thank you for your contribution to Astropy! 🌌 This checklist is meant to remind the package maintainers who will review this pull request of some common things to look for.

  - [ ] Do the proposed changes actually accomplish desired goals?
  - [ ] Do the proposed changes follow the [Astropy coding guidelines](https://docs.astropy.org/en/latest/development/codeguide.html)?
  - [ ] Are tests added/updated as required? If so, do they follow the [Astropy testing guidelines](https://docs.astropy.org/en/latest/development/testguide.html)?
  - [ ] Are docs added/updated as required? If so, do they follow the [Astropy documentation guidelines](https://docs.astropy.org/en/latest/development/docguide.html)?
  - [ ] Is rebase and/or squash necessary? If so, please provide the author with appropriate instructions. Also see instructions for [rebase](https://docs.astropy.org/en/latest/development/development_details.html#rebase-if-necessary) and [squash](https://docs.astropy.org/en/latest/development/development_details.html#squash-if-necessary).
  - [ ] Did the CI pass? If no, are the failures related? If you need to run daily and weekly cron jobs as part of the PR, please apply the "Extra CI" label. Codestyle issues can be fixed by the [bot](https://docs.astropy.org/en/latest/development/development_details.html#pre-commit).
  - [ ] Is a change log needed? If yes, did the change log check pass? If no, add the "no-changelog-entry-needed" label. If this is a manual backport, use the "skip-changelog-checks" label unless special changelog handling is necessary.
  - [ ] Is this a big PR that makes a "What's new?" entry worthwhile and if so, is (1) a "what's new" entry included in this PR and (2) the "whatsnew-needed" label applied?
  - [ ] At the time of adding the milestone, if the milestone set requires a backport to release branch(es), apply the appropriate "backport-X.Y.x" label(s) *before* merge.
astrofrog: The Polynomial2D thing was a red herring, fitters will just call ``fit_deriv`` with 1D arrays even when there are multiple dimensions, so this won't matter. We are down to three test failures: one because there is a difference of NaN vs Inf for the numerical vs analytical derivatives, and the other one is a slight numerical difference in the final model, which might just be a numerical issue, so might be safe to increase the relative tolerance?

Patch:
diff --git a/astropy/modeling/core.py b/astropy/modeling/core.py
index 653f043da8..8f9fd76e62 100644
--- a/astropy/modeling/core.py
+++ b/astropy/modeling/core.py
@@ -3163,8 +3163,6 @@ class CompoundModel(Model):
             raise ModelDefinitionError("Illegal operator: ", self.op)
         self.name = name
         self._fittable = None
-        self.fit_deriv = None
-        self.col_fit_deriv = None
         if op in ("|", "+", "-"):
             self.linear = left.linear and right.linear
         else:
@@ -3294,6 +3292,86 @@ class CompoundModel(Model):
 
         return self._apply_operators_to_value_lists(leftval, rightval, **kw)
 
+    @property
+    def fit_deriv(self):
+        # If either side of the model is missing analytical derivative then we can't compute one
+        if self.left.fit_deriv is None or self.right.fit_deriv is None:
+            return None
+
+        # Only the following operators are supported
+        op = self.op
+        if op not in ["-", "+", "*", "/"]:
+            return None
+
+        def _calc_compound_deriv(*args, **kwargs):
+            args, kw = self._get_kwarg_model_parameters_as_positional(args, kwargs)
+            left_inputs = self._get_left_inputs_from_args(args)
+            left_params = self._get_left_params_from_args(args)
+
+            right_inputs = self._get_right_inputs_from_args(args)
+            right_params = self._get_right_params_from_args(args)
+
+            left_deriv = self.left.fit_deriv(*left_inputs, *left_params)
+            if not self.left.col_fit_deriv:
+                left_deriv = (np.asanyarray(left_deriv).T).tolist()
+            right_deriv = self.right.fit_deriv(*right_inputs, *right_params)
+            if not self.right.col_fit_deriv:
+                right_deriv = (np.asanyarray(right_deriv).T).tolist()
+
+            # We now have to use various differentiation rules to apply the
+            # arithmetic operators to the derivatives.
+            # If we consider an example of a compound model
+            # h(x, a, b, c) made up of two models g(x, a)
+            # and h(x, b, c), one with one parameter and
+            # the other with two parameters, the derivatives
+            # are evaluated as follows:
+
+            # Addition
+            # h(x, a, b, c) = f(x, a) + g(x, b, c)
+            # fit_deriv = [df/da, dg/db, dg/dc]
+
+            # Subtraction
+            # h(x, a, b, c) = f(x, a) - g(x, b, c)
+            # fit_deriv = [df/da, -dg/db, -dg/dc]
+
+            # Multiplication
+            # h(x, a, b, c) = f(x, a) * g(x, b, c)
+            # fit_deriv = [g(x, b, c) * df/da,
+            #              f(x, a) * dg/db,
+            #              f(x, a) * dg/dc]
+
+            # Division - Quotient rule
+            # h(x, a, b, c) = f(x, a) / g(x, b, c)
+            # fit_deriv = [df/da / g(x, b, c),
+            #              -f(x, a) * dg/db / g(x, b, c)**2,
+            #              -f(x, a) * dg/dc / g(x, b, c)**2]
+
+            if op in ["+", "-"]:
+                if op == "-":
+                    right_deriv = [-x for x in right_deriv]
+
+                return left_deriv + right_deriv
+
+            leftval = self.left.evaluate(*left_inputs, *left_params)
+            rightval = self.right.evaluate(*right_inputs, *right_params)
+
+            if op == "*":
+                return (
+                    [rightval * dparam for dparam in left_deriv] +
+                    [leftval * dparam for dparam in right_deriv]
+                )  # fmt: skip
+            if op == "/":
+                return (
+                    [dparam / rightval for dparam in left_deriv] +
+                    [-leftval * (dparam / rightval**2) for dparam in right_deriv]
+                )  # fmt: skip
+
+        return _calc_compound_deriv
+
+    @property
+    def col_fit_deriv(self):
+        return True
+
     @property
     def n_submodels(self):
         if self._leaflist is None:
diff --git a/docs/changes/modeling/17034.perf.rst b/docs/changes/modeling/17034.perf.rst
new file mode 100644
index 0000000000..090e3b633c
--- /dev/null
+++ b/docs/changes/modeling/17034.perf.rst
@@ -0,0 +1,1 @@
+``CompoundModel`` now implements numerical derivatives of parameters when using the +, -, * or / operators. This improves the speed of fitting these models because numerical derivatives of the parameters are not calculated.


Test patch:
diff --git a/astropy/modeling/tests/test_compound.py b/astropy/modeling/tests/test_compound.py
index 8b33ddc407..5f657d60b8 100644
--- a/astropy/modeling/tests/test_compound.py
+++ b/astropy/modeling/tests/test_compound.py
@@ -1216,3 +1216,61 @@ def test_fit_mixed_recursive_compound_model_with_mixed_units():
     assert_quantity_allclose(fit.intercept_2, 10 * u.m)
     assert_quantity_allclose(fit.slope_3, 0 * u.kg / u.s)
     assert_quantity_allclose(fit.intercept_3, 5 * u.kg)
+
+
+def numerical_partial_deriv(model, *inputs, param_idx, delta=1e-5):
+    """
+    Evaluate the central difference approximation of the derivative for param_idx.
+
+    Parameters
+    ----------
+    model
+        The model to evaluate
+    inputs
+        The inputs to the model
+    param_idx
+        The index of the parameter to compute the partial derivative for.
+    delta
+        The step size with which to compute the central difference.
+    """
+    param = model.parameters
+
+    param_down = param.copy()
+    param_down[param_idx] = param[param_idx] - delta
+    param_up = param.copy()
+    param_up[param_idx] = param[param_idx] + delta
+
+    up = model.evaluate(*inputs, *param_up)
+    down = model.evaluate(*inputs, *param_down)
+
+    return (up - down) / (2 * delta)
+
+
+@pytest.mark.parametrize(
+    "model",
+    [
+        pytest.param(
+            m, id=m._format_expression(format_leaf=lambda i, l: type(l).__name__)
+        )
+        for m in [
+            Gaussian1D(5, 2, 3) + Linear1D(2, 3),
+            Gaussian1D(5, 2, 3) - Linear1D(2, 3),
+            Polynomial1D(2) * Gaussian1D(),
+            Polynomial1D(2) / Gaussian1D(),
+            Polynomial1D(2) + Gaussian1D(),
+        ]
+    ],
+)
+def test_compound_fit_deriv(model):
+    """
+    Given some compound models compare the numerical derivatives to analytical ones.
+    """
+    x = np.linspace(1, 5, num=10)
+    numerical = [
+        numerical_partial_deriv(model, x, param_idx=i)
+        for i in range(len(model.parameters))
+    ]
+    analytical = model.fit_deriv(x, *model.parameters)
+
+    for n, a in zip(numerical, analytical):
+        assert np.allclose(n, a)
diff --git a/astropy/modeling/tests/test_fitting_parallel.py b/astropy/modeling/tests/test_fitting_parallel.py
index 92c6e7aec2..863518467c 100644
--- a/astropy/modeling/tests/test_fitting_parallel.py
+++ b/astropy/modeling/tests/test_fitting_parallel.py
@@ -637,10 +637,10 @@ def test_compound_model():
         scheduler="synchronous",
     )
 
-    assert_allclose(model_fit.amplitude_0.value, [2, 1.63349282])
-    assert_allclose(model_fit.mean_0.value, [1.1, 1.145231])
-    assert_allclose(model_fit.stddev_0.value, [0.1, 0.73632987])
-    assert_allclose(model_fit.amplitude_1.value, [2, 2])
+    assert_allclose(model_fit.amplitude_0.value, [2, 1.633], atol=0.001)
+    assert_allclose(model_fit.mean_0.value, [1.1, 1.145], atol=0.001)
+    assert_allclose(model_fit.stddev_0.value, [0.1, 0.736], atol=0.001)
+    assert_allclose(model_fit.amplitude_1.value, [2, 2], atol=0.001)
 
 
 def test_model_dimension_mismatch():