Approximating Partial Derivatives

Contents

Approximating Partial Derivatives#

OpenMDAO allows you to specify analytic derivatives for your models, but it is not a requirement. If certain partial derivatives are not available, you can ask the framework to approximate the derivatives by using the declare_partials method inside setup_partials, and give it a method that is either ‘fd’ for finite difference or ‘cs’ for complex step.

Component.declare_partials(of, wrt, dependent=True, rows=None, cols=None, val=None, method='exact', step=None, form=None, step_calc=None, minimum_step=None)[source]

Declare information about this component’s subjacobians.

Parameters:
ofstr or iter of str

The name of the residual(s) that derivatives are being computed for. May also contain a glob pattern.

wrtstr or iter of str

The name of the variables that derivatives are taken with respect to. This can contain the name of any input or output variable. May also contain a glob pattern.

dependentbool(True)

If False, specifies no dependence between the output(s) and the input(s). This is only necessary in the case of a sparse global jacobian, because if ‘dependent=False’ is not specified and declare_partials is not called for a given pair, then a dense matrix of zeros will be allocated in the sparse global jacobian for that pair. In the case of a dense global jacobian it doesn’t matter because the space for a dense subjac will always be allocated for every pair.

rowsndarray of int or None

Row indices for each nonzero entry. For sparse subjacobians only.

colsndarray of int or None

Column indices for each nonzero entry. For sparse subjacobians only.

valfloat or ndarray of float or scipy.sparse

Value of subjacobian. If rows and cols are not None, this will contain the values found at each (row, col) location in the subjac.

methodstr

The type of approximation that should be used. Valid options include: ‘fd’: Finite Difference, ‘cs’: Complex Step, ‘exact’: use the component defined analytic derivatives. Default is ‘exact’.

stepfloat

Step size for approximation. Defaults to None, in which case the approximation method provides its default value.

formstr

Form for finite difference, can be ‘forward’, ‘backward’, or ‘central’. Defaults to None, in which case the approximation method provides its default value.

step_calcstr

Step type for computing the size of the finite difference step. It can be ‘abs’ for absolute, ‘rel_avg’ for a size relative to the absolute value of the vector input, or ‘rel_element’ for a size relative to each value in the vector input. In addition, it can be ‘rel_legacy’ for a size relative to the norm of the vector. For backwards compatibilty, it can be ‘rel’, which is now equivalent to ‘rel_avg’. Defaults to None, in which case the approximation method provides its default value.

minimum_stepfloat

Minimum step size allowed when using one of the relative step_calc options.

Returns:
dict

Metadata dict for the specified partial(s).

Usage#

You may use glob patterns as arguments to to and wrt.

import numpy as np

import openmdao.api as om

class FDPartialComp(om.ExplicitComponent):
    def setup(self):
        self.add_input('x', shape=(4,))
        self.add_input('y', shape=(2,))
        self.add_input('y2', shape=(2,))
        self.add_output('f', shape=(2,))

    def setup_partials(self):
        self.declare_partials('f', 'y*', method='fd')
        self.declare_partials('f', 'x', method='fd')

    def compute(self, inputs, outputs):
        f = outputs['f']

        x = inputs['x']
        y = inputs['y']

        f[0] = x[0] + y[0]
        f[1] = np.dot([0, 2, 3, 4], x) + y[1]

model = om.Group()
model.add_subsystem('example', FDPartialComp())

problem = om.Problem(model=model)
problem.setup()
problem.run_model()
totals = problem.compute_totals(['example.f'], ['example.x', 'example.y'])
print(totals['example.f', 'example.x'])
[[ 1. -0. -0. -0.]
 [-0.  2.  3.  4.]]
print(totals['example.f', 'example.y'])
[[ 1. -0.]
 [-0.  1.]]

For finite difference approximations (method=’fd’), we have three (optional) parameters: the form, step size, and the step_calc. The form should be one of the following:

  • form=’forward’ (default): Approximates the derivative as \(\frac{\partial f}{\partial x} \approx \frac{f(x+\delta, y) - f(x,y)}{||\delta||}\). Error scales like \(||\delta||\).

  • form=’backward’: Approximates the derivative as \(\frac{\partial f}{\partial x} \approx \frac{f(x,y) - f(x-\delta, y) }{||\delta||}\). Error scales like \(||\delta||\).

  • form=’central’: Approximates the derivative as \(\frac{\partial f}{\partial x} \approx \frac{f(x+\delta, y) - f(x-\delta,y)}{2||\delta||}\). Error scales like \(||\delta||^2\), but requires an extra function evaluation.

The step size can be any nonzero number, but should be positive (one can change the form to perform backwards finite difference formulas), small enough to reduce truncation error, but large enough to avoid round-off error. Choosing a step size can be highly problem dependent, but for double precision floating point numbers and reasonably bounded derivatives, 10−6 can be a good place to start.

The step_calc can be either ‘abs’ for absolute or one of (‘rel’, ‘rel_avg’, ‘rel_legacy’, ‘rel_element’) for different forms of relative stepping. This determines whether the stepsize ie absolute or a percentage of the input value. The following table details how the relative step is calculated:

step_calc

Step size is scaled by

“rel_avg”

Average absolute value of the vector.

