Robust Estimation with RANSAC#

cuNLS ships two RANSAC minimizers, RansacGaussNewtonMinimizer and RansacLevenbergMarquardtMinimizer, with the common base RansacMinimizer (Python pycunls; C++ header cunls/minimizer/ransac_minimizer.h). They solve the same Problem as the regular minimizers, but they assume that some measurements are gross outliers: wrong data association, not just noise. They find the estimate that most measurements agree with, refine it on those measurements, and report which measurements were inliers.

This page explains when to use them, the theory behind them, exactly what the implementation does, how to call it from Python and C++, and how to tune it. How to write custom factors and states that work with RANSAC is covered in Custom Factors and States (Python and C++). Complete runnable programs are in python/examples/ransac_pnp.py (Python) and examples/ransac_pnp (C++); python/examples/tartan_vio.py combines RANSAC with an IMU factor and priors (always-on batches) in a visual-inertial odometry on real data.

When to use RANSAC #

A least-squares solver weighs every residual by its square. A single measurement that is wrong by 100 noise standard deviations contributes as much as 10,000 good ones, so a few outliers move the solution arbitrarily far. There are two families of remedies:

Robust losses (Huber, Cauchy, Tukey, …; see Robustifier API): down-weight large residuals inside the ordinary solver. They are cheap and work well when the outlier ratio is moderate and the initial guess is already close: the loss decides what “large” means from the current estimate, so a poor start can make the outliers look like the inliers.
RANSAC (RANdom SAmple Consensus) does not trust the initial estimate: to tell inliers from outliers. It generates many candidate solutions (“hypotheses”) from tiny random subsets of the measurements, keeps the one that most measurements agree with, and only then refines it. It tolerates very high outlier ratios (80–90% in the PnP benchmarks) and returns an explicit inlier/outlier classification.

Use the RANSAC minimizers when:

a significant fraction of the measurements can be completely wrong (feature mismatches, wrong loop closures, spurious returns);
you need the inlier set, not just the estimate;
the free state is small: the sum of the tangent dimensions of all non-constant states must be at most kMaxRansacTangentDim = 64 (a camera pose is 6, a pose plus a focal length is 7, a rig of 10 poses is 60). Constant states (e.g. known 3D points) do not count, however many there are, and the number of factors is unlimited.

Typical problems: PnP (pose from 3D-2D matches), point-cloud registration (pose from 3D-3D matches), relative pose / extrinsic calibration, fitting a low-dimensional model (line, plane, homography-like) to data, a multi-camera rig pose with known extrinsics.

Theory #

The problem #

The regular minimizers solve

\[\min_x \; \frac{1}{2} \sum_i \rho_i\!\left(\|r_i(x)\|^2\right)\]

over the free state \(x\) (dimension \(D\)) for residual blocks \(r_i\) of dimension \(m_i\). RANSAC assumes the residual blocks are of two kinds:

inliers: \(r_i(x^\star)\) at the true state \(x^\star\) is small, of the size of the measurement noise;
outliers: \(r_i(x^\star)\) is arbitrary.

Which is which is unknown. The goal is the state that explains as many measurements as possible within the noise level, and that partition.

In cuNLS every residual batch has a role:

RansacRole.sampled (C++ RansacRole::kSampled): Data factors that may be outliers. RANSAC samples from them and classifies every one of them. Each sampled batch has an inlier threshold \(\tau\): factor \(i\) is an inlier of state \(x\) iff \(\|r_i(x)\| \le \tau\).
RansacRole.always_on (C++ RansacRole::kAlwaysOn): Trusted factors: priors, motion models, known extrinsic constraints. They are part of every solve and never classified.

Hypotheses from minimal samples #

A minimal sample is a set of \(s\) sampled factors that, together with the always-on factors, determines the state. Each factor contributes \(m\) equations, so \(s = \lceil D / m_{\min} \rceil\) factors are enough in general (\(m_{\min}\) is the smallest residual dimension among the sampled batches). For PnP, \(D = 6\) and \(m = 2\), so \(s = 3\) correspondences.

If all \(s\) factors of a sample are inliers, the state fitted to them is close to \(x^\star\), and most other inliers will agree with it. If any is an outlier, the fitted state is essentially random and few factors agree with it. RANSAC therefore draws many samples and keeps the hypothesis with the most support.

