\section{Introduction}
\label{sec:intro}

A direct optimization on the test-time evaluation metric has shown to be beneficial in deep learning for many vision tasks \cite{arulkumaran2017brief,lee2018context,balle2018variational,berrada2018smooth,lapin2016loss}. 
Training the model directly on the evaluation metric is infeasible when the metric is non-differentiable. 
In such cases, training with a surrogate of the metric has shown to be effective. Examples include surrogate losses of average precision and recall@k for image retrieval \cite{brown2020smooth,patel2022recall,cakir2019deep,ramzi2021robust}, perceptual loss for image compression \cite{patel2021saliency,patel2019deep}, intersection-over-union loss for object detection \cite{rezatofighi2019generalized,yu2016unitbox,patel2020learning}, edit distance loss for scene text recognition \cite{patel2020learning,patel2021feds}, etc. 
However, a differentiable version of RANSAC~\cite{RANSAC} still does not exist, prohibiting a direct optimization on the evaluation metrics for several tasks that are utmost important in computer vision and robots.

Robust estimation is a fundamental component in vision pipelines, including relative pose estimation~\cite{chum2004epipolar}, wide baseline matching~\cite{pritchett1998wide,matas2004robust,mishkin2015mods}, multi-model fitting~\cite{isack2012energy,pham2014interacting}, image-based localization~\cite{Brachmann_2017_CVPR}, motion segmentation~\cite{torr1993outlier}, and pose graph initialization of Structure-from-Motion (SfM) algorithms~\cite{schoenberger2016sfm, schoenberger2016mvs}.
While several robust estimators have been proposed throughout the years~\cite{holland1977robust,illingworth1988survey,li2009consensus,xiao2016hypergraph}, randomized hypothesize-and-verify approaches, like RANSAC~\cite{RANSAC} and its recent variants~\cite{barath2018graph,barath2019magsac,barath2020magsac++,ivashechkin2021vsac}, have become the most widely used methods due to their robustness, simplicity, and efficiency. 
RANSAC repeatedly selects random minimal subsets of the input data and fits a model, \eg, a 3D plane to three points or fundamental matrix to seven point correspondences.
The model quality is then computed as the cardinality of its support, \ie, the number of points that fall closer than a manually set inlier threshold. 
The so-far-the-best model is updated if a new model is found with higher quality. 
RANSAC terminates when the probability of finding an all-inlier sample falls below a threshold, controlled by a confidence parameter typically set to 0.95 or 0.99. 
Finally, the model with the highest quality, polished, \eg, by least-squares fitting on all inliers, is returned.

\begin{figure}
    \centering
   
    \includegraphics[width=1\columnwidth]{Figures/fig1new.pdf}
   
   
   
   
    
   
   
   
   
   
    \caption{$\nabla$-RANSAC gets a set of tentative correspondences as input, predicts their inlier probabilities, runs a differentiable RANSAC procedure and, finally, estimates the model parameters, \eg, fundamental matrix (shown with its inliers here).
   
    }
    \label{fig:example}
\end{figure}

\begin{figure*}[t]
	\centering
    \includegraphics[width = 2.1\columnwidth, trim=0mm 0mm 0mm 0mm, clip]{Figures/pipeline.pdf}
    \caption{\textbf{Pipeline of $\nabla$-RANSAC.}
    Given an image pair, a feature matcher and detector is applied to obtain a set of tentative correspondences. 
    The point coordinates and other side information, like SNN ratio, feature orientation and scale, are fed into the consensus learning network from \cite{clnet2021} to predict inlier probabilities.
    The predictions go into the fully differentiable RANSAC module, $\nabla$-RANSAC.
    In each iteration, the differentiable and randomized Gumbel sampler (Section~\ref{sec:sampler}) selects a minimal sample $\mathcal{S}$ of $m$ correspondences. 
    Model $\hat{\theta}$ is estimated by a differentiable minimal solver (Section~\ref{sec:diff_solver}) and its score is calculated (Section~\ref{sec:quality}) from the set of tentative correspondences.}
   \label{fig:flowchart}
\end{figure*}

Since the publication of RANSAC, there have been a number of important modifications proposed to improve specific components of the original algorithm.
Improvements to the classical framework include local optimization (LO-RANSAC, LO$^+$-RANSAC, GC-RANSAC)~\cite{chum2003locally,lebeda2012fixing,barath2018graph}, 
advanced scoring functions (MLESAC, MSAC, MAGSAC, and MAGSAC++)~\cite{torr2000mlesac,barath2019magsac,barath2020magsac++,barath2022learning}, 
speed-ups using probabilistic sampling of hypotheses (PROSAC, NAPSAC, and P-NAPSAC)~\cite{torr2002napsac,chum2005matching,barath2020magsac++} or by preemptive verification strategies (SPRT, SP-RANSAC)~\cite{chum2008optimal,barath2021space}, 
degeneracy checks (DEGENSAC, QDEGSAC, NeFSAC)~\cite{chum2005two,frahm2006ransac,cavalli2022nefsac}, and strategies to auto-tune the inlier-outlier threshold (MINPRAN, a contrario RANSAC)~\cite{stewart1995minpran,moisan2012automatic}.
According to the recent CVPR 2022 Image Matching Challenge~\cite{imc2022}, the most accurate method for relative pose estimation currently is MAGSAC++~\cite{barath2020magsac++} with PROSAC sampling~\cite{chum2005matching} and DEGENSAC-based degeneracy~\cite{chum2005two} testing. 
Regardless of their good performance, none of the RANSAC variants is differentiable, keeping the randomized robust estimation one of the only components in geometric image matching pipelines that can not be trained end-to-end.





In the recent years, several algorithms were proposed approaching the robust relative pose estimation problem using neural networks. 
Context normalization networks (CNe)~\cite{cne2018} is the first paper on the topic which proposes to use PointNet (MLP)~\cite{qi2017pointnet}, \ie a multi-layer perceptron (MLP), with batch normalization~\cite{batchnorm2015} as context mechanism. 
Attentive context normalization networks~\cite{acne2020} introduces a special architectural block for the task. 
Deep Fundamental matrix estimation~\cite{dfe2018} uses differentiable iteratively re-weighted least-squares to predict inliers weights for the input correspondences. 
The OANet algorithm~\cite{oanet2019} introduces several architectural blocks for the correspondence filtering.
Neural Guided RANSAC~\cite{brachmann2019neural} uses a CNe-like architecture with a different training objective.
It assumes the inlier probabilities to follow the categorical distribution and approximates the gradients.
The predicted correspondence scores are exploited inside RANSAC by using a guided sampling algorithm that helps in finding accurate models early. 
CLNet~\cite{clnet2021} introduces several algorithmic and architectural improvements  to remove gross outliers with iterative pruning. 
It propagates the gradients between the input data and estimated models by using the weighted eight-point algorithm~\cite{hartley1997defense}.
Exploiting the differentiable pooling layers in OANet~\cite{oanet2019}, MS2DG-Net~\cite{Dai_2022_CVPR} uses additional semantic features based on sparse similarities while maintaining its permutation equivariant property. 
DSAC~\cite{Brachmann_2017_CVPR} is the first method to propose a differentiable RANSAC-like robust estimator, however, only for a highly specific task, \ie, one-shot camera re-localization from an RGB image.
Instead of finding directly the best model, DSAC employs soft probabilistic hypothesis selection. 

Even though most of these algorithms are trainable in an end-to-end manner, \eg via employing iteratively re-weighted least-squares fitting with the predicted weights, they can be considered elaborate inlier pre-filtering approaches -- a final RANSAC on the top-\textit{k} correspondences \textit{always} improves the results compared to LS fitting~\cite{cvpr2020ransactutorial}.

In this paper, we focus on making RANSAC end-to-end differentiable and propose $\nabla$-RANSAC (see Fig.~\ref{fig:example}). 
The main contribution is as follows:
we investigate all components of the RANSAC pipeline and propose differentiable alternatives that allow learning inlier probabilities in an end-to-end manner.
This includes a new random sampling approach based on a clever re-parametrization strategy, \ie Gumbel Softmax sampler, that allows for propagating the gradients through the entire random procedure.
Moreover, we implement differentiable versions of the widely used minimal solvers, \eg five and seven-point algorithms~\cite{nister2004efficient}.
Finally, we propose a model quality function, avoiding the non-differentiable best-model-selection of RANSAC, via marginalizing over all models estimated within the procedure.
$\nabla$-RANSAC allows robust estimators to use test-time evaluation metrics to optimize end-to-end training.

As a proof of its potential to unlock the end-to-end training of geometric pipelines, we train $\nabla$-RANSAC together with a recent detector-free feature matcher, LoFTR \cite{Sun_2021_CVPR}, with which we achieve improved accuracy. 
We keep it as future work to find the best way, \eg losses and network architecture, that maximizes performance when the entire pipeline is trained end-to-end. 

\section{Fully Differentiable RANSAC}
\label{sec:differentiable_RANSAC}

In this section, we discuss the algorithmic components of the proposed $\nabla$-RANSAC making it capable of end-to-end training.
A pipeline employing $\nabla$-RANSAC is in Fig~\ref{fig:flowchart}.

The input is a set of tentative point correspondences, possibly equipped with extra information from the detector and matcher, \eg, feature orientation, scale, affine shape, and matching confidence. 
An initial step of the algorithm is based on consensus learning via exploiting the pruning block from the recent \cite{clnet2021}.
Instead of using the weights to iteratively remove outliers as in \cite{clnet2021}, we consider them initial inlier probabilities that are fed into the subsequent fully differentiable RANSAC module, $\nabla$-RANSAC. 
Similarly to plain RANSAC, the proposed $\nabla$-RANSAC is an iterative random sampling of $m$ data points, model estimation via a minimal solver, and model scoring.
The devil is, however, in the details. 
The sampler is the \textit{Gumbel Softmax Sampler} that exploits the input probabilities as guidance in a randomized manner while propagating the gradients. 
The differentiable minimal solver estimates the model parameters from the drawn sample. 
The model quality is computed in a supervised way using the ground truth, \eg, it is the relative pose error calculated from the estimated essential matrix and the known pose from a reference reconstruction.\footnote{Note that the algorithm should work with self-supervised losses, like inlier ratio, too.
Exploring such approaches we keep as future work.}
While, in RANSAC, the best-model-selection is a non-differentiable step,
we aggregate the model qualities calculated in the $k$ iterations performed throughout the entire procedure.
The minimized loss comes from the average of these qualities, guiding the network to predict inlier probabilities that \textit{on average} lead to accurate models.

\subsection{Gumbel Softmax Sampler}
\label{sec:sampler}


The first step of $\nabla$-RANSAC is hypothesis generation which requires sampling $m$ data points from a set of $n$ total samples $\{x_{i}\}$. 
The sampling distribution is either governed by the importance scores for each sample $\{s_{i}\}$ or follows a uniform distribution. 
The standard sampling operation is either non-differentiable or has sparse gradients that are not useful for training \cite{doersch2016tutorial}. 
Further, even the deterministic selection of top-$k$ samples with highest importance scores is non-differentiable as ranking and counting operations require a Heaviside step function \cite{patel2022recall}. Therefore, to make randomized importance score sampling differentiable for RANSAC, we extend the Gumbel-Softmax \cite{jang2016categorical,maddison2016concrete} with the straight through trick \cite{van2017neural} to multiple samples.

The sampling is performed on a distribution that is governed by the importance scores for each samples. Let this distribution be $\mathbf{A}$ as follows:
\begin{equation}
\mathbf{A}_{i} = \frac{\exp(s_{i} + \gamma_{i})}{\sum_{j=1}^{n} \exp(s_{j} + \gamma_{j})},
\end{equation}
where ${\gamma_{i}}$ are i.i.d random samples drawn from a Gumbel distribution $\text{G}(0,1)$. 
Parameters $\{s_{i}\}$ are the importance scores for each sample in the set. 
Note that if the scores are not available, $\{s_{i}\}$ are set to a constant value to perform uniform random sampling. 
\begin{figure}[t]
	\centering
    \includegraphics[width = 1\columnwidth, trim=0mm 0mm 0mm 0mm, clip]{Figures/sampler.pdf}
    \caption{Overview of \textbf{Gumbel Softmax Sampler} 
 used in $\nabla$-RANSAC. The input to the sampler is the importance scores $s_1, s_2, ..., s_n$. The process starts by i.i.d random sampling $\gamma_1, \gamma_2, ..., \gamma_n$ from a Gumbel(0,1) distribution followed by a Softmax operation. 
 Based on the distribution obtained from the Softmax, $m$ data points from a total of $n$ are drawn by using the top-k operation. 
 Since top-k is non-differentiable, the straight through trick~\cite{van2017neural} is used and the Softmax output is added and subtracted with stop-gradient (sg) to allow backward flow of the gradients.
 The figure is best viewed in color. 
 The green arrows show the connections through which gradient flow is possible. 
 Red arrows show connections without any gradients.}
    \label{fig:gumbel_sampler}
