\section{Combinatorial methods}
\label{sec:combinatorial}

While the features obtained through neural network architectures based on functional maps are remarkably descriptive, the differentiability requirement on the matching algorithm often yields undesirable matches on disparate but geometrically similar regions of the shape.
Another line of work seeks to exploit the manifold structure and geometric coherence of feature matching through powerful combinatorial optimisation techniques.

\subsection{Background}
Combinatorial methods provide powerful tools for shape matching by formulating correspondence problems as discrete optimisation tasks over finite (but potentially very large) solution spaces.
As such, combinatorial methods are particularly well-suited for handling discrete shape representations such as meshes and point clouds~\cite{schrijver2003combinatorial,korte2011combinatorial} and usually seek an optimal assignment between discrete elements of shapes, such as vertex-to-vertex correspondences.


Formally, in the shape matching context, combinatorial methods aim to find an optimal solution $x^*\in\mathcal{F}\subset\{0,1\}^n$ 
(which usually resembles an assignment between surface elements) within a set of feasible solutions $\mathcal{F}$ while minimizing some (matching cost) function $\Energy(x) :\;\mathcal{F} \rightarrow \mathbb{R}$ such that 
\begin{equation} \label{eq:combinatial-general}
    x^*=\underset{x \in\mathcal{F}}{\min} \;\Energy(x).
\end{equation}
Consequently, the advantage of combinatorial formalisms is that they allow (i) to enforce that the solution is guaranteed to have certain properties by constraining the feasible set $\mathcal{F}$ adequately, and (ii) to add soft constraints and similarity measures by designing appropriate cost functions $\Energy(x)$.
In addition to the binary variables $x\in\mathcal{F}$ used in \eqref{eq:combinatial-general}, some methods also introduce additional continuous variables to model deformation alongside matching~\cite{bernard2020mina, gao_sigma_2023}.

Among the most prominent examples of combinatorial matching methods are the linear assignment problem (LAP) and the quadratic assignment problem (QAP).
The LAP was introduced to the shape matching community in~\cite{belongie2001matching}.
Yet, the LAP does not take into account relations between surface elements, such as their geodesic distances, and is therefore not suitable for shape matching when geometric relations must be preserved.
On the contrary, the QAP allows for incorporating such relations, but due to the quadratic objective function and resulting NP-hardness~\cite{rendl1994quadratic} it does not scale to practically relevant resolutions.
Thus, shape matching approaches based on QAP-formalisms rely on finding solutions via approximations, including convex relaxations~\cite{schellewald2005probabilistic,dym2017relax,bernard2018ds,kushinsky2019sinkhorn,droge_kissing_2023}, and heuristics~\cite{dubrovina2011approximately, le2017alternating,vestner2017efficient, haller2022comparative}.
To obtain more scalable solutions, other combinatorial shape matching methods exploit the discrete representation of 3D shapes and incorporate neighborhood relations using linear constraints to restrict the feasible set.
Among these, one of the most prominent examples is the work by Windheuser et al.~\cite{windheuser2011geometrically} which sparked many follow-up approaches building on similar concepts. 

In the following, we summarise recent combinatorial methods. 
We further categorise approaches into landmark-based formalism, surface-to-surface formalisms, multi-shape matching approaches as well as methods based upon quantum computing. 
We summarise all methods in \cref{tab:combinatorial}.