Classic RANSAC fits each sample with a problem-specific closed-form “minimal solver” (e.g. P3P for PnP). cuNLS instead fits each sample with a few Gauss-Newton or Levenberg-Marquardt iterations from the initial guess, using the problem’s own factors. This needs no problem-specific code, so any factor (built-in or custom) works, at the price of needing an initial guess inside the convergence basin of the sample problem. In the PnP benchmarks this basin is wide: rotations off by 0.3 rad and translations off by 10% of the depth converge reliably.

Scoring #

Every hypothesis \(x_k\) is scored against all sampled factors. The default rule is MSAC (M-estimator SAmple Consensus), a truncated quadratic:

\[\mathrm{score}(x_k) = \sum_{i \in \text{sampled}} \min\!\left(\|r_i(x_k)\|^2,\ \tau_i^2\right) \quad (+\ 2 \cdot \mathrm{cost}_{\text{always-on}}(x_k))\]

Lower is better. Inliers contribute their (small) squared residual, outliers a constant \(\tau^2\). Compared with counting inliers, MSAC also prefers the hypothesis whose inliers fit more tightly, which breaks ties between hypotheses with the same support. RansacScoring::kInlierCount counts inliers instead (ties broken by MSAC). When score_always_on is set (the default), the always-on cost is added so that hypotheses violating a trusted prior are penalized.

Informative inliers. Some factors report a zero residual for configurations they cannot evaluate. PnPFactorBatch, for example, returns zero residual and zero Jacobian for points behind the camera. A hypothesis that puts every point behind the camera would then look perfect. With require_informative_inliers (default on), a factor counts as an inlier only if its Jacobian has a non-zero entry on a free state.

How many hypotheses?#

Let \(w\) be the inlier ratio. A random sample of \(s\) factors is all-inlier with probability \(w^s\). To draw at least one all-inlier sample with probability \(p\) (confidence), one needs

\[K \ \ge\ \frac{\log(1 - p)}{\log(1 - w^s)}\]

hypotheses. For \(p = 0.999\):

inlier ratio \(w\)	\(s = 2\)	\(s = 3\)	\(s = 4\)	\(s = 6\)
0.9	5	6	7	10
0.5	25	52	108	439
0.3	74	253	850	9,473
0.2	170	861	4,314	107,931
0.1	688	6,905	69,075	\(6.9 \cdot 10^6\)

\(w\) is unknown in advance, so cuNLS works in rounds of hypotheses_per_round (default 256) hypotheses and stops adaptively: after each round it sets \(w\) to the best inlier ratio found so far and stops as soon as the hypotheses drawn reach the bound above (or the inlier ratio reaches early_stop_inlier_ratio, or max_rounds is hit). The table shows why small samples matter: prefer factors with larger residual dimension (fewer factors per sample) and keep \(D\) small.

Choosing the inlier threshold #

\(\tau\) is in the units of the raw residual (before any loss function), and it is the single most important parameter. If the residual of an inlier is Gaussian with per-component standard deviation \(\sigma\), then \(\|r\|^2 / \sigma^2\) follows a \(\chi^2\) distribution with \(m\) degrees of freedom, and

\[\tau = \sigma \sqrt{\chi^2_m(q)}\]

keeps a fraction \(q\) of the true inliers:

residual dim \(m\)	\(q = 0.95\)	\(q = 0.99\)	\(q = 0.999\)
1	\(1.96\sigma\)	\(2.58\sigma\)	\(3.29\sigma\)
2	\(2.45\sigma\)	\(3.03\sigma\)	\(3.72\sigma\)
3	\(2.80\sigma\)	\(3.37\sigma\)	\(4.03\sigma\)
6	\(3.55\sigma\)	\(4.10\sigma\)	\(4.74\sigma\)

Too small a threshold rejects good measurements and makes the estimate noisy; too large a threshold accepts outliers that lie close to the model. If the noise is not isotropic, wrap the factor batch in an InformationFactorBatch with the square-root information matrix: the residual is then whitened (\(\sigma = 1\)) and \(\tau = \sqrt{\chi^2_m(q)}\) directly.

What the implementation does #

One call to minimize(stream, problem) (C++ Minimize):