\end{figure}

The feed-forward through the sampler is performed such that the output of the selected indices are binary and gradients can backpropagate for training. Such feed-forward is achieved via the straight through trick \cite{van2017neural} as follows:
\begin{equation}
\mathbf{\hat{A}} = \text{top-k}(\textbf{A},m) + \textbf{A} - \text{sg}(\textbf{A}),
\end{equation}
where sg is the stop-gradient operation and top-k selects the indices of the $m$ highest values. During the forward pass, $\hat{\textbf{A}}$ contains only ones and zeros where one represents the samples that are selected. During back-propagation, the gradients only flow through $\textbf{A}$ as $\text{top-k}(\mathbf{A},m)$ and $\text{sg}(\mathbf{A})$ do not have any gradients. An overview of the proposed differentiable sampler is shown in \cref{fig:gumbel_sampler}. 

\subsection{Differentiable Minimal Solver}
\label{sec:diff_solver}

Minimal solvers~\cite{kukelova2008automatic} are a fundamental part of RANSAC-like hypothesize-and-verify approaches, allowing to progressively explore the input data in a randomized manner. 
Such methods estimate the sought model parameters from a minimal set $\mathcal{S} \subseteq \mathcal{P}$ of $m \in \mathbb{N}$ data points, where $\mathcal{P}$ is the set all points. 
For example, an essential matrix can be estimated from five point correspondences or a 3D plane from a minimum of three 3D points. 
This allows keeping the complexity of the robust estimation problem low via sampling from the possible $\binom{|\mathcal{P}|}{m}$ minimal configurations. 

Most of the minimal solvers commonly used in computer vision are actually differentiable after taking each algorithmic component into careful considerations. 
While we will consider two utmost important problems in this paper, \ie fundamental and essential matrix estimation, the previous statement holds for all geometric model estimation problems that we are aware of in vision. 

Recent learning-based approaches\cite{clnet2021} propagate the gradients through the well-known normalized eight-point (8PC) algorithm~\cite{hartley1997defense} when estimating both essential and fundamental matrices.
There are two main issues with the 8PC algorithm in practice.
First, and most importantly, the 8PC and 7PC solvers have a degeneracy when the points stem from a close-to-planar underlying 3D structure.
In this case, a degenerate model is estimated that, while often having a large number of inliers~\cite{chum2005two}, is incorrect w.r.t.\ the actual scene geometry.  
This misguides the learning in scenes dominated by planar structures and deteriorates the performance. 
Second, using eight correspondences instead of the minimal five (E) or seven (F) largely increases the complexity of the problem. 
This prevents the network from learning the maximal information available in a given correspondence set, and makes the learning sub-optimal. 

For essential matrix estimation, there are a number of 5PC solvers~\cite{nister2004efficient,li2006five,batra2007alternative,stewenius2008minimal,bujnak2012making} proposed throughout the years.
Most frequently, in practical applications (\eg, SfM~\cite{schonberger2016structure,theia-manual}), one of the following methods are used:
the method of Stewenius et al.~\cite{stewenius2008minimal} due to its stability, and that of David Nister~\cite{nister2004efficient} due to its effectiveness. 
The steps of~\cite{stewenius2008minimal} are as follows:
creating the coefficient matrix from the input points, 
decomposing it by SVD, solving a linear system of equations by, \eg, the Moore-Penrose pseudo-inverse, and finally, recovering the linear combination of the null vectors via the Eigen decomposition of the action matrix. 
Even though both SVD~\cite{wang2021robust} and Eigen decompositions~\cite{song2022batch} have differentiable solutions, we found the Eigen decomposition to be unstable for gradient propagation. 

The steps of the solver of Nister~\cite{nister2004efficient} are as follows: creating the coefficient matrix from the input correspondences, 
decomposing it by SVD, solving a linear system of equations, solving a set of polynomial equations, and finally, basic arithmetic operations. 
We implemented a differentiable polynomial solver based on Sturm sequences~\cite{gonzalez1989sturm} and, also, one using companion matrices.
While both algorithms work well, the companion matrix-based solution is the fastest. 
We apply this solver in the proposed $\nabla$-RANSAC.

Fundamental matrix estimation is an easier problem, where the 7PC~\cite{hartley2003multiple} algorithm, besides running an SVD decomposition on the coefficient matrix, only solves a single cubic polynomial.
As we are given a closed-form solution for the cubic problem, we can straightforwardly make the entire algorithm differentiable. 
Note that, in our experiments, we have not observed improvements by replaced 8PC with 7PC.
This might be due to the shared degeneracy.

Most of the minimal solvers return multiple solutions.
Even though they are consistent with the constraints, most of them are inconsistent with the actual underlying scene geometry.
At inference, from the multiple solutions the best one is selected based on the evaluation metric. To mimic the test-time evaluation in training, for each iteration the best solution is selected and then used for the loss computation. Note that the gradients for the selected solution will be used for training, those for not selected solutions will be zero.

We are not aware of public implementations of such solvers, neither in open-source libraries~\cite{van2014scikit,riba2020kornia} nor in standalone public repositories.
Therefore, we consider this as a technical contribution to the community and will make the source codes publicly available.

\subsection{Trainable Quality Function}
\label{sec:quality}
In RANSAC, the quality of an estimated model is calculated as the cardinality of its support, \ie, the inlier number.
Since RANSAC, a number of algorithms~\cite{torr2000mlesac,barath2019magsac,barath2020magsac++,barath2022learning} have improved the performance by better modelling the noise both in the inliers and outliers.
However, all such methods perform a best model selection step based on the maximal quality which renders the procedure non-differentiable. 
Some works~\cite{Brachmann_2017_CVPR} tackle with this problem by employing soft probabilistic hypothesis selection. 
Other methods~\cite{cne2018, clnet2021} combine classification loss with regression and geometry-induced losses~\cite{brachmann2019neural} to reason about the quality of the estimated model and inlier assignments.

Instead of the above solutions, we exploit all models $\{ \hat{\theta}_i \}_{i = 1}^t$ estimated during the fixed $t \in \mathbb{N}^{>0}$ iterations.
In each iteration, the estimated model is compared to the ground truth and its implied loss is calculated, \eg, as the relative pose error in case of essential matrix estimation. 
The minimized loss, in case of relative pose estimation, is $L_\text{pose}$ calculated as the average over all models as follows:
\begin{equation}
    L_\text{pose} = \frac{1}{2t} \sum_{i = 1}^t \epsilon_\mathbf{R}(\hat{\theta}_i^\mathbf{R}, \theta_i^\mathbf{R}) + \epsilon_\mathbf{t}(\hat{\theta}_i^\mathbf{t}, \theta_i^\mathbf{t}),
    \label{eqn:loss_pose}
\end{equation}
where the upper index refers to the rotation ($\mathbf{R} \in \text{SO}(3)$) and translation components  ($\mathbf{t} \in \mathbb{R}^3$) of the pose, $\theta$ is the ground truth, and functions $\epsilon_\mathbf{R}$ and $\epsilon_\mathbf{t}$ return the rotation and translation errors, respectively. 
Functions 
$\epsilon_\textbf{R}(\hat{\theta}_i^\mathbf{R}, \theta^\mathbf{R}) = \cos^{-1} ( ( \text{tr} ( \widehat{\mathbf{R}} \mathbf{R}^\text{T} ) - 1 ) / 2)$ 
and 
$\epsilon_\mathbf{t}(\hat{\theta}_i^\mathbf{t}, \theta^\mathbf{t}) = \cos^{-1} ( \mathbf{\hat{t}}^\text{T} \mathbf{t} / (|\mathbf{\hat{t}}| |\mathbf{t}|))$.

Besides $L_\text{pose}$ we employ other widely used objectives too.
Given the probabilities $\mathcal{\hat P}_i = \{ \hat p_j^i \}_{j = 1}^n$ that model $\theta_i$ implies in the $i$th iteration ($\hat p_j^i \in [0, 1]$), we calculate classification loss $L_\text{inlier}$ w.r.t.\ the inlier probabilities $\mathcal{P} = \{ p_j \}_{j}^n$ determined by the ground truth model.
Probabilities $\{ p_j \}_{j = 1}^n$ can be calculated from the point-to-model residuals making certain assumptions about the noise.
In our algorithm, we follow the RANSAC assumption and $p_i \in \{0, 1\}$ is a binary mask. 
The implied classification loss is as follows:
\begin{equation}
    L_\text{inlier} = \frac{1}{t} \sum_{i = 1}^t {\text{BCE}}( \mathcal{\hat P}_i, \mathcal{P}). 
    \label{eqn:loss_inlier}
\end{equation}
Also, we compute a loss $L_\text{epi}$ summing the average symmetric epipolar errors as follows:
\begin{equation}
    L_\text{epi} = \frac{1}{t |I|} \sum_{i = 1}^t \sum_{p \in \mathcal{I}} \epsilon_\text{epi} (\theta_i, p),
    \label{eqn:loss_epi}
\end{equation}
where $\mathcal{I}$ is the inlier set selected by the ground truth model.


The overall loss minimized in our training procedure is the linear combination of the above ones as follows:
\begin{equation}
    L_\text{overall}  = w_\alpha L_\text{inlier} + w_\beta L_\text{pose} + w_\gamma L_\text{epi},
    \label{eqn:overall}
\end{equation}
where $w_\alpha, w_\beta, w_\gamma \in \mathbb{R}$ are the weighting parameters.
Note that these test-time evaluation metrics cannot be used in the case the input parameter is non-differentiable, \ie, robust estimation results from standard RANSAC.
Thus, it is meaningful to utilize test-time evaluation metrics in the fully differentiable RANSAC we propose.





\subsection{Training and Testing Details}