\begin{table}[t]
    \centering
    \setlength{\tabcolsep}{3pt}
    \begin{tabular}{l@{}cccccc}
    \toprule
    \small
     \textbf{Method} & \textbf{LB} & \textbf{S2S} & \textbf{MS} & \textbf{Q} & \textbf{LO} & \textbf{RF} \\ 
     \midrule
     Maron et al.~\cite{maron2016point} & \cmark & \cmark & \xmark & \xmark & \xmark & \cmark \\
     MINA~\cite{bernard2020mina} & \cmark & \xmark & \xmark & \xmark & \xmark & \xmark \\
     SIGMA~\cite{gao_sigma_2023} & \cmark & \xmark & \xmark & \xmark & \xmark & \xmark  \\
     Ren et al.~\cite{ren2021discrete} & \xmark & \cmark & \xmark & \xmark & \xmark & \xmark  \\
     Pai et al.~\cite{pai2021fast} & \xmark & \cmark & \xmark & \xmark & \xmark & \xmark \\
     SM-Comb~\cite{roetzer_scalable_2022} & \xmark & \cmark & \xmark & \xmark &\cmark & \xmark  \\
     DiscoMatch~\cite{roetzer2024discomatch} & \xmark & \cmark & \xmark & \xmark &\cmark & \xmark  \\
     SpiderMatch~\cite{roetzer_spidermatch_2024} & \xmark & \cmark & \xmark & \xmark &\cmark & \xmark \\
     Roetzer et al.~\cite{roetzer2025higherorder} & \xmark & \cmark & \xmark & \xmark &\cmark & \xmark \\
     GeCo~\cite{roetzer2025geco} & \xmark & \cmark & \xmark & \xmark &\cmark & \xmark \\
     SuPaMatch~\cite{amrani2025highres} & \xmark & \cmark & \xmark & \xmark &\cmark & \xmark \\
     Gao et al.~\cite{gao2021isometric} & \xmark & \xmark & \cmark & \xmark & \xmark & \xmark \\
     Kahl et al.~\cite{kahl2025towards} & \xmark & \xmark  & \cmark & \xmark & \xmark & \xmark \\
     Bhatia et al.~\cite{bhatia_ccuantumm_2023} & \xmark & \xmark & \cmark & \cmark  & \xmark & \xmark \\
     Benkner et al.~\cite{benkner2020adiabatic} & \xmark & \xmark & \xmark & \cmark & \xmark & \xmark \\
     Seelbach et al.~\cite{seelbach2021qmatch} & \xmark & \xmark & \xmark & \cmark & \xmark & \xmark \\
     Meli et al.~\cite{meli_qucoop_2025} & \xmark & \xmark & \xmark & \cmark & \xmark & \xmark  \\
     \bottomrule
    \end{tabular}
    \caption{Overview over \textbf{combinatorial} matching methods and classification whether methods are \textbf{(LB)} landmark-based, are \textbf{(S2S)} Surface-to-Surface methods, are applicable to \textbf{MS} multiple shapes, are \textbf{(Q)} Quantum, have a \textbf{(LO)} linear objective $f(x)$, or have a \textbf{(RF)} relax the feasible set $\mathcal{F}$.
    }
    \label{tab:combinatorial}
\end{table}







\subsection{Landmark-based approaches}
Landmark-based approaches, sometimes also referred to as \emph{sparse} shape matching approaches, focus on a small subset of vertices on each shape, 
denoted by $\Tilde{X}_\mathcal{M}\subset X_\mathcal{M}$ and $\Tilde{X}_\mathcal{N}\subset X_\mathcal{N}$, see \cref{fig:key-points}. 
These vertices, often referred to as \emph{landmarks} or sometimes also as \emph{keypoints}, are intended to capture 
salient or distinctive locations on the respective shapes $\mathcal{M}$ and $\mathcal{N}$.
Generally, such landmarks are not available and are often defined by human annotators.
Landmark-based formalisms usually aim to find correspondences only among the landmarks with the underlying assumption that correspondences between shapes can be inferred or approximated effectively from the relationships between their landmarks.
Consequently, using landmarks often allows for efficient solvability while still enabling the recovery of dense correspondences via post-processing.
Yet, landmark-based approaches assume consistent landmark positions across shapes and thus the solution quality of these methods is heavily based on their position.
\begin{figure}
    \centering
    \begin{tabular}{cc}
        \small
         \includegraphics[height=3cm]{figures/paul/keypoints/shrec20_bison_elephant_a_M}&
         \includegraphics[height=3cm]{figures/paul/keypoints/shrec20_bison_elephant_a_N}
    \end{tabular} 
    \caption{Visualisation of correspondence between \textbf{landmarks} (in colour) lying on the surface (grey) of two 3D shapes. Image source: \cite{gao_sigma_2023}}
    \label{fig:key-points}