Validate and lay out the problem (see Limits and requirements): find the free states and their tangent columns (\(D\) total), the role of every residual batch, the sample size \(s\).
Rounds (at most max_rounds), each with \(K\) = hypotheses_per_round hypotheses solved in parallel on the GPU:
1. Every hypothesis gets its own copy (“replica”) of the free state batches, initialized from the problem’s current state values.
2. Hypothesis \(k\) draws its sample: the first \(s\) entries of a keyed pseudo-random permutation of all sampled factors (key = seed, round, \(k\)). No sorting, no RNG state; the same seed gives the same samples.
3. hypothesis_iterations (default 5) GN or LM iterations per hypothesis on its sample plus the always-on factors. Each iteration builds the small dense normal equations \(J^\top J\,\delta = -J^\top r\) per hypothesis, solves them with a pivoted LDLᵀ (default) or Cholesky, and applies the step through the state batch’s Plus. A hypothesis whose first solve fails (degenerate sample) is marked invalid.
4. Every hypothesis is scored on all sampled factors (MSAC). For large problems, two-stage scoring (when there are more than 2 × scoring_subset_size sampled factors) first scores all hypotheses on a random subset of 16,384 factors and then scores only the best scoring_finalists (default 4) on all factors.
5. The best hypothesis so far is kept on the device. The adaptive stopping rule decides whether another round is needed.
Refine the best hypothesis on its inliers (final_iterations, default 20), re-classify, write the state back into the problem’s state batches, and fill the summary.

All hypotheses of a round are evaluated together. The minimizer does not call your factor once per hypothesis. It calls FactorBatch::Evaluate once per residual batch for all hypotheses, using the item parameters factor_ids / num_factor_ids to say which factor each output row belongs to and which hypothesis’s state it reads. Likewise it calls StateBatch::Plus once per state batch for all hypotheses, using num_replicas. This is why custom factors and states must honor those parameters (Custom Factors and States (Python and C++)). Everything else (normal equations, solves, scoring, selection) runs in a handful of fused kernels without atomics.

Deterministic. Every reduction runs in a fixed order and sampling is counter-based, so the same problem and the same seed give bitwise identical results, run after run.

Synchronization. The host reads a few scalars once per round (to decide whether to stop) and at the end. Everything else is asynchronous on the given stream.

Usage #

Important

Capacity vs. active count. Factor and state batches are constructed with their capacity (how many factors / states their buffers hold) and start with zero active entries: call set_num_active_factors(n) / set_num_active_states(n) (C++: SetNumActiveFactors / SetNumActiveStates) before solving, and again whenever the problem size changes. See Capacity and active count.

Build the problem exactly as for the regular minimizers. Then choose the roles and thresholds, run the minimizer, and read the inlier mask. The example below is the PnP problem from python/examples/ransac_pnp.py (C++: examples/ransac_pnp): one SE(3) pose state and one PnPFactorBatch with a factor per 3D-2D correspondence, all pointing to the pose.

Python #

import pycunls

# ... states, factors and problem built as usual:
#   pose_state.set_num_active_states(1)
#   pnp.set_num_active_factors(num_matches)
#   problem.add_state_batch(pose_state)
#   problem.add_factor_batch(pnp, pointers)          # residual batch 0

options = pycunls.RansacLevenbergMarquardtMinimizerOptions()
ransac = options.base_options          # a reference: edits change `options`
ransac.factor_batches = [               # assign a whole list, one entry per batch
    pycunls.RansacFactorBatchOptions(pycunls.RansacRole.sampled, 0.01),
]
ransac.seed = 1

minimizer = pycunls.RansacLevenbergMarquardtMinimizer(options)
summary = minimizer.minimize(stream, problem)     # writes the estimate back
mask = minimizer.inlier_mask(0)                   # numpy uint8, 1 = inlier
print(summary)                                     # RansacSummary(rounds=..., inliers=...)

Note

options.factor_batches returns a copy of the list: options.factor_batches.append(...) has no effect. Always assign a complete list.

C++#

#include "cunls/cunls.h"

// ... states, factors and problem built as usual:
//   pose_state.SetNumActiveStates(1);
//   pnp.SetNumActiveFactors(num_matches);
//   problem.AddStateBatch(&pose_state);
//   problem.AddFactorBatch(&pnp, pointers);          // residual batch 0