The input of $\nabla$-RANSAC is a set of correspondences obtained by any feature detector and matcher. 
The number of  matches is fixed to 2000.
We choose the best 2000 ones, based on the matching score, if we are given more.
In case of having fewer correspondences, we fill up the missing values with zeros.
We extract local and global features from the correspondences by a consensus learning block as backbone~\cite{clnet2021}.
We integrate over the extracted geometric information, \eg, scale, and orientation of the local descriptors to help learn the qualities of correspondences.
\begin{table*}[t]
    \centering
  \resizebox{2.1\columnwidth}{!}{\begin{tabular}{lcccccccccc}
      \toprule
  Dataset / Method  &\tabincell{c}{LMEDS \\\cite{rousseeuw2005robust}} & \tabincell{c}{RSC \\\cite{RANSAC}} 
  &  \tabincell{c}{ GC-RSC\\ \cite{barath2018graph}}& \tabincell{c}{MSC\\  \cite{barath2019magsac}}& \tabincell{c}{MSC++ \\
  \cite{barath2020magsac++}} & \tabincell{c}{EAS \\ \cite{fan2021efficient}} & \tabincell{c}{OANet\\ \cite{oanet2019}} &\tabincell{c}{CLNet \\ \cite{clnet2021}} &\tabincell{c}{NG-RSC \\\cite{brachmann2019neural}} & \tabincell{c}{$\nabla$-RANSAC} \\ \midrule
  Avg. time (ms) & $\textbf{21.00}$  & $\underline{30.67}$  & $73.25$  & $281.3$ &$ 318.13$ & $325.83$ &$\textbf{21.00}$ &$34.83$& $25.85$  & $35.08$  \\
  \midrule
  
  Buckingham Palace & $23.26$ & $24.95$ & $26.48$& $27.23$ & $26.28$  & $\underline{28.28} $& $24.70$ & $27.46$ &$28.06$&$\textbf{30.94} $ \\
  
  Brandenburg Gate & $31.74$ & $39.70$ & $43.01$& $42.02$ & $42.43$ &$34.84$ & $42.49$ & $39.69$ & $\underline{43.19}$& $\textbf{46.34}$\\ 
  
  Colosseum Exterior & $43.32$ & $48.25$ & $50.86$& $51.76$ & $51.56$ &$50.49$ & $40.50$& $47.12$ & $\underline{52.42}$ & $\textbf{54.99}$\\ 
  
  Grand Place Brussels & $27.45$ & $31.42$ & $\underline{33.31}$& $32.60$& $33.08 $& $32.40$&$29.10$ & $31.80$& $32.13$& $\textbf{35.49}$\\
  
  Notre Dame Front Facade & $30.39$ & $37.20$ & $40.48$& $39.35$ & $39.00$ &$39.98$& $33.17$ & $38.20$ & $\underline{40.59}$  & $\textbf{44.26}$\\
  
  Palace of Westminster & $22.20$ & $28.29$ & $33.15$& $31.94$ & $31.56$ & $32.42$ & $31.33$ & $32.96$ & $\underline{33.54}$ & $\textbf{38.88}$\\
  
  Pantheon Exterior & $54.22$ & $59.23$ & $\underline{62.26}$& $61.86$ & $61.60$ & $60.54$&$49.89$ &$56.81$& $61.31$ & $\textbf{63.81}$ \\
  
  Prague Old Town Square & $26.61$ & $30.14$ & $32.48$& $32.39$ & $31.30$ & $34.35$& $34.05$ & $\textbf{37.58}$ & $33.13$& $\underline{36.95}$\\
  
  Sacre Coeur & $41.58$ & $49.03$& $56.36$ & $53.23$ &$53.06$ & $45.10$ & $41.30$ &$45.09$ & $\underline{56.61}$ &$\textbf{60.15}$\\ 
  
  Taj Mahal & $38.43$ & $48.44$ & $51.51$& $50.71$ & $50.43$ &$51.63$ & $50.04$ &$52.17$ &$\underline{54.71}$ & $\textbf{57.79}$\\
  
  Trevi Fountain & $29.85$ & $31.67$ & $34.99$& $33.94$ & $34.28$ &$33.75$& $27.69$ &$30.82$ & $\underline{34.61}$& $\textbf{37.66}$\\
  
  Westminster Abbey & $50.97$ & $52.27$ & $\underline{55.99}$& $55.15$ & $54.91$ & $53.07$ &$42.10$ & $48.37$ & $53.61$ &$\textbf{56.28}$\\ 
  \cline{0-0}
  
  Avg.\ over all scenes & $35.00$ & $40.10$ & $43.41$ & $42.68$ & $42.46$ & $42.91$& $36.91$ & $40.67$ & $\underline{43.66}$& $\textbf{46.96}$ \\ \bottomrule
    \end{tabular} }
    \caption
  Fundamental matrix estimation on 12000 image pairs from the PhotoTourism~\cite{IMC2020}. 
  We report the the average run-times (ms) in the first row and F1 scores on each scene. 
  The score averaged over all scenes is shown in the last row. 
  We use the threshold tuned in~\cite{barath2021marginalizing} for RANSAC, GC-RANSAC, MAGSAC, and MAGSAC++. 
  We tuned the parameters of EAS, and retrained OANet, CLNet, NG-RANSAC on the same datasets that is used for training the proposed $\nabla$-RANSAC.
  The results with the pre-trained models provided by the authors are in Table~\ref{tab:pretrained_clnet_compare_f}. 
  The best results are \textbf{bold}.
  The second best results are \underline{underlined}.
 
    }
  \label{tab:photo_result_f}
\end{table*}
\begin{figure*}[t]
    \centering
   
    \includegraphics[width=2.1\columnwidth]{Figures/Layouttest-eps-converted-to.pdf}
   
   
   
   
   
   
   
   
   
    \caption{ The cumulative distribution functions (CDF) of the epipolar errors (left; in pixels), run-times (middle; in seconds), and F1 score (right; in percentages) for $\matr{F}$ estimation.
   
   
    We use the thresholds tuned in \cite{barath2021marginalizing} for the traditional algorithms. 
    We trained OANet, CLNet, NG-RANSAC, and $\nabla$-RANSAC on the same datasets.
    In the left two plots, being close to the top-left corner indicates accurate results. 
    In the right one (F1 score), the bottom-right corner is preferable. 
    }
    \label{fig:cdf}
\end{figure*}

\noindent{\textbf{Training.}}
First, we apply a 1000 epoch-long weight initialization procedure that minimizes the Kullback–Leibler divergence~\cite{van2014renyi} of the predicted and target probabilities~\cite{brachmann2019neural}. 
Along with the initialized weights, the gradient clipping~\cite{chen2020understanding} technique is used to avoid exploding gradients, make the training stable and accelerate the convergence. 
\begin{table}[t]
    \centering
    \resizebox{0.95\columnwidth}{!}{\begin{tabular}{lcccc}
    \toprule 
    Method &  \tabincell{c}{OANet \\ \cite{oanet2019}}& \tabincell{c}{CLNet \\ \cite{clnet2021}} & \tabincell{c}{NG-RSC \\ \cite{brachmann2019neural}} & $\nabla$-RANSAC \\ \midrule
    F1 Score (\%) & $42.29$ & $38.61$ & $45.80$ & $\textbf{46.96}$\\
    Avg. Med.\ epi.\ error (px) & \phantom{1}$2.50$ & \phantom{1}$8.73$ & \phantom{1}$1.51$ &\phantom{1}$\textbf{1.00}$ \\
    Time (ms) & $\textbf{21.75}$ &$ 27.83$ & $25.90$ &$35.08$ \\
    \bottomrule
    \end{tabular}}
        \caption{Comparison of \textbf{F} estimation on the PhotoTourism dataset with the models given by the authors of NG-RANSAC, OANet and CLNet. CLNet and OANet were trained on more than 541K image pairs from YFCC~\cite{thomee2016yfcc100m}. Note that we train on \num{4950} image pairs of one specific scene and show good generalization ability to real-world data. NG-RANSAC uses twice as many images as $\nabla$-RANSAC.
    }
    \label{tab:pretrained_clnet_compare_f}
\end{table}

In the \textbf{F} matrix case, we normalize the point coordinates for consensus learning, and use the original unnormalized ones in minimal solvers. 
We train the pipeline using the 8PC algorithm for \textbf{F} estimation with a fixed 1000 iterations.
For \textbf{E} matrices, we use the 5PC algorithm and 100 iterations.
The entire training pipeline is implemented in PyTorch. 

\noindent
\textbf{Inference.}
Since in test time, the back-propagation of the gradients is not required, we equip $\nabla$-RANSAC with state-of-the-art components.
The end-to-end trained model provides importance scores to the Gumbel Softmax Sampler.
For selecting the best model, we use the MAGSAC++~\cite{barath2020magsac++} model quality function, marginalizing over a range of noise scales. 
We also apply an inner RANSAC-based local optimization~\cite{lebeda2012fixing} to further improve the accuracy.
Also, we perform the Levenberg-Marquardt~\cite{more1978levenberg} numerical optimization minimizing the pose error on all inliers as a final step.
The testing algorithm is implemented in C++ to make it run fast. 




\section{Experimental Results}
\label{sec:results}

We test epipolar geometry estimation on 13 scenes from the training set of the CVPR IMW 2020 PhotoTourism benchmark \cite{IMC2020} that provides images, intrinsic and extrinsic camera parameters from a reference COLMAP reconstruction, and pre-detected RootSIFT features~\cite{RootSIFT2012}.  
We train and validate the method on scene St.\ Peter's Square consisting a total of 4950 image pairs. 
We split it $3$ to $1$ into training and validation sets, respectively.
The other 12 scenes are used for testing.
Note that using a bigger training set would likely boost the performance even further.

We compare $\nabla$-RANSAC on fundamental (Section~\ref{sec:fundamental}) and essential matrix (Section~\ref{sec:essential}) estimation, to classical robust estimators, \ie, RANSAC~\cite{RANSAC}, LMEDS~\cite{rousseeuw2005robust}, or their recent alternatives, like GC-RANSAC~\cite{barath2018graph}, EAS~\cite{fan2021efficient}, MAGSAC~\cite{barath2019magsac}, and MAGSAC++~\cite{barath2020magsac++}. 
Also, we run state-of-the-art learning-based methods, OANet~\cite{oanet2019}, CLNet~\cite{clnet2021} and NG-RANSAC~\cite{brachmann2019neural}.
To make fair comparison, we retrained NG-RANSAC~\cite{brachmann2019neural}, CLNet~\cite{clnet2021}, and OANet~\cite{oanet2019} on the same data as what we use for training. 


\noindent{\textbf{Technical details.}}
When training, we provide the consensus learning module with the SNN ratio~\cite{lowe1999object} coming from RootSIFT descriptor matching, and the underlying feature scales and orientations as learnable side-information. 
Also, we pre-filter the correspondences by an SNN ratio threshold set to $0.8$.
We use $0.75$ pixels as the inlier-outlier threshold for robust estimators, normalized by the focal lengths in case of essential matrix estimation. 

All the experiments were conducted on Ubuntu 20.04 with GTX 3090Ti, OpenCV 4.5.5, and PyTorch 1.11.1 with Cuda 11.3.1. 
Due to the batch processing of PyTorch, we re-implemented RANSAC in PyTorch and made it easy to connect with other modules in an end-to-end fashion. 

\subsection{Fundamental Matrix Estimation}
\label{sec:fundamental}

For fundamental matrix estimation, 
the coordinates are normalized by the image sizes before training.
Benefiting from the fully differentiable $\nabla$-RANSAC, we trained the model parameters end-to-end jointly with the prediction network to learn the statistical features of tentative matches and predict inlier probabilities. 
We adopt the consensus learning block of~\cite{clnet2021} but do not filter the points with their predicted probabilities. 
The model is trained for 10 epochs, optimized by Adam~\cite{kingma2014adam}, with the weights of the layers initialized by 1000-epoch procedure.
For testing, we use \num{1000} randomly chosen image pairs from each of the remaining \num{12} scenes -- these will be made public to maximize reproducibility. 
Thus, the methods are tested on a total of \num{12000} image pairs.
To measure the quality of the estimated fundamental matrices, we use F1 score, in percentages, and median epipolar errors in pixels.
We use the normalized 8PC algorithm when training on fundamental matrix estimation as it leads to more stable solutions than the 7PC solver in our experiments.
The iteration number is fixed to 1000. 

The average F1-scores and the run-time of the robust estimation (in milliseconds) are reported in Table~\ref{tab:photo_result_f} on each scene of the dataset and, also, averaged over all in the last row. 
The proposed $\nabla$-RANSAC achieves the most accurate results on all but one scene, where it is the second best by a small margin. 
On average, it improves by ${\sim}3$\% compared to the second best method that is, interestingly, NG-RANSAC in this case.
We consider NG-RANSAC being the second best an interesting result given that it is based on simple MLPs with batch normalization. 
$\nabla$-RANSAC runs at comparable speed to its fastest, but less accurate alternatives.
It is far real time. 
The most efficient method is LMeDS achiving 11.61\% lower F1 score than $\nabla$-RANSAC.
Note that these timings do not include the inference time which is around $2$ milliseconds on average. 