\end{figure}



Among the landmark-based approaches, one line of work seeks to uncover synergies between shape correspondences and shape reconstruction.
\def\assignmentMatrix{\Pi}
In~\cite{bernard2020mina}, Bernard et al.~propose a mixed integer programming framework which optimises for a $m\times m$ permutation matrix $\assignmentMatrix$ between landmarks $\Tilde{X}_\mathcal{M}\subset X_\mathcal{M}$, $\Tilde{X}_\mathcal{N}\subset X_\mathcal{N}$, and  a set of reconstructed landmark vertex-positions $\Tilde{X}_\mathcal{M}^\text{rec}$.
While $\Tilde{X}_\mathcal{M}$ represents the initial positions of the selected landmarks on shape $\mathcal{M}$, $\Tilde{X}_\mathcal{M}^\text{rec}$ denotes their optimized coordinates in the target space that best align with the landmarks $\Tilde{X}_\mathcal{M}$.
On a high-level the goal is to minimise
\begin{equation}
    \begin{aligned}
    \underset{
        \assignmentMatrix \in \mathbb{P}, \;
        \Tilde{X}_\mathcal{M}^\text{rec}\in \mathbb{R}^{m\times r}
    }
    {\min} 
    ||\Tilde{X}_\mathcal{M}^\text{rec} - 
    \assignmentMatrix 
    \Tilde{X}_\mathcal{N}||  + 
    \Energy_\text{rec}(\Tilde{X}_\mathcal{M}^\text{rec}, \mathcal{M})
    \end{aligned}
\end{equation}
Here, $\Energy_\text{rec}(\Tilde{X}_\mathcal{M}^\text{rec}, \mathcal{M})$ is a reconstruction function that aims for a smooth reconstruction of shape $\mathcal{M}$ using the optimised landmark positions $\Tilde{X}_\mathcal{M}^\text{rec}$.

Building on the ideas presented in~\cite{bernard2020mina},
Gao et al.~\cite{gao_sigma_2023} propose a similar mixed integer programming formalism
and introduce the 
\emph{projected Laplace-Beltrami} (LBO) operator for the reconstruction, which reads
\begin{equation}
    L_\mathcal{M}^\text{proj} \coloneqq {M}^T \Delta^\text{stiff}_\mathcal{M} {M}.
\end{equation}
Here, $\Delta^\text{stiff}_\mathcal{M}$ is the stiffness matrix component of the Laplace-Beltrami operator $L_\mathcal{M}$,
${M}\coloneqq I -
\Tilde{X}_\mathcal{M}
(
\Tilde{X}_\mathcal{M}^T  
\Tilde{X}_\mathcal{M}
)^{-1} \Tilde{X}_\mathcal{M}^T$
is a projection matrix
and $ \Tilde{X}_\mathcal{M}\coloneqq( {X}_\mathcal{M}, \boldsymbol{1})$.
As shown in \cite{gao_sigma_2023}, the projected LBO $L_\mathcal{M}^\text{proj}$ is scale-invariant and invariant under rigid body transformations.
Furthermore, it retains more high-frequency geometric details of the shapes compared to the Laplace-Beltrami operator.
This helps to increase reconstruction quality and, in turn, leads to better matching results.

While landmark-based approaches offer improved scalability due to smaller optimisation problems and, with that, greater mathematical modelling flexibility, many argue that full surface-to-surface maps must be considered to compute faithful matchings between shapes.

