Debugging an Indeterminate Linear System

This guide shows how to diagnose IndeterminateSystemException, which indicates that a linearized factor graph is underdetermined, indefinite, or too poorly conditioned to solve reliably. We will reproduce a missing-gauge failure, inspect the Jacobian null space, and repair the model.

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

POSE2_DIMENSION = 3

What the exception tells you¶

In C++, catch gtsam::IndeterminateSystemException and call nearbyVariable(). Python currently surfaces the failure as a RuntimeError whose message names the nearby key. The reported key is where elimination detected the problem; it is not necessarily where the modeling error originated. Its identity can change with the elimination ordering.

The exception is not limited to exactly singular systems. GTSAM also raises it for nearly indeterminate systems that are mathematically full rank but too poorly conditioned to solve reliably. Cholesky elimination checks every pivot relative to its original diagonal entry; this normalization is invariant to diagonal changes of variable units but still depends on elimination ordering. For example, a prior with sigmas near 10^-6 combined with measurement models whose sigmas are near 10² produces raw Jacobian weights spanning roughly eight orders of magnitude (and a raw Hessian condition number that can span roughly sixteen). That disparity can expose a weakly observed direction, but it may also reflect nothing more than different variable units; the raw condition number alone is therefore not sufficient evidence.

The former C++ spelling, IndeterminantLinearSystemException, is available only when the GTSAM 4.3 deprecated API is enabled.

Reproduce and preserve the failing graph¶

This graph constrains the relative pose between x0 and x1, but nothing anchors the global translation or rotation. Preserve both the graph and the exact values that failed; linearization depends on those values.

odometry_noise = gtsam.noiseModel.Diagonal.Sigmas(
    np.array([0.1, 0.1, 0.1])
)

graph = gtsam.NonlinearFactorGraph()
graph.add(
    gtsam.BetweenFactorPose2(
        X(0), X(1), gtsam.Pose2(1.0, 0.0, 0.0), odometry_noise
    )
)

values = gtsam.Values()
values.insert(X(0), gtsam.Pose2(0.2, 0.4, 0.1))
values.insert(X(1), gtsam.Pose2(1.5, -0.2, -0.1))

try:
    gtsam.GaussNewtonOptimizer(graph, values).optimize()
except RuntimeError as exception:
    first_line = next(line for line in str(exception).splitlines() if line)
    print(first_line)
    print("Graph keys:", [gtsam.DefaultKeyFormatter(k) for k in graph.keyVector()])

Inspect the Jacobian before the Hessian¶

Linearize at the failing values and choose an explicit ordering so the matrix columns have a known key order. The Jacobian $A$ is preferable for rank tests because forming $A^T A$ squares the condition number. A rank smaller than the number of columns means the system has a null space.

ordering = gtsam.Ordering()
for key in (X(0), X(1)):
    ordering.push_back(key)

linear_graph = graph.linearize(values)
A, b = linear_graph.jacobian(ordering)
_, singular_values, Vt = np.linalg.svd(A, full_matrices=True)
tolerance = singular_values[0] * max(A.shape) * np.finfo(A.dtype).eps
rank = int(np.sum(singular_values > tolerance))
nullity = A.shape[1] - rank

print(f"Jacobian shape: {A.shape}")
print(f"rank: {rank}, nullity: {nullity}")
print("singular values:", singular_values)

Here the Jacobian has six columns (two Pose2 variables with three tangent dimensions each) but only rank three. The three null directions are the unconstrained global $x$ , $y$ , and heading gauges.

For a mixed graph, map column ranges using each variable’s tangent dimension. Large block norms in a null vector identify variables participating in an unobservable direction; they do not by themselves identify a faulty factor.

null_space = Vt[rank:]
for direction_index, direction in enumerate(null_space):
    block_norms = {}
    for block_index, key in enumerate((X(0), X(1))):
        start = block_index * POSE2_DIMENSION
        stop = start + POSE2_DIMENSION
        block_norms[gtsam.DefaultKeyFormatter(key)] = np.linalg.norm(
            direction[start:stop]
        )
    print(f"null direction {direction_index}: {block_norms}")

Repair the model and verify it¶

A prior removes the global gauge in this example. A diagnostic finite prior can confirm the missing degrees of freedom, but making it arbitrarily strong is not a robust substitute for an exact constraint: combined with much looser factor noise, it can expose a weak direction that produces a very small normalized elimination pivot.

When the modeling intent is to fix the gauge exactly, as it often is for a gauge-fixing prior, use a hard constraint such as noiseModel.Constrained.All(dimension). It states that intent directly and does not introduce an extreme finite weight, so the gauge constraint itself will not trigger the conditioning check. The exception can still reveal a separate underconstraint or ill-conditioned direction elsewhere in the graph.

fixed_graph = gtsam.NonlinearFactorGraph(graph)
prior_noise = gtsam.noiseModel.Constrained.All(POSE2_DIMENSION)
fixed_graph.add(gtsam.PriorFactorPose2(X(0), gtsam.Pose2(), prior_noise))

fixed_linear_graph = fixed_graph.linearize(values)
fixed_A, _ = fixed_linear_graph.jacobian(ordering)
print(f"fixed Jacobian rank: {np.linalg.matrix_rank(fixed_A)} / {fixed_A.shape[1]}")

result = gtsam.GaussNewtonOptimizer(fixed_graph, values).optimize()
print("optimized x0:", result.atPose2(X(0)))
print("optimized x1:", result.atPose2(X(1)))

Checklist for real systems¶

Save the exact graph and values at failure. Also record the elimination ordering and incremental or fixed-lag update that triggered it.
Inspect graph structure. Look for values absent from every factor, accidental new keys, disconnected components, and missing priors or gauge constraints. Start near the reported key, but trace the whole connected component.
Linearize at the failing values. Compute the Jacobian singular spectrum and null space with an explicit ordering. Prefer the Jacobian over a determinant or the Hessian for numerical rank diagnosis.
Map weak directions back to variable blocks. Check whether they represent an expected gauge or a missing observable degree of freedom. Degenerate triangulation, small camera baselines, straight-line or zero-excitation IMU motion, and aggressive marginalization are common causes.
Check numerical scale and curvature. Very different units or noise sigmas can make the raw matrix appear ill-conditioned; distinguish removable diagonal scaling from weak directions that produce a small normalized elimination pivot. Remember that the pivot test depends on elimination ordering. A very strong finite prior together with loose measurement noise is a common example to investigate. If custom HessianFactors are present, inspect the symmetric Hessian eigenvalues for unintended negative curvature.
Test one modeling hypothesis at a time. Temporarily add a physically meaningful prior or remove a suspect measurement block, then recompute rank. If the prior is intended to fix a gauge exactly, replace it with a hard constrained noise model rather than an extreme finite weight. Damping or dense optimization may make one solve succeed without restoring observability, so do not treat that alone as a fix.
For incremental estimators, replay updates. Find the first update that loses rank and inspect factors removed or marginalized at that step. Observability can be lost over time even when every individual measurement is valid.