Gradients - GTSAM Docs

Gradients helps differentiate through GTSAM, for example to learn factor parameters from a loss evaluated on optimized estimates. The application provides loss derivatives on the estimates, and parameterized factors provide derivatives of their raw factor residuals. Gradients combines these to return a loss gradient for a requested parameter vector, including the noise weighting.

Two noisy observations provide a worked example of learning an observation parameter. The Python tutorial extends the workflow to odometry calibration and a training update. The mathematical derivation is in the final advanced section.

The PyTorch tutorial trains a small MLP through several factor graphs using raw factor VJPs.

import numpy as np
import gtsam
from gtsam.symbol_shorthand import X

assert hasattr(gtsam, "ParameterTag"), "Use a GTSAM development build with ParameterTag."

Identifying a parameter¶

A ParameterTag identifies one application-owned parameter vector. Create the tag once and use it in both the factors and the gradient query. The tag holds no parameter values, dimensions, or factor registrations; the values can be ordinary NumPy arrays or belong to an automatic-differentiation library.

Each new tag has a distinct identity, even if its display name matches another tag’s name. Copying a tag preserves its identity and keeps it alive. Python equality and dictionary lookup compare that shared identity.

observation = gtsam.ParameterTag("observation")
same_parameter = gtsam.ParameterTag(observation)
another_parameter = gtsam.ParameterTag("observation")
assert observation == same_parameter
assert observation != another_parameter
parameters = {observation: np.array([1.0])}
np.testing.assert_array_equal(parameters[same_parameter], [1.0])

A factor with a parameter¶

Two noisy observations constrain a scalar state $x$ : the learnable observation $z_0=1$ and the fixed observation $z_1=3$ . Their standard deviations are $\sigma_0=1$ and $\sigma_1=2$ . The raw factor residual for observation $i$ is $e_i=x-z_i$ ; GTSAM applies the noise model to obtain the whitened factor residual $r_i=e_i/\sigma_i$ .

A Python CustomFactor supplies the raw factor residual and its state Jacobian in its usual error callback. Its optional fourth constructor argument supplies a vector–Jacobian product (VJP) for a requested parameter: multiply the transpose of the raw factor residual’s parameter Jacobian by the supplied vector. That supplied vector is called a factor residual cotangent; it describes how the loss changes with a raw factor residual change.

Here $\partial e_0/\partial z_0=-1$ , so the VJP returns the negative of its input. Gradients has already accounted for noise weighting and the effect of re-solving the graph. Return None for any tag the factor does not recognize. The solved state and supplied cotangent stay fixed during this calculation.

key = X(0)
observation_value = parameters[observation].copy()
sigmas = np.array([1.0, 2.0])


def observation_error(this, values, jacobians):
    if jacobians is not None:
        jacobians[0] = np.ones((1, 1))  # de/dx
    return np.array([values.atDouble(key) - observation_value[0]])


def observation_vjp(this, values, parameter, factor_residual_cotangent):
    if parameter != observation:
        return None
    return -factor_residual_cotangent  # (de/dz0).T @ cotangent


graph = gtsam.NonlinearFactorGraph()
graph.add(gtsam.CustomFactor(
    gtsam.noiseModel.Isotropic.Sigma(1, sigmas[0]), [key],
    observation_error, observation_vjp,
))
graph.resize(2)  # An empty slot is harmless; parameter queries need no indices.
graph.add(gtsam.PriorFactorDouble(
    key, 3.0, gtsam.noiseModel.Isotropic.Sigma(1, sigmas[1]),
))
initial = gtsam.Values()
initial.insert(key, 0.0)
solution = gtsam.GaussNewtonOptimizer(graph, initial).optimize()
linearized = graph.linearize(solution)
assert linearized.gradientAtZero().norm() < 1e-12
np.testing.assert_allclose(solution.atDouble(key), 1.4, atol=1e-12)
print(f"Solved x: {solution.atDouble(key):.6f}")

Solved x: 1.400000

The parameter loss gradient¶

The loss $L(x)=x$ has derivative $dL/dx=1$ at the optimized value $x^\star=1.4$ . Supply that derivative in a VectorValues, which associates each vector with its state key. These loss derivatives are called state cotangents in the API. Omitted state keys have zero loss derivative; a squared-error loss $\tfrac12\|x-\mathrm{target}\|^2$ would instead supply $x-\mathrm{target}$ .

parameterGradient asks every active factor for its contribution to the requested tag and sums the returned vectors. Ordinary factors and custom factors without a parameter callback declare no parameter dependencies. They still influence the optimized state and its derivative.