\subsection{Surface-to-surface approaches}
Surface-to-surface means finding a dense matching for all vertices $X_\mathcal{M}$ of a shape $\mathcal{M}$, rather than only for a subset of vertices, as in landmark-based approaches.
To find such a matching, some lines of work build on the \emph{Procrustes matching} problem where the goal is to simultaneously find correspondences between two sets of points $X_\mathcal{M}, X_\mathcal{N} \in \mathbb{R}^{n\times d}$ (in form of a assignment matrix $\Pi_{\mathcal{M}\mathcal{N}}$), as well as a rigid alignment $T \in SE(d)$. 
This is formulated as the following optimization problem
\begin{equation}
    \underset{
        \Pi_{\mathcal{M}\mathcal{N}} \in\mathbb{P}, T \in SE(d),
                  }{\min} ||\bar{X}_\mathcal{M} T - \Pi_{\mathcal{M}\mathcal{N}} \bar{X}_\mathcal{N} ||_{F}^2,
\end{equation}
with $\bar{X}_\mathcal{M} := [X_\mathcal{M}, \mathds{1}_{n}] \in \mathbb{R}^{n \times (d + 1)}$ denoting the homogeneous coordinates of shape $\mathcal{M}$, where $\mathds{1}_{n} \in \mathbb{R}^{n}$ is the all-ones vector. 
The homogeneous coordinates $\bar{X}_\mathcal{N}$ are defined analogously.
An iterative approximation to this formulation was popularised by the seminal iterative-closest point (ICP) algorithm~\cite{besl1992method}, see also \cref{sec:deformation}.
Other works~\cite{maron2016point} propose convex relaxations for the Procrustes matching problem.

Ren et al.~\cite{ren2021discrete} 
propose a discrete optimisation approach which computes surface maps between two shapes $\mathcal{M}$ and $\mathcal{N}$ by constraining a functional map $C_{\mathcal{M}\mathcal{N}}$ such that it can be related to a point-wise map $\Pi_{\mathcal{M}\mathcal{N}}$.
These are then also called \emph{proper functional maps}, see~\cref{eq:fmap-to-pmap}.
As such, proper functional maps have been further explored 
in \cite{pai2021fast} using convex relaxations of point-wise maps, as well as in
deep shape matching methods as a loss function~\cite{cao_unsupervised_2023, attaiki_understanding_2023}.


In contrast, Windheuser et al.~\cite{windheuser2011geometrically} enforce consistent neighbourhood relations of surface elements by utilising linear constraints in an integer linear programming formalism based upon a so-called \emph{product space}.
While mathematically elegant, this formalism is hard to solve in practice and thus Windheuser et al.~\cite{windheuser2011geometrically} propose to solve it using linear programming (LP) relaxations and iterative variable fixations. 
Yet, even with this approximate scheme, it only scales to shapes with around $50-100$ triangles on modern hardware and with modern LP solvers; see~\cite{roetzer2025geco} for runtime comparisons.
This inspired follow-up works~\cite{roetzer_scalable_2022, roetzer2024discomatch} to develop approximate solvers for this formalism.
Both solvers~\cite{roetzer_scalable_2022, roetzer2024discomatch}
build on dual decomposition schemes presented in~\cite{lange2021efficient,abbas2022fastdog} and run on CPU~\cite{roetzer_scalable_2022} and GPU~\cite{roetzer2024discomatch} respectively.
These solvers show significant improvements in scalability, i.e.~\cite{roetzer_scalable_2022} scales to shapes with up to a $1000$ triangles and \cite{roetzer2024discomatch} to $750$ triangles.
Yet, due to runtime and memory limitations, scalability is nevertheless limited. 

This motivated follow-up works to develop new formalisms based on \emph{product graphs}. Product graphs are denoted 
\begin{equation}
    \mathcal{P}=\mathcal{G}_1\boxtimes\mathcal{G}_2,
