Check the Jacobian of a discipline

In this example, the Jacobian of an MDODiscipline is checked by derivative approximation.

from __future__ import annotations

from typing import TYPE_CHECKING

from numpy import array
from numpy import exp

from gemseo import configure_logger
from gemseo.core.discipline import MDODiscipline

if TYPE_CHECKING:
    from collections.abc import Iterable

configure_logger()

<RootLogger root (INFO)>

First, we create a discipline computing \(f(x,y)=e^{-(1-x)^2-(1-y)^2}\) and \(g(x,y)=x^2+y^2-1\) and introduce an error in the implementation of \(\frac{\partial f(x,y)}{\partial x}\).

class BuggedDiscipline(MDODiscipline):
    def __init__(self) -> None:
        super().__init__()
        self.input_grammar.update_from_names(["x", "y"])
        self.output_grammar.update_from_names(["f", "g"])
        self.default_inputs = {"x": array([0.0]), "y": array([0.0])}

    def _run(self) -> None:
        x, y = self.get_inputs_by_name(["x", "y"])
        self.local_data["f"] = exp(-((1 - x) ** 2) - (1 - y) ** 2)
        self.local_data["g"] = x**2 + y**2 - 1

    def _compute_jacobian(
        self,
        inputs: Iterable[str] | None = None,
        outputs: Iterable[str] | None = None,
    ) -> None:
        x, y = self.get_inputs_by_name(["x", "y"])
        self._init_jacobian()
        g_jac = self.jac["g"]
        g_jac["x"][:] = 2 * x
        g_jac["y"][:] = 2 * y
        f_jac = self.jac["f"]
        aux = 2 * exp(-((1 - x) ** 2) - (1 - y) ** 2)
        f_jac["x"][:] = aux  # this is wrong.
        f_jac["y"][:] = aux * (1 - y)

We want to check if the implemented Jacobian is correct. For practical applications where Jacobians are needed, this is not a simple task. GEMSEO automates such tests thanks to the MDODiscipline.check_jacobian() method.

Finite differences (default)

discipline = BuggedDiscipline()
discipline.check_jacobian(
    input_data={"x": array([0.0]), "y": array([1.0])},
    show=True,
    plot_result=True,
    step=1e-1,
)

   ERROR - 09:00:21: BuggedDiscipline Jacobian: dp f/d x is wrong by 1.9226823669921054%.
    INFO - 09:00:21: Approximate jacobian =
    INFO - 09:00:21: [[0.76978625]]
    INFO - 09:00:21: Provided by linearize method =
    INFO - 09:00:21: [[0.73575888]]
    INFO - 09:00:21: Difference of jacobians =
    INFO - 09:00:21: [[0.03402737]]
   ERROR - 09:00:21: BuggedDiscipline Jacobian: dp f/d y is wrong by 3.5312032605514854%.
    INFO - 09:00:21: Approximate jacobian =
    INFO - 09:00:21: [[-0.03660462]]
    INFO - 09:00:21: Provided by linearize method =
    INFO - 09:00:21: [[0.]]
    INFO - 09:00:21: Difference of jacobians =
    INFO - 09:00:21: [[-0.03660462]]
   ERROR - 09:00:21: BuggedDiscipline Jacobian: dp g/d x is wrong by 9.090909090909099%.
    INFO - 09:00:21: Approximate jacobian =
    INFO - 09:00:21: [[0.1]]
    INFO - 09:00:21: Provided by linearize method =
    INFO - 09:00:21: [[0.]]
    INFO - 09:00:21: Difference of jacobians =
    INFO - 09:00:21: [[0.1]]
   ERROR - 09:00:21: BuggedDiscipline Jacobian: dp g/d y is wrong by 3.2258064516129616%.
    INFO - 09:00:21: Approximate jacobian =
    INFO - 09:00:21: [[2.1]]
    INFO - 09:00:21: Provided by linearize method =
    INFO - 09:00:21: [[2.]]
    INFO - 09:00:21: Difference of jacobians =
    INFO - 09:00:21: [[0.1]]
    INFO - 09:00:21: Linearization of MDODiscipline: BuggedDiscipline is wrong.
/home/docs/checkouts/readthedocs.org/user_builds/gemseo/envs/develop/lib/python3.9/site-packages/gemseo/utils/derivatives/derivatives_approx.py:659: UserWarning: set_ticklabels() should only be used with a fixed number of ticks, i.e. after set_ticks() or using a FixedLocator.
  axe.set_xticklabels(x_labels, fontsize=14)

False

The step here is chosen big enough to underline the truncation error. From this graph, we can see that almost all the provided components of the Jacobians (blue dots) are close but distinct from the approximated by finite differences using a step of 0.1 (red dots). This kind of graph can be used to spot implementation mistakes in fact we can already spot a large mistake in the wrong components.

# The ``derr_approx`` argument can be either ``finite_differences``, ``centered_differences`` or
# ``complex_step``.

# Centered differences
# --------------------
discipline.check_jacobian(
    input_data={"x": array([0.0]), "y": array([1.0])},
    derr_approx=discipline.ApproximationMode.CENTERED_DIFFERENCES,
    show=True,
    plot_result=True,
    step=1e-1,
)

# With the same step the truncation error is in this case much smaller.


# Complex step
# ------------
discipline.check_jacobian(
    input_data={"x": array([0.0]), "y": array([1.0])},
    derr_approx=discipline.ApproximationMode.COMPLEX_STEP,
    show=True,
    plot_result=True,
    step=1e-1,
)

