Robustifier API

The robustifier module defines GPU-batched robust loss functions used by the minimizer to reduce the influence of outliers in non-linear least squares.

C++ — cunls/robustifier

Python — pycunls

What robustifier functions are for

In least squares, a few bad measurements (outliers) can pull the solution away from the true optimum. A robust loss function \(\rho(s)\) replaces the squared residual term \(s = \|f_i\|^2\) with \(\rho(s)\), so that large residuals contribute less to the cost. The solver then minimizes \(\frac{1}{2}\sum_i \rho(\|f_i\|^2)\) instead of \(\frac{1}{2}\sum_i \|f_i\|^2\). Typical choices of \(\rho\) behave like \(s\) for small \(s\) (inliers) and grow sublinearly for large \(s\) (outliers), so inliers are fitted normally and outliers are down-weighted.

Inputs and outputs

Each robustifier evaluates at squared residual values \(s = \|f\|^2\) (one \(s\) per residual vector). For each \(s\) it returns a float3 \((\rho(s), \rho'(s), \rho''(s))\):

• First component \(\rho(s)\): robustified cost value; the term contributed to the total cost is \(\frac{1}{2}\rho(s)\).

• Second component \(\rho'(s)\): first derivative of \(\rho\) with respect to \(s\); used to scale residuals and Jacobians in the robustified Gauss-Newton step.

• Third component \(\rho''(s)\): second derivative of \(\rho\) with respect to \(s\); used in the Triggs correction when forming the robustified normal equations.

Calling the evaluator with negative \(s\) is invalid; implementations need not handle that case. Common choices of \(\rho\) satisfy \(\rho(0)=0\), \(\rho'(0)=1\), and in the outlier region \(\rho'(s) < 1\) and \(\rho''(s) < 0\).

LossFunctionBatch

Abstract base (cunls/robustifier/loss_function_batch.h).

bool Evaluate(

float *s,

float3 *out,

int num_losses,

cudaStream_t stream

) const

Evaluates the loss for a batch of squared residuals.

Parameters:

• s – [in] Device pointer to squared residual values \(s = \|f\|^2\).

• out – [out] Device pointer to float3 tuples \((\rho(s), \rho'(s), \rho''(s))\) for each input.

• num_losses – [in] Number of residual values to process.

• stream – [in] CUDA stream for asynchronous execution.

Returns:

true on success.

TrivialLossFunctionBatch

Header: cunls/robustifier/trivial_loss_function_batch.h

TrivialLossFunctionBatch()

Returns:

Constructor has no return value.

Formula (unscaled)

\[\rho(s) = s,\qquad \rho'(s) = 1,\qquad \rho''(s) = 0.\]

Identity loss: no robustification; equivalent to standard least squares.

HuberLossFunctionBatch

Header: cunls/robustifier/huber_loss_function_batch.h

HuberLossFunctionBatch(float delta)

Parameters:

delta – [in] Inlier/outlier threshold (scale); quadratic for \(s \le \delta^2\), linear for \(s > \delta^2\).

Returns:

Constructor has no return value.

Formula (scaled with :math:`delta`)

\[\begin{split}\rho(s) = \begin{cases} s & s \le \delta^2 \\ 2\delta\sqrt{s} - \delta^2 & s > \delta^2 \end{cases}\end{split}\]

\[\begin{split}\rho'(s) = \begin{cases} 1 & s \le \delta^2 \\ \delta/\sqrt{s} & s > \delta^2 \end{cases} ,\qquad \rho''(s) = \begin{cases} 0 & s \le \delta^2 \\ -\rho'(s)/(2s) & s > \delta^2 \end{cases}.\end{split}\]

CauchyLossFunctionBatch

Header: cunls/robustifier/cauchy_loss_function_batch.h

CauchyLossFunctionBatch(float b, float c)

Parameters:

• b – [in] Output scale parameter.

• c – [in] Shape parameter (larger \(c\) makes the loss grow more slowly).

Returns:

Constructor has no return value.

Formula

\[\rho(s) = b\,\ln(1 + c\,s),\qquad \rho'(s) = \frac{b\,c}{1 + c\,s},\qquad \rho''(s) = -\frac{c^2 b}{(1+c\,s)^2}.\]

Unscaled case: \(\rho(s) = \ln(1+s)\) (e.g. \(b=1\), \(c=1\)).

ArctanLossFunctionBatch

Header: cunls/robustifier/arctan_loss_function_batch.h

ArctanLossFunctionBatch(float a, float b)

Parameters:

• a – [in] Scale parameter (argument scale in \(\arctan(s/a)\)).

• b – [in] Shape parameter, typically \(1/a^2\) for derivative scaling.

Returns:

Constructor has no return value.

Formula

With \(s\) the squared residual, the implementation uses \(\rho(s) = a\,\arctan(s/a)\) and \(\rho'(s) = 1/(1 + s^2 b)\) with \(b = 1/a^2\):

\[\rho(s) = a\,\arctan\frac{s}{a},\qquad \rho'(s) = \frac{1}{1 + (s/a)^2},\qquad \rho''(s) = -\frac{2s/a^2}{(1+(s/a)^2)^2}.\]

Unscaled case: \(\rho(s) = \arctan(s)\) (e.g. \(a=1\), \(b=1\)).

SoftLOneLossFunctionBatch

Header: cunls/robustifier/soft_lone_loss_function_batch.h

SoftLOneLossFunctionBatch(float b, float c)

Parameters:

• b – [in] Scale parameter.

• c – [in] Shape parameter (larger \(c\) makes the loss grow more slowly).

Returns:

Constructor has no return value.

Formula

\[\rho(s) = 2b\left(\sqrt{1 + c\,s} - 1\right),\qquad \rho'(s) = \frac{b\,c}{\sqrt{1+c\,s}},\qquad \rho''(s) = -\frac{c^2 b}{2(1+c\,s)^{3/2}}.\]

Unscaled case: \(\rho(s) = 2(\sqrt{1+s}-1)\) (e.g. \(b=1\), \(c=1\)).

TolerantLossFunctionBatch

Header: cunls/robustifier/tolerant_loss_function_batch.h

TolerantLossFunctionBatch(float a, float b)

Parameters:

• a – [in] Offset parameter (soft threshold).

• b – [in] Scale parameter (smoothing).

Returns:

Constructor has no return value.

Formula

With \(c = b\,\ln(1 + e^{-a/b})\) so that \(\rho(0)=0\):

\[\rho(s) = b\,\ln\left(1 + e^{(s-a)/b}\right) - c,\qquad \rho'(s) = \frac{e^{(s-a)/b}}{1 + e^{(s-a)/b}},\qquad \rho''(s) = \frac{1}{4b\,\cosh^2\bigl((s-a)/(2b)\bigr)}.\]

TukeyLossFunctionBatch

Header: cunls/robustifier/tukey_loss_function_batch.h

TukeyLossFunctionBatch(float a)

Parameters:

a – [in] Cutoff threshold; residuals with \(s > a^2\) get zero weight.

Returns:

Constructor has no return value.

Formula

With \(s\) the squared residual and \(a^2\) the squared cutoff:

\[\begin{split}\rho(s) = \begin{cases} \displaystyle\frac{a^2}{3}\left(1 - \left(1 - \frac{s}{a^2}\right)^3\right) & s \le a^2 \\[0.5em] \displaystyle\frac{a^2}{3} & s > a^2 \end{cases}\end{split}\]

\[\begin{split}\rho'(s) = \begin{cases} \displaystyle\left(1 - \frac{s}{a^2}\right)^2 & s \le a^2 \\ 0 & s > a^2 \end{cases} ,\qquad \rho''(s) = \begin{cases} \displaystyle -\frac{2}{a^2}\left(1 - \frac{s}{a^2}\right) & s \le a^2 \\ 0 & s > a^2 \end{cases}.\end{split}\]

ScaledLossFunctionBatch

Header: cunls/robustifier/scaled_loss_function_batch.h

template<class T>
ScaledLossFunctionBatch(

float a,

Args&&... loss_args

)

Parameters:

• a – [in] Positive scale factor applied to all loss outputs.

• loss_args – [in] Arguments forwarded to the wrapped loss function constructor (e.g. delta when T is HuberLossFunctionBatch).

Returns:

Constructor has no return value.

Throws:

std::invalid_argument – if a <= 0.

T must derive from LossFunctionBatch. The inner loss is owned by value and constructed from the forwarded arguments, following the same decorator pattern as InformationFactorBatch<T>.

Formula

Given an inner loss \(f(s)\) and a positive scalar \(a\):

\[\rho(s) = a\,f(s),\qquad \rho'(s) = a\,f'(s),\qquad \rho''(s) = a\,f''(s).\]

C++ example

// Scale Huber loss by 0.5 — the delta=1.0 argument is forwarded
// to the HuberLossFunctionBatch constructor.
cunls::ScaledLossFunctionBatch<cunls::HuberLossFunctionBatch> loss(0.5f, 1.0f);

Theory — How robustifier outputs are used in optimization

The non-linear least squares problem with robustification is

\[\min_{\mathbf{x}} \quad \frac{1}{2}\sum_i \rho_i\bigl(\|f_i(x_{i_1},\ldots,x_{i_k})\|^2\bigr),\]

where \(f_i\) are residual vectors and \(\rho_i\) are loss functions. Let \(s = \|f\|^2\) and write \(\rho\), \(\rho'\), \(\rho''\) for the loss and its derivatives at \(s\). The contribution of one term to the total cost is \(\frac{1}{2}\rho(s)\); the robustifier API returns \((\rho(s), \rho'(s), \rho''(s))\) so the solver can form the robustified gradient and Gauss-Newton system without recomputing \(\rho\).

Robustified gradient

For a single residual block, the gradient of \(\frac{1}{2}\rho(\|f(x)\|^2)\) with respect to the parameters \(x\) is

\[g = \rho'\, J^\top f,\]

where \(J\) is the Jacobian of \(f\) with respect to \(x\). So \(\rho'\) (the second component of the robustifier output) scales the gradient and thus down-weights the contribution of large residuals.

Robustified Gauss-Newton Hessian

The Gauss-Newton approximation to the Hessian of \(\frac{1}{2}\rho(\|f\|^2)\) involves both \(\rho'\) and \(\rho''\). With \(r = f(x)\) and \(s = \|f\|^2 = r^\top r\), the Hessian contribution (ignoring second derivatives of \(f\)) is

\[H = J^\top \left( \rho'\, I + 2\rho''\, r r^\top \right) J.\]

When \(\rho'' < 0\) (typical for robust losses in the outlier region), \(H\) can be indefinite. To keep a positive-definite approximation and still use a Jacobian-based solver, the implementation rescales the residual and Jacobian so that the robustified problem looks like a standard least squares problem in the rescaled variables.

Rescaling (Triggs correction)

Let \(\alpha\) be a root of

\[\frac{1}{2}\alpha^2 - \alpha - \frac{\rho''}{\rho'}\|f\|^2 = 0.\]

Then the rescaled residual and Jacobian

\[\tilde{f} = \frac{\sqrt{\rho'}}{1-\alpha}\, f,\qquad \tilde{J} = \sqrt{\rho'}\,\left(I - \alpha\, \frac{f f^\top}{\|f\|^2}\right) J\]

yield a Gauss-Newton step equivalent to the robustified problem. When \(2\rho''\|f\|^2 + \rho' \lesssim 0\), \(\alpha\) is capped (e.g. \(\alpha \le 1-\epsilon\)) to avoid numerical issues. The robustifier output \((\rho(s), \rho'(s), \rho''(s))\) is used to compute \(\alpha\) and the scaling factors \(\sqrt{\rho'}\) and \((1-\alpha)^{-1}\) applied to residuals and Jacobians in the solver. For more detail, see the Ceres Solver documentation on LossFunction and the references therein (e.g. Triggs).

Python API (pycunls)

All Python loss function batches share the same base class pycunls.LossFunctionBatch. The formulas and parameters are identical to the C++ versions above. Pass a loss function instance to Problem.add_factor_batch to apply robustification:

loss = pycunls.HuberLossFunctionBatch(delta=1.0)
problem.add_factor_batch(factor_batch, loss, state_pointers)

Python class	Constructor
TrivialLossFunctionBatch	TrivialLossFunctionBatch()
HuberLossFunctionBatch	HuberLossFunctionBatch(delta: float)
CauchyLossFunctionBatch	CauchyLossFunctionBatch(b: float, c: float)
ArctanLossFunctionBatch	ArctanLossFunctionBatch(a: float, b: float)
SoftLOneLossFunctionBatch	SoftLOneLossFunctionBatch(b: float, c: float)
TolerantLossFunctionBatch	TolerantLossFunctionBatch(a: float, b: float)
TukeyLossFunctionBatch	TukeyLossFunctionBatch(a: float)
ScaledLossFunctionBatch	ScaledLossFunctionBatch(loss_function: LossFunctionBatch, a: float)