“rel_element”

Absolute value of each vector element.

“rel_legacy”

Norm of the vector.

“rel”

Same as “rel_avg”.

import numpy as np

import openmdao.api as om

class FDPartialComp(om.ExplicitComponent):

    def setup(self):
        self.add_input('x', shape=(4,))
        self.add_input('y', shape=(2,))
        self.add_input('y2', shape=(2,))
        self.add_output('f', shape=(2,))

    def setup_partials(self):
        self.declare_partials('f', 'y*', method='fd', form='backward', step=1e-6)
        self.declare_partials('f', 'x', method='fd', form='central', step=1e-4)

    def compute(self, inputs, outputs):
        f = outputs['f']

        x = inputs['x']
        y = inputs['y']

        f[0] = x[0] + y[0]
        f[1] = np.dot([0, 2, 3, 4], x) + y[1]

model = om.Group()
model.add_subsystem('example', FDPartialComp())

problem = om.Problem(model=model)
problem.setup()
problem.run_model()
totals = problem.compute_totals(['example.f'], ['example.x', 'example.y'])
print(totals['example.f', 'example.x'])
[[ 1. -0. -0. -0.]
 [-0.  2.  3.  4.]]

Complex Step#

If you have a pure python component (or an external code that can support complex inputs and outputs), then you can also choose to use complex step to calculate the Jacobian of that component. This will give more accurate derivatives that are less sensitive to the step size. Like finite difference, complex step runs your component using the apply_nonlinear or solve_nonlinear functions, but it applies a step in the complex direction. You can activate it using the declare_partials method inside setup_partials and giving it a method of ‘cs’. In many cases, this will require no other changes to your code, as long as all of the calculation in your solve_nonlinear and apply_nonlinear support complex numbers. During a complex step, the incoming inputs vector will return a complex number when a variable is being stepped. Likewise, the outputs and residuals vectors will accept complex values. If you are allocating temporary numpy arrays, remember to conditionally set their dtype based on the dtype in the outputs vector.

Here is how to turn on complex step for all input/output pairs in the Sellar problem:

class SellarDis1CS(om.ExplicitComponent):
    """
    Component containing Discipline 1 -- no derivatives version.
    Uses Complex Step
    """

    def setup(self):

        # Global Design Variable
        self.add_input('z', val=np.zeros(2))

        # Local Design Variable
        self.add_input('x', val=0.)

        # Coupling parameter
        self.add_input('y2', val=1.0)

        # Coupling output
        self.add_output('y1', val=1.0)

    def setup_partials(self):
        # Finite difference all partials.
        self.declare_partials('*', '*', method='cs')

    def compute(self, inputs, outputs):
        """
        Evaluates the equation
        y1 = z1**2 + z2 + x1 - 0.2*y2
        """
        z1 = inputs['z'][0]
        z2 = inputs['z'][1]
        x1 = inputs['x']
        y2 = inputs['y2']

        outputs['y1'] = z1**2 + z2 + x1 - 0.2*y2
class SellarDis2CS(om.ExplicitComponent):
    """
    Component containing Discipline 2 -- no derivatives version.
    Uses Complex Step
    """

    def setup(self):
        # Global Design Variable
        self.add_input('z', val=np.zeros(2))

        # Coupling parameter
        self.add_input('y1', val=1.0)

        # Coupling output
        self.add_output('y2', val=1.0)

    def setup_partials(self):
        # Finite difference all partials.
        self.declare_partials('*', '*', method='cs')

    def compute(self, inputs, outputs):
        """
        Evaluates the equation
        y2 = y1**(.5) + z1 + z2
        """

        z1 = inputs['z'][0]
        z2 = inputs['z'][1]
        y1 = inputs['y1']

        # Note: this may cause some issues. However, y1 is constrained to be
        # above 3.16, so lets just let it converge, and the optimizer will
        # throw it out
        if y1.real < 0.0:
            y1 *= -1

        outputs['y2'] = y1**.5 + z1 + z2

Sometimes you need to know when you are under a complex step so that your component can correctly handle complex inputs (e.g, in case you need to allocate intermediate arrays as complex.) All Components and Groups provide the attribute under_complex_step that you can use to tell if you are under a complex step. In the following example, we print out the incoming complex value when the “compute” method is called while computing this component’s derivatives under complex step.

import openmdao.api as om

class SimpleComp(om.ExplicitComponent):

    def setup(self):
        self.add_input('x', val=1.0)
        self.add_output('y', val=1.0)

        self.declare_partials(of='y', wrt='x', method='cs')

    def compute(self, inputs, outputs):
        outputs['y'] = 3.0*inputs['x']

        if self.under_complex_step:
            print("Under complex step")
            print("x", inputs['x'])
            print("y", outputs['y'])

prob = om.Problem()
prob.model.add_subsystem('comp', SimpleComp())

prob.model.add_design_var('comp.x', lower=-100, upper=100)
prob.model.add_objective('comp.y')

prob.setup(force_alloc_complex=True)

prob.run_model()

prob.compute_totals(of=['comp.y'], wrt=['comp.x'])
Under complex step
x [1.+1.e-40j]
y [3.+3.e-40j]
{('comp.y', 'comp.x'): array([[3.]])}