This guide explains GTSAM’s reusable imperative multifrontal solver for Gaussian factor graphs. It covers the solver lifecycle, supported inputs, structural and numerical choices, parallel traversal, diagnostics, and every MultifrontalParameters option with its default.
import gtsam
import numpy as np
from gtsam.symbol_shorthand import X
SCALAR_IDENTITY = np.eye(1)
UNIT_NOISE = gtsam.noiseModel.Unit.Create(1)Where the solver fits¶
A Gaussian factor graph represents a least-squares objective of the form
Multifrontal elimination uses a variable ordering to construct a symbolic junction tree. Each numerical clique gathers its local factors and child updates, eliminates its frontal variables, and passes a Schur-complement update to its parent. A top-down pass then back-substitutes through the resulting conditionals.
MultifrontalSolver is the reusable C++ implementation of this process. It separates symbolic setup and fixed solve allocation from repeated numerical loads, factorizations, and solves. On the first load, each clique builds cached load plans and lazily allocates packed Jacobian fallback storage; later compatible loads reuse it. The class itself is not currently wrapped in Python, so the executable cells below use the public Python Gaussian-factor APIs to illustrate the same problem and reference solution. The C++ sections show the MultifrontalSolver-specific API.
A small Gaussian chain¶
The first factor anchors at zero. Each binary factor then constrains . Eliminating the graph multifrontally should recover the scalar sequence 0, 1, 2, 3.
graph = gtsam.GaussianFactorGraph()
graph.add(X(0), SCALAR_IDENTITY, np.array([0.0]), UNIT_NOISE)
for index in range(3):
graph.add(
X(index),
-SCALAR_IDENTITY,
X(index + 1),
SCALAR_IDENTITY,
np.array([1.0]),
UNIT_NOISE,
)
ordering = gtsam.Ordering()
for index in range(4):
ordering.push_back(X(index))
bayes_tree = graph.eliminateMultifrontal(ordering)
solution = bayes_tree.optimize()
estimated = np.array([solution.at(X(index))[0] for index in range(4)])
np.testing.assert_allclose(estimated, np.arange(4.0), atol=1e-12)
estimatedC++ lifecycle¶
A solver instance is reusable while the graph’s keys, factor types, factor-to-key structure, dimensions, row layout, and ordering remain unchanged. Only the numerical factor values may change between calls to load. QR cliques materialize every factor row and reserve frontal damping rows. Cholesky cliques materialize only Jacobian fallback rows; eligible batch factors update the Hessian directly and consume no packed Jacobian rows.
#include <gtsam/linear/MultifrontalSolver.h>
using namespace gtsam;
MultifrontalSolver::Parameters parameters;
MultifrontalSolver solver(graph, ordering, parameters);
// Separate load and elimination: useful when timing or inspecting phases.
solver.load(graph);
solver.eliminateInPlace();
const VectorValues& delta = solver.updateSolution();
// A later graph may reuse the symbolic structure and allocated storage.
solver.load(graphWithUpdatedNumbers);
solver.eliminateInPlace();
const VectorValues& updatedDelta = solver.updateSolution();The overload eliminateInPlace(graph) fuses loading and elimination into one post-order traversal. Use computeBayesTree() after elimination when a persistent GaussianBayesTree is needed. After updateSolution(), deltaError() returns the cached change in linearized error and can optionally return the old and new errors.
Reusing symbolic precomputation¶
Precompute() separates symbolic analysis from solver construction. This is useful when several solver instances share one graph structure and ordering.
auto data = MultifrontalSolver::Precompute(graph, ordering);
MultifrontalSolver solver(std::move(data), ordering, parameters);
solver.load(graph);
solver.eliminateInPlace();
const VectorValues& delta = solver.updateSolution();Calling eliminateInPlace() before load(), computeBayesTree() before elimination, or deltaError() before updateSolution() is an error. The solution returned by updateSolution() is owned by the solver and remains valid only as a reference to its cached storage.
Partial elimination and retained factors¶
A partial solver eliminates only the leading firstPhaseSize keys and leaves the remaining variables assembled in clique-level Hessian factors.
MultifrontalSolver solver(graph, ordering, firstPhaseSize, parameters);
solver.eliminatePartialInPlace(graph);
GaussianFactorGraph remaining = solver.remainingFactorGraph();Each exported HessianFactor owns a compact retained information matrix, so it remains valid after clique views change or the solver is destroyed. Export performs one unavoidable copy of the active upper triangle directly into that final owner, then transfers ownership through the rvalue SymmetricBlockMatrix constructor; it does not copy the complete clique matrix or create an intermediate dense self-adjoint matrix. The existing const-reference HessianFactor constructor remains available when callers need to retain their matrix. For end-to-end comparisons, include remainingFactorGraph() immediately after partial elimination.
Complete parameter reference¶
The following defaults come directly from MultifrontalParameters.h.
| Option | Type | Default | Meaning |
|---|---|---|---|
leafMergeDimCap | size_t | 256 | Algebraically merges compatible multi-factor sibling leaves with identical separators into dimension-capped cliques. One-factor leaves remain separate. 0 disables this pass. |
leafAggregationProblemSize | size_t | 2048 | Total problem-size budget for grouping independent sibling leaves into one scheduled bottom-up task. 0 disables leaf task aggregation. |
leafMode | LeafMode | LeafMode::Bounded | Bounded aggregates only scheduled work. SameSeparator also accumulates eligible Cholesky leaf updates sharing a separator before one parent scatter. Neither mode merges the algebraic cliques. |
mergeDimCap | size_t | 32 | Bottom-up general clique-merge threshold. A child is merged into its parent when the child’s frontal dimension plus the parent’s total dimension is strictly less than this value. 0 disables this pass. |
qrMode | QRMode | QRMode::Allow | Chooses whether eligible leaf cliques use QR or Cholesky: Off, Allow, or Force. Details appear below. |
qrAspectRatio | double | 2.0 | In Allow mode, an eligible leaf uses QR when its total frontal-plus-separator dimension is greater than this value times its frontal dimension. |
compactCholeskySeparatorDimThreshold | size_t | 256 | Minimum separator dimension for eligible non-star Cholesky leaves to avoid materializing the separator Hessian. Fused one-frontal-block star leaves bypass this cutoff. |
reportStream | std::ostream* | nullptr | Receives symbolic structure summaries during construction. A null pointer disables reporting. The caller owns the stream. |
eliminationParallelThreshold | int | 10 | Problem-size threshold used when creating tasks in the bottom-up elimination traversal. Smaller values expose more task parallelism; larger values keep more work in recursive tasks. |
solutionParallelThreshold | int | 4096 | Problem-size threshold used when creating tasks in the top-down solution traversal. |
numThreads | size_t | 0 | Worker count. Zero selects three quarters of std::thread::hardware_concurrency(), rounded down and clamped to at least one worker. A positive value requests exactly that many workers. |
Leaf scheduling and structural merge options¶
leafMergeDimCap controls the first symbolic pass. It groups multi-factor leaf siblings with identical separators in stable order, packs them under the dimension cap, and replaces each group with one algebraic clique. One-factor leaves are excluded so a sole direct batch factor can still select the fused leaf path after load planning.
mergeDimCap runs next as a general symbolic bottom-up pass that can merge non-leaf children into parents. After both symbolic passes, leafAggregationProblemSize bounds groups of remaining independent sibling leaves processed by one scheduled bottom-up task. LeafMode::Bounded changes scheduling only. LeafMode::SameSeparator also accumulates eligible Cholesky updates with identical separators in separator-local information matrices before one parent scatter; it does not merge those cliques. Set each cap to zero independently to disable that operation.
QR and Cholesky selection¶
QR selection applies only to leaf cliques with at least as many factor rows as frontal scalar dimensions.
QRMode::Off always uses Cholesky.
QRMode::Allow, the default, uses QR for an eligible leaf when frontalDim + separatorDim > qrAspectRatio * frontalDim, except that a one-frontal leaf containing exactly one directly updatable batch factor uses the fused Cholesky path.
QRMode::Force uses QR for every eligible leaf. A clique that is not an eligible leaf still uses Cholesky.
QR avoids explicitly forming normal equations for the selected leaves and can be preferable for tall local Jacobians or when numerical conditioning matters. It materializes every factor row and reserves frontal damping capacity in the clique’s Ab matrix. Cholesky is usually cheaper when the clique is not tall. Its Ab matrix contains only fallback rows: ordinary Jacobian factors and batch factors that cannot use direct Hessian updates. Direct batch factors consume no Ab rows. qrAspectRatio has no effect in Off or Force mode.
A fully eliminated non-root Cholesky leaf with one frontal block, directly updatable batch factors, and optional frontal-only Jacobian factors uses the fused star path at every separator size. It forms and factors only [Hff Hfs gf], adds original retained terms directly to the destination, and subtracts [S d]ᵀ[S d] without a temporary separator Hessian. Other eligible leaves use the generic compact path when the separator is larger than the frontal block and reaches compactCholeskySeparatorDimThreshold. Both paths retain only necessary fallback rows in Ab. Leaf kernels are single-threaded; concurrency remains at the parent child-gather level.
Parallel traversal and reporting¶
The two parallel thresholds are scheduling controls, not matrix dimensions or thread counts. Each clique reports a problem size equal to its frontal dimension plus separator dimension. The traversal uses the applicable threshold to decide where additional tasks are worthwhile. Their best values depend on tree shape, clique size, scheduler overhead, and hardware; tune them with representative graphs rather than assuming that more tasks are always faster.
numThreads controls the solver’s worker pool independently of those thresholds. Set it to 1 for deterministic single-worker performance measurements. The default 0 derives the worker count from available hardware.
When reportStream is non-null, construction prints summaries of the original symbolic cluster structure and the structures after each enabled merge pass. This is useful for confirming whether a parameter actually changed clique count, dimensions, or fan-out.
Configuring every option¶
This example writes every field explicitly. The assigned values are the defaults, so removing any line preserves the same behavior.
MultifrontalSolver::Parameters parameters;
parameters.leafMergeDimCap = 256;
parameters.leafAggregationProblemSize = 2048;
parameters.leafMode = MultifrontalParameters::LeafMode::Bounded;
parameters.mergeDimCap = 32;
parameters.qrMode = MultifrontalParameters::QRMode::Allow;
parameters.qrAspectRatio = 2.0;
parameters.compactCholeskySeparatorDimThreshold = 256;
parameters.reportStream = nullptr;
parameters.eliminationParallelThreshold = 10;
parameters.solutionParallelThreshold = 4096;
parameters.numThreads = 0;
MultifrontalSolver solver(graph, ordering, parameters);Supported factors and constraints¶
MultifrontalSolver accepts JacobianFactor and BatchJacobianFactor inputs. Other Gaussian factor types, including HessianFactor, are rejected during precomputation or loading.
A constrained JacobianFactor is supported only when it is unary, fully constrained, and has a zero residual. Such a factor fixes its key to zero during the linear solve. Mixed-key, partially constrained, or infeasible constrained Jacobian factors are rejected. Constrained BatchJacobianFactor inputs are not supported.
The graph passed to load() must preserve the structure used during construction: factor positions, factor types, keys, dimensions, and row layout must match. Use a new precomputation and solver if the structure changes.
Using it from a nonlinear optimizer¶
In C++, Gauss-Newton and Levenberg-Marquardt can select the reusable solver through NonlinearOptimizerParams::MULTIFRONTAL_SOLVER. The multifrontal parameters are stored in the optimizer parameters.
GaussNewtonParams parameters;
parameters.linearSolverType =
NonlinearOptimizerParams::MULTIFRONTAL_SOLVER;
parameters.multifrontalParams.qrMode =
MultifrontalParameters::QRMode::Allow;
parameters.multifrontalParams.numThreads = 0;
GaussNewtonOptimizer optimizer(nonlinearGraph, initialValues, parameters);
const Values result = optimizer.optimize();LM damping policy is configured separately through LevenbergMarquardtParams. MultifrontalParameters controls symbolic structure, numerical leaf selection, traversal, worker count, and reporting; it does not choose the LM damping formula.
Practical tuning workflow¶
Start with the defaults and a representative graph. If performance matters, record end-to-end runtime and peak memory, not elimination time alone. Then change one family of controls at a time:
Set reportStream to std::cout and inspect the symbolic structure.
Compare leafMergeDimCap while keeping general merging and task aggregation fixed.
Compare leaf aggregation disabled, Bounded, and SameSeparator while keeping both symbolic merge caps fixed.
Compare mergeDimCap and compactCholeskySeparatorDimThreshold with one change at a time.
Compare QRMode::Off, Allow, and Force; tune qrAspectRatio only for Allow.
Fix numThreads while tuning traversal thresholds so worker-count changes do not confound the result.
For partial elimination, make elimination plus immediate retained-factor export the primary timing. Keep elimination-only and export-only timings as diagnostics. The timeSfmPartialElimination target reports these alongside compact complete assembly.
Validate every tuned configuration against the untuned solution and linearized error.
The defaults are general-purpose heuristics, not universal optima. Chains, balanced trees, and bundle-adjustment graphs have very different clique shapes.
Summary¶
MultifrontalSolver is most useful when the same symbolic Gaussian graph structure is solved repeatedly with new numerical values. Keep the ordering and graph structure fixed, reuse precomputation and lazily packed numerical storage, choose QR only where its leaf eligibility rules apply, and treat leaf scheduling, compact Cholesky, merge, and traversal settings as measurable performance heuristics. For partial elimination, compare complete pipelines including retained-factor export. For failures caused by rank deficiency or poor conditioning, see Debugging an Indeterminate Linear System.