\end{equation}
where the \emph{strong product} 
between two graphs $\mathcal{G}_1$ and $\mathcal{G}_2$, see~\cite{hammack2011handbook} for more details and \cref{fig:prodgraph} for a visualization.
\begin{figure}
    \centering
    \includegraphics[width=0.5\linewidth]{figures/paul/productgraph_new.drawio.pdf}
    \caption{Visualisation of the (in-layers-structured) \textbf{product graph} $\mathcal{P}$ between a cyclic chain graph $\mathcal{C}$ and the graph of a triangle mesh of a 3D shape $\mathcal{N}$. Correspondences between $\mathcal{C}$ and $\mathcal{N}$ can be computed by finding the shortest cyclic path that visits all layers. Image source: \cite{roetzer_spidermatch_2024}.}
    \label{fig:prodgraph}
\end{figure}
In the past, these product graphs were successfully used for dynamic time-warping~\cite{sakoe1978dtw},
i.e.~finding point-to-point correspondence between two contours with start and end,
matching a shape to an image~\cite{schoenemann2009combinatorial}, 2D-to-2D shape matching~\cite{schmidt2009planar}, and 2D-to-3D shape matching~\cite{lahner2016efficient}.
To obtain an efficient-to-solve formalism for 3D-3D shape matching, SpiderMatch~\cite{roetzer_spidermatch_2024} changes the representation of shape $\mathcal{M}$, i.e.~the \emph{source} shape:
it uses a single cyclic, self-intersecting curve $\mathcal{C}$ on the surface of the shape $\mathcal{M}$ to represent the respective shape.
This allows casting the 3D shape matching problem as finding a minimum cost cyclic path in a product graph $\mathcal{P}$ with additional constraints.
As shown in~\cite{roetzer_spidermatch_2024}, this problem is globally optimal and efficiently solvable in practice, scaling to shape resolutions of up to $3500$ triangles.

In~\cite{roetzer2025higherorder}, authors utilise the curve representation proposed in SpiderMatch~\cite{roetzer_spidermatch_2024}, drop the additional constraints, and show that good shape matching results can be computed when using minimum ratio cycles with local rigidity regularisers~\cite{sorkine2007rigid,roetzer_conjugate_2023}, achieving global optimality in practice.

Contrary to these methods, \cite{roetzer2025geco} argues that representing a 3D shape with a single cyclic curve $\mathcal{C}$ and matching this curve to the target shape $\mathcal{N}$ does not lead to a faithful surface-to-surface matching.
Therefore, \cite{roetzer2025geco} proposes to use multiple cyclic curves $\mathcal{C}_1,\dots,\mathcal{C}_n$ to represent each
triangle of the source shape $\mathcal{M}$.
This allows formulating 3D shape matching via multiple product graphs $\mathcal{P}_1,\dots,\mathcal{P}_n$ as finding multiple shortest cyclic paths which have to adhere to additional coupling constraints. 
In turn, this leads to a faithful surface-to-surface matching that is globally optimal in practice and scales to shapes with up to $3000$ triangles.
Building on this formalism, \cite{amrani2025highres} propose a two-stage approach in which they first transfer patches from source to target shape 
and subsequently find dense surface-to-surface matchings for each pair of corresponding patches.
This allows scaling to shapes with up to $10$k triangles without the need for coarse-to-fine strategies.
This has the advantage of fewer discretisation artifacts which usually occur in approaches relying on coarse-to-fine strategies.

While \cite{roetzer2025geco} has shown that multiple coupled problems lead to high-quality matching results, these coupled problems only involve two 3D shapes.
As a result, their formalism is not directly applicable to problems that involve collections of shapes, which we discuss next.


\subsection{Multi-shape matching approaches}
Multi-shape matching approaches simultaneously optimise maps between a collection of shapes.
Assuming bijective mappings between shapes, \emph{cycle consistency} emerges as a fundamental property that should be fulfilled in most cases~\cite{schmidt2007intrinsic,huang2013consistent,cosmo2017consistent}.
This means that the matching from a vertex on one shape via multiple other shapes back to the first shape must end in the same vertex again.
For non-bijective mappings between shapes, the definition of cycle consistency is more intricate~\cite{bernard2019synchronisation}.
In particular, considering the bijective case, for a collection of $N$-many shapes $\mathcal{M}_1,\dots,\mathcal{M}_N$, it must hold that
\begin{equation} \label{eq:cycle-consistency}
    \forall i,j,k \in \{1,\dots,N\}:\; \Pi_{\mathcal{M}_i\mathcal{M}_j} \Pi_{\mathcal{M}_j\mathcal{M}_k} = \Pi_{\mathcal{M}_i\mathcal{M}_k}.