\begin{table}[t]
    \centering
   
    \resizebox{0.95\columnwidth}{!}{\begin{tabular}{lccccc}
    \toprule 
    Method & AUC@$5^\circ$ & AUC@$10^\circ$ & AUC@$20^\circ$& Time (ms)\\ \midrule
    LMEDS\cite{rousseeuw2005robust} & 0.24 & 0.30 & 0.37 & \textbf{\phantom{1}27}\\
    RSC\cite{RANSAC}& 0.26 & 0.32 & 0.40 & \phantom{1}88\\
    GC-RSC\cite{barath2018graph} & 0.33 & 0.37 & 0.42 & 175\\
    MSC\cite{barath2019magsac} & 0.37 & 0.42 &0.47 & 239\\
    MSC++\cite{barath2020magsac++}& 0.37 & 0.42 & 0.47 & 113\\
    EAS \cite{fan2021efficient} & 0.24 & 0.28 &0.34 & 325\\
    OANet~\cite{oanet2019}&0.29 & 0.33 & 0.39 & \phantom{1}49\\
    CLNet~\cite{clnet2021}&0.34 &0.40 & 0.47 & \phantom{1}58\\
    NG-RSC~\cite{brachmann2019neural} &0.35&0.41& 0.47& \phantom{1}80\\
    $\nabla$-RANSAC & \textbf{0.41} & \textbf{0.45} & \textbf{0.50} & 122\\
    \bottomrule
    \end{tabular}}
    \caption{The average AUC scores of $\nabla$-RANSAC and comparison methods over 12 scenes from PhotoTourism, under different thresholds.  We are the most accurate method for \textbf{E} estimation
    }
    \label{tab:e}
\end{table}

Furthermore, in Fig.~\ref{fig:cdf}, we show the cumulative distribution functions (CDF) of the epipolar errors, run-times and F1 scores calculated from the 12000 tested image pairs.
In the two left plots, the epipolar error and the run-times CDFs are shown. Being close to the top-left corner indicates better performance.
The epipolar error curve of the proposed $\nabla$-RANSAC is above other methods on the entire plot.
While it is not the fastest, it finishes in $0.1$ seconds in $100$\% of the test cases.
This confirms that $\nabla$-RANSAC is applicable in time-sensitive applications. 
The right figure shows the CDFs of the F1-scores. 
In contrast to the other plots, here, being close to the bottom-right corner is preferred. 
$\nabla$-RANSAC has better score on the entire range.
It coincides with other methods only at the end of the range, where all methods are supposed to have $1$. 

The avg.\ results on all scenes when $\nabla$-RANSAC is trained with the 5PC, 7PC and norm.\ 8PC solvers are reported in Table~\ref{tab:F_solvers}.
For \textbf{F} estimation, the norm.\ 8PC solver leads to the most accurate results while also being the fastest.
The 5PC solver is close to it in terms of accuracy.

\subsection{Essential Matrix Estimation}
\label{sec:essential}

We evaluate our method for \textbf{E} estimation with the same train, validation, and test scenes as \textbf{F} estimation.
The correspondences are normalized by the intrinsic matrices. 
We use the differentiable 5PC algorithm when training on essential matrix estimation as it leads to more stable solutions than the 8PC solver in our experiments.
We train the method end-to-end for 10 epochs with weight initialization and gradient clipping.
The iteration number is fixed to 100.

To evaluate the estimated essential matrix, we decompose the \textbf{E} matrix to rotation and translation, calculate their errors $\epsilon_\mathbf{R}$, $\epsilon_\mathbf{t}$ and report the maximum $\max(\epsilon_\mathbf{R}, \epsilon_\mathbf{t})$.
We calculate the Area Under the Recall curve (AUC) thresholded at $5^\circ$, $10^\circ$ and $20^\circ$~\cite{brachmann2019neural}. 

The AUC scores, averaged over all scenes, are reported in Table~\ref{tab:e}.
The highest AUC scores, at all thresholds, are achieved by $\nabla$-RANSAC. For example, its AUC@$5^\circ$ score is higher than that of the second best methods (\ie, MAGSAC and MAGSAC++) by AUC $3$ points. 
It has comparable run-times to other less accurate alternatives. 

The results of minimal solvers are reported in Table~\ref{tab:E_solvers}.
The highest AUC scores, at all thresholds, are achieved by the five-point algorithm, confirming the necessity of use better minimal solvers instead of the 8PC algorithm.


\subsection{Ablation Studies}
\label{subsec:ablation}
\begin{figure}[t]
	\centering
 \begin{subfigure}[b]{0.495\columnwidth}
    \includegraphics[width =1\columnwidth, trim=0.2mm 0mm 0mm 0mm, clip]{Figures/technical_imA-eps-converted-to.pdf}
    \end{subfigure}
    \begin{subfigure}[b]{0.495\columnwidth}
    \includegraphics[width = 1\columnwidth, trim=0mm 0mm 0mm 0mm, clip]{Figures/Graph36A-eps-converted-to.pdf}
    \end{subfigure}
    \caption{\textbf{Ablation study}: AUC@10$^\circ$ score of \textbf{E} estimation on scene Sacre Coeur from~\cite{IMC2020}.
    \textit{Left:} $\nabla$-RANSAC trained with the 5PC and 8PC solvers with and without weight initialization.
    \textit{Right:} Score when using uniform sampler and Gumbel sampler with SNN ratios or model trained by minimized different loss functions.}
    \label{fig:technical}
\end{figure}

The effect of the weight initialization with the Kullback–Leibler divergence~\cite{van2014renyi} with different minimal solvers on essential matrix estimation is shown in the left plot of Fig.~\ref{fig:technical}.
For both minimal solvers, the initialization clearly improves the final accuracy.
Training with the 5PC solver without performing initialization leads to the same results as training with 8PC with initialization.
The best results are obtained by the 5PC solver with pre-trained weights. 
This again confirms the importance of using other solvers than the 8PC algorithm. 

The right plot of Fig.~\ref{fig:technical} shows different sampling weights.
For this plot, we only used scene Sacre Coeur as test set. 
Uniform sampler is basically the results of plain RANSAC.
Gumbel Sampler with SNN ratios requires no training.
It uses the pre-calculated SNN ratios to guide the sampling.
Finally, the effect of training with different loss functions is shown. 
On this dataset and problem, the epipolar error leads to the best results. 


\begin{table}
    \centering
    \resizebox{0.75\columnwidth}{!}{\begin{tabular}{l | ccc}
    \toprule
    \textbf{E} estimation &  5PC~\cite{nister2004efficient} & 7PC~\cite{hartley2003multiple} & 8PC~\cite{hartley2003multiple} \\ 
    \midrule 
         AUC@$5^\circ$  &\textbf{0.41}&0.37 &0.38\\
         AUC@$10^\circ$ & \textbf{0.45}& 0.42 & 0.43\\
         AUC@$20^\circ$ &\textbf{0.50}& 0.46 & 0.48\\
        \bottomrule
    \end{tabular}}
    \caption{
    Performance of \textbf{differentiable minimal solvers} on \textbf{E} estimation used for training $\nabla$-RANSAC.
    The AUC scores are reported on 12000 image pairs from PhotoTourism~\cite{IMC2020}.}
    \label{tab:E_solvers}
\end{table}

\begin{table}
    \centering
    \resizebox{0.85\columnwidth}{!}{\begin{tabular}{lccc}
    \toprule
    \textbf{F} estimation &  5PC~\cite{nister2004efficient} & 7PC~\cite{hartley2003multiple} & 8PC~\cite{hartley2003multiple} \\  \midrule 
         F1 score (\%) & 45.6 & 43.7 & \textbf{46.9}\\ 
         Med.\ epi.\ error (px) & \phantom{1}1.5 & \phantom{1}2.1 & \phantom{1}\textbf{1.0} \\
       
        \bottomrule
    \end{tabular}}
    \caption{Performance of \textbf{differentiable minimal solvers} on \textbf{F} estimation used for training $\nabla$-RANSAC.
    The AUC scores are reported on 12000 image pairs from PhotoTourism~\cite{IMC2020}.}
    \label{tab:F_solvers}
\end{table}

\subsection{Learning Feature Matching with $\nabla$-RANSAC}

An advantage of the fully differentiable RANSAC is that it allows the RANSAC pipeline to be optimized for test-time evaluation metrics such as pose error, inlier classification, and epipolar error (see~\cref{sec:quality}). Note that these metrics are differentiable in themselves. However, their input comes from the RANSAC pipeline, so if the RANSAC is non-differentiable, these metrics cannot be directly used as a loss function. Additionally, if the RANSAC is differentiable, the entire feature computation and matching module can be directly optimized on the above-mentioned metrics. In this section, we tune an end-to-end feature matching approach, LoFTR~\cite{Sun_2021_CVPR} on pose error, inlier classification, and epipolar error using the proposed $\nabla$-RANSAC.

Following the work of Patel {\em et al.}\cite{patel2020learning,patel2021feds} on learning surrogate loss functions, the LoFTR and the $\nabla$-RANSAC models are initialized with the pre-trained weights and post-tuned. To train LoFTR jointly with RANSAC and optimize the entire pipeline on the evaluation metric directly, we consider three training settings: LoFTR training, joint training, and alternating training. The {\em LoFTR training} setup refers to the setting where the pre-trained $\nabla$-RANSAC model is kept frozen, and only the LoFTR model is trained. The {\em joint training} setup refers to the setting where both the LoFTR and $\nabla$-RANSAC models are trained jointly. This setup involves adding a differentiable module on top of the LoFTR to obtain final predictions, and the loss is computed on these final predictions. The third setup of {\em alternating training} is similar to the training setup in \cite{patel2020learning} where at first the parameters of the LoFTR model are kept frozen, and only the $\nabla$-RANSAC is trained, and then only the LoFTR model is trained when $\nabla$-RANSAC is kept frozen. This alternating training is repeated for each epoch. In the last setup, the $\nabla$-RANSAC is part of the loss computation.

The evaluation is performed on the $12$ scenes from the PhotoTourism dataset, which contains $12$K image pairs in total. For the three training setups mentioned above, the results are presented in Table \ref{tab:loftr}. 
Note that for inference, the maximum number of RANSAC iterations is set to $5$K, while the actual number of iterations adapts dynamically based on the number of inliers in each iteration. The results show that the proposed $\nabla$-RANSAC can be used to improve end-to-end feature matching approaches and is robust to different training setups.
The best results were observed with the setup where only the LoFTR model is trained, and the $\nabla$-RANSAC is kept frozen with the pre-trained weights. This observation clearly shows the robustness of $\nabla$-RANSAC as a pre-trained model without any adaptation can be used in a plug-and-play manner.


\begin{table}[t]
    \centering
    \resizebox{0.99\columnwidth}{!}{\begin{tabular}{lccccc}
    \toprule 
    LoFTR ~\cite{Sun_2021_CVPR} & Robust Estimator & Training Setup & F1 score (\%) & run-time (ms) \\ \midrule
    pre-trained & RANSAC~\cite{RANSAC} & -  & 26.70 & \phantom{1}59.06 \\
    pre-trained & pre-trained $\nabla$-RSC & - & 29.38 & 138.67 \\
   
    post-tuned & pre-trained $\nabla$-RSC & LoFTR & \textbf{29.66} & \phantom{1}\textbf{58.85}\\
    post-tuned & post-tuned $\nabla$-RSC& joint & 29.44 & \phantom{1}63.87\\
    post-tuned & post-tuned $\nabla$-RSC  & alternating & 29.52 & \phantom{1}61.37\\
    \bottomrule
    \end{tabular}}
    \caption{The performance of training feature matcher with $\nabla$-RANSAC compared with the pre-trained models and non-differentiable robust estimator. The testing is performed on $12$K image pairs from PhotoTourism dataset with varying image sizes. All the setups, \eg, LoFTR training, joint training, and alternating training outperforms pre-trained LoFTR and $\nabla$-RANSAC models.
   
    }
    \label{tab:loftr}
\end{table}

\section{Conclusion}
\label{sec:conclusion}

In this paper, we propose the fully differentiable $\nabla$-RANSAC, \ie, the first end-to-end trainable randomized robust estimator. 
It predicts the inlier probabilities of the input data points, exploits the predictions in a guided sampler, and estimates the model parameters (e.g., fundamental matrix) and its quality while propagating the gradients through the entire procedure.    
$\nabla$-RANSAC leads to the most accurate fundamental and essential matrices compared to state-of-the-art robust estimators, on 12K image pairs from the PhotoTourism dataset.
Also, it is one of the fastest algorithms. 
To demonstrate its potential in unlocking the end-to-end training of geometric pipelines, we train $\nabla$-RANSAC together with a recent detector-free feature matcher, LoFTR \cite{Sun_2021_CVPR}, with which we achieve improved accuracy. 
We believe that the community will benefit from the algorithm.
We will make the code repository and the trained model public. 