cunls::RansacLevenbergMarquardtMinimizerOptions options;
cunls::RansacMinimizerOptions &ransac = options.base_options;
// One entry per residual batch, in the order they were added to the problem.
ransac.factor_batches = {{cunls::RansacRole::kSampled, /*inlier_threshold=*/0.01f}};
ransac.seed = 1;

cunls::RansacLevenbergMarquardtMinimizer minimizer(options);
cunls::RansacSummary summary = minimizer.Minimize(stream, problem);

// The estimate is now in the problem's state batches (here: pose_state).
// Inlier mask of residual batch 0: device memory, one byte per factor.
std::vector<uint8_t> mask(minimizer.InlierMaskSize(0));
cudaMemcpy(mask.data(), minimizer.InlierMask(0), mask.size(), cudaMemcpyDeviceToHost);

printf("%zu rounds, %zu inliers (%.1f%%)\n", summary.num_rounds, summary.num_inliers,
       100.f * summary.inlier_ratio);

RansacGaussNewtonMinimizer takes a RansacMinimizerOptions directly (Python: pycunls.RansacGaussNewtonMinimizer(options) with a pycunls.RansacMinimizerOptions):

cunls::RansacMinimizerOptions options;
options.default_inlier_threshold = 0.01f;   // every batch sampled with this threshold
cunls::RansacGaussNewtonMinimizer minimizer(options);

Roles in practice #

One data batch, nothing else (PnP, registration, model fitting): leave factor_batches empty and set default_inlier_threshold; every batch is then sampled with that threshold.
Data plus priors (e.g. a pose prior from odometry): mark the prior batches RansacRole.always_on (C++ kAlwaysOn). They take part in every hypothesis solve and the refinement and, with score_always_on, in the score. If the prior may be wrong by much more than its stated uncertainty, set score_always_on = False (C++ false) so it does not veto the correct hypothesis.
Several data batches (e.g. one per camera of a rig): mark each RansacRole.sampled (C++ kSampled) with its own threshold. Samples are drawn from the union of all sampled factors, and inlier_mask(i) (C++ InlierMask(i)) gives the mask of batch i.
Known quantities (3D landmarks in PnP, rig extrinsics): put them in state batches with constant states (const_state_ids); they cost nothing towards \(D\).

Gauss-Newton or Levenberg-Marquardt?#

Both variants share every option and step (they are the two RansacMinimizer subclasses; a RansacMinimizer& accepts either). RansacGaussNewtonMinimizer takes full Gauss-Newton steps and is the fastest. RansacLevenbergMarquardt Minimizer damps each hypothesis’s steps with its own \(\lambda\), which widens the convergence basin when the initial guess is far off or the problem is strongly nonlinear, at some extra cost per iteration. Start with GN; switch to LM if hypotheses fail to converge from your initial guesses.

Options reference and tuning #

Defaults are tuned on PnP and work for most small problems.