\end{equation}
\cite{bhatia_ccuantumm_2023} enforce cycle consistency constraints by considering triplets of shapes and by proposing a quantum-hybrid approach building on ideas of the famous alpha-expansion algorithm~\cite{boykov2002fast}.
\cite{eisenberger2023g} learns maps between multiple shapes, which are used as edge weights for all pairs $i,j \in \{1,\dots,N\}$ in the shape graph, to then propagate maps along shortest paths in the shape graph to enforce cycle consistency.
\cite{kahl2025towards} propose an approximate algorithm for cycle-consistent multi-graph matching, which exploits parallelism, performs local graph matching and maintains cycle consistency at all times.


To avoid enforcing the cycle consistency constraints \eqref{eq:cycle-consistency} explicitly, other works revert to shape-to-universe matching where cycle consistency is enforced implicitly~\cite{pachauri2013solving}.
Among these, \cite{gao2021isometric} introduces a multi-shape matching formalism for isometric deformed shapes and further proposes an algorithm based upon convex relaxations.
Later, the shape-to-universe matching scheme was also employed in unsupervised learning approaches~\cite{cao2022unsupervised}, see also \cref{tab:spectral_methods}.

Multi-shape matching approaches and many other combinatorial methods are most often hard to solve and as a consequence, on classical hardware, their scalability is limited. 
Motivated by that, other lines of work aim for formalisms which can be employed on quantum hardware
which can find global optima naturally in some configurations, 
such that formalisms might become much more scalable with advances in quantum computing.

\subsection{Quantum annealing}

Quantum annealers can efficiently solve quadratic unconstrained binary optimisation problems in the form of
\begin{equation}
 \min_{x \in \{0,1\}^n} x^\top Q x + q^\top x. \label{eq:qubo}
\end{equation}
Here $Q\in\mathbb{R}^{n\times n}$ is a symmetric (possibly dense) cost matrix containing pairwise costs and $q \in \mathbb{R}^{n}$ a linear bias. 
Current available systems, for example, the D-Wave line~\cite{dwave_website}, contain several thousand qubits (\emph{quantum bits}).
In theory, the number of qubits corresponds directly to $n$ in \eqref{eq:qubo}, however, due to hardware limitations in the connectivity of the annealer, this only holds for certain forms of $Q$ and in reality, many qubits are blocked when representing a dense $Q$ matrix~\cite{seelbach2025compensating}. 
The natural use of binary variables makes quantum annealers primed for solving combinatorial versions like \cref{eq:combinatial-general}, but the unconstrained setting cannot model many of the desired properties of the solution, such as bijectivity.
The first approach encoded the permutation constraint in $Q$, but this occupies most of the qubits on fulfilling these constraints and could only correspondence between up to 4 points~\cite{benkner2020adiabatic}. 
However, by solving the problem iteratively in subproblems based on a version of $\alpha$-expansion, quantum annealers can be scaled to several thousand vertices~\cite{seelbach2021qmatch} and even solve the hard multi-shape case with a linear runtime scaling in the number of shapes~\cite{bhatia_ccuantumm_2023}.
The permutation set can also be parametrised by the binary variables directly, which leads to an efficient approximation problem without the need for an iterative approach~\cite{meli_qucoop_2025}.
For an introduction to quantum annealing in the context of computer vision and more details on correspondence estimation approaches relying on quantum computational paradigms, interested readers can refer to the recent survey by Kuete Meli et al.~\cite{meli2025quantum}. 
