{\small
\bibliographystyle{ieee_fullname}

\section{Introduction}
\label{sec:intro}

A direct optimization on the test-time evaluation metric has shown to be beneficial in deep learning for many vision tasks \cite{arulkumaran2017brief,lee2018context,balle2018variational,berrada2018smooth,lapin2016loss}. 
Training the model directly on the evaluation metric is infeasible when the metric is non-differentiable. 
In such cases, training with a surrogate of the metric has shown to be effective. Examples include surrogate losses of average precision and recall@k for image retrieval \cite{brown2020smooth,patel2022recall,cakir2019deep,ramzi2021robust}, perceptual loss for image compression \cite{patel2021saliency,patel2019deep}, intersection-over-union loss for object detection \cite{rezatofighi2019generalized,yu2016unitbox,patel2020learning}, edit distance loss for scene text recognition \cite{patel2020learning,patel2021feds}, etc. 
However, a differentiable version of RANSAC~\cite{RANSAC} still does not exist, prohibiting a direct optimization on the evaluation metrics for several tasks that are utmost important in computer vision and robots.

Robust estimation is a fundamental component in vision pipelines, including relative pose estimation~\cite{chum2004epipolar}, wide baseline matching~\cite{pritchett1998wide,matas2004robust,mishkin2015mods}, multi-model fitting~\cite{isack2012energy,pham2014interacting}, image-based localization~\cite{Brachmann_2017_CVPR}, motion segmentation~\cite{torr1993outlier}, and pose graph initialization of Structure-from-Motion (SfM) algorithms~\cite{schoenberger2016sfm, schoenberger2016mvs}.
While several robust estimators have been proposed throughout the years~\cite{holland1977robust,illingworth1988survey,li2009consensus,xiao2016hypergraph}, randomized hypothesize-and-verify approaches, like RANSAC~\cite{RANSAC} and its recent variants~\cite{barath2018graph,barath2019magsac,barath2020magsac++,ivashechkin2021vsac}, have become the most widely used methods due to their robustness, simplicity, and efficiency. 
RANSAC repeatedly selects random minimal subsets of the input data and fits a model, \eg, a 3D plane to three points or fundamental matrix to seven point correspondences.
The model quality is then computed as the cardinality of its support, \ie, the number of points that fall closer than a manually set inlier threshold. 
The so-far-the-best model is updated if a new model is found with higher quality. 
RANSAC terminates when the probability of finding an all-inlier sample falls below a threshold, controlled by a confidence parameter typically set to 0.95 or 0.99. 
Finally, the model with the highest quality, polished, \eg, by least-squares fitting on all inliers, is returned.

\begin{figure}
    \centering
   
    \includegraphics[width=1\columnwidth]{Figures/fig1new.pdf}
   
   
   
   
    
   
   
   
   
   
    \caption{$\nabla$-RANSAC gets a set of tentative correspondences as input, predicts their inlier probabilities, runs a differentiable RANSAC procedure and, finally, estimates the model parameters, \eg, fundamental matrix (shown with its inliers here).
   
    }
    \label{fig:example}
\end{figure}

\begin{figure*}[t]
	\centering
    \includegraphics[width = 2.1\columnwidth, trim=0mm 0mm 0mm 0mm, clip]{Figures/pipeline.pdf}
    \caption{\textbf{Pipeline of $\nabla$-RANSAC.}
    Given an image pair, a feature matcher and detector is applied to obtain a set of tentative correspondences. 
    The point coordinates and other side information, like SNN ratio, feature orientation and scale, are fed into the consensus learning network from \cite{clnet2021} to predict inlier probabilities.
    The predictions go into the fully differentiable RANSAC module, $\nabla$-RANSAC.
    In each iteration, the differentiable and randomized Gumbel sampler (Section~\ref{sec:sampler}) selects a minimal sample $\mathcal{S}$ of $m$ correspondences. 
    Model $\hat{\theta}$ is estimated by a differentiable minimal solver (Section~\ref{sec:diff_solver}) and its score is calculated (Section~\ref{sec:quality}) from the set of tentative correspondences.}
   \label{fig:flowchart}
\end{figure*}

Since the publication of RANSAC, there have been a number of important modifications proposed to improve specific components of the original algorithm.
Improvements to the classical framework include local optimization (LO-RANSAC, LO$^+$-RANSAC, GC-RANSAC)~\cite{chum2003locally,lebeda2012fixing,barath2018graph}, 
advanced scoring functions (MLESAC, MSAC, MAGSAC, and MAGSAC++)~\cite{torr2000mlesac,barath2019magsac,barath2020magsac++,barath2022learning}, 
speed-ups using probabilistic sampling of hypotheses (PROSAC, NAPSAC, and P-NAPSAC)~\cite{torr2002napsac,chum2005matching,barath2020magsac++} or by preemptive verification strategies (SPRT, SP-RANSAC)~\cite{chum2008optimal,barath2021space}, 
degeneracy checks (DEGENSAC, QDEGSAC, NeFSAC)~\cite{chum2005two,frahm2006ransac,cavalli2022nefsac}, and strategies to auto-tune the inlier-outlier threshold (MINPRAN, a contrario RANSAC)~\cite{stewart1995minpran,moisan2012automatic}.
According to the recent CVPR 2022 Image Matching Challenge~\cite{imc2022}, the most accurate method for relative pose estimation currently is MAGSAC++~\cite{barath2020magsac++} with PROSAC sampling~\cite{chum2005matching} and DEGENSAC-based degeneracy~\cite{chum2005two} testing. 
Regardless of their good performance, none of the RANSAC variants is differentiable, keeping the randomized robust estimation one of the only components in geometric image matching pipelines that can not be trained end-to-end.





In the recent years, several algorithms were proposed approaching the robust relative pose estimation problem using neural networks. 
Context normalization networks (CNe)~\cite{cne2018} is the first paper on the topic which proposes to use PointNet (MLP)~\cite{qi2017pointnet}, \ie a multi-layer perceptron (MLP), with batch normalization~\cite{batchnorm2015} as context mechanism. 
Attentive context normalization networks~\cite{acne2020} introduces a special architectural block for the task. 
Deep Fundamental matrix estimation~\cite{dfe2018} uses differentiable iteratively re-weighted least-squares to predict inliers weights for the input correspondences. 
The OANet algorithm~\cite{oanet2019} introduces several architectural blocks for the correspondence filtering.
Neural Guided RANSAC~\cite{brachmann2019neural} uses a CNe-like architecture with a different training objective.
It assumes the inlier probabilities to follow the categorical distribution and approximates the gradients.
The predicted correspondence scores are exploited inside RANSAC by using a guided sampling algorithm that helps in finding accurate models early. 
CLNet~\cite{clnet2021} introduces several algorithmic and architectural improvements  to remove gross outliers with iterative pruning. 
It propagates the gradients between the input data and estimated models by using the weighted eight-point algorithm~\cite{hartley1997defense}.
Exploiting the differentiable pooling layers in OANet~\cite{oanet2019}, MS2DG-Net~\cite{Dai_2022_CVPR} uses additional semantic features based on sparse similarities while maintaining its permutation equivariant property. 
DSAC~\cite{Brachmann_2017_CVPR} is the first method to propose a differentiable RANSAC-like robust estimator, however, only for a highly specific task, \ie, one-shot camera re-localization from an RGB image.
Instead of finding directly the best model, DSAC employs soft probabilistic hypothesis selection. 

Even though most of these algorithms are trainable in an end-to-end manner, \eg via employing iteratively re-weighted least-squares fitting with the predicted weights, they can be considered elaborate inlier pre-filtering approaches -- a final RANSAC on the top-\textit{k} correspondences \textit{always} improves the results compared to LS fitting~\cite{cvpr2020ransactutorial}.

In this paper, we focus on making RANSAC end-to-end differentiable and propose $\nabla$-RANSAC (see Fig.~\ref{fig:example}). 
The main contribution is as follows:
we investigate all components of the RANSAC pipeline and propose differentiable alternatives that allow learning inlier probabilities in an end-to-end manner.
This includes a new random sampling approach based on a clever re-parametrization strategy, \ie Gumbel Softmax sampler, that allows for propagating the gradients through the entire random procedure.
Moreover, we implement differentiable versions of the widely used minimal solvers, \eg five and seven-point algorithms~\cite{nister2004efficient}.
Finally, we propose a model quality function, avoiding the non-differentiable best-model-selection of RANSAC, via marginalizing over all models estimated within the procedure.
$\nabla$-RANSAC allows robust estimators to use test-time evaluation metrics to optimize end-to-end training.

As a proof of its potential to unlock the end-to-end training of geometric pipelines, we train $\nabla$-RANSAC together with a recent detector-free feature matcher, LoFTR \cite{Sun_2021_CVPR}, with which we achieve improved accuracy. 
We keep it as future work to find the best way, \eg losses and network architecture, that maximizes performance when the entire pipeline is trained end-to-end. 

\section{Fully Differentiable RANSAC}
\label{sec:differentiable_RANSAC}

In this section, we discuss the algorithmic components of the proposed $\nabla$-RANSAC making it capable of end-to-end training.
A pipeline employing $\nabla$-RANSAC is in Fig~\ref{fig:flowchart}.

The input is a set of tentative point correspondences, possibly equipped with extra information from the detector and matcher, \eg, feature orientation, scale, affine shape, and matching confidence. 
An initial step of the algorithm is based on consensus learning via exploiting the pruning block from the recent \cite{clnet2021}.
Instead of using the weights to iteratively remove outliers as in \cite{clnet2021}, we consider them initial inlier probabilities that are fed into the subsequent fully differentiable RANSAC module, $\nabla$-RANSAC. 
Similarly to plain RANSAC, the proposed $\nabla$-RANSAC is an iterative random sampling of $m$ data points, model estimation via a minimal solver, and model scoring.
The devil is, however, in the details. 
The sampler is the \textit{Gumbel Softmax Sampler} that exploits the input probabilities as guidance in a randomized manner while propagating the gradients. 
The differentiable minimal solver estimates the model parameters from the drawn sample. 
The model quality is computed in a supervised way using the ground truth, \eg, it is the relative pose error calculated from the estimated essential matrix and the known pose from a reference reconstruction.\footnote{Note that the algorithm should work with self-supervised losses, like inlier ratio, too.
Exploring such approaches we keep as future work.}
While, in RANSAC, the best-model-selection is a non-differentiable step,
we aggregate the model qualities calculated in the $k$ iterations performed throughout the entire procedure.
The minimized loss comes from the average of these qualities, guiding the network to predict inlier probabilities that \textit{on average} lead to accurate models.

\subsection{Gumbel Softmax Sampler}
\label{sec:sampler}


The first step of $\nabla$-RANSAC is hypothesis generation which requires sampling $m$ data points from a set of $n$ total samples $\{x_{i}\}$. 
The sampling distribution is either governed by the importance scores for each sample $\{s_{i}\}$ or follows a uniform distribution. 
The standard sampling operation is either non-differentiable or has sparse gradients that are not useful for training \cite{doersch2016tutorial}. 
Further, even the deterministic selection of top-$k$ samples with highest importance scores is non-differentiable as ranking and counting operations require a Heaviside step function \cite{patel2022recall}. Therefore, to make randomized importance score sampling differentiable for RANSAC, we extend the Gumbel-Softmax \cite{jang2016categorical,maddison2016concrete} with the straight through trick \cite{van2017neural} to multiple samples.

The sampling is performed on a distribution that is governed by the importance scores for each samples. Let this distribution be $\mathbf{A}$ as follows:
\begin{equation}
\mathbf{A}_{i} = \frac{\exp(s_{i} + \gamma_{i})}{\sum_{j=1}^{n} \exp(s_{j} + \gamma_{j})},
\end{equation}
where ${\gamma_{i}}$ are i.i.d random samples drawn from a Gumbel distribution $\text{G}(0,1)$. 
Parameters $\{s_{i}\}$ are the importance scores for each sample in the set. 
Note that if the scores are not available, $\{s_{i}\}$ are set to a constant value to perform uniform random sampling. 
\begin{figure}[t]
	\centering
    \includegraphics[width = 1\columnwidth, trim=0mm 0mm 0mm 0mm, clip]{Figures/sampler.pdf}
    \caption{Overview of \textbf{Gumbel Softmax Sampler} 
 used in $\nabla$-RANSAC. The input to the sampler is the importance scores $s_1, s_2, ..., s_n$. The process starts by i.i.d random sampling $\gamma_1, \gamma_2, ..., \gamma_n$ from a Gumbel(0,1) distribution followed by a Softmax operation. 
 Based on the distribution obtained from the Softmax, $m$ data points from a total of $n$ are drawn by using the top-k operation. 
 Since top-k is non-differentiable, the straight through trick~\cite{van2017neural} is used and the Softmax output is added and subtracted with stop-gradient (sg) to allow backward flow of the gradients.
 The figure is best viewed in color. 
 The green arrows show the connections through which gradient flow is possible. 
 Red arrows show connections without any gradients.}
    \label{fig:gumbel_sampler}
\end{figure}

The feed-forward through the sampler is performed such that the output of the selected indices are binary and gradients can backpropagate for training. Such feed-forward is achieved via the straight through trick \cite{van2017neural} as follows:
\begin{equation}
\mathbf{\hat{A}} = \text{top-k}(\textbf{A},m) + \textbf{A} - \text{sg}(\textbf{A}),
\end{equation}
where sg is the stop-gradient operation and top-k selects the indices of the $m$ highest values. During the forward pass, $\hat{\textbf{A}}$ contains only ones and zeros where one represents the samples that are selected. During back-propagation, the gradients only flow through $\textbf{A}$ as $\text{top-k}(\mathbf{A},m)$ and $\text{sg}(\mathbf{A})$ do not have any gradients. An overview of the proposed differentiable sampler is shown in \cref{fig:gumbel_sampler}. 

\subsection{Differentiable Minimal Solver}
\label{sec:diff_solver}

Minimal solvers~\cite{kukelova2008automatic} are a fundamental part of RANSAC-like hypothesize-and-verify approaches, allowing to progressively explore the input data in a randomized manner. 
Such methods estimate the sought model parameters from a minimal set $\mathcal{S} \subseteq \mathcal{P}$ of $m \in \mathbb{N}$ data points, where $\mathcal{P}$ is the set all points. 
For example, an essential matrix can be estimated from five point correspondences or a 3D plane from a minimum of three 3D points. 
This allows keeping the complexity of the robust estimation problem low via sampling from the possible $\binom{|\mathcal{P}|}{m}$ minimal configurations. 

Most of the minimal solvers commonly used in computer vision are actually differentiable after taking each algorithmic component into careful considerations. 
While we will consider two utmost important problems in this paper, \ie fundamental and essential matrix estimation, the previous statement holds for all geometric model estimation problems that we are aware of in vision. 

Recent learning-based approaches\cite{clnet2021} propagate the gradients through the well-known normalized eight-point (8PC) algorithm~\cite{hartley1997defense} when estimating both essential and fundamental matrices.
There are two main issues with the 8PC algorithm in practice.
First, and most importantly, the 8PC and 7PC solvers have a degeneracy when the points stem from a close-to-planar underlying 3D structure.
In this case, a degenerate model is estimated that, while often having a large number of inliers~\cite{chum2005two}, is incorrect w.r.t.\ the actual scene geometry.  
This misguides the learning in scenes dominated by planar structures and deteriorates the performance. 
Second, using eight correspondences instead of the minimal five (E) or seven (F) largely increases the complexity of the problem. 
This prevents the network from learning the maximal information available in a given correspondence set, and makes the learning sub-optimal. 

For essential matrix estimation, there are a number of 5PC solvers~\cite{nister2004efficient,li2006five,batra2007alternative,stewenius2008minimal,bujnak2012making} proposed throughout the years.
Most frequently, in practical applications (\eg, SfM~\cite{schonberger2016structure,theia-manual}), one of the following methods are used:
the method of Stewenius et al.~\cite{stewenius2008minimal} due to its stability, and that of David Nister~\cite{nister2004efficient} due to its effectiveness. 
The steps of~\cite{stewenius2008minimal} are as follows:
creating the coefficient matrix from the input points, 
decomposing it by SVD, solving a linear system of equations by, \eg, the Moore-Penrose pseudo-inverse, and finally, recovering the linear combination of the null vectors via the Eigen decomposition of the action matrix. 
Even though both SVD~\cite{wang2021robust} and Eigen decompositions~\cite{song2022batch} have differentiable solutions, we found the Eigen decomposition to be unstable for gradient propagation. 

The steps of the solver of Nister~\cite{nister2004efficient} are as follows: creating the coefficient matrix from the input correspondences, 
decomposing it by SVD, solving a linear system of equations, solving a set of polynomial equations, and finally, basic arithmetic operations. 
We implemented a differentiable polynomial solver based on Sturm sequences~\cite{gonzalez1989sturm} and, also, one using companion matrices.
While both algorithms work well, the companion matrix-based solution is the fastest. 
We apply this solver in the proposed $\nabla$-RANSAC.

Fundamental matrix estimation is an easier problem, where the 7PC~\cite{hartley2003multiple} algorithm, besides running an SVD decomposition on the coefficient matrix, only solves a single cubic polynomial.
As we are given a closed-form solution for the cubic problem, we can straightforwardly make the entire algorithm differentiable. 
Note that, in our experiments, we have not observed improvements by replaced 8PC with 7PC.
This might be due to the shared degeneracy.

Most of the minimal solvers return multiple solutions.
Even though they are consistent with the constraints, most of them are inconsistent with the actual underlying scene geometry.
At inference, from the multiple solutions the best one is selected based on the evaluation metric. To mimic the test-time evaluation in training, for each iteration the best solution is selected and then used for the loss computation. Note that the gradients for the selected solution will be used for training, those for not selected solutions will be zero.

We are not aware of public implementations of such solvers, neither in open-source libraries~\cite{van2014scikit,riba2020kornia} nor in standalone public repositories.
Therefore, we consider this as a technical contribution to the community and will make the source codes publicly available.

\subsection{Trainable Quality Function}
\label{sec:quality}
In RANSAC, the quality of an estimated model is calculated as the cardinality of its support, \ie, the inlier number.
Since RANSAC, a number of algorithms~\cite{torr2000mlesac,barath2019magsac,barath2020magsac++,barath2022learning} have improved the performance by better modelling the noise both in the inliers and outliers.
However, all such methods perform a best model selection step based on the maximal quality which renders the procedure non-differentiable. 
Some works~\cite{Brachmann_2017_CVPR} tackle with this problem by employing soft probabilistic hypothesis selection. 
Other methods~\cite{cne2018, clnet2021} combine classification loss with regression and geometry-induced losses~\cite{brachmann2019neural} to reason about the quality of the estimated model and inlier assignments.

Instead of the above solutions, we exploit all models $\{ \hat{\theta}_i \}_{i = 1}^t$ estimated during the fixed $t \in \mathbb{N}^{>0}$ iterations.
In each iteration, the estimated model is compared to the ground truth and its implied loss is calculated, \eg, as the relative pose error in case of essential matrix estimation. 
The minimized loss, in case of relative pose estimation, is $L_\text{pose}$ calculated as the average over all models as follows:
\begin{equation}
    L_\text{pose} = \frac{1}{2t} \sum_{i = 1}^t \epsilon_\mathbf{R}(\hat{\theta}_i^\mathbf{R}, \theta_i^\mathbf{R}) + \epsilon_\mathbf{t}(\hat{\theta}_i^\mathbf{t}, \theta_i^\mathbf{t}),
    \label{eqn:loss_pose}
\end{equation}
where the upper index refers to the rotation ($\mathbf{R} \in \text{SO}(3)$) and translation components  ($\mathbf{t} \in \mathbb{R}^3$) of the pose, $\theta$ is the ground truth, and functions $\epsilon_\mathbf{R}$ and $\epsilon_\mathbf{t}$ return the rotation and translation errors, respectively. 
Functions 
$\epsilon_\textbf{R}(\hat{\theta}_i^\mathbf{R}, \theta^\mathbf{R}) = \cos^{-1} ( ( \text{tr} ( \widehat{\mathbf{R}} \mathbf{R}^\text{T} ) - 1 ) / 2)$ 
and 
$\epsilon_\mathbf{t}(\hat{\theta}_i^\mathbf{t}, \theta^\mathbf{t}) = \cos^{-1} ( \mathbf{\hat{t}}^\text{T} \mathbf{t} / (|\mathbf{\hat{t}}| |\mathbf{t}|))$.

Besides $L_\text{pose}$ we employ other widely used objectives too.
Given the probabilities $\mathcal{\hat P}_i = \{ \hat p_j^i \}_{j = 1}^n$ that model $\theta_i$ implies in the $i$th iteration ($\hat p_j^i \in [0, 1]$), we calculate classification loss $L_\text{inlier}$ w.r.t.\ the inlier probabilities $\mathcal{P} = \{ p_j \}_{j}^n$ determined by the ground truth model.
Probabilities $\{ p_j \}_{j = 1}^n$ can be calculated from the point-to-model residuals making certain assumptions about the noise.
In our algorithm, we follow the RANSAC assumption and $p_i \in \{0, 1\}$ is a binary mask. 
The implied classification loss is as follows:
\begin{equation}
    L_\text{inlier} = \frac{1}{t} \sum_{i = 1}^t {\text{BCE}}( \mathcal{\hat P}_i, \mathcal{P}). 
    \label{eqn:loss_inlier}
\end{equation}
Also, we compute a loss $L_\text{epi}$ summing the average symmetric epipolar errors as follows:
\begin{equation}
    L_\text{epi} = \frac{1}{t |I|} \sum_{i = 1}^t \sum_{p \in \mathcal{I}} \epsilon_\text{epi} (\theta_i, p),
    \label{eqn:loss_epi}
\end{equation}
where $\mathcal{I}$ is the inlier set selected by the ground truth model.


The overall loss minimized in our training procedure is the linear combination of the above ones as follows:
\begin{equation}
    L_\text{overall}  = w_\alpha L_\text{inlier} + w_\beta L_\text{pose} + w_\gamma L_\text{epi},
    \label{eqn:overall}
\end{equation}
where $w_\alpha, w_\beta, w_\gamma \in \mathbb{R}$ are the weighting parameters.
Note that these test-time evaluation metrics cannot be used in the case the input parameter is non-differentiable, \ie, robust estimation results from standard RANSAC.
Thus, it is meaningful to utilize test-time evaluation metrics in the fully differentiable RANSAC we propose.





\subsection{Training and Testing Details}

The input of $\nabla$-RANSAC is a set of correspondences obtained by any feature detector and matcher. 
The number of  matches is fixed to 2000.
We choose the best 2000 ones, based on the matching score, if we are given more.
In case of having fewer correspondences, we fill up the missing values with zeros.
We extract local and global features from the correspondences by a consensus learning block as backbone~\cite{clnet2021}.
We integrate over the extracted geometric information, \eg, scale, and orientation of the local descriptors to help learn the qualities of correspondences.
\begin{table*}[t]
    \centering
  \resizebox{2.1\columnwidth}{!}{\begin{tabular}{lcccccccccc}
      \toprule
  Dataset / Method  &\tabincell{c}{LMEDS \\\cite{rousseeuw2005robust}} & \tabincell{c}{RSC \\\cite{RANSAC}} 
  &  \tabincell{c}{ GC-RSC\\ \cite{barath2018graph}}& \tabincell{c}{MSC\\  \cite{barath2019magsac}}& \tabincell{c}{MSC++ \\
  \cite{barath2020magsac++}} & \tabincell{c}{EAS \\ \cite{fan2021efficient}} & \tabincell{c}{OANet\\ \cite{oanet2019}} &\tabincell{c}{CLNet \\ \cite{clnet2021}} &\tabincell{c}{NG-RSC \\\cite{brachmann2019neural}} & \tabincell{c}{$\nabla$-RANSAC} \\ \midrule
  Avg. time (ms) & $\textbf{21.00}$  & $\underline{30.67}$  & $73.25$  & $281.3$ &$ 318.13$ & $325.83$ &$\textbf{21.00}$ &$34.83$& $25.85$  & $35.08$  \\
  \midrule
  
  Buckingham Palace & $23.26$ & $24.95$ & $26.48$& $27.23$ & $26.28$  & $\underline{28.28} $& $24.70$ & $27.46$ &$28.06$&$\textbf{30.94} $ \\
  
  Brandenburg Gate & $31.74$ & $39.70$ & $43.01$& $42.02$ & $42.43$ &$34.84$ & $42.49$ & $39.69$ & $\underline{43.19}$& $\textbf{46.34}$\\ 
  
  Colosseum Exterior & $43.32$ & $48.25$ & $50.86$& $51.76$ & $51.56$ &$50.49$ & $40.50$& $47.12$ & $\underline{52.42}$ & $\textbf{54.99}$\\ 
  
  Grand Place Brussels & $27.45$ & $31.42$ & $\underline{33.31}$& $32.60$& $33.08 $& $32.40$&$29.10$ & $31.80$& $32.13$& $\textbf{35.49}$\\
  
  Notre Dame Front Facade & $30.39$ & $37.20$ & $40.48$& $39.35$ & $39.00$ &$39.98$& $33.17$ & $38.20$ & $\underline{40.59}$  & $\textbf{44.26}$\\
  
  Palace of Westminster & $22.20$ & $28.29$ & $33.15$& $31.94$ & $31.56$ & $32.42$ & $31.33$ & $32.96$ & $\underline{33.54}$ & $\textbf{38.88}$\\
  
  Pantheon Exterior & $54.22$ & $59.23$ & $\underline{62.26}$& $61.86$ & $61.60$ & $60.54$&$49.89$ &$56.81$& $61.31$ & $\textbf{63.81}$ \\
  
  Prague Old Town Square & $26.61$ & $30.14$ & $32.48$& $32.39$ & $31.30$ & $34.35$& $34.05$ & $\textbf{37.58}$ & $33.13$& $\underline{36.95}$\\
  
  Sacre Coeur & $41.58$ & $49.03$& $56.36$ & $53.23$ &$53.06$ & $45.10$ & $41.30$ &$45.09$ & $\underline{56.61}$ &$\textbf{60.15}$\\ 
  
  Taj Mahal & $38.43$ & $48.44$ & $51.51$& $50.71$ & $50.43$ &$51.63$ & $50.04$ &$52.17$ &$\underline{54.71}$ & $\textbf{57.79}$\\
  
  Trevi Fountain & $29.85$ & $31.67$ & $34.99$& $33.94$ & $34.28$ &$33.75$& $27.69$ &$30.82$ & $\underline{34.61}$& $\textbf{37.66}$\\
  
  Westminster Abbey & $50.97$ & $52.27$ & $\underline{55.99}$& $55.15$ & $54.91$ & $53.07$ &$42.10$ & $48.37$ & $53.61$ &$\textbf{56.28}$\\ 
  \cline{0-0}
  
  Avg.\ over all scenes & $35.00$ & $40.10$ & $43.41$ & $42.68$ & $42.46$ & $42.91$& $36.91$ & $40.67$ & $\underline{43.66}$& $\textbf{46.96}$ \\ \bottomrule
    \end{tabular} }
    \caption
  Fundamental matrix estimation on 12000 image pairs from the PhotoTourism~\cite{IMC2020}. 
  We report the the average run-times (ms) in the first row and F1 scores on each scene. 
  The score averaged over all scenes is shown in the last row. 
  We use the threshold tuned in~\cite{barath2021marginalizing} for RANSAC, GC-RANSAC, MAGSAC, and MAGSAC++. 
  We tuned the parameters of EAS, and retrained OANet, CLNet, NG-RANSAC on the same datasets that is used for training the proposed $\nabla$-RANSAC.
  The results with the pre-trained models provided by the authors are in Table~\ref{tab:pretrained_clnet_compare_f}. 
  The best results are \textbf{bold}.
  The second best results are \underline{underlined}.
 
    }
  \label{tab:photo_result_f}