Option	Default	Meaning and advice
`factor_batches`	empty	Role and inlier threshold per residual batch (same order as `Problem::GetResidualBatches()`). Empty: all sampled with `default_inlier_threshold`.
`default_inlier_threshold`	1.0	Threshold \(\tau\) used when `factor_batches` is empty. Set it from the noise level (see above); the default is only a placeholder.
`hypotheses_per_round`	256	Hypotheses solved in parallel per round. Larger rounds use the GPU better but can overshoot the number actually needed; 256–1024 is a good range.
`max_rounds`	8	Upper bound on rounds. The maximum number of hypotheses is `hypotheses_per_round * max_rounds`; raise it for very high outlier ratios (see the table above).
`confidence`	0.999	Target probability of having drawn at least one all-inlier sample.
`early_stop_inlier_ratio`	1.0	Stop as soon as the best inlier ratio reaches this value (1 disables).
`sample_size`	0	Factors per minimal sample; 0 = \(\lceil D / m_{\min} \rceil\). Increase it if minimal samples are often degenerate (e.g. collinear points), at the price of more hypotheses.
`seed`	0	Sampler seed. Same seed and problem, bitwise identical result.
`scoring`	MSAC	`RansacScoring.msac` or `RansacScoring.inlier_count` (C++ `kMSAC` or `kInlierCount`).
`score_always_on`	true	Add the always-on cost to the score.
`require_informative_inliers`	true	Count a factor as an inlier only if its Jacobian is non-zero on a free state. Disable only for factors that never report zero residuals for invalid configurations; it saves one Jacobian evaluation per scored factor.
`scoring_subset_size`	16384	Two-stage scoring for problems with more than twice this many sampled factors. 0 always scores every hypothesis on every factor.
`scoring_finalists`	4	Hypotheses re-scored on all factors in two-stage scoring (at most 64).
`scoring_memory_budget_bytes`	64 MiB	Device memory for scoring buffers; bounds how many hypotheses are scored at once. Results do not depend on it.
`hypothesis_iterations`	5	GN / LM iterations per hypothesis. More helps a poor initial guess.
`final_iterations`	20	Iterations of the final refinement (stops early on convergence).
`state_tolerance`, `cost_tolerance`	1e-10, 1e-7	Per-hypothesis convergence on the squared step norm and the relative cost decrease.
`linear_solver`	LDLT	`RansacLinearSolverType.ldlt` (C++ `kLDLT`) gives rank-deficient directions a zero step, which suits the near-singular systems minimal samples produce. `cholesky` (C++ `kCholesky`) marks such hypotheses invalid.

RansacLevenbergMarquardtMinimizerOptions adds the damping parameters of the regular LM minimizer (initial_lambda = 1e-3, lambda_upscale = 2, lambda_downscale = 0.5, lambda_min/lambda_max, step_accept_threshold, lambda_downscale_threshold), applied per hypothesis.

Reading the summary. RansacSummary extends MinimizerSummary (initial_cost over all factors, final_cost of the refinement over the inliers and always-on factors, num_iterations and iteration_costs of the refinement) with num_rounds, num_hypotheses, num_valid_hypotheses, num_inliers, inlier_ratio, best_score and refinement_reverted.

Troubleshooting.

Too few inliers / wrong estimate: check the threshold against the actual residual noise (evaluate the residuals at a known good state); raise max_rounds for low inlier ratios; use the LM variant or more hypothesis_iterations if the initial guess is far.
Many invalid hypotheses (num_valid_hypotheses much smaller than num_hypotheses): samples are often degenerate. Increase sample_size by one or two.
Inliers that are clearly wrong: the threshold is too large, or a factor returns zero residuals for invalid configurations without require_informative_inliers.

Limits and requirements #

minimize raises ValueError (C++: Minimize throws std::invalid_argument) with an explanatory message when:

the free tangent dimension \(D\) is 0 or exceeds 64;
no residual batch is sampled, or the sample size exceeds the number of sampled factors;
factor_batches is non-empty but its length differs from the number of residual batches, or a sampled batch has a non-positive threshold;
a factor references more than 8 states, or a state that belongs to no registered state batch;
a residual batch uses numeric Jacobians (JacobianMode::kNumeric), which RANSAC does not support yet;
a residual batch is a constraint (solved by AugmentedLagrangianMinimizer), a state batch has box bounds, or the problem has a subproblem partition or state stages: RANSAC does not implement these and rejects them rather than ignoring them;
Problem.check_consistency() (C++ Problem::CheckConsistency()) fails, or an option is out of range (hypotheses_per_round, max_rounds or hypothesis_iterations of 0, confidence outside (0, 1)).

Every factor and state batch must support the item / replica parameters of Evaluate and Plus. All built-in batches do; for your own types see Custom Factors and States (Python and C++).

Performance #

PnP, 30% outliers, one round of 256 hypotheses, median wall time (NVIDIA RTX A6000):

correspondences	LM + Cauchy loss	RANSAC-GN	RANSAC-LM
1,000	0.7 ms	1.4 ms	1.4 ms
100,000	3.7 ms	7.3 ms	7.5 ms
1,000,000	31.6 ms	14.9 ms	47.3 ms

In the same benchmarks RANSAC succeeds in 100% of trials from 0% to 90% uniformly random outliers, from both small and large initial errors, and against coherent outliers (a second, competing pose); plain GN/LM fail from 10% outliers and LM with a Huber loss from 60%.