Compute the Jacobian of a discipline with finite differences

In this example, we will compute the Jacobians of some outputs of an Discipline with respect to some inputs using the finite difference method.

from __future__ import annotations

from numpy import array

from gemseo.disciplines.auto_py import AutoPyDiscipline

First, we create a discipline, e.g. an AutoPyDiscipline:

def f(a=0.0, b=0.0):
    y = a**2 + b
    z = a**3 + b**2
    return y, z


discipline = AutoPyDiscipline(f)

We can execute it with its default input values:

discipline.execute()
discipline.local_data

{'a': array([0.]), 'b': array([0.]), 'y': array([0.]), 'z': array([0.])}

or with custom ones:

discipline.execute({"a": array([1.0])})
discipline.local_data

{'a': array([1.]), 'b': array([0.]), 'y': array([1.]), 'z': array([1.])}

Then, we use the method Discipline.linearize() to compute the derivatives:

jacobian_data = discipline.linearize()
jacobian_data

{}

There is no Jacobian data because we need to set the input variables with respect to which to compute the Jacobian of the output ones. For that, we use the method add_differentiated_inputs(). We also need to set these output variables: with the method add_differentiated_outputs(). For instance, we may want to compute the derivative of "z" with respect to "a" only:

discipline.add_differentiated_inputs(["a"])
discipline.add_differentiated_outputs(["z"])
jacobian_data = discipline.linearize()
jacobian_data

{'z': {'a': array([[1.e-14]])}}

We can have a quick look at the values of these derivatives and verify that they are equal to the analytical results, up to the numerical precision.

By default, GEMSEO uses Discipline.default_input_data as input data for which to compute the Jacobian one. We can change them with input_data:

jacobian_data = discipline.linearize({"a": array([1.0])})
jacobian_data

{'z': {'a': array([[3.0000003]])}}

We can also force the discipline to compute the derivatives of all the outputs with respect to all the inputs:

jacobian_data = discipline.linearize(compute_all_jacobians=True)
jacobian_data

{'y': {'a': array([[1.e-07]]), 'b': array([[1.]])}, 'z': {'a': array([[1.e-14]]), 'b': array([[1.e-07]])}}

we can change the approximation type to complex step and compare the results:

discipline.set_jacobian_approximation(
    jac_approx_type=discipline.ApproximationMode.COMPLEX_STEP
)
jacobian_data_complex_step = discipline.linearize(compute_all_jacobians=True)
jacobian_data_complex_step

{'y': {'a': array([[0.]]), 'b': array([[0.]])}, 'z': {'a': array([[0.]]), 'b': array([[0.]])}}

Lastly, We can change the approximation type to centered differences and compare the results:

discipline.set_jacobian_approximation(
    jac_approx_type=discipline.ApproximationMode.CENTERED_DIFFERENCES
)
jacobian_data_centered_differences = discipline.linearize(compute_all_jacobians=True)
jacobian_data_centered_differences

{'y': {'a': array([[0.]]), 'b': array([[1.]])}, 'z': {'a': array([[1.e-14]]), 'b': array([[0.]])}}