\end{table*}
\begin{figure*}[t]
    \centering
   
    \includegraphics[width=2.1\columnwidth]{Figures/Layouttest-eps-converted-to.pdf}
   
   
   
   
   
   
   
   
   
    \caption{ The cumulative distribution functions (CDF) of the epipolar errors (left; in pixels), run-times (middle; in seconds), and F1 score (right; in percentages) for $\matr{F}$ estimation.
   
   
    We use the thresholds tuned in \cite{barath2021marginalizing} for the traditional algorithms. 
    We trained OANet, CLNet, NG-RANSAC, and $\nabla$-RANSAC on the same datasets.
    In the left two plots, being close to the top-left corner indicates accurate results. 
    In the right one (F1 score), the bottom-right corner is preferable. 
    }
    \label{fig:cdf}
\end{figure*}

\noindent{\textbf{Training.}}
First, we apply a 1000 epoch-long weight initialization procedure that minimizes the Kullback–Leibler divergence~\cite{van2014renyi} of the predicted and target probabilities~\cite{brachmann2019neural}. 
Along with the initialized weights, the gradient clipping~\cite{chen2020understanding} technique is used to avoid exploding gradients, make the training stable and accelerate the convergence. 
\begin{table}[t]
    \centering
    \resizebox{0.95\columnwidth}{!}{\begin{tabular}{lcccc}
    \toprule 
    Method &  \tabincell{c}{OANet \\ \cite{oanet2019}}& \tabincell{c}{CLNet \\ \cite{clnet2021}} & \tabincell{c}{NG-RSC \\ \cite{brachmann2019neural}} & $\nabla$-RANSAC \\ \midrule
    F1 Score (\%) & $42.29$ & $38.61$ & $45.80$ & $\textbf{46.96}$\\
    Avg. Med.\ epi.\ error (px) & \phantom{1}$2.50$ & \phantom{1}$8.73$ & \phantom{1}$1.51$ &\phantom{1}$\textbf{1.00}$ \\
    Time (ms) & $\textbf{21.75}$ &$ 27.83$ & $25.90$ &$35.08$ \\
    \bottomrule
    \end{tabular}}
        \caption{Comparison of \textbf{F} estimation on the PhotoTourism dataset with the models given by the authors of NG-RANSAC, OANet and CLNet. CLNet and OANet were trained on more than 541K image pairs from YFCC~\cite{thomee2016yfcc100m}. Note that we train on \num{4950} image pairs of one specific scene and show good generalization ability to real-world data. NG-RANSAC uses twice as many images as $\nabla$-RANSAC.
    }
    \label{tab:pretrained_clnet_compare_f}