gradients = gtsam.Gradients(graph, solution)
state_cotangents = gtsam.VectorValues()
state_cotangents.insert(key, np.array([1.0]))  # dL/dx for L(x) = x
observation_gradient = gradients.parameterGradient(state_cotangents, observation)
np.testing.assert_allclose(observation_gradient, [0.8], atol=1e-12)
print("Loss gradient for the observation:", observation_gradient)

assert gradients.parameterGradient(state_cotangents, another_parameter) is None
np.testing.assert_array_equal(
    gradients.parameterGradient(gtsam.VectorValues(), observation), [0.0],
)

Loss gradient for the observation: [0.8]

Increasing the first observation increases the estimate, and hence this loss, by approximately 0.8 times the observation change. An unrecognized tag returns None, whereas a recognized parameter with zero gradient returns an explicit zero vector. All participating factors must return finite, nonempty vectors of the same dimension for a given tag.

Several factors can recognize the same tag, and one factor can recognize several tags by branching in its callback. Query each parameter separately to choose which blocks to train. Each query reuses the factorization and performs one adjoint calculation; no factor list or parameter registration is needed.

Parameter updates and automatic differentiation¶

A training step adds any direct parameter dependence of the loss to parameterGradient’s result. For a parameter vector $\theta$ and penalty weight $\alpha$ , a penalty $\tfrac12\alpha\|\theta\|^2$ contributes $\alpha\theta$ in addition to the derivative through the optimized states. The odometry tutorial computes both contributions and performs an update.

An automatic-differentiation library can implement parameterVjp by differentiating the dot product of the supplied cotangent and the raw factor residual with respect to the requested parameter. Treat the solved states and the cotangent as constants, or detach them from the library’s derivative graph. The callback returns a vector in that parameter’s coordinates. Noise models remain fixed; the callback does not differentiate or apply them.

Changing the parameters requires a new optimized solution and a new Gradients object. Reuse the parameter tags across these training steps. For repeated queries at one solution, keep the factors, their parameter values, active sets, and noise models unchanged.

Factor residual queries¶

pullback and pushforward expose lower-level operations on whitened factor residuals. They use Errors, GTSAM’s collection of factor-associated vectors, which appears as a list of NumPy arrays in Python. Each entry corresponds to an original graph slot; null or inactive slots have empty vectors. Use these operations when the application already works in factor residual coordinates.

pullback returns the loss gradient for additive changes in each whitened factor residual. The scalar example has one coordinate per observation and an empty block for the unused slot.

factor_residual_cotangents = gradients.pullback(state_cotangents)
print("First factor:", factor_residual_cotangents[0])
print("Empty factor slot:", factor_residual_cotangents[1])
print("Second factor:", factor_residual_cotangents[2])
np.testing.assert_allclose(factor_residual_cotangents[0], [-0.8], atol=1e-12)
assert factor_residual_cotangents[1].size == 0
np.testing.assert_allclose(factor_residual_cotangents[2], [-0.4], atol=1e-12)

First factor: [-0.8]
Empty factor slot: []
Second factor: [-0.4]

A factor residual direction specifies how the whitened factor residuals change for a chosen parameter direction, holding the solved state fixed. Increasing $z_0$ at unit rate gives the first factor direction $-1/\sigma_0$ and leaves the second factor unchanged. pushforward returns the corresponding change in the optimized state. Pairing the factor residual direction with the pullback result gives the same loss derivative as parameterGradient.

factor_residual_direction = [np.array([-1.0 / sigmas[0]]), np.empty(0), np.zeros(1)]
state_direction = gradients.pushforward(factor_residual_direction)
np.testing.assert_allclose(state_direction.at(key), observation_gradient, atol=1e-12)
np.testing.assert_allclose(
    gtsam.dot(factor_residual_cotangents, factor_residual_direction),
    observation_gradient[0], atol=1e-12,
)
print("State change per unit observation change:", state_direction.at(key))

State change per unit observation change: [0.8]

Whitening distinguishes the low-level queries from the parameter callback. For raw factor residuals $e$ and fixed whitening matrix $W$ , the whitened factor residuals are $r=We$ . A pullback result $\bar r$ becomes $W^T\bar r$ before it can be passed to a raw factor VJP. parameterGradient performs this conversion automatically, including for correlated Gaussian noise. Signs and coordinate order must match the factor’s raw error callback.

The query methods differ in the quantities they accept and return:

Method	Supply	Receive
`parameterGradient(g, tag)`	Keyed state loss gradients and a parameter tag	A parameter loss gradient vector, or `None` if no factor recognizes the tag
`pullback(g)`	Keyed state loss gradients	Loss gradients on whitened factor residuals, grouped by factor
`pushforward(q)`	Whitened factor residual directions, grouped by factor	Keyed changes in optimized states
`adjoint(g)`	Keyed state loss gradients	An intermediate state vector; see the advanced section

State results contain every active key. Unknown keys, incorrect dimensions, nonfinite inputs, and incorrect factor block counts raise errors. A factor with no state variables retains its factor residual block but has zero state response.

Constructors and factorization reuse¶

The default constructor factors the linearized system using Cholesky elimination and an automatic key ordering called COLAMD. QR is an alternative factorization, selected with gtsam.Gradients.Factorization.QR. An explicit Ordering lists each active state key exactly once; extra entries in the solution Values are ignored.

The Gaussian constructors accept explicit factor residual rows or an existing Gaussian Bayes tree, which stores the factorization used for sparse solves. For the ordinary factor residuals here, graph.linearize(solution) produces the same factor residual rows as the nonlinear constructor.

ordering = gtsam.Ordering()
ordering.push_back(key)
qr = gtsam.Gradients.Factorization.QR
ordered = gtsam.Gradients(graph, solution, ordering, qr)
from_rows = gtsam.Gradients(linearized, ordering, qr)
tree = linearized.eliminateMultifrontal(ordering)
from_tree = gtsam.Gradients(linearized, tree)

The nonlinear constructors retain both the factor residual Jacobians and the context needed for parameter VJPs. They obtain raw factor residual Jacobians from each active NoiseModelFactor, apply Gaussian whitening, and retain those rows even if a factor has a reduced linearize() override. The sum of squared factor residuals must match the optimized objective.

The Gaussian constructors support the lower-level queries using explicit JacobianFactor rows. They lack nonlinear factor callbacks and therefore reject parameterGradient. If the supplied rows combine or project original factor residuals, the application supplies any mapping back to those errors. Hessian-only factors do not define factor residual coordinates and are unsupported.

Construction requires finite, positive-definite Gaussian noise and a nonsingular state system. Robust losses and hard constraints are unsupported. The caller optimizes the graph and checks stationarity before construction.

Gradients owns its numerical factor residual rows and normally computes its own Bayes tree. The existing-tree constructor copies tree topology and shares conditional and noise-model pointers, following native GaussianBayesTree ownership. Shared numerical data must remain unchanged. A reused tree must factor the same undamped system; matching keys and dimensions cannot establish numerical freshness, especially for a partially relinearized iSAM2 tree.

Nonlinear construction also retains the original factors by shared ownership and copies the solved Values for VJP calls. The source containers can go out of scope, but the factors, parameter values captured by callbacks, and noise models must remain unchanged during queries. Query workspaces are independent; concurrent parameter queries additionally require thread-safe callbacks.

for operator in (gradients, ordered, from_rows, from_tree):
    for _ in range(3):
        for actual, expected in zip(operator.pullback(state_cotangents), factor_residual_cotangents):
            np.testing.assert_allclose(actual, expected, atol=1e-12)
print("All constructor choices give the same factor residual loss gradients.")

All constructor choices give the same factor residual loss gradients.

Advanced: implicit differentiation and the Gauss-Newton approximation¶

The implicit-function theorem (IFT) relates changes in optimized variables to changes in factor parameters by differentiating the equation satisfied at a stationary solution. Stationarity means that the estimation objective’s gradient is zero.

State derivatives use the local coordinates in which GTSAM perturbs the supplied estimate, called its tangent coordinates. Hold that coordinate basis fixed throughout one derivative calculation. Let $X$ collect the state variables, $\theta$ the factor parameters, $X^\star$ the supplied solution, and $r(X,\theta)$ the whitened factor residual vector. The matrices $J=\partial_Xr$ and $P=\partial_\theta r$ are the factor residual Jacobians with respect to states and parameters, respectively. The equations collect factor residual coordinates into mathematical vectors; the API keeps them grouped by factor in Errors.

The gradient of the estimation objective provides the equation to differentiate. Write $\phi(X,\theta)=\tfrac12\|r(X,\theta)\|^2$ for that objective and $F=\partial_X\phi$ for its state gradient. At the solution,

F(X^\star,\theta)=J^T r=0.

(1)