# With the same step the truncation error is also smaller than finite differences.
# This confirms again that an implementation mistake was done.

# Advantages and drawbacks of each method
# ---------------------------------------
# Finite differnces and complex are first-order methods, they use one
# sampling point per input and the truncation error goes down linearly with the step.
# Centered differences are second-order methods which use twice as many points as finite
# differences and complex step. Complex step derivatives are less prone to numerical
# cancellation errors so that a tiny step can be used. On the other hand complex step is
# not compatible with discipline not supporting complex inputs.

discipline.check_jacobian(
    input_data={"x": array([0.0]), "y": array([1.0])},
    derr_approx=discipline.ApproximationMode.COMPLEX_STEP,
    show=True,
    plot_result=True,
    step=1e-10,
)

   ERROR - 09:00:22: BuggedDiscipline Jacobian: dp f/d x is wrong by 0.14163403944972192%.
    INFO - 09:00:22: Approximate jacobian =
    INFO - 09:00:22: [[0.73330393]]
    INFO - 09:00:22: Provided by linearize method =
    INFO - 09:00:22: [[0.73575888]]
    INFO - 09:00:22: Difference of jacobians =
    INFO - 09:00:22: [[-0.00245495]]
    INFO - 09:00:22: Jacobian: dp f/dp y succeeded.
    INFO - 09:00:22: Jacobian: dp g/dp x succeeded.
    INFO - 09:00:22: Jacobian: dp g/dp y succeeded.
    INFO - 09:00:22: Linearization of MDODiscipline: BuggedDiscipline is wrong.
/home/docs/checkouts/readthedocs.org/user_builds/gemseo/envs/develop/lib/python3.9/site-packages/gemseo/utils/derivatives/derivatives_approx.py:659: UserWarning: set_ticklabels() should only be used with a fixed number of ticks, i.e. after set_ticks() or using a FixedLocator.
  axe.set_xticklabels(x_labels, fontsize=14)
   ERROR - 09:00:22: BuggedDiscipline Jacobian: dp f/d x is wrong by 0.1409521643331837%.
    INFO - 09:00:22: Approximate jacobian =
    INFO - 09:00:22: [[0.73820893]]
    INFO - 09:00:22: Provided by linearize method =
    INFO - 09:00:22: [[0.73575888]]
    INFO - 09:00:22: Difference of jacobians =
    INFO - 09:00:22: [[0.00245004]]
    INFO - 09:00:22: Jacobian: dp f/dp y succeeded.
    INFO - 09:00:22: Jacobian: dp g/dp x succeeded.
    INFO - 09:00:22: Jacobian: dp g/dp y succeeded.
    INFO - 09:00:22: Linearization of MDODiscipline: BuggedDiscipline is wrong.
/home/docs/checkouts/readthedocs.org/user_builds/gemseo/envs/develop/lib/python3.9/site-packages/gemseo/utils/derivatives/derivatives_approx.py:659: UserWarning: set_ticklabels() should only be used with a fixed number of ticks, i.e. after set_ticks() or using a FixedLocator.
  axe.set_xticklabels(x_labels, fontsize=14)
    INFO - 09:00:23: Jacobian: dp f/dp x succeeded.
    INFO - 09:00:23: Jacobian: dp f/dp y succeeded.
    INFO - 09:00:23: Jacobian: dp g/dp x succeeded.
    INFO - 09:00:23: Jacobian: dp g/dp y succeeded.
    INFO - 09:00:23: Linearization of MDODiscipline: BuggedDiscipline is correct.
/home/docs/checkouts/readthedocs.org/user_builds/gemseo/envs/develop/lib/python3.9/site-packages/gemseo/utils/derivatives/derivatives_approx.py:659: UserWarning: set_ticklabels() should only be used with a fixed number of ticks, i.e. after set_ticks() or using a FixedLocator.
  axe.set_xticklabels(x_labels, fontsize=14)

True

Automatic time step

Finite differences and centered differences steps need to be chosen as a trade between truncation and numerical errors. For this reason, the auto_set_step option can be used to automatically compute the step where the total error is minimized.

discipline.check_jacobian(
    input_data={"x": array([0.0]), "y": array([1.0])},
    derr_approx=discipline.ApproximationMode.CENTERED_DIFFERENCES,
    show=True,
    plot_result=True,
    auto_set_step=True,
)

    INFO - 09:00:23: Set optimal step for finite differences. Estimated approximation errors =
    INFO - 09:00:23: [1.20985073e-07 1.20985073e-07]
    INFO - 09:00:23: Jacobian: dp f/dp x succeeded.
    INFO - 09:00:23: Jacobian: dp f/dp y succeeded.
    INFO - 09:00:23: Jacobian: dp g/dp x succeeded.
    INFO - 09:00:23: Jacobian: dp g/dp y succeeded.
    INFO - 09:00:23: Linearization of MDODiscipline: BuggedDiscipline is correct.
/home/docs/checkouts/readthedocs.org/user_builds/gemseo/envs/develop/lib/python3.9/site-packages/gemseo/utils/derivatives/derivatives_approx.py:659: UserWarning: set_ticklabels() should only be used with a fixed number of ticks, i.e. after set_ticks() or using a FixedLocator.
  axe.set_xticklabels(x_labels, fontsize=14)

True