\end{table}

In the \textbf{F} matrix case, we normalize the point coordinates for consensus learning, and use the original unnormalized ones in minimal solvers. 
We train the pipeline using the 8PC algorithm for \textbf{F} estimation with a fixed 1000 iterations.
For \textbf{E} matrices, we use the 5PC algorithm and 100 iterations.
The entire training pipeline is implemented in PyTorch. 

\noindent
\textbf{Inference.}
Since in test time, the back-propagation of the gradients is not required, we equip $\nabla$-RANSAC with state-of-the-art components.
The end-to-end trained model provides importance scores to the Gumbel Softmax Sampler.
For selecting the best model, we use the MAGSAC++~\cite{barath2020magsac++} model quality function, marginalizing over a range of noise scales. 
We also apply an inner RANSAC-based local optimization~\cite{lebeda2012fixing} to further improve the accuracy.
Also, we perform the Levenberg-Marquardt~\cite{more1978levenberg} numerical optimization minimizing the pose error on all inliers as a final step.
The testing algorithm is implemented in C++ to make it run fast. 




\section{Experimental Results}
\label{sec:results}

We test epipolar geometry estimation on 13 scenes from the training set of the CVPR IMW 2020 PhotoTourism benchmark \cite{IMC2020} that provides images, intrinsic and extrinsic camera parameters from a reference COLMAP reconstruction, and pre-detected RootSIFT features~\cite{RootSIFT2012}.  
We train and validate the method on scene St.\ Peter's Square consisting a total of 4950 image pairs. 
We split it $3$ to $1$ into training and validation sets, respectively.
The other 12 scenes are used for testing.
Note that using a bigger training set would likely boost the performance even further.

We compare $\nabla$-RANSAC on fundamental (Section~\ref{sec:fundamental}) and essential matrix (Section~\ref{sec:essential}) estimation, to classical robust estimators, \ie, RANSAC~\cite{RANSAC}, LMEDS~\cite{rousseeuw2005robust}, or their recent alternatives, like GC-RANSAC~\cite{barath2018graph}, EAS~\cite{fan2021efficient}, MAGSAC~\cite{barath2019magsac}, and MAGSAC++~\cite{barath2020magsac++}. 
Also, we run state-of-the-art learning-based methods, OANet~\cite{oanet2019}, CLNet~\cite{clnet2021} and NG-RANSAC~\cite{brachmann2019neural}.
To make fair comparison, we retrained NG-RANSAC~\cite{brachmann2019neural}, CLNet~\cite{clnet2021}, and OANet~\cite{oanet2019} on the same data as what we use for training. 


\noindent{\textbf{Technical details.}}
When training, we provide the consensus learning module with the SNN ratio~\cite{lowe1999object} coming from RootSIFT descriptor matching, and the underlying feature scales and orientations as learnable side-information. 
Also, we pre-filter the correspondences by an SNN ratio threshold set to $0.8$.
We use $0.75$ pixels as the inlier-outlier threshold for robust estimators, normalized by the focal lengths in case of essential matrix estimation. 

All the experiments were conducted on Ubuntu 20.04 with GTX 3090Ti, OpenCV 4.5.5, and PyTorch 1.11.1 with Cuda 11.3.1. 
Due to the batch processing of PyTorch, we re-implemented RANSAC in PyTorch and made it easy to connect with other modules in an end-to-end fashion. 

\subsection{Fundamental Matrix Estimation}
\label{sec:fundamental}

For fundamental matrix estimation, 
the coordinates are normalized by the image sizes before training.
Benefiting from the fully differentiable $\nabla$-RANSAC, we trained the model parameters end-to-end jointly with the prediction network to learn the statistical features of tentative matches and predict inlier probabilities. 
We adopt the consensus learning block of~\cite{clnet2021} but do not filter the points with their predicted probabilities. 
The model is trained for 10 epochs, optimized by Adam~\cite{kingma2014adam}, with the weights of the layers initialized by 1000-epoch procedure.
For testing, we use \num{1000} randomly chosen image pairs from each of the remaining \num{12} scenes -- these will be made public to maximize reproducibility. 
Thus, the methods are tested on a total of \num{12000} image pairs.
To measure the quality of the estimated fundamental matrices, we use F1 score, in percentages, and median epipolar errors in pixels.
We use the normalized 8PC algorithm when training on fundamental matrix estimation as it leads to more stable solutions than the 7PC solver in our experiments.
The iteration number is fixed to 1000. 

The average F1-scores and the run-time of the robust estimation (in milliseconds) are reported in Table~\ref{tab:photo_result_f} on each scene of the dataset and, also, averaged over all in the last row. 
The proposed $\nabla$-RANSAC achieves the most accurate results on all but one scene, where it is the second best by a small margin. 
On average, it improves by ${\sim}3$\% compared to the second best method that is, interestingly, NG-RANSAC in this case.
We consider NG-RANSAC being the second best an interesting result given that it is based on simple MLPs with batch normalization. 
$\nabla$-RANSAC runs at comparable speed to its fastest, but less accurate alternatives.
It is far real time. 
The most efficient method is LMeDS achiving 11.61\% lower F1 score than $\nabla$-RANSAC.
Note that these timings do not include the inference time which is around $2$ milliseconds on average. 