Differentiating this zero-gradient equation relates a parameter change to the corresponding change in the optimized states. Let $\dot\theta$ be a parameter change direction and $\dot X$ the corresponding state change. With $A=\partial_XF$ and $B=\partial_\theta F$ , the equation is $A\dot X+B\dot\theta=0$ . When $A$ is nonsingular, the IFT gives

\dot X=-A^{-1}B\dot\theta.

(2)

The derivative terms used by Gradients¶

The exact state derivative $A$ and mixed state-parameter derivative $B$ each contain a term weighted by the factor residual values. Writing $r_i$ for factor residual coordinate $i$ , their expressions are

A=J^TJ+\sum_i r_i\,\partial^2_{XX}r_i,\qquad B=J^TP+\sum_i r_i\,\partial^2_{X\theta}r_i.

(3)

Gradients uses $H=J^TJ$ in place of $A$ and $J^TP$ in place of $B$ . These are the Gauss-Newton terms. The approximation therefore applies to both the state Hessian, the second state derivative of the estimation objective, and the mixed state-parameter derivative. An exact state Hessian alone does not make the parameter response exact if the mixed derivative still has a term weighted by factor residuals.

Forward mode applies the approximate state derivative to a parameter direction. Define the induced factor residual direction as $q=P\dot\theta$ ; then pushforward(q) computes

\dot X=-H^{-1}J^Tq.

(4)

Reverse mode propagates the supplied state loss gradient through the same linear system. Let $g$ be that gradient, and let the intermediate state vector $\lambda$ , called the state adjoint, solve $H\lambda=g$ . adjoint(g) returns $\lambda$ , and pullback(g) returns the factor residual loss gradient $\bar r$ :

H\lambda=g,\qquad \bar r=-J\lambda.

(5)

The total parameter loss gradient adds the direct loss derivative to the contribution through the states. With the direct term evaluated at fixed $X$ ,

\nabla_\theta L= \left.\partial_\theta L\right|_{X\ \mathrm{fixed}}+P^T\bar r.

(6)

Two triangular solves through the Bayes-tree factorization apply $H^{-1}$ . Queries reuse that factorization without assembling a global Jacobian, Hessian, or inverse. parameterGradient evaluates $P^T\bar r$ through the factors’ VJP methods; the application adds any direct parameter loss derivative. For a raw factor residual $e$ with fixed whitening $r=We$ , it passes $\bar e=W^T\bar r$ to parameterVjp, which returns $(\partial_\theta e)^T\bar e$ . This also handles correlated Gaussian noise, where $W$ need not be symmetric.

Agreement with exact differentiation¶

Gauss-Newton and exact implicit differentiation agree when the omitted terms weighted by factor residuals vanish. This happens when all factor residuals are zero or when every factor residual is a linear function plus a constant in the combined state and parameter coordinates, called a jointly affine factor residual. The scalar prior example and odometry tutorial have this form.

For nonlinear problems with nonzero factor residuals, exact reoptimization can give a different derivative from the Gauss-Newton calculation. Convergence establishes stationarity; it does not remove the second derivatives weighted by factor residuals.

Independent columns of $J$ , also called full column rank, ensure that $H=J^TJ$ is nonsingular. The operations also require consistent Gaussian whitening, fixed tangent coordinates, and a matching undamped factorization. They apply no damping or pseudoinverse fallback.

state_adjoint = gradients.adjoint(state_cotangents)
np.testing.assert_allclose(state_adjoint.at(key), [0.8], atol=1e-12)
tree_adjoint = tree.backSubstitute(tree.backSubstituteTranspose(state_cotangents))
np.testing.assert_allclose(tree_adjoint.at(key), state_adjoint.at(key), atol=1e-12)
print("Intermediate state adjoint:", state_adjoint.at(key))

Intermediate state adjoint: [0.8]

Loss gradients for manifold-valued states¶

Pose loss gradients must be expressed in the local update coordinates used by the optimizer at the solved pose. For example, a Pose2 update has two translation coordinates and one rotation coordinate. A loss written in other coordinates needs a derivative that converts it to those local coordinates.

A comparison with a reference pose contributes its own Jacobian to that conversion. Let $e_i=\operatorname{Local}_{x_i^{\mathrm{ref}}}(x_i)$ be the comparison error, $V_i$ a fixed symmetric loss-weight matrix, and $J_{e_i}$ the Jacobian of that error with respect to local updates at the solved $x_i$ . For loss $\tfrac12e_i^TV_i e_i$ , supply $J_{e_i}^TV_i e_i$ as the state cotangent. This comparison Jacobian is generally not identity. The scalar example’s ordinary coordinates need no such conversion.