\begin{table}[t]
    \centering
   
    \resizebox{0.95\columnwidth}{!}{\begin{tabular}{lccccc}
    \toprule 
    Method & AUC@$5^\circ$ & AUC@$10^\circ$ & AUC@$20^\circ$& Time (ms)\\ \midrule
    LMEDS\cite{rousseeuw2005robust} & 0.24 & 0.30 & 0.37 & \textbf{\phantom{1}27}\\
    RSC\cite{RANSAC}& 0.26 & 0.32 & 0.40 & \phantom{1}88\\
    GC-RSC\cite{barath2018graph} & 0.33 & 0.37 & 0.42 & 175\\
    MSC\cite{barath2019magsac} & 0.37 & 0.42 &0.47 & 239\\
    MSC++\cite{barath2020magsac++}& 0.37 & 0.42 & 0.47 & 113\\
    EAS \cite{fan2021efficient} & 0.24 & 0.28 &0.34 & 325\\
    OANet~\cite{oanet2019}&0.29 & 0.33 & 0.39 & \phantom{1}49\\
    CLNet~\cite{clnet2021}&0.34 &0.40 & 0.47 & \phantom{1}58\\
    NG-RSC~\cite{brachmann2019neural} &0.35&0.41& 0.47& \phantom{1}80\\
    $\nabla$-RANSAC & \textbf{0.41} & \textbf{0.45} & \textbf{0.50} & 122\\
    \bottomrule
    \end{tabular}}
    \caption{The average AUC scores of $\nabla$-RANSAC and comparison methods over 12 scenes from PhotoTourism, under different thresholds.  We are the most accurate method for \textbf{E} estimation
    }
    \label{tab:e}
\end{table}

Furthermore, in Fig.~\ref{fig:cdf}, we show the cumulative distribution functions (CDF) of the epipolar errors, run-times and F1 scores calculated from the 12000 tested image pairs.
In the two left plots, the epipolar error and the run-times CDFs are shown. Being close to the top-left corner indicates better performance.
The epipolar error curve of the proposed $\nabla$-RANSAC is above other methods on the entire plot.
While it is not the fastest, it finishes in $0.1$ seconds in $100$\% of the test cases.
This confirms that $\nabla$-RANSAC is applicable in time-sensitive applications. 
The right figure shows the CDFs of the F1-scores. 
In contrast to the other plots, here, being close to the bottom-right corner is preferred. 
$\nabla$-RANSAC has better score on the entire range.
It coincides with other methods only at the end of the range, where all methods are supposed to have $1$. 

The avg.\ results on all scenes when $\nabla$-RANSAC is trained with the 5PC, 7PC and norm.\ 8PC solvers are reported in Table~\ref{tab:F_solvers}.
For \textbf{F} estimation, the norm.\ 8PC solver leads to the most accurate results while also being the fastest.
The 5PC solver is close to it in terms of accuracy.

\subsection{Essential Matrix Estimation}
\label{sec:essential}

We evaluate our method for \textbf{E} estimation with the same train, validation, and test scenes as \textbf{F} estimation.
The correspondences are normalized by the intrinsic matrices. 
We use the differentiable 5PC algorithm when training on essential matrix estimation as it leads to more stable solutions than the 8PC solver in our experiments.
We train the method end-to-end for 10 epochs with weight initialization and gradient clipping.
The iteration number is fixed to 100.

To evaluate the estimated essential matrix, we decompose the \textbf{E} matrix to rotation and translation, calculate their errors $\epsilon_\mathbf{R}$, $\epsilon_\mathbf{t}$ and report the maximum $\max(\epsilon_\mathbf{R}, \epsilon_\mathbf{t})$.
We calculate the Area Under the Recall curve (AUC) thresholded at $5^\circ$, $10^\circ$ and $20^\circ$~\cite{brachmann2019neural}. 

The AUC scores, averaged over all scenes, are reported in Table~\ref{tab:e}.
The highest AUC scores, at all thresholds, are achieved by $\nabla$-RANSAC. For example, its AUC@$5^\circ$ score is higher than that of the second best methods (\ie, MAGSAC and MAGSAC++) by AUC $3$ points. 
It has comparable run-times to other less accurate alternatives. 

The results of minimal solvers are reported in Table~\ref{tab:E_solvers}.
The highest AUC scores, at all thresholds, are achieved by the five-point algorithm, confirming the necessity of use better minimal solvers instead of the 8PC algorithm.


\subsection{Ablation Studies}
\label{subsec:ablation}
\begin{figure}[t]
	\centering
 \begin{subfigure}[b]{0.495\columnwidth}
    \includegraphics[width =1\columnwidth, trim=0.2mm 0mm 0mm 0mm, clip]{Figures/technical_imA-eps-converted-to.pdf}
    \end{subfigure}
    \begin{subfigure}[b]{0.495\columnwidth}
    \includegraphics[width = 1\columnwidth, trim=0mm 0mm 0mm 0mm, clip]{Figures/Graph36A-eps-converted-to.pdf}
    \end{subfigure}
    \caption{\textbf{Ablation study}: AUC@10$^\circ$ score of \textbf{E} estimation on scene Sacre Coeur from~\cite{IMC2020}.
    \textit{Left:} $\nabla$-RANSAC trained with the 5PC and 8PC solvers with and without weight initialization.
    \textit{Right:} Score when using uniform sampler and Gumbel sampler with SNN ratios or model trained by minimized different loss functions.}
    \label{fig:technical}
\end{figure}

The effect of the weight initialization with the Kullback–Leibler divergence~\cite{van2014renyi} with different minimal solvers on essential matrix estimation is shown in the left plot of Fig.~\ref{fig:technical}.
For both minimal solvers, the initialization clearly improves the final accuracy.
Training with the 5PC solver without performing initialization leads to the same results as training with 8PC with initialization.
The best results are obtained by the 5PC solver with pre-trained weights. 
This again confirms the importance of using other solvers than the 8PC algorithm. 

The right plot of Fig.~\ref{fig:technical} shows different sampling weights.
For this plot, we only used scene Sacre Coeur as test set. 
Uniform sampler is basically the results of plain RANSAC.
Gumbel Sampler with SNN ratios requires no training.
It uses the pre-calculated SNN ratios to guide the sampling.
Finally, the effect of training with different loss functions is shown. 
On this dataset and problem, the epipolar error leads to the best results. 


\begin{table}
    \centering
    \resizebox{0.75\columnwidth}{!}{\begin{tabular}{l | ccc}
    \toprule
    \textbf{E} estimation &  5PC~\cite{nister2004efficient} & 7PC~\cite{hartley2003multiple} & 8PC~\cite{hartley2003multiple} \\ 
    \midrule 
         AUC@$5^\circ$  &\textbf{0.41}&0.37 &0.38\\
         AUC@$10^\circ$ & \textbf{0.45}& 0.42 & 0.43\\
         AUC@$20^\circ$ &\textbf{0.50}& 0.46 & 0.48\\
        \bottomrule
    \end{tabular}}
    \caption{
    Performance of \textbf{differentiable minimal solvers} on \textbf{E} estimation used for training $\nabla$-RANSAC.
    The AUC scores are reported on 12000 image pairs from PhotoTourism~\cite{IMC2020}.}
    \label{tab:E_solvers}
\end{table}

\begin{table}
    \centering
    \resizebox{0.85\columnwidth}{!}{\begin{tabular}{lccc}
    \toprule
    \textbf{F} estimation &  5PC~\cite{nister2004efficient} & 7PC~\cite{hartley2003multiple} & 8PC~\cite{hartley2003multiple} \\  \midrule 
         F1 score (\%) & 45.6 & 43.7 & \textbf{46.9}\\ 
         Med.\ epi.\ error (px) & \phantom{1}1.5 & \phantom{1}2.1 & \phantom{1}\textbf{1.0} \\
       
        \bottomrule
    \end{tabular}}
    \caption{Performance of \textbf{differentiable minimal solvers} on \textbf{F} estimation used for training $\nabla$-RANSAC.
    The AUC scores are reported on 12000 image pairs from PhotoTourism~\cite{IMC2020}.}
    \label{tab:F_solvers}
\end{table}

\subsection{Learning Feature Matching with $\nabla$-RANSAC}

An advantage of the fully differentiable RANSAC is that it allows the RANSAC pipeline to be optimized for test-time evaluation metrics such as pose error, inlier classification, and epipolar error (see~\cref{sec:quality}). Note that these metrics are differentiable in themselves. However, their input comes from the RANSAC pipeline, so if the RANSAC is non-differentiable, these metrics cannot be directly used as a loss function. Additionally, if the RANSAC is differentiable, the entire feature computation and matching module can be directly optimized on the above-mentioned metrics. In this section, we tune an end-to-end feature matching approach, LoFTR~\cite{Sun_2021_CVPR} on pose error, inlier classification, and epipolar error using the proposed $\nabla$-RANSAC.

Following the work of Patel {\em et al.}\cite{patel2020learning,patel2021feds} on learning surrogate loss functions, the LoFTR and the $\nabla$-RANSAC models are initialized with the pre-trained weights and post-tuned. To train LoFTR jointly with RANSAC and optimize the entire pipeline on the evaluation metric directly, we consider three training settings: LoFTR training, joint training, and alternating training. The {\em LoFTR training} setup refers to the setting where the pre-trained $\nabla$-RANSAC model is kept frozen, and only the LoFTR model is trained. The {\em joint training} setup refers to the setting where both the LoFTR and $\nabla$-RANSAC models are trained jointly. This setup involves adding a differentiable module on top of the LoFTR to obtain final predictions, and the loss is computed on these final predictions. The third setup of {\em alternating training} is similar to the training setup in \cite{patel2020learning} where at first the parameters of the LoFTR model are kept frozen, and only the $\nabla$-RANSAC is trained, and then only the LoFTR model is trained when $\nabla$-RANSAC is kept frozen. This alternating training is repeated for each epoch. In the last setup, the $\nabla$-RANSAC is part of the loss computation.

The evaluation is performed on the $12$ scenes from the PhotoTourism dataset, which contains $12$K image pairs in total. For the three training setups mentioned above, the results are presented in Table \ref{tab:loftr}. 
Note that for inference, the maximum number of RANSAC iterations is set to $5$K, while the actual number of iterations adapts dynamically based on the number of inliers in each iteration. The results show that the proposed $\nabla$-RANSAC can be used to improve end-to-end feature matching approaches and is robust to different training setups.
The best results were observed with the setup where only the LoFTR model is trained, and the $\nabla$-RANSAC is kept frozen with the pre-trained weights. This observation clearly shows the robustness of $\nabla$-RANSAC as a pre-trained model without any adaptation can be used in a plug-and-play manner.


\begin{table}[t]
    \centering
    \resizebox{0.99\columnwidth}{!}{\begin{tabular}{lccccc}
    \toprule 
    LoFTR ~\cite{Sun_2021_CVPR} & Robust Estimator & Training Setup & F1 score (\%) & run-time (ms) \\ \midrule
    pre-trained & RANSAC~\cite{RANSAC} & -  & 26.70 & \phantom{1}59.06 \\
    pre-trained & pre-trained $\nabla$-RSC & - & 29.38 & 138.67 \\
   
    post-tuned & pre-trained $\nabla$-RSC & LoFTR & \textbf{29.66} & \phantom{1}\textbf{58.85}\\
    post-tuned & post-tuned $\nabla$-RSC& joint & 29.44 & \phantom{1}63.87\\
    post-tuned & post-tuned $\nabla$-RSC  & alternating & 29.52 & \phantom{1}61.37\\
    \bottomrule
    \end{tabular}}
    \caption{The performance of training feature matcher with $\nabla$-RANSAC compared with the pre-trained models and non-differentiable robust estimator. The testing is performed on $12$K image pairs from PhotoTourism dataset with varying image sizes. All the setups, \eg, LoFTR training, joint training, and alternating training outperforms pre-trained LoFTR and $\nabla$-RANSAC models.
   
    }
    \label{tab:loftr}
\end{table}

\section{Conclusion}
\label{sec:conclusion}

In this paper, we propose the fully differentiable $\nabla$-RANSAC, \ie, the first end-to-end trainable randomized robust estimator. 
It predicts the inlier probabilities of the input data points, exploits the predictions in a guided sampler, and estimates the model parameters (e.g., fundamental matrix) and its quality while propagating the gradients through the entire procedure.    
$\nabla$-RANSAC leads to the most accurate fundamental and essential matrices compared to state-of-the-art robust estimators, on 12K image pairs from the PhotoTourism dataset.
Also, it is one of the fastest algorithms. 
To demonstrate its potential in unlocking the end-to-end training of geometric pipelines, we train $\nabla$-RANSAC together with a recent detector-free feature matcher, LoFTR \cite{Sun_2021_CVPR}, with which we achieve improved accuracy. 
We believe that the community will benefit from the algorithm.
We will make the code repository and the trained model public. 

{\small
\bibliographystyle{ieee_fullname}
