# An Exact Hypergraph Matching Algorithm for Nuclear Identification in
Embryonic Caenorhabditis elegans

Andrew Lauziere
Department of Mathematics
University of Maryland, College Park
College Park, MD 20742
<EMAIL_ADDRESS>
&Ryan Christensen
Laboratory of High Resolution Optical Imaging
National Institutes of Health
Bethesda, MD 20892
<EMAIL_ADDRESS>
&Hari Shroff
Laboratory of High Resolution Optical Imaging
National Institutes of Health
Bethesda, MD 20892
<EMAIL_ADDRESS>
&Radu Balan
Center for Scientific Computation and Mathematical Modeling
University of Maryland, College Park
College Park, MD 20742
<EMAIL_ADDRESS>

###### Abstract

Point-set matching is a common task in computer vision. Many matching
applications feature affine point transformations which can be adequately
modeled via lower order objective functions [1, 2]. However, point-set
matching tasks may require richer detail in order to characterize underlying
relationships. Hypergraphs, an extension of traditional graphs, have emerged
to more intricately model relationships between points. Existing
hypergraphical point-set matching methods are limited to heuristic algorithms
which do not easily scale to handle higher degree hypergraphs [3, 4, 5]. Our
proposed algorithm, Exact Hypergraph Matching (EHGM), adapts the classical
branch-and-bound paradigm to dynamically identify a globally optimal
correspondence between point-sets under an arbitrarily intricate
hypergraphical model. EHGM is applicable to conservatively sized ($n\leq 20$)
point-set matching problems in which relationships between points require
increased context to adequately characterize. The methodology is motivated by
Caenorhabditis elegans, a model organism used frequently in developmental
biology and neurobiology [6, 7, 8, 9]. The C. elegans embryo can be used for
cell tracking studies to understand how cell movement drives the development
of specific embryonic tissues. However, twitching due to muscular activity in
late-stage embryos invalidates traditional cell tracking approaches. The
embryo possesses a small set of cells which together act as fiducial markers
to approximate the coiled embryo’s posture, serving as a frame of reference to
track cells of various tissues during late-stage embryogenesis [10]. Current
approaches to posture identification rely on time-consuming manual efforts by
trained users which limits the efficiency of subsequent cell tracking. EHGM
with biologically inspired hypergraphical models identifies posture more
accurately than established point-set matching methods, correctly identifying
twice as many sampled postures as a heuristic graphical approach.

## Introduction

Point-set matching describes the task of finding an alignment between two sets
of points. The problem appears in computer vision applications such as point-
set registration [11], object recognition [12], and multiple object tracking
[13]. Often the point-sets are modeled via graphs, abstract mathematical
objects in which points are represented as vertices and edges define
relationships between pairs of vertices.

User defined attributes characterize the vertices and edges, such as
coordinate positions or shape descriptions and lengths of chords connecting
vertices, respectively. Specified attributes give insight to observable
relationships between vertices and allow for structural analyses of graphs.
Graph matching is the optimization problem defined by the search for a
correspondence of vertices between a pair of attributed graphs. The
optimization problem uses binary variables $x_{ij}$ to specify a matching
between vertex i in the first graph to vertex j of the second. The graph
matching domain consists of assignment matrices of size $n_{1}\times n_{2}$,
for matching graphs of size $n_{1}$ and $n_{2}$.

$\mathcal{X}=\\{X\in\\{0,1\\}^{n_{1}\times n_{2}}:\forall
j,\sum_{i=1}^{n_{1}}x_{ij}\leq 1,\forall i\sum_{j=1}^{n_{2}}x_{ij}=1\\}$ (1)

The space $\mathcal{X}$ (Eq 1) comprises assignment matrices which each
describe a one-to-one alignment between nodes of the two graphs. The
specification of the graph matching optimization objective function allows for
joint assignment costs: i.e., how the assignment of a pair of vertex-to-vertex
assignments changes the quality of the match. Let C be an $n_{1}\times n_{2}$
matrix and D be a $n_{1}\times n_{2}\times n_{1}\times n_{2}$ tensor storing
the vertex-to-vertex and edge-to-edge dissimilarities, respectively. The graph
matching optimization problem is expressed in Eq 2, which takes the form of
the quadratic assignment problem (QAP).

$\displaystyle\underset{X\in\mathcal{X}}{\text{minimize}}$
$\displaystyle\sum_{i=1}^{n_{1}}\sum_{j=1}^{n_{2}}\sum_{k=1}^{n_{1}}\sum_{l=1}^{n_{2}}d_{ijkl}x_{ij}x_{kl}+\sum_{i=1}^{n_{1}}\sum_{j=1}^{n_{2}}c_{ij}x_{ij}$
(2)

Graphs are limited in their expressive power as edges can only relate pairs of
vertices; hypergraphs extend the definition of a graph to include hyperedges
which can specify relationships among an arbitrary number of vertices.
Hypergraph matching then concerns finding an optimal vertex correspondence
between pairs of attributed hypergraphs. The number of vertices aligned by the
most comprehensive hyperedge defines the degree of a hypergraph.

Maximum degree hypergraphs with hyperedges composed of all $n_{1}$ vertices
yield the most comprehensive point-set matching function possible. The
optimization objective function captures the dissimilarity arising between the
matching:
$(l_{1},l_{2},\dots,l_{n_{1}})\mapsto(l^{\prime}_{1},l^{\prime}_{2},l^{\prime}_{3},\dots,l^{\prime}_{n_{1}})$.
Then, for a given assignment matrix $X\in\mathcal{X}$, the hypergraph matching
objective can be expressed using $n_{1}$ dissimilarity tensors of dimension
$2,4,\dots,2d,\dots,2n_{1}$, each measuring dissimilarity between degree $d$
hyperedges, respectively. Define $\mathbf{Z}^{(d)}$ as the tensor mapping the
dissimilarity for the degree d hyperedges. The hypergraph matching objective
is expressed in Eq 7.

$\mathit{f}(X|\mathbf{Z}^{(1)},\mathbf{Z}^{(2)},\dots,\mathbf{Z}^{(n_{1})})=\sum_{l_{1}=1}^{n_{1}}\sum_{l^{\prime}_{1}=1}^{n_{2}}\mathbf{Z}^{(1)}_{l_{1}l^{\prime}_{1}}x_{l_{1}l^{\prime}_{1}}+\sum_{l_{1}=1}^{n_{1}}\sum_{l^{\prime}_{1}=1}^{n_{2}}\sum_{l_{2}=l_{1}+1}^{n_{1}}\sum_{l^{\prime}_{2}=1}^{n_{2}}\mathbf{Z}^{(2)}_{l_{1}l^{\prime}_{1}l_{2}l^{\prime}_{2}}x_{l_{1}l^{\prime}_{1}}x_{l_{2}l^{\prime}_{2}}\\\
+\sum_{l_{1}=1}^{n_{1}}\sum_{l^{\prime}_{1}=1}^{n_{2}}\sum_{l_{2}=l_{1}+1}^{n_{1}}\sum_{l^{\prime}_{2}=1}^{n_{2}}\sum_{l_{3}=l_{2}+1}^{n_{1}}\sum_{l^{\prime}_{3}=1}^{n_{2}}\mathbf{Z}^{(3)}_{l_{1}l^{\prime}_{1}l_{2}l^{\prime}_{2}l_{3}l^{\prime}_{3}}x_{l_{1}l^{\prime}_{1}}x_{l_{2}l^{\prime}_{2}}x_{l_{3}l^{\prime}_{3}}+...\\\
+\sum_{l_{1}=1}^{n_{1}}\sum_{l^{\prime}_{1}=1}^{n_{2}}...\sum_{l_{n_{1}}=l_{n_{1}-1}+1}^{n_{1}}\sum_{l^{\prime}_{n_{1}}=1}^{n_{2}}\mathbf{Z}^{(n_{1})}_{l_{1}l^{\prime}_{1}\dots
l_{n_{1}}l^{\prime}_{n_{1}}}x_{l_{1}l^{\prime}_{1}}\dots
x_{l_{n_{1}}l^{\prime}_{n_{1}}}$ (3)

Hypergraph matching allows for the modeling of intricate point-set matching
problems through high multiplicity assignment objective function formulations.
The $\mathbf{Z}^{(d)}$ dissimilarity terms measure degree d hyperedge
dissimilarity comprising d simultaneous vertex assignments. The range in
assignment problem objective complexity from d=1 to d=$n_{1}$ trades off model
capacity for increased computation. The traditional linear assignment problem
(d=1) is solvable in polynomial time [14], but treats points between sets
independently. Existing graphical methods (d=2) and hypergraphical methods
(d$>$2) rely on approximate searches and do not generalize to high degree
formulations of Eq 7. Exact Hypergraph Matching (EHGM) is able to find
globally optimal solutions to hypergraph matching problems of arbitrary
degree, allowing for the modeling of intricate point-set matching tasks.

### Related Research

Finding an exact solution to the QAP is an $\mathcal{NP}$-hard problem. That
is, unless P=NP, there does not exist a polynomial time solution to exactly
solve the QAP [15]. Higher order assignment problems (i.e. hypergraph
matching) are also $\mathcal{NP}$-hard as they are at least as hard as the QAP
[16]. As a result, recent methods for graph matching and lower-degree
hypergraph matching focus on heuristic solutions which offer no guarantee on
performance [1, 2, 3, 5, 4]. Heuristic hypergraph matching methods are adapted
from existing graph matching algorithms. In particular, spectral methods for
solving graph matching (Eq 2) have been extended to solve hypergraph matching.
Duchenne et al. [3] adapt Leordeanu’s [11] work to obtain a rank-1
approximation of the affinity tensor via higher order power iteration.
However, calculating affinity tensors ($\mathbf{Z}^{(d)}$ terms) is
computationally prohibitive, especially for higher degree hypergraphs due to
the exponentially growing number of entries in the tensors. Simplifying
assumptions such as super-symmetry and sparseness are used with sampling
methods to build large affinity tensors [3, 17]. Chertok and Keller propose
similar methodology to [3], but instead unfold the affinity tensor and use the
leading left singular vector to approximate the adjacency matrix [4]. All such
methods operate outside the permutation matrix space. The Hungarian algorithm
or similar binarization step is used to yield a valid assignment, e.g. as in
[11].

Exactness allows for a more rigorous analysis of a hypergraphical point-set
matching model than is possible using heuristic techniques. The guarantee of a
globally optimal correspondence allows an iterative tuning of the underlying
model in pursuit of accurate characterization, whereas the output of a
heuristic algorithm could be incorrect due either to the stochasticity of the
search or to inadequacy of the optimization objective. Branch-and-bound is a
paradigm originally developed to exactly solve the the travelling salesman
problem, a type of QAP [18, 19]. Branch-and-bound methods recursively commit
partial assignments and solve successive subproblems within $\mathcal{X}$. The
paradigm iteratively partitions the search space while bounding the optimum at
each branch. At each step the method prunes branches which cannot contain lead
to the optimum. Convergence occurs when only feasible assignments achieving a
global optimum remain. The $\mathcal{NP}-$hardness of the QAP implies
convergence occurs only after implicit enumeration of $\mathcal{X}$.

###  0.1 Overview of EHGM & Application to C. elegans

EHGM deviates from recent graph matching and hypergraph matching methodology
as an exact method, guaranteeing convergence to a globally optimal solution
(S1:Convergence of EHGM). Heuristic hypergraph matching methods approximate
the assignment matrix using the dissimilarity tensor [3, 4] whereas EHGM
builds upon the seminal branch-and-bound algorithm [18]. EHGM extends the
methodology to branch and prune based upon a given hypergraphical model. A
k-tuple of nodes at branch m are greedily selected while another step
encapsulates the full hypergraphical objective upon selection. These changes
enable flexibility in altering the hypergraph matching objective, particularly
in allowing for high degree hypergraphical modeling.

EHGM is applied to model posture in embryonic Caenorhabditis elegans (C.
elegans), a small, free-living roundworm. The nematode features approximately
550 cells upon hatching, including a set of twenty seam cells and two
associated neuroblasts. The seam cells and neuroblasts form in lateral pairs
along the left and right sides of the worm, resulting in eleven pairs upon
hatching [7]. The neuroblasts appear in the final hours of development, just
prior to hatching. The pairs of cells are named, posterior to anterior: T, V6,
V5, Q (neuroblasts), V4, V3, V2, V1, H2, H1, and H0. Each pair’s left and
right cell is named accordingly; for example, H1L and H1R comprise the H1
pair. Fig 1-A depicts center points of seam cell nuclei located in an example
image volume as imaged in the eggshell (left) and straightened to reveal the
bilateral symmetry in seam cell locations (right). Fig 1-B shows four
sequential images of an embryo, five minutes between images.

Figure 1: High spatial resolution, low temporal resolution imaging
necessitates posture identification. A: Manually identified and seam cell
nuclei from an imaged C. elegans embryo. The cells form in pairs; they are
labelled posterior to anterior: T, V6, …, H0. The identification of all seam
cells reveals the embryo’s posture. Natural cubic splines through the left and
right-side seam cells estimate the coiled body. The left image depicts
identified nuclei connected to outline the embryonic worm. The fit splines are
used to untwist the worm, generating the remapped straightened points in the
diagram on the right. B: Labelled nuclear coordinates from a sequence of four
images. The embryo repositions in the five minute intervals between images,
causing failure of traditional tracking approaches.

We define posture as the identification of all seam cells and neuroblasts,
which together reveal the shape of the coiled embryo. Posture identification
allows for traditional frame-to-frame tracking of imaged cells belonging to
various tissues such as the gut, nerve ring, and bands of muscle [10]. Images
are captured in five minute intervals (Fig 1-B) in order to achieve necessary
resolution to track cells of other tissues without disturbing embryo
development. Fig 2-A highlights muscle cell nuclei (red dots) with the
identified seam cells to contextualize the embryo’s positioning. The posture
is used to remap the muscle cells such that traditional cell tracking
approaches can be applied in the late-stage embryo (Fig 2-B). Fig 2-C depicts
the cell remapping process [10]. The muscle cells are remapped according to
splines fitted to the posture. The untwisted cell positions are then tracked
frame-to-frame (Fig 2-D).

Figure 2: Posture identification allows the tracking of other cells during
late-stage embrygenesis. A: Seam cell nuclei coordinates (black) and muscle
nuclei coordinates (red) in a sequence of three sequential volumetric images.
The untwisting process (green arrows) uses the seam cells to remap muscle
coordinates to a common frame of reference. B: The remapped muscle nuclei are
tracked frame-to-frame (blue arrows). C: A higher magnification view from the
right coordinate plot of A. The left, right, and midpoint splines are used to
create a change of basis defined by the tangent (black), normal (blue), and
binormal (cyan) vectors. Ellipses are inscribed along the tangent of the
midpoint spline, approximating the skin of the coiled embryo. D: A portion of
the left (red) and center (blue) remapped muscle coordinates. Black lines
connect the coordinates, frame-to-frame.

Current methods for posture identification rely on trained users to manually
annotate the imaged nuclei using a 3D rendering tool [20]. The process takes
several minutes per image volume and must be performed on approximately 100
image volumes per embryo [10]. Manual annotation strategies motivated us to
develop EHGM, as established methods for point-set matching fail to adequately
capture the relationships between seam cells throughout myriad twists and
deformations of the developing embryo. Fig 3 depicts manually identified
postures in the first two successive image volumes of Fig 1-B. Manual
identification is performed in Medical Imaging, Processing, Analysis and
Visualization (MIPAV), a 3D rendering program used for manual annotation [20].

Figure 3: Manual posture identification in two successive image volumes of Fig
1-B using MIPAV. The 20 fluorescently imaged seam cell nuclei rendered in two
successive image volumes. Scale bar: 10 $\mu m$. A & B: Seam cell nuclei
appearing in two successive image volumes visualized in MIPAV. The five minute
interval allows the embryo to reposition between images, yielding entirely
different postures. C & D: Manual seam cell identification by trained users
reveals the posture. The curved lines are cubic splines as described in Fig
2-C.

EHGM uses hypergraphical models comprising biologically driven geometric
features to more accurately identify posture than established graphical
methods. The limited expressive power of graphical models hinders accurate
seam cell identification; graphical models accurately identify posture in 27%
of samples compared to 56% using a hypergraphical model. User labelling of the
posterior-most seam cell nuclei improves the success of hypergraph matching to
correctly identifying all nuclei in 77% of samples. The improved accuracy in
posture identification attributed to high-degree hypergraphical modeling
solved via EHGM paves a path toward automatic posture identification while
presenting a general framework for approaching similarly challenging point-set
matching tasks.

## Results

###  0.2 Posture Identification Models

Posture was predicted via EHGM according to three models: a graphical model,
denoted Sides, and two hypergraphical models. The two hypergraphical models,
Pairs and Posture, showcase EHGM as existing algorithms cannot find solutions
under such high degree hypergraphs. Each of the three models incrementally use
higher degree terms to describe posture. Sides follows the form of Eq 2 and
leverages pairwise assignments to calculate lengths and widths of portions of
the embryo. Pairs uses degrees four and six hyperedges to better model local
regions of the embryo than is possible with graphical methods which rely on
pairwise relationships. Posture further demonstrates the capabilities of EHGM
by including a degree $n_{1}$ hyperedge to maximize context in evaluating a
hypothesized posture. Geometric features such as pair-to-pair rotation angles
and left-right flexion angles were developed to more accurately measure and
compare posture hypotheses. The calculation of each angle or distance requires
identification of multiple seam cells in tandem to calculate, necessitating
the use of hyperedges.

Fig 4 demonstrates four types of models applied to perform posture
identification on the first two sampled images in Fig 1-B. Linear models (Fig
4-A & Fig 4-B) are ill-equipped to identify posture due to the repositioning
of the embryo between successive images, so linear models are not evaluated on
sampled data. The graphical model Sides (Fig 4-C & Fig 4-D) associates local
seam cells via edges (purple). Edge-wise features such as lengths and widths
vary if the embryo coils tightly, but are otherwise approximately static
frame-to-frame. However, the similarity in these measurements throughout the
embryo yields a model incapable of differentiating portions of the embryo.
Hypergraphical models Pairs (Fig 4-E & Fig 4-F) and Posture (Fig 4-G & Fig
4-H) use aforementioned hyperedges to more strongly characterize embryonic
posture.

Figure 4: Posture identification applied to the two successive images in Fig 3
according to a series of increasingly intricate models. The embryo repositions
between images. A & B: Linear models (LAP) cannot quantify relationships
between seam cells; posture identification is impossible without context of
neighboring cell identities. C & D: A graphical model (Sides) specifies edges
(purple) between pairs of seam cell nuclei. Edge lengths are relatively static
frame-to-frame, but the similarity of edge lengths throughout the embryo
causes the edges to have a weak signal in identifying seam cells. E & F: The
Pairs model uses degrees four (red) and six (blue) hyperedges to model a
greater local context than is possible in a graphical model. G & H: The
Posture model extends the Pairs model to use a degree $n_{1}$ (black)
hyperedge to evaluate all seam cell assignments jointly.

###  0.3 Posture Identification Accuracy

Annotators curated a dataset of seam cell nuclei center coordinates from 16
imaged embryos. Each imaged embryo yielded approximately 80 image volumes for
a total of N=1264 labelled seam cell nuclei coordinate sets. Homogeneity in C.
elegans embryo development allowed use of samples spanning multiple embryos to
fit models via a leave-one-out approach (S1:Model Fitting, S1: Posture
Modeling). EHGM allows for known correspondences, henceforth referred to as
seeds, to be given as input prior to search initialization. The algorithm was
evaluated both in a traditional point-set matching scenario given no a priori
information, and in a series of seeded simulations. Seeded trials assumed
incrementally more pairs given sequentially from the tail pair, T, to the
fourth pair, V4 (or Q for $n_{1}$=22 samples). KerGM [2], a leading algorithm
for heuristic graph matching, was applied to posture identification. The
algorithm used the same connectivity matrix as Sides, but processed results
frame-to-frame serially, relying on the correct posture identification at the
prior image as input to search.

EHGM is able to store complete assignments encountered during the search as it
compares against the current solution at the final branch. This allowed for an
analysis of the similarity between cost minimizing posture hypotheses and
progressively higher cost solutions encountered during search. The top x
accuracy describes the percentage of all N samples in which EHGM returned the
correct posture in the x lowest cost solutions; i.e. the top 1 accuracy
describes the percentage of samples in which the correct posture was returned
as the cost minimizing posture, and the top 3 accuracy is the percentage of
samples in which the correct assignment was among 3 lowest cost posture
hypotheses returned by the search. Top x accuracies are reported alongside the
median runtime and the median cost ratio. The cost ratio is defined as the
ratio of the correct posture’s objective to the cost minimizing posture’s
objective. A cost ratio greater than one implies the objective of the
hypothesized posture is lower than that of the correct posture, suggesting the
model is not aptly characterizing posture as an incorrect posture hypothesis
was preferred by the model.

Table 1 shows the percentage of all N samples in which the correct posture
(correct identification of all seam cells) was returned as the minimizer
according to KerGM and each of the models solved via EHGM: Sides, Pairs, and
Posture. KerGM identified 27% of sampled postures correctly, outperforming
Sides (10%). Pairs and Posture more effectively identified posture with 52%
and 56% top 1 accuracies, respectively. Both hypergraphical models also
reported a median cost ratio of 1.00, compared to 1.28 of Sides, suggesting
the hypergraphical representations of coiled posture provided enhanced
discriminatory power across samples. The hypergraphical models demonstrated
small trade-offs between accuracy and runtime. The Posture model’s $n_{1}$
degree hypergraphical features improved accuracy over Pairs, 56% to 52%, in
exchange for longer median runtime, 60 minutes to 43 minutes. Differences
between the top 1 and top 3 accuracies reflect the challenge in posture
identification. The optimums under the Pairs and Posture models were often
similar to those of similar posture hypotheses. Notably, the Posture model
returned the correct posture in the top 3 hypotheses in approximately 67% of
samples, an approximate 20% increase in relative accuracy over the top 1
percentage, 56%.

| Top 1 (%) | Top 2 (%) | Top 3 (%) | Top 5 (%) | Top 10 (%) | R (minutes) | CR
---|---|---|---|---|---|---|---
KerGM | 27 | 27 | 27 | 27 | 27 | .01 |
Sides | 10 | 14 | 15 | 16 | 16 | 5.97 | 1.28
Pairs | 52 | 60 | 63 | 65 | 65 | 43.22 | 1.00
Posture | 56 | 65 | 67 | 68 | 68 | 60.35 | 1.00
Table 1: Hypergraphical model Posture achieves highest accuracy. Posture
identification accuracies across all N=1264 samples. KerGM is compared to
proposed models. The first columns list the top x accuracy as a percentage of
samples. The column titled R shows the median runtime of each model in
minutes. CR reports the median cost ratio, defined as the ratio of the correct
posture cost to the returned posture cost.

Posture identification results were stratified by the presence of the Q
neuroblasts; 875 of the 1264 samples contain only the seam cells while the
remaining 389 samples are mature enough to have the Q neuroblasts. Table 2
depicts the findings presented in Table 1 split by Q neuroblast presence.
KerGM and all models solved via EHGM achieved a higher accuracy on Q samples.
Notably, the Posture model’s top 3 accuracy is higher on the Q samples (82%)
than the pre-Q samples (60%). The extra pair of coordinates provided
substantial context, further defining the coiled shape and helping to penalize
incorrect postures.

| Top 1 (%) | Top 2 (%) | Top 3 (%) | Top 5 (%) | Top 10 (%) | R (minutes) | CR
---|---|---|---|---|---|---|---
KerGM | 25 | 25 | 25 | 25 | 25 | .01 |
Sides | 7 | 10 | 11 | 12 | 12 | 4.81 | 1.36
Pairs | 44 | 51 | 55 | 57 | 58 | 34.25 | 1.04
Posture | 48 | 57 | 60 | 61 | 62 | 51.12 | 1.00

| Top 1 (%) | Top 2 (%) | Top 3 (%) | Top 5 (%) | Top 10 (%) | R (minutes) | CR
---|---|---|---|---|---|---|---
KerGM | 35 | 35 | 35 | 35 | 35 | .01 |
Sides | 19 | 25 | 26 | 26 | 26 | 9.66 | 1.16
Pairs | 71 | 80 | 82 | 82 | 82 | 56.58 | 1.00
Posture | 72 | 81 | 82 | 83 | 83 | 72.60 | 1.00

Table 2: Hypergraphical models leverage Q neuroblasts to identify posture. The
samples are split according to the absence (top) or presence (bottom) of the Q
neuroblasts, which form in the last two hours of development. There are $875$
$n_{1}$=20 cell samples and $389$ $n_{1}$=22 Q samples. Reported methods more
accurately identify embryonic posture in the Q samples, suggesting the
increased continuity along the body of the embryo allows for more consistent
posture identification.

Seeded experiments specifying nuclear identities provided a priori information
starting with the tail pair, and incrementally identified more pairs in the
posterior region. Each experiment was given five minutes of maximum runtime; a
semi-automated solution requiring more runtime was deemed infeasible. Top 1
and top 3 accuracy percentages are reported by EHGM models and number of
seeded pairs in Table 3. Seeding yielded decreasing marginal improvements to
accuracy and runtime. Fig 5 depicts top 1 accuracies and median runtimes
across seeded experiments for the Pairs and Posture models split by Q pair
labelling. Particularly, seeding the first two pairs, T and V6, greatly
reduces the median runtime while also netting the largest gains in top 1
accuracy, partially attributable to EHGM converging in the given timeframe.

| Top 1 (%) | Top 3 (%)
---|---|---
| None | T | T-V6 | T-V5 | T-V4 | None | T | T-V6 | T-V5 | T-V4
Sides | 9 | 10 | 22 | 29 | 37 | 13 | 15 | 27 | 35 | 43
Pairs | 34 | 49 | 72 | 79 | 84 | 38 | 54 | 77 | 83 | 87
Posture | 25 | 36 | 68 | 79 | 84 | 27 | 39 | 73 | 84 | 87
Table 3: Seeding posterior pair identities promotes accurate posture
identification and reduces runtime. Top 1 and top 3 seeded posture
identification accuracies across all samples. All trials had a five-minute
maximum runtime. The rows again correspond to each model. Columns specify
which pairs were given as seeds prior to search. The None columns recreate the
original no information task. The subsequent columns specify which pairs are
correctly identified prior to search. Figure 5: Evaluating the Pairs and
Posture models as seam cell identities were seeded. The Pairs and Posture
models top 1 accuracies and median runtimes by Q pair labelling. Posterior
pair seeding drastically improved top 1 accuracy and reduced runtime when
applying both models. Q pair samples required more runtime ($n_{1}$=22 as
opposed to $n_{1}$=20), but the added context improved posture identification
accuracy. The majority of samples converged within 5 minutes when seeded with
the T and V6 pairs of nuclei.

## Discussion

We have presented EHGM as a dynamic and effective tool for intricate point-set
matching tasks. The hypergraph matching algorithm provides a method in which
to gauge the efficacy of modeling point correspondences in conservatively-
sized problems; problems featuring larger numbers of points likely contain the
context required to match adequately via lower degree models. For example,
postures in samples containing Q nuclei were more accurately identified across
models, but the largest marginal gain in accuracy came from Sides (d=2) to
Pairs (d=4,6). The results suggest that added context throughout the embryo
would further improve posture identification accuracy, reducing the reliance
on higher degree (and thus more computationally expensive) hypergraphical
objective function formulations. EHGM specifically addresses a gap in
literature concerning challenging point-set matching applications in which
domain-specific features lead to rigorously testable models. Seeding allows a
wider range of problems to be approached, and mitigates the computational
expense of the algorithm for scenarios featuring larger point-sets.

Posture identification in embryonic C. elegans is a challenging problem
benefiting from high degree hypergraphical modeling. EHGM equipped with
biologically inspired hypergraphical models led to substantial improvement in
posture identification. The top 1 accuracy doubled from 27% with a graphical
model to 56% via the Posture model (Table 1). The top 3 accuracy rate improved
to 67%, highlighting the challenge in precisely specifying the coiled embryo
due to the similarity of competing posture hypotheses. The presence of Q
neuroblasts further contributed to accurate posture identification. The added
context empowered the Posture model to identify the correct posture in 82% of
Q samples (Table 2.

The top x percentage accuracy metric reflects the need to correctly identify
all seam cells in order to recover the underlying posture, but does not
distinguish between hypotheses that are incorrect due to one cell identity
swap or a more systemic modeling inadequacy. A qualitative analysis
highlighted a few themes among incorrectly predicted postures. The foremost
errors concern the tail pair cells, TL and TR; spurious identifications
occurred when the tail pair coiled against another the body of the embryo,
causing one tail cell identity to be interchanged with a cell of a nearby body
pair. The variance of feature measurements in the posterior region resulted in
similar costs for postures with minor differences about the posterior region.

Pair seeding allows for the strengths of EHGM to compensate for the most
challenging aspect of posture identification. The posterior region of the
embryonic worm is especially flexible and contributes to the majority of
reported errors. Feature engineering stands to create hypergraphical models
more capable of reliable posture identification, particularly in
contextualizing the posterior region. The method and application outline a
protocol for challenging point-set matching tasks.

## Methods

### Exact Hypergraph Matching

EHGM extends the branch-and-bound paradigm to exactly solve hypergraph
matching. The algorithm performs the search in the permutation space
$\mathcal{X}$ subject to a given branch size k which specifies the number of
vertices assigned at each branch. A size $n_{1}$ hypergraph will require
$M\vcentcolon=\frac{n_{1}}{k}$ branch steps, where branch m concerns the
assignment of vertices $((m-1)k+1,(m-1)k+2,\dots,mk)$; vertices $1,2,\dots,mk$
have been assigned upon completion of the $m^{th}$ branch. The set
$\mathbf{P}$ contains all possible permutations of the indices of the
unordered point set, $|\mathbf{P}|=\frac{n_{2}!}{(n_{2}-k)!}$. $\mathbf{P}$ is
incrementally subset into queues $\mathbf{Q}_{m}\subseteq\mathbf{P}$ at
branches $m=1,2,\dots,M$ at each branching. The queue $\mathbf{Q}_{m}$ is
subset according to both a pruning rule which eliminates permutations leading
to a suboptimal solution as well as the one-to-one constraints of
$\mathcal{X}$. The search converges to a global optimum upon the implicit
enumeration of $\mathbf{Q}_{1}=\mathbf{P}$.

The objective function f is further stratified according to the branch size k.
Lower degree ($d\leq 2k$) hyperedge dissimilarity tensors are computed prior
to search. Branches comprising k-tuples of vertices are partially assigned in
a greedy manner according these lower degree hyperedge dissimilarities via the
selection rule H. Later branches accrue higher degree ($d>2k$) hyperedge
dissimilarities which are calculated at time of branching; the intent of the
method is to rely on lower degree terms to steer the search towards an optimum
in effort to minimize the number of branches explored. The aggregation rule I
accrues higher degree hyperedge dissimilarity terms upon branching, further
guiding the pruning step and ensuring the complete specification of the
objective $f$.

The branching and selection rules are designed to reduce computation performed
throughout the search. A partial assignment at branch $m$:
$\mathbf{K}_{m}=(l^{\prime}_{(m-1)k+1},l^{\prime}_{(m-1)k+2},\dots,l^{\prime}_{mk})\in\mathbf{Q}_{m}$
is selected via precomputed lower degree hyperedge dissimilarity tensors
$\mathbf{Z}^{(1)},\dots,\mathbf{Z}^{(2k)}$. A larger branch size k results in
a selection rule with larger scope of the optimization landscape, better
equipped to place optimal branches earlier in each queue $\mathbf{Q}_{m}$ at
time of branching. However, computing the lower degree dissimilarity tensors
prior to search can be prohibitively expensive for larger point-sets.

### Selection & Aggregation

The first branch permutation
$\mathbf{K}_{1}=(l^{\prime}_{1},l^{\prime}_{2},\dots,l^{\prime}_{k})\in\mathbf{Q}_{1}=\mathbf{P}$
assigns vertices $(l_{1},l_{2},\dots,l_{k})$ to points
$(l^{\prime}_{1},l^{\prime}_{2},\dots,l^{\prime}_{k})$ according to the
initial branch selection rule $H_{1}$. Eq 4 defines a cost given dissimilarity
tensors $\mathbf{Z}^{(1)},\mathbf{Z}^{(2)},\dots\mathbf{Z}^{(k)}$ according to
a permutation $\mathbf{K}_{1}$. The k pairs of constraints given by the branch
m and permutation of point indices $\mathbf{K}_{m}$:
$\\{(l_{1},l^{\prime}_{1}),\dots,(l_{k},l^{\prime}_{k})\\}$ enables a
simplification in the objective formulation.

$H_{1}(\mathbf{K}_{1}|\mathbf{Z}^{(1)},\mathbf{Z}^{(2)},...,\mathbf{Z}^{(k)})\vcentcolon=\\\
\sum_{i_{1}=1}^{k}\mathbf{Z}^{(1)}_{l_{i_{1}}l^{\prime}_{i_{1}}}+\sum_{i_{1}=1}^{k}\sum_{i_{2}=i_{1}+1}^{k}\mathbf{Z}^{(2)}_{l_{i_{1}}l^{\prime}_{i_{1}}l_{i_{2}}l^{\prime}_{i_{2}}}+...+\sum_{i_{1}=1}^{k}\sum_{i_{2}=i_{1}+1}^{k}\dots\sum_{i_{k}=i_{k-1}+1}^{k}\mathbf{Z}^{(k)}_{l_{1}l^{\prime}_{i_{1}}l_{i_{2}}l^{\prime}_{i_{2}}\dots
l_{i_{k}}l^{\prime}_{i_{k}}}$ (4)

Subsequent branches $m=2,3,\dots M$ then use the general selection rule
$H_{m}$ to order the permutations of the $m^{th}$ branch:
$\mathbf{K}_{m}=(l^{\prime}_{(m-1)k+1},l^{\prime}_{(m-1)k+2},\dots
l^{\prime}_{mk})\in\mathbf{Q}_{m}$. Branch $\mathbf{K}_{m}$ incurs a selection
rule cost $H_{m}$ according to Eq 5 comprising lower degree hyperedge
dissimilarities for assignments both within branch $m$ and the assignments
between branches $1,2,\dots,m-1$ and branch $m$. The partial assignment
constraints $\mathbf{K}_{m}$ allow further simplification of notation; the
reversed order of summation indices satisfies the criteria that only hyperedge
dissimilarities pertaining to branch $m$ assignments are considered via
$H_{m}$.

$H_{m}(\mathbf{K}_{m}|\mathbf{K}_{1},...,\mathbf{K}_{m-1},\mathbf{Z}^{(1)},...,\mathbf{Z}^{(2k)})\vcentcolon=\\\
\sum_{i_{1}=(m-1)k+1}^{mk}\mathbf{Z}^{(1)}_{l_{i_{1}}l^{\prime}_{i_{1}}}+\sum_{i_{2}=(m-1)k+1}^{mk}\sum_{i_{1}=1}^{i_{2}-1}\mathbf{Z}^{(2)}_{l_{i_{1}}l^{\prime}_{i_{1}}l_{i_{2}}l^{\prime}_{i_{2}}}\\\
+\sum_{i_{3}=(m-1)k+1}^{mk}\sum_{i_{2}=1}^{i_{3}-1}\sum_{i_{1}=1}^{i_{2}-1}\mathbf{Z}^{(3)}_{l_{i_{1}}l^{\prime}_{i_{1}}l_{i_{2}}l^{\prime}_{i_{2}}l_{i_{3}}l^{\prime}_{i_{3}}}+...+\sum_{i_{2k}=(m-1)k+1}^{mk}\sum_{i_{2k-1}=1}^{i_{2k}-1}\dots\sum_{i_{1}=1}^{i_{2}-1}\mathbf{Z}^{(2k)}_{l_{i_{1}}l^{\prime}_{i_{1}}\dots
l_{i_{2k}}l^{\prime}_{i_{2k}}}$ (5)

The greedy selection rule orders queues $\mathbf{Q}_{m}$, but does not account
for higher degree ($2k<d\leq n_{1}$) hyperedge dissimilarities. Precomputing
higher degree dissimilarity tensors can be both computationally expensive, and
inefficient as ideally only a small percentage of combinations are queried
throughout the search. The aggregation rule $I_{m}$, $m=3,4,\dots,M$ measures
the dissimilarity attributable to higher degree ($2k<d\leq mk$) hyperedges
accessible due to branch $m$ partial assignments. The aggregation rule updates
the cost of branch $\mathbf{K}_{m}$ assignments, further informing the pruning
step to subset the next queue $\mathbf{Q}_{m+1}$. The greedy selection rule
$H_{m}$ in tandem with the aggregation rule $I_{m}$ aim to minimize the total
computation performed in finding an optimum. The definition $I_{m}$ follows
from the general selection rule $H_{m}$, but is applied to the higher degree
hyperedge dissimilarities. The aggregation rule $I_{m}$ (Eq 6) can be
expressed as the degree $d$ dissimilarities calculable upon assignments of
branch $m$ assignments for degrees $2k<d\leq mk$.

$I_{m}(\mathbf{K}_{m}|\mathbf{K}_{1},\mathbf{K}_{2},\dots,\mathbf{K}_{m-1},\mathbf{Z}^{(2k+1)},\dots,\mathbf{Z}^{(mk)})\vcentcolon=\sum_{d=2k+1}^{mk}\sum_{i_{d}=(m-1)k+1}^{mk}\sum_{i_{d-1}=1}^{i_{d}-1}...\sum_{i_{1}=1}^{i_{2}-1}\mathbf{Z}^{(d)}_{l_{i_{1}}l^{\prime}_{i_{1}}...l_{i_{d}}l^{\prime}_{i_{d}}}$
(6)

The $m^{th}$ branch allows for hyperedge dissimilarities up to degree mk
concerning the first mk assignments. The $M^{th}$ branch yields a complete
assignment, allowing the evaluation of maximum degree $n_{1}$ hyperedge
dissimilarities. The partitioning and further regrouping of each $H_{m}$ and
$I_{m}$ as defined fully accounts for the objective f while allowing efficient
computation during the search (S1:Hypergraphical Objective Decomposition,
S1:Convergence of EHGM).

### Posture Identification in Embryonic C. elegans

Caenorhabditis elegans (C. elegans) is a small, free-living nematode found
across the world. The worm is often studied as a model of nervous system
development due to its relative simplicity [6, 9]. The adult worm features
only 302 neurons, the morphology and synaptic patterning of which have been
determined via electron microscopy [6]. The complete embryonic cell lineage
has also been determined [7]; methods and technology have been developed to
allow study of cell position and tissue development in the embryo [21, 22, 23,
24, 25, 26]. Systems-level studies of these processes may be able to discover
larger-scale principles underlying developmental events.

The embryo features a set of twenty seam cells and two associated neuroblasts.
The seam cells and neuroblasts together describe anatomical structure in the
coiled embryo, acting as a type of “skeleton” outlining its body.
Identification of the seam cells and neuroblasts defines the embryo’s posture.
Fluorescent proteins are used to label cell nuclei, including the seam cell
nuclei so that they may be visualized during imaging, e.g. with light sheet
microscopy [27]. Volumetric images are captured at five minute intervals in
order to capture subcellular resolution without damaging the worm’s
development [10]. Seam cell nuclei appear in the fluorescent images as
homogeneous spheroids. Their positions relative to other nuclei and other
salient cues present in the image volumes comprise the information that
trained users employ to manually identify seam cells. Fig 6 shows the two
rendered fluorescent images from Fig 1-A in Medical Image Processing, Analysis
and Visualization (MIPAV), a 3D rendering tool [20]. The interface is used to
annotate both seam cell nuclei and track remapped nuclei, as in Fig 3 [10].

Figure 6: Rendered image volumes in the MIPAV GUI. The imaged twisted embryo
(left) and imaged straightened embryo (right) rendered in Medical Image
Processing, Analysis and Visualization (MIPAV) [20]. The fluorescent images
are those depicted in Fig 1-A. Trained users navigate the MIPAV GUI to
identify seam cells based upon relative positioning and other salient features
such as specks of fluorescence on the skin. Correct identification of all
imaged nuclei reveals the coiled embryonic posture. Green (left), red
(center), and purple (right) splines yield an approximation of the coiled
embryo’s posture. Yellow lines connect seam cell nuclei laterally. The splines
are used with the image volume to sweep planes orthogonal to the center
spline, yielding the straightened embryo image.

We cast posture identification as hypergraph matching and use EHGM to solve
the resulting optimization problem. The proposed models: Sides, Pairs, and
Posture trade off modeling capacity for increased computation to identify
optimal solutions. Sides expresses posture identification as graph matching;
edge-wise (degree d=2) features take the form of standardized chord lengths
between nuclei laterally and sequentially along each side. The first
hypergraphical model, Pairs, employs a greater local context than Sides using
degrees four and six hyperedges to describe relationships between seam cells.
Hyperedges formed by two or three sequential pairs (d=4,6) better detail local
regions throughout the embryo than is capable of a graphical model. Fig 7-A
presents the hyperedge connectivity among nodes in the Pairs model [28]. The
Posture model extends the Pairs model by leveraging complete posture
(d=$n_{1}$) features in effort to further discriminate between posture
hypotheses that appear similar in sequential regions of the embryo. Geometric
features help contextualize the coiled posture. Fig 8 illustrates three of the
features used in the Pairs and Posture models. The angle $\Theta$ measures the
angle between three successive pair midpoints. The angles $\Theta$ decrease
throughout development as the worm elongates. Pair-to-pair twist angles
$\varphi$ and $\tau$ penalize posture hypotheses in which posterior to
anterior transitions are jagged and unnatural in appearance. See S1:Posture
Modeling for further details and specification of model features.

Figure 7: The Pairs hypergraphical model uses expansive local contexts about
each portion of the embryo. A: The Pairs hyperedges connect local seam cell
nuclei in sets of four and six. B: Degree four hyperedges connect sequential
pairs of seam cells while degree six hyperedges connect sequential triplets of
pairs. The posterior-most degree four hyperedge and a central degree six
hyperedge are bolded. Figure 8: Hypergraphical geometric features
contextualize seam cell assignments. Anatomically inspired geometric features
describe bend and twist of a posture assignment. A: Three pairs of sequential
nuclei: red, green, blue. Rectangles represent pair midpoints. The angle
$\Theta$ in red is used as a degree six feature given six point to nuclei
assignments. B, C: Degree four hypergraphical features measuring twist angles
$\varphi$ and $\tau$. These angles measure posterior to anterior twist pair-
to-pair and left-right twist, respectively.

The traditional point-set matching task requires a labelled point-set and a
second unidentified point-set. Higher order features such as bend and twist
angles may vary largely frame-to-frame depending on the posture at moment of
imaging. However, elongation throughout late-stage development causes
macroscopic trends in these geometric features. We estimate a template posture
as a composite of feature measurements from a corpus of manually annotated
postures. The templates are time dependent to reflect the elongation from the
first point of imaging throughout development until hatching. See S1:Model
Fitting for details on template estimation.

Together, the fitted models are used with EHGM to identify posture in imaged
C. elegans embryos. The branch size k=2 is set for all models, i.e. a lateral
pair of seam cell identities are assigned at each branch starting with the
tail pair cells TL and TR. The successive pair cells, V6L and V6R, are
assigned given the established cells and hypergraphical relationships
accessible with the hypothesized identities. Fig 9 depicts EHGM applied to the
sample image depicted in Fig 1-A. The initial pair (TL and TR) is selected,
instantiating a search tree (Fig 9-A). Successive seam cell identities are
partially assigned according to the given hypergraphical model in a pair-wise
fashion. Each branch greedily queues hypothesized point-pair assignments
conditioned on the previous branch assignments (black arrows within a branch).
The next leading V6 pair (Fig 9-E) is chosen upon exhaustion of the leading
hypothesized V6 pair (Fig 9-B). EHGM continues the recursion to implicitly
identify a globally optimal posture under the given hypergraphical model; each
possible initial pair will follow this illustrated process subject to pruning
of the minimizing posture accessed via the hypothesized tail pair in Fig 9-A.

Figure 9: EHGM applied to the sample image depicted in Fig 1-A. A: Two points
are selected at the initial branch for TL and TR, respectively. Candidates for
the successive pair, V6L and V6R, are queued based on hypergraphical
relationships between the established cell identities TL and TR and each
hypothesized V6 pair (lower costs are green to higher costs in red). B: The
leading hypothesis at branch m=2 given the initial branch pair is chosen. The
recursion continues to queue V5 pair choices at branch m=3. Black arrows
within branch m specify the ordering of the branch given established cell
assignments. Each branch creates a new subproblem of completing the posture
given partially assigned identities. C: The tree continuing from the V5 pair
hypothesis is fully explored according to the established recursion. D: The
next leading V5 hypothesis is initiated upon exhaustion of the subtree formed
at panel C. E: Implicit enumeration of the subtree formed at panel B causes
the search to progress to the second leading V6 hypothesis.

## Acknowledgments

This work utilized the computational resources of the NIH HPC Biowulf cluster.
(http://hpc.nih.gov). Dr. Evan Ardiel was instrumental in developing
descriptive features for identifying worm posture. Post-Baccelaureate research
fellows Brandon Harvey and Nensi Karaj were supportive in providing data and
discussions concerning the modeling. Dr. Zhen Zhang and Dr. Arye Nehorai
provided assistance in using KerGM [2]. Dr. Vincent Lyzinski also provided
insight on the methods. We also thank Dr. Hank Eden and Dr. Matthew Guay for
their careful readings and suggestions. The code and data are available at
https://github.com/lauziere/EHGM.

## References

  * [1] Zhou F, Torre FDl.  Factorized Graph Matching.  IEEE Transactions on Pattern Analysis and Machine Intelligence. 2016;38(9):1774–1789.  doi:10.1109/TPAMI.2015.2501802.
  * [2] Zhang Z, Xiang Y, Wu L, Xue B, Nehorai A.  KerGM: Kernelized Graph Matching.  In: Wallach H, Larochelle H, Beygelzimer A, Alché-Buc Fd, Fox E, Garnett R, editors. Advances in Neural Information Processing Systems. vol. 32. Curran Associates, Inc.; 2019. p. 1.  Available from: https://proceedings.neurips.cc/paper/2019/file/cd63a3eec3319fd9c84c942a08316e00-Paper.pdf.
  * [3] Duchenne O, Bach F, Kweon IS, Ponce J.  A Tensor-Based Algorithm for High-Order Graph Matching.  IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE. 2010; p. 13.
  * [4] Chertok M, Keller Y.  Efficient High Order Matching.  IEEE Transactions on Pattern Analysis and Machine Intelligence. 2010;32(12):2205–2215.  doi:10.1109/TPAMI.2010.51.
  * [5] Lee J, Cho M, Lee KM.  Hyper-graph matching via reweighted random walks.  In: CVPR 2011; 2011. p. 1633–1640.
  * [6] White JG, Southgate E, Thomson JN, Brenner S.  The structure of the nervous system of the nematode Caenorhabditis elegans.  Philosophical Transactions of the Royal Society of London B, Biological Sciences. 1986;314(1165):1–340.  doi:10.1098/rstb.1986.0056.
  * [7] Sulston JE, Schierenberg E, White JG, Thomson JN.  The embryonic cell lineage of the nematode Caenorhabditis elegans.  Developmental Biology. 1983;100(1):64–119.  doi:10.1016/0012-1606(83)90201-4.
  * [8] Chisholm AD, Hutter H, Jin Y, Wadsworth WG.  The Genetics of Axon Guidance and Axon Regeneration in Caenorhabditis elegans.  Genetics. 2016;204(3):849–882.  doi:10.1534/genetics.115.186262.
  * [9] Rapti G.  A perspective on C. elegans neurodevelopment: from early visionaries to a booming neuroscience research.  Journal of Neurogenetics. 2020;34(3-4):259–272.  doi:10.1080/01677063.2020.1837799.
  * [10] Christensen RP, Bokinsky A, Santella A, Wu Y, Marquina-Solis J, Guo M, et al.  Untwisting the Caenorhabditis elegans embryo.  eLife. 2015;4:e10070.  doi:10.7554/eLife.10070.
  * [11] Leordeanu M, Hebert M.  A spectral technique for correspondence problems using pairwise constraints.  In: Tenth IEEE International Conference on Computer Vision (ICCV’05) Volume 1. Beijing, China: IEEE; 2005. p. 1482–1489 Vol. 2.  Available from: http://ieeexplore.ieee.org/document/1544893/.
  * [12] Berg AC, Berg TL, Malik J.  Shape matching and object recognition using low distortion correspondences.  In: 2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’05). vol. 1; 2005. p. 26–33 vol. 1.
  * [13] Wen L, Li W, Yan J, Lei Z, Yi D, Li SZ.  Multiple Target Tracking Based on Undirected Hierarchical Relation Hypergraph.  In: 2014 IEEE Conference on Computer Vision and Pattern Recognition; 2014. p. 1282–1289.
  * [14] Kuhn HW.  The Hungarian method for the assignment problem.  Naval Research Logistics Quarterly. 1955;2(1-2):83–97.  doi:https://doi.org/10.1002/nav.3800020109.
  * [15] Sahni S, Gonzales T.  P-complete problems and approximate solutions.  In: 15th Annual Symposium on Switching and Automata Theory (swat 1974). USA: IEEE; 1974. p. 28–32.  Available from: http://ieeexplore.ieee.org/document/4569755/.
  * [16] Pardalos P, Du DZ, Graham R.  Handbook of Combinatorial Optimization.  Springer Reference; 2013.
  * [17] Zaslavskiy M, Bach F, Vert JP.  A Path Following Algorithm for the Graph Matching Problem.  IEEE Transactions on Pattern Analysis and Machine Intelligence. 2009;31(12):2227–2242.  doi:10.1109/TPAMI.2008.245.
  * [18] Land AH, Doig AG.  An Automatic Method of Solving Discrete Programming Problems.  Econometrica. 1960;28(3):497–520.  doi:10.2307/1910129.
  * [19] Little JDC, Murty KG, Sweeney DW, Karel C.  An Algorithm for the Traveling Salesman Problem.  Operations Research. 1963;11(6):972–989.  doi:10.1287/opre.11.6.972.
  * [20] McAuliffe MJ, Lalonde FM, McGarry D, Gandler W, Csaky K, Trus BL.  Medical Image Processing, Analysis and Visualization in clinical research.  In: Proceedings 14th IEEE Symposium on Computer-Based Medical Systems. CBMS 2001; 2001. p. 381–386.
  * [21] Bao Z, Murray JI, Boyle T, Ooi SL, Sandel MJ, Waterston RH.  Automated cell lineage tracing in Caenorhabditis elegans.  Proceedings of the National Academy of Sciences of the United States of America. 2006;103(8):2707–2712.  doi:10.1073/pnas.0511111103.
  * [22] Boyle TJ, Bao Z, Murray JI, Araya CL, Waterston RH.  AceTree: a tool for visual analysis of Caenorhabditis elegans embryogenesis.  BMC Bioinformatics. 2006;7(1):275.  doi:10.1186/1471-2105-7-275.
  * [23] Santella A, Catena R, Kovacevic I, Shah P, Yu Z, Marquina-Solis J, et al.  WormGUIDES: an interactive single cell developmental atlas and tool for collaborative multidimensional data exploration.  BMC Bioinformatics. 2015;16(1):189.  doi:10.1186/s12859-015-0627-8.
  * [24] Mace DL, Weisdepp P, Gevirtzman L, Boyle T, Waterston RH.  A High-Fidelity Cell Lineage Tracing Method for Obtaining Systematic Spatiotemporal Gene Expression Patterns in Caenorhabditis elegans.  G3: Genes, Genomes, Genetics. 2013;3(5):851–863.  doi:10.1534/g3.113.005918.
  * [25] Cao J, Guan G, Ho VWS, Wong MK, Chan LY, Tang C, et al.  Establishment of a morphological atlas of the Caenorhabditis elegans embryo using deep-learning-based 4D segmentation.  Nature Communications. 2020;11(1):6254.  doi:10.1038/s41467-020-19863-x.
  * [26] Wang S, Ochoa SD, Khaliullin RN, Gerson-Gurwitz A, Hendel JM, Zhao Z, et al.  A high-content imaging approach to profile C. elegans embryonic development.  Development. 2019;146(7).  doi:10.1242/dev.174029.
  * [27] Wu Y, Wawrzusin P, Senseney J, Fischer RS, Christensen R, Santella A, et al.  Spatially isotropic four-dimensional imaging with dual-view plane illumination microscopy.  Nature Biotechnology. 2013;31(11):1032–1038.  doi:10.1038/nbt.2713.
  * [28] Valdivia P, Buono P, Plaisant C, Dufournaud N, Fekete JD.  Analyzing Dynamic Hypergraphs with Parallel Aggregated Ordered Hypergraph Visualization.  IEEE Transactions on Visualization and Computer Graphics. 2021;27(1):1–13.  doi:10.1109/TVCG.2019.2933196.

## Supporting Information

###  EHGM Pseudocode

EHGM requires as input the branching step $k\in\\{1,2,\dots,n_{1}\\}$, the
dissimilarity tensors $\mathbf{Z}^{(d)}$, $d=1,2,\dots,n_{1}$, the size $k$
permutation set $\mathbf{P}$, and optionally an initial upper bound $C_{0}$ on
the global minimum $C^{*}$. The dissimilarity tensors are calculated given the
reference hypergraph of size $n_{1}$, and either a previous frame hypergraph
or a template hypergraph as described in Model Fitting. The lower degree
dissimilarity tensors $\mathbf{Z}^{(d)},d\leq 2k$ are calculated prior to the
search and used to select branches. The higher degree dissimilarity terms
$d>2k$ are calculated during the search as required. Algorithm 1 initializes
the search from the first candidate set $\mathbf{Q}_{1}=\mathbf{P}$. The
search is parallelized via initializing several first branches. Each explores
a disjoint section of the domain $\mathcal{X}$.

Algorithm 1 initializes arrays and variables to start the recursive branch
search (Algorithm 2). Eligible branch candidates are subset from the general
queue $\mathbf{P}$ into $\mathbf{Q}_{m}$ via the Enqueue procedure (Algorithm
3). Each $\mathbf{Q}_{m}$ contains the potential assignments for the next k
terms that satisfy both the pruning constraints and assignment constraints
specified by $\mathcal{X}$. The current assignment cost $\tilde{C}$ is checked
against the current minimum $C^{*}$ upon reaching a complete assignment. The
Backtrack procedure (Algorithm 4) removes $\mathbf{k}_{m-1}$ from
$\mathbf{Q}_{m-1}$ when the path from $\mathbf{k}_{m-1}$ is exhausted, which
occurs when $\mathbf{Q}_{m}=\emptyset$. The recursion will continue until
$\mathbf{Q}_{1}$ is empty, signaling the complete enumeration of the search
space $S_{n}$.

Input: $k,C_{0},\mathbf{P},\mathbf{Z}^{(1)},\dots,\mathbf{Z}^{(2k)}$

Output: $\mathbf{x}^{*},C^{*}=\mathit{f}(\mathbf{x}^{*})$

Initialization

$C^{*}\leftarrow C_{0}$

$\tilde{\mathbf{H}}=[]$

$\tilde{\mathbf{I}}=[]$

$\tilde{\mathbf{x}}\leftarrow\emptyset$

$\tilde{C}\leftarrow 0$

$m\leftarrow 1$

$\mathbf{Q}_{1}\leftarrow
Enqueue(\tilde{\mathbf{x}},\mathbf{P},\tilde{C},C^{*},1)$

while  _$\mathbf{Q}_{1}\neq\emptyset$_ do

$\mathbf{k}_{1}\leftarrow\mathbf{Q}_{1}.pop()$

$\hat{H_{1}}\leftarrow
H_{1}(\mathbf{k}_{1}|\mathbf{Z}^{(1)},\dots,\mathbf{Z}^{(k)})$

$\tilde{C}\leftarrow\hat{H_{1}}$

$\tilde{\mathbf{H}}[1]=\hat{H_{1}}$

$\tilde{\mathbf{x}}\leftarrow\tilde{\mathbf{x}}\cup\\{\mathbf{k}_{1}\\}$

$m\mathrel{+}=1$

$Visit(\mathbf{P},\tilde{\mathbf{x}},\tilde{C},C^{*},m)$

end while

Return: $\mathbf{x}^{*},C^{*}$

Algorithm 1 EHGM

Input: $\mathbf{P},\tilde{\mathbf{x}},\tilde{C},C^{*},m$

$\mathbf{Q}_{m}\leftarrow
Enqueue(\tilde{\mathbf{x}},\mathbf{P},\tilde{C},C^{*},m)$

while  _$\mathbf{Q}_{m}\neq\emptyset$_ do

$\hat{\mathbf{k}}_{m}\leftarrow\mathbf{Q}_{m}.pop()$

$\tilde{\mathbf{x}}\leftarrow\tilde{\mathbf{x}}\cup\\{\hat{\mathbf{k}}_{m}\\}$

$m\mathrel{+}=1$

$\hat{H}_{m}\leftarrow
H_{m}(\mathbf{k}_{m}|\tilde{\mathbf{x}},\mathbf{Z}^{(1)},\dots,\mathbf{Z}^{(2k)})$

$\tilde{\mathbf{H}}[m]=\hat{H}_{m}$

$\tilde{C}\mathrel{+}=\hat{H}_{m}$

if  _$m\geq 3$_ then

$\hat{I}_{m}\leftarrow
I_{m}(\mathbf{k}_{m}|\tilde{\mathbf{x}},\mathbf{Z}^{(2k+1)},\dots,\mathbf{Z}^{(mk)})$

$\tilde{\mathbf{I}}[m]=\hat{I}_{m}$

$\tilde{C}\mathrel{+}=\hat{I}_{m}$

end if

if  _$m <M$_ then

$Visit(\mathbf{P},\tilde{\mathbf{x}},\tilde{C},C^{*},m)$

else if  _$m=M$_ then

if  _$\tilde{C}\leq C^{*}$_ then

$\mathbf{x}^{*}\leftarrow\tilde{\mathbf{x}}$

$C^{*}\leftarrow\tilde{C}$

end if

$Backtrack(\tilde{\mathbf{x}},\mathbf{Q}_{m-1},m)$

end while

Algorithm 2 Visit

Input: $\tilde{\mathbf{x}},\mathbf{P},\tilde{C},C^{*},m$

Output: $\mathbf{Q}_{m}$

$\mathbf{Q}_{m}\leftarrow\emptyset$

for  _$\mathbf{k}\in\mathbf{P}$_ do

if
_$(\mathbf{k}\cap\tilde{\mathbf{x}}=\emptyset)\wedge(\tilde{C}+H_{m}(\mathbf{k}|\tilde{\mathbf{x}})
<C^{*})$_ then

$\mathbf{Q}_{m}\leftarrow\mathbf{Q}_{m}\cup\mathbf{k}$

end if

end for

Algorithm 3 Enqueue

Input: $\tilde{\mathbf{x}},\mathbf{Q}_{m-1},m$

$\mathbf{Q}_{m-1}\leftarrow\mathbf{Q}_{m-1}\setminus\tilde{\mathbf{x}}_{m}$

$\tilde{\mathbf{x}}\leftarrow\\{\hat{\mathbf{x}}_{1},\hat{\mathbf{x}}_{2},\dots,\hat{\mathbf{x}}_{m-1}\\}$

$\tilde{C}\leftarrow\Sigma_{j=1}^{m-1}(\mathbf{C}_{m}+\mathbf{I}_{m})$

Algorithm 4 Backtrack

### Hypergraphical Objective Decomposition

The hypergraphical optimization objective can be decomposed according to
hyperedge multiplicity and branching step. The stratification enables
efficient search via EHGM.

######  Theorem 1.

Assume an assignment problem objective f is in the form:

$\mathit{f}(X|\mathbf{Z}^{(1)},\mathbf{Z}^{(2)},\dots,\mathbf{Z}^{(n_{1})})=\sum_{l_{1}=1}^{n_{1}}\sum_{l^{\prime}_{1}=1}^{n_{2}}\mathbf{Z}^{(1)}_{l_{1}l^{\prime}_{1}}x_{l_{1}l^{\prime}_{1}}+\sum_{l_{1}=1}^{n_{1}}\sum_{l^{\prime}_{1}=1}^{n_{2}}\sum_{l_{2}=l_{1}+1}^{n_{1}}\sum_{l^{\prime}_{2}=1}^{n_{2}}\mathbf{Z}^{(2)}_{l_{1}l^{\prime}_{1}l_{2}l^{\prime}_{2}}x_{l_{1}l^{\prime}_{1}}x_{l_{2}l^{\prime}_{2}}\\\
+\sum_{l_{1}=1}^{n_{1}}\sum_{l^{\prime}_{1}=1}^{n_{2}}\sum_{l_{2}=l_{1}+1}^{n_{1}}\sum_{l^{\prime}_{2}=1}^{n_{2}}\sum_{l_{3}=l_{2}+1}^{n_{1}}\sum_{l^{\prime}_{3}=1}^{n_{2}}\mathbf{Z}^{(3)}_{l_{1}l^{\prime}_{1}l_{2}l^{\prime}_{2}l_{3}l^{\prime}_{3}}x_{l_{1}l^{\prime}_{1}}x_{l_{2}l^{\prime}_{2}}x_{l_{3}l^{\prime}_{3}}+...\\\
+\sum_{l_{1}=1}^{n_{1}}\sum_{l^{\prime}_{1}=1}^{n_{2}}...\sum_{l_{n_{1}}=l_{n_{1}-1}+1}^{n_{1}}\sum_{l^{\prime}_{n_{1}}=1}^{n_{2}}\mathbf{Z}^{(n_{1})}_{l_{1}l^{\prime}_{1}\dots
l_{n_{1}}l^{\prime}_{n_{1}}}x_{l_{1}l^{\prime}_{1}}\dots
x_{l_{n_{1}}l^{\prime}_{n_{1}}}$ (7)

Then, for $k\in\\{1,2,\dots,n_{1}$}, the stratification fully describes the
objective f after $M=\frac{n_{1}}{k}$ branches. Define of $H_{1}$, $H_{m}$,
and $I_{m}$:

$H_{1}(\mathbf{K}_{1}|\mathbf{Z}^{(1)},\mathbf{Z}^{(2)},...,\mathbf{Z}^{(k)})\vcentcolon=\\\
\sum_{i_{1}=1}^{k}\mathbf{Z}^{(1)}_{l_{i_{1}}l^{\prime}_{i_{1}}}+\sum_{i_{1}=1}^{k}\sum_{i_{2}=i_{1}+1}^{k}\mathbf{Z}^{(2)}_{l_{i_{1}}l^{\prime}_{i_{1}}l_{i_{2}}l^{\prime}_{i_{2}}}+...+\sum_{i_{1}=1}^{k}\sum_{i_{2}=i_{1}+1}^{k}\dots\sum_{i_{k}=i_{k-1}+1}^{k}\mathbf{Z}^{(k)}_{l_{1}l^{\prime}_{i_{1}}l_{i_{2}}l^{\prime}_{i_{2}}\dots
l_{i_{k}}l^{\prime}_{i_{k}}}$

$H_{m}(\mathbf{K}_{m}|\mathbf{K}_{1},...,\mathbf{K}_{m-1},\mathbf{Z}^{(1)},...,\mathbf{Z}^{(2k)})\vcentcolon=\\\
\sum_{i_{1}=(m-1)k+1}^{mk}\mathbf{Z}^{(1)}_{l_{i_{1}}l^{\prime}_{i_{1}}}+\sum_{i_{2}=(m-1)k+1}^{mk}\sum_{i_{1}=1}^{i_{2}-1}\mathbf{Z}^{(2)}_{l_{i_{1}}l^{\prime}_{i_{1}}l_{i_{2}}l^{\prime}_{i_{2}}}\\\
+\sum_{i_{3}=(m-1)k+1}^{mk}\sum_{i_{2}=1}^{i_{3}-1}\sum_{i_{1}=1}^{i_{2}-1}\mathbf{Z}^{(3)}_{l_{i_{1}}l^{\prime}_{i_{1}}l_{i_{2}}l^{\prime}_{i_{2}}l_{i_{3}}l^{\prime}_{i_{3}}}+...+\sum_{i_{2k}=(m-1)k+1}^{mk}\sum_{i_{2k-1}=1}^{i_{2k}}\dots\sum_{i_{1}=1}^{i_{2}-1}\mathbf{Z}^{(2k)}_{l_{i_{1}}l^{\prime}_{i_{1}}\dots
l_{i_{2k}}l^{\prime}_{i_{2k}}}$

$I_{m}(\mathbf{K}_{m}|\mathbf{K}_{1},\mathbf{K}_{2},\dots,\mathbf{K}_{m-1},\mathbf{Z}^{(2k+1)},\dots,\mathbf{Z}^{(mk)})\vcentcolon=\sum_{d=2k+1}^{mk}\Xi^{(d)}_{m}$

where

$\Xi^{(d)}_{m}(\mathbf{K}_{m}|\mathbf{K}_{1},\mathbf{K}_{2},\dots,\mathbf{K}_{m-1},\mathbf{Z}^{(2k+1)},\dots,\mathbf{Z}^{(mk)})\vcentcolon=\\\
\sum_{i_{d}=(m-1)k+1}^{mk}\sum_{i_{d-1}=1}^{i_{d}-1}...\sum_{i_{1}=1}^{i_{2}-1}\mathbf{Z}^{(d)}_{l_{i_{1}}l^{\prime}_{i_{1}}...l_{i_{d}}l^{\prime}_{i_{d}}}$

Then, the degree $n_{1}$ hypergraph matching objective f can be expressed

$\mathit{f}(X|\mathbf{Z}^{(1)},\mathbf{Z}^{(2)},\dots,\mathbf{Z}^{(n_{1})})=\sum_{m=1}^{m}\mathit{H}_{m}+\sum_{m=3}^{m}I_{m}$

###### Proof.

First consider the single branching case $k=1$. This yields
$M=\frac{n_{1}}{k}=\frac{n_{1}}{1}=n_{1}$ branches. Each branch yields one
assignment; i.e. $K_{m}=l^{\prime}_{m}$ is assigned to the vertex $l_{m}$. The
initial branch selection rule $H_{1}$ can only utilize the first order term:

$H_{1}(K_{1}|\mathbf{Z}^{(1)})=\mathbf{Z}^{(1)}_{l_{1}l^{\prime}_{1}}$

Then the general selection rule for the second branch will: gather the first
order costs for the second assignment as well as the quadratic (second order)
costs between the first two assignments:

$H_{2}(K_{2}|K_{1},\mathbf{Z}^{(1)},\mathbf{Z}^{(2)})=\mathbf{Z}^{(1)}_{l_{2}l^{\prime}_{2}}+\sum_{i_{2}=2}^{2}\sum_{i_{1}=1}^{i_{2}}\mathbf{Z}^{(2)}_{l_{i_{1}}l^{\prime}_{i_{1}}l_{i_{2}}l^{\prime}_{i_{2}}}=\mathbf{Z}^{(1)}_{l_{2}}+\mathbf{Z}^{(2)}_{l_{1}l^{\prime}_{1}l_{2}l^{\prime}_{2}}$

The third branching step will include $H_{3}$ and $I_{3}$. $H_{3}$ follows
from $H_{2}$:

$H_{3}(K_{3}|K_{1},K_{2},\mathbf{Z}^{(1)},\mathbf{Z}^{(2)})=\mathbf{Z}^{(1)}_{l_{3}l^{\prime}_{3}}+\sum_{i_{2}=3}^{3}\sum_{i_{1}=1}^{i_{2}-1}\mathbf{Z}^{(2)}_{l_{i_{1}}l^{\prime}_{i_{1}}l_{i_{2}}l^{\prime}_{i_{2}}}=\mathbf{Z}^{(1)}_{l_{3}l^{\prime}_{3}}+\mathbf{Z}^{(2)}_{l_{1}l^{\prime}_{1}l_{3}l^{\prime}_{3}}+\mathbf{Z}^{(2)}_{l_{2}l^{\prime}_{2}l_{3}l^{\prime}_{3}}$

$I_{3}(K_{3}|K_{1},K_{2},\mathbf{Z}^{(3)})=\Xi^{(3)}_{m}=\mathbf{Z}^{(3)}_{l_{1}l^{\prime}_{1}l_{2}l^{\prime}_{2}l_{3}l^{\prime}_{3}}$

Note that if $n_{1}=3$, then $H_{1}+H_{2}+H_{3}+I_{3}$ fully describes the
third order assignment problem:

$H_{1}+H_{2}+H_{3}+I_{3}\\\
=\underbrace{\mathbf{Z}^{(1)}_{l_{1}l^{\prime}_{1}}}_{H_{1}}+\underbrace{\mathbf{Z}^{(1)}_{l_{2}l^{\prime}_{2}}+\mathbf{Z}^{(2)}_{l_{1}l^{\prime}_{1}l_{2}l^{\prime}_{2}}}_{H_{2}}+\underbrace{\mathbf{Z}^{(1)}_{l_{3}l^{\prime}_{3}}+\mathbf{Z}^{(2)}_{l_{1}l^{\prime}_{1}l_{3}l^{\prime}_{3}}+\mathbf{Z}^{(2)}_{l_{2}l^{\prime}_{2}l_{3}l^{\prime}_{3}}}_{H_{3}}+\underbrace{\mathbf{Z}^{(3)}_{l_{1}l^{\prime}_{1}l_{2}l^{\prime}_{2}l_{3}l^{\prime}_{3}}}_{I_{3}}\\\
=\sum_{i_{1}=1}^{3}\mathbf{Z}^{(1)}_{l_{i_{1}}l^{\prime}_{i_{1}}}+\sum_{i_{1}=1}^{3}\sum_{i_{2}=i_{1}+1}^{3}\mathbf{Z}^{(2)}_{l_{i_{1}}l^{\prime}_{i_{1}}l_{i_{2}}l^{\prime}_{i_{2}}}+\sum_{i_{1}=1}^{3}\sum_{i_{2}=i_{1}+1}^{3}\sum_{i_{3}=i_{2}+1}^{3}\mathbf{Z}^{(3)}_{l_{i_{1}}l^{\prime}_{i_{1}}l_{i_{2}}l^{\prime}_{i_{2}}l_{i_{3}}l^{\prime}_{i_{3}}}\\\
=f(X|\mathbf{Z}^{(1)},\mathbf{Z}^{(2)},\mathbf{Z}^{(3)})$ (8)

Now consider the extension to $n_{1}=4$, yielding a fourth degree assignment
problem. The fourth branch will assign the next term, $K_{4}=l^{\prime}_{4}$.
The terms $H_{4}$ and $I_{4}$ will then fully specify the fourth degree
problem:

$H_{4}(K_{4}|K_{1},K_{2},K_{3},\mathbf{Z}^{(1)},\mathbf{Z}^{(2)})=\mathbf{Z}^{(1)}_{l_{4}l^{\prime}_{4}}+\mathbf{Z}^{(2)}_{l_{1}l^{\prime}_{1}l_{4}l^{\prime}_{4}}+\mathbf{Z}^{(2)}_{l_{2}l^{\prime}_{2}l_{4}l^{\prime}_{4}}+\mathbf{Z}^{(2)}_{l_{3}l^{\prime}_{3}l_{4}l^{\prime}_{4}}$

The second aggregation rule $I_{4}$ will consider third order terms between
branches $1,2$ and $4$ as well as the fourth order term using all four
assignments:

$I_{4}(K_{4}|K_{1},K_{2},K_{3},\mathbf{Z}^{(3)},\mathbf{Z}^{(4)})=\Xi^{(3)}_{4}+\Xi^{(4)}_{4}=\mathbf{Z}^{(3)}_{l_{1}l^{\prime}_{1}l_{2}l^{\prime}_{2}l_{4}l^{\prime}_{4}}+\mathbf{Z}^{(3)}_{l_{2}l^{\prime}_{2}l_{3}l^{\prime}_{3}l_{4}l^{\prime}_{4}}+\mathbf{Z}^{(4)}_{l_{1}l^{\prime}_{1}l_{2}l^{\prime}_{2}l_{3}l^{\prime}_{3}l_{4}l^{\prime}_{4}}$

Joining the fourth branch:

$H_{1}+H_{2}+H_{3}+I_{3}+H_{4}+I_{4}=\\\
\underbrace{\sum_{i_{1}=1}^{3}\mathbf{Z}^{(1)}_{l_{i_{1}}l^{\prime}_{i_{1}}}+\sum_{i_{1}=1}^{3}\sum_{i_{2}=i_{1}+1}^{3}\mathbf{Z}^{(2)}_{l_{i_{1}}l^{\prime}_{i_{1}}l_{i_{2}}l^{\prime}_{i_{2}}}+\sum_{i_{1}=1}^{3}\sum_{i_{2}=i_{1}+1}^{3}\sum_{i_{3}=i_{2}+1}^{3}\mathbf{Z}^{(3)}_{l_{i_{1}}l^{\prime}_{i_{1}}l_{i_{2}}l^{\prime}_{i_{2}}l_{i_{3}}l^{\prime}_{i_{3}}}}_{H_{1}+H_{2}+H_{3}+I_{3}}\\\
+\underbrace{\mathbf{Z}^{(1)}_{l_{4}l^{\prime}_{4}}+\mathbf{Z}^{(2)}_{l_{1}l^{\prime}_{1}l_{4}l^{\prime}_{4}}+\mathbf{Z}^{(2)}_{l_{2}l^{\prime}_{2}l_{4}l^{\prime}_{4}}+\mathbf{Z}^{(2)}_{l_{3}l^{\prime}_{3}l_{4}l^{\prime}_{4}}}_{H_{4}}+\underbrace{\mathbf{Z}^{(3)}_{l_{1}l^{\prime}_{1}l_{2}l^{\prime}_{2}l_{4}l^{\prime}_{4}}+\mathbf{Z}^{(3)}_{l_{2}l^{\prime}_{2}l_{3}l^{\prime}_{3}l_{4}l^{\prime}_{4}}+\mathbf{Z}^{(4)}_{l_{1}l^{\prime}_{1}l_{2}l^{\prime}_{2}l_{3}l^{\prime}_{3}l_{4}l^{\prime}_{4}}}_{I_{4}}\\\
=\sum_{i_{1}=1}^{4}\mathbf{Z}^{(1)}_{l_{i_{1}}l^{\prime}_{i_{1}}}+\sum_{i_{1}=1}^{4}\sum_{i_{2}=i_{1}+1}^{4}\mathbf{Z}^{(2)}_{l_{i_{1}}l^{\prime}_{i_{1}}l_{i_{2}}l^{\prime}_{i_{2}}}+\sum_{i_{1}=1}^{4}\sum_{i_{2}=i_{1}+1}^{4}\sum_{i_{3}=i_{2}+1}^{4}\mathbf{Z}^{(3)}_{l_{i_{1}}l^{\prime}_{i_{1}}l_{i_{2}}l^{\prime}_{i_{2}}l_{i_{3}}l^{\prime}_{i_{3}}}\\\
+\sum_{i_{1}=1}^{4}\sum_{i_{2}=i_{1}+1}^{4}\sum_{i_{3}=i_{2}+1}^{4}\sum_{i_{4}=i_{3}+1}^{4}\mathbf{Z}^{(4)}_{l_{i_{1}}l^{\prime}_{i_{1}}l_{i_{2}}l^{\prime}_{i_{2}}l_{i_{3}}l^{\prime}_{i_{3}}l_{i_{4}}l^{\prime}_{i_{4}}}\\\
=f(X|\mathbf{Z}^{(1)},\mathbf{Z}^{(2)},\mathbf{Z}^{(3)},\mathbf{Z}^{(4)})$ (9)

Now consider the arbitrary $(m+1)^{st}$ branch. This will yield the full
objective for an assignment problem of size $m+1$ up to degree $m+1$.

$\sum_{p=1}^{m+1}H_{p}+\sum_{p=3}^{m+1}I_{p}=\sum_{p=1}^{m}H_{p}+\sum_{p=3}^{m}I_{p}+H_{m+1}+I_{m+1}=\\\
\underbrace{\sum_{i_{1}=1}^{m}\mathbf{Z}^{(1)}_{l_{i_{1}}l^{\prime}_{i_{1}}}+\sum_{i_{1}=1}^{m}\sum_{i_{2}=i_{1}+1}^{m}\mathbf{Z}^{(2)}_{l_{i_{1}}l^{\prime}_{i_{1}}l_{i_{2}}l^{\prime}_{i_{2}}}+\dots+\sum_{i_{1}=1}^{m}\sum_{i_{2}=i_{1}+1}^{m}...\sum_{i_{m}=i_{m-1}+1}^{m}\mathbf{Z}^{(m)}_{l_{i_{1}}l^{\prime}_{i_{1}}l_{i_{2}}l^{\prime}_{i_{2}}...l_{i_{m}}l^{\prime}_{i_{m}}}}_{\sum_{p=1}^{m}H_{p}+\sum_{p=3}^{m}I_{p}}\\\
+\underbrace{\mathbf{Z}^{(1)}_{l_{i_{m+1}}l^{\prime}_{i_{m+1}}}+\sum_{i_{1}=1}^{m}\mathbf{Z}^{(2)}_{l_{i_{1}}l^{\prime}_{i_{1}}l_{m+1}l^{\prime}_{m+1}}}_{H_{m+1}}+\underbrace{\sum_{d=3}^{m+1}\Xi^{(d)}_{m+1}}_{I_{m+1}}$
(10)

It is sufficient to show each degree $d\in\\{1,2,\dots,m+1\\}$ hyperedge is
fully accounted for across all $m+1$ points to prove the $(m+1)^{st}$ branch
satisfies the objective f. The hyperedge costs across all points will be
decomposed into three disjoint sets, and each set considered at a time:

$\\{1\\},\\{2\\},\\{3,\dots,m\\},\\{m+1\\}$

The first and final of the four cases are trivial. The first degree terms are
enumerated via the first term in $H_{m+1}$, while $Xi^{(m+1)}_{m+1}$
explicitly addresses the degree $m+1$ hyperedge comprising all assignments:
$\mathbf{Z}^{(m+1)}_{l_{1}l^{\prime}_{1}l_{2}l^{\prime}_{2}...l_{m+1}l^{\prime}_{m+1}}$.
We will focus on the second and third cases. The degree $d=2$ terms are formed
by the addition of branch $m+1$ are considered in term $H_{m+1}$:

$\sum_{i_{1}=1}^{m}\sum_{i_{2}=i_{1}+1}^{m}\mathbf{Z}^{(2)}_{l_{i_{1}}l^{\prime}_{i_{1}}l_{i_{2}}l^{\prime}_{i_{2}}}+\sum_{i_{1}=1}^{m}\mathbf{Z}^{(2)}_{l_{i_{1}}l^{\prime}_{i_{1}}l_{i_{m+1}}l^{\prime}_{i_{m+1}}}=\sum_{i_{1}=1}^{m+1}\sum_{i_{2}=i_{1}+1}^{m+1}\mathbf{Z}^{(2)}_{l_{i_{1}}l^{\prime}_{i_{1}}l_{i_{2}}l^{\prime}_{i_{2}}}$

Let $d\in\\{3,\dots,m\\}$. The completion is similar to the $d=2$ degree case;
however, the term $\Xi^{(d)}_{m+1}$ in $I_{m+1}$ address higher degree
hyperedges up to and including degree $m$ concerning branch $m+1$:

$\sum_{i_{1}=1}^{m}...\sum_{i_{d}=i_{d-1}+1}^{m}\mathbf{Z}^{(d)}_{l_{i_{1}}l^{\prime}_{i_{1}}...l_{i_{d}}l^{\prime}_{i_{d}}}+\Xi^{(d)}_{m+1}=\sum_{i_{1}=1}^{m+1}...\sum_{i_{d}=i_{d-1}+1}^{m+1}\mathbf{Z}^{(d)}_{l_{i_{1}}l^{\prime}_{i_{1}}...l_{i_{d}}l^{\prime}_{i_{d}}}$

Therefore, the $(m+1)^{st}$ step fully accrues the objective f:

$\sum_{p=1}^{m+1}H_{p}+\sum_{p=3}^{m+1}I_{p}=f(X|\mathbf{Z}^{(1)},\mathbf{Z}^{(2)},\dots,\mathbf{Z}^{(m+1)})$

Then inductively, the stratification holds such that:

$\sum_{m=1}^{n_{1}}H_{m}+\sum_{m=3}^{n_{1}}I_{m}=f(X|\mathbf{Z}^{(1)},\mathbf{Z}^{(2)},\dots\mathbf{Z}^{(n_{1})})$

Now consider the plural branching rule $k>1$. The proof will follow from the
single assignment branching case. The base case at the fourth branch will be
established, followed by the induction hypothesis demonstrating the branching
from $m$ to $m+1$. First, define the terms $H_{1},H_{2},H_{3},I_{3},H_{4}$,
and $I_{4}$:

$H_{1}=\sum_{i_{1}=1}^{k}\mathbf{Z}^{(1)}_{l_{i_{1}}l^{\prime}_{i_{1}}}+\sum_{i_{1}=1}^{k}\sum_{i_{2}=i_{1}+1}^{k}\mathbf{Z}^{(2)}_{l_{i_{1}}l^{\prime}_{i_{1}}l_{i_{2}}l^{\prime}_{i_{2}}}+...+\sum_{i_{1}=1}^{k}\sum_{i_{2}=i_{1}+1}^{k}...\sum_{i_{k}=i_{k-1}+1}^{k}\mathbf{Z}^{(k)}_{l_{i_{1}}l^{\prime}_{i_{1}}l_{i_{2}}l^{\prime}_{i_{2}}...l_{i_{k}}l^{\prime}_{i_{k}}}$
(11)

$H_{2}=\sum_{i_{1}=k+1}^{2k}\mathbf{Z}^{(1)}_{l_{i_{1}}l^{\prime}_{i_{1}}}+\sum_{i_{2}=k+1}^{2k}\sum_{i_{1}=1}^{i_{2}-1}\mathbf{Z}^{(2)}_{l_{i_{1}}l^{\prime}_{i_{1}}l_{i_{2}}l^{\prime}_{i_{2}}}+...+\sum_{i_{2k}=k+1}^{2k}...\sum_{i_{2}=1}^{i_{3}-1}\sum_{i_{1}=1}^{i_{2}-1}\mathbf{Z}^{(2k)}_{l_{i_{1}}l^{\prime}_{i_{1}}...l_{i_{2k}}l^{\prime}_{i_{2k}}}$
(12)

$H_{3}=\sum_{i_{1}=2k+1}^{3k}\mathbf{Z}^{(1)}_{l_{i_{1}}l^{\prime}_{i_{1}}}+\sum_{i_{2}=2k+1}^{3k}\sum_{i_{1}=1}^{i_{2}-1}\mathbf{Z}^{(2)}_{l_{i_{1}}l^{\prime}_{i_{1}}l_{i_{2}}l^{\prime}_{i_{2}}}+...+\sum_{i_{2k}=2k+1}^{3k}...\sum_{i_{1}=1}^{i_{2}-1}\mathbf{Z}^{(2k)}_{l_{i_{1}}l^{\prime}_{i_{1}}...l_{i_{2k}}l^{\prime}_{i_{2k}}}$
(13)

$I_{3}=\sum_{d=2k+1}^{3k}\Xi^{(d)}_{3}$ (14)

$H_{4}=\sum_{i_{1}=3k+1}^{4k}\mathbf{Z}^{(1)}_{l_{i_{1}}l^{\prime}_{i_{1}}}+\sum_{i_{2}=3k+1}^{4k}\sum_{i_{1}=1}^{i_{2}-1}\mathbf{Z}^{(2)}_{l_{i_{1}}l^{\prime}_{i_{1}}l_{i_{2}}l^{\prime}_{i_{2}}}+...+\sum_{i_{2k}=3k+1}^{4k}...\sum_{i_{1}=1}^{i_{2}-1}\mathbf{Z}^{(2k)}_{l_{i_{1}}l^{\prime}_{i_{1}}...l_{i_{2k}}l^{\prime}_{i_{2k}}}$
(15)

$I_{4}=\sum_{d=3k+1}^{4k}\Xi^{(d)}_{4}$ (16)

The terms presented thus far for the general $k>1$ case fully describe all
terms concerning assignments $1,2,\dots 4k$ up to degree $4k$. The hyperedge
multiplicities will again be partitioned into disjoint groups:

$\\{1\\},\\{2,\dots,k\\},\\{k+1,\dots,2k\\},\\{2k+1,\dots,3k\\},\\{3k+1,\dots,4k\\}$

The first case is trivial, just as in the single assignment branching ($k=1$)
proof. Unary terms are accounted for in the first summand of each $H_{m}$.
Then, consider $d\in\\{2,\dots,k\\}$:

$\underbrace{\sum_{i_{1}=1}^{k}...\sum_{i_{d}=i_{d-1}+1}^{k}\mathbf{Z}^{(d)}_{l_{i_{1}}l^{\prime}_{i_{1}}...l_{i_{d}}l^{\prime}_{i_{d}}}}_{H_{1}}+\underbrace{\sum_{i_{d}=k+1}^{2k}\sum_{i_{d-1}=1}^{i_{d}-1}...\sum_{i_{1}=1}^{i_{2}-1}\mathbf{Z}^{(d)}_{l_{i_{1}}l^{\prime}_{i_{1}}...l_{i_{d}}l^{\prime}_{i_{d}}}}_{H_{2}}\\\
+\underbrace{\sum_{i_{d}=2k+1}^{3k}\sum_{i_{d-1}=1}^{i_{d}-1}...\sum_{i_{1}=1}^{i_{2}-1}\mathbf{Z}^{(d)}_{l_{i_{1}}l^{\prime}_{i_{1}}...l_{i_{d}}l^{\prime}_{i_{d}}}}_{H_{3}}+\underbrace{\sum_{i_{d}=3k+1}^{4k}\sum_{i_{d-1}=1}^{i_{d}-1}...\sum_{i_{1}=1}^{i_{2}-1}\mathbf{Z}^{(d)}_{l_{i_{1}}l^{\prime}_{i_{1}}...l_{i_{d}}l^{\prime}_{i_{d}}}}_{H_{4}}\\\
=\sum_{i_{1}=1}^{4k}\sum_{i_{2}=i_{1}+1}^{4k}...\sum_{i_{d}=i_{d-1}+1}^{4k}\mathbf{Z}^{(d)}_{l_{i_{1}}l^{\prime}_{i_{1}}...l_{i_{d}}l^{\prime}_{i_{d}}}$

The proof for degree $d\in\\{k+1,\dots,2k\\}$ follows immediately from the
grouping presented above, but without the initial branch selection rule term
$H_{1}$. Next, assume $d\in\\{2k+1,\dots,(m-1)k\\}$. Degree $d$ hyperedge
dissimilarities will be contained in both $I_{3}$ and $I_{4}$ terms:

$\underbrace{\sum_{i_{d}=2k+1}^{3k}\sum_{i_{d-1}=1}^{i_{d}-1}...\sum_{i_{1}=1}^{i_{2}-1}\mathbf{Z}^{(d)}_{l_{i_{1}}l^{\prime}_{i_{1}}...l_{i_{d}}l^{\prime}_{i_{d}}}}_{I_{3}}+\underbrace{\sum_{i_{d}=3k+1}^{4k}\sum_{i_{d-1}=1}^{i_{d}-1}...\sum_{i_{1}=1}^{i_{2}-1}\mathbf{Z}^{(d)}_{l_{i_{1}}l^{\prime}_{i_{1}}...l_{i_{d}}l^{\prime}_{i_{d}}}}_{I_{4}}\\\
=\sum_{i_{1}=1}^{4k}...\sum_{i_{d-1}=i_{d-2}+1}^{4k}\sum_{i_{d}=i_{d-1}+1}^{4k}\mathbf{Z}^{(d)}_{l_{i_{1}}l^{\prime}_{i_{1}}...l_{i_{d}}l^{\prime}_{i_{d}}}$

Since $d\leq 2k+1$, the terms only appear in the third branch term $I_{3}$
when the assignment $2k+1$ is committed. The final set arises from the
definition of $I_{4}$ which accrues hyperedges of degree
$d\in\\{3k+1,\dots,4k\\}$ across assignments in branches $m=1,2,3,4$. The base
case is fully established for the arbitrary $k>1$ case. The final step of the
proof is to establish the extension of the $(m+1)^{st}$ branch:

$\sum_{p=1}^{m+1}H_{p}+\sum_{p=3}^{m+1}I_{p}=\sum_{p=1}^{m}H_{p}+\sum_{p=3}^{m}I_{p}+H_{m+1}+I_{m+1}=\\\
\underbrace{\sum_{i_{1}=1}^{mk}\mathbf{Z}^{(1)}_{l_{i_{1}}l^{\prime}_{i_{1}}}+\sum_{i_{1}=1}^{mk}\sum_{i_{2}=i_{1}+1}^{mk}\mathbf{Z}^{(2)}_{l_{i_{1}}l^{\prime}_{i_{1}}l_{i_{2}}l^{\prime}_{i_{2}}}+\dots+\sum_{i_{1}=1}^{mk}\sum_{i_{2}=i_{1}+1}^{mk}...\sum_{i_{mk}=i_{mk-1}+1}^{mk}\mathbf{Z}^{(mk)}_{l_{i_{1}}l^{\prime}_{i_{1}}...l_{i_{mk}}l^{\prime}_{i_{mk}}}}_{\sum_{p=1}^{m}H_{p}+\sum_{p=3}^{m}I_{p}}\\\
+\underbrace{\sum_{i_{1}=mk+1}^{(m+1)k}\mathbf{Z}^{(1)}_{l_{i_{1}}l^{\prime}_{i_{1}}}+...+\sum_{i_{2k}=mk+1}^{(m+1)k}\sum_{i_{2k-1}=1}^{i_{2k}-1}...\sum_{i_{1}=1}^{i_{2}-1}\mathbf{Z}^{(2k)}_{l_{i_{1}}l^{\prime}_{i_{1}}...l_{i_{2k}}l^{\prime}_{i_{2k}}}}_{H_{m+1}}+\underbrace{\sum_{d=2k+1}^{(m+1)k}\Xi^{(d)}_{m+1}}_{I_{m+1}}$
(17)

The $(m+1)k$ hyperedge multiplicities will be stratified into four groups:

$\\{1\\},\\{2,\dots,2k\\},\\{2k+1,\dots,mk\\},\\{mk+1,\dots,(m+1)k\\}$

Just as in the singular $k=1$ case, the proof for the first and last groups
are trivial. The unary terms are again evident from the first term in
$H_{m+1}$, while the $mk+1\leq d\leq(m+1)k$ terms in $I_{m+1}$ fully
encapsulates the fourth group. The steps in the remaining two cases will
follow that of the $k=1$ case.

First, assume $d\in\\{2,\dots,2k\\}$. The extension of the $(m+1)^{st}$ branch
uses exclusively the selection rule $H_{m+1}$:

$\underbrace{\sum_{i_{1}=1}^{mk}...\sum_{i_{d-1}=i_{d-2}+1}^{mk}\sum_{i_{d}=i_{d-1}+1}^{mk}\mathbf{Z}^{(d)}_{l_{i_{1}}l^{\prime}_{i_{1}}...l_{i_{d}}l^{\prime}_{i_{d}}}}_{\sum_{p=1}^{m}H_{m}}+\underbrace{\sum_{i_{d}=mk+1}^{(m+1)k}\sum_{i_{d-1}=1}^{i_{d}-1}...\sum_{i_{1}=1}^{i_{2}-1}\mathbf{Z}^{(d)}_{l_{i_{1}}l^{\prime}_{i_{1}}...l_{i_{d}}l^{\prime}_{i_{d}}}}_{H_{m+1}}=\\\
\sum_{i_{1}=1}^{(m+1)k}...\sum_{i_{d-1}=i_{d-2}+1}^{(m+1)k}\sum_{i_{d}=i_{d-1}+1}^{(m+1)k}\mathbf{Z}^{(d)}_{l_{i_{1}}l^{\prime}_{i_{1}}...l_{i_{d}}l^{\prime}_{i_{d}}}$
(18)

Next, assume $d\in\\{2k+1,\dots,mk\\}$. These terms are captured in $I_{m+1}$
using each definition of $\Xi^{(d)}_{m+1}$:

$\sum_{i_{1}=1}^{mk}...\sum_{i_{d-1}=i_{d-2}+1}^{mk}\sum_{i_{d}=i_{d-1}+1}^{mk}\mathbf{Z}^{(d)}_{l_{i_{1}}l^{\prime}_{i_{1}}...l_{i_{d}}l^{\prime}_{i_{d}}}+\underbrace{\sum_{i_{d}=mk+1}^{(m+1)k}\sum_{i_{d-1}=1}^{i_{d}-1}...\sum_{i_{1}=1}^{i_{2}-1}\mathbf{Z}^{(d)}_{l_{i_{1}}l^{\prime}_{i_{1}}...l_{i_{d}}l^{\prime}_{i_{d}}}}_{I_{m+1}}=\\\
\sum_{i_{1}=1}^{(m+1)k}...\sum_{i_{d-1}=i_{d-2}+1}^{(m+1)k}\sum_{i_{d}=i_{d-1}+1}^{(m+1)k}\mathbf{Z}^{(d)}_{l_{i_{1}}l^{\prime}_{i_{1}}...l_{i_{d}}l^{\prime}_{i_{d}}}$
(19)

All four results together show that every degree hyperedge $1,\dots,(m+1)k$ is
accounted for in the $(m+1)^{st}$ branch, thus proving the induction
hypothesis:

$\sum_{p=1}^{m+1}H_{p}+\sum_{p=3}^{m+1}I_{p}=f(X|\mathbf{Z}^{(1)},\mathbf{Z}^{(2)},\dots\mathbf{Z}^{((m+1)k)})$

The $M^{th}$ branch completes the degree $n_{1}$ assignment problem of size
$n_{1}$. For any $k\in\\{1,2,\dots n_{1}\\}$, the selection and aggregation
rules yield the full degree $n_{1}$ assignment problem objective:

$\sum_{p=1}^{M}H_{p}+\sum_{p=3}^{M}I_{p}=f(X|\mathbf{Z}^{(1)},\mathbf{Z}^{(2)},\dots\mathbf{Z}^{(n_{1})})$

∎

### Convergence & Exactness of EHGM

######  Theorem 2.

EHGM (algorithm 1) will converge to a globally optimal solution of the
following hypergraph matching optimization problem given input
$k\in\\{1,2,\dots,n_{1}\\}$:

$\underset{X\in\mathcal{X}}{\text{min}}\sum_{l_{1}=1}^{n_{1}}\sum_{l^{\prime}_{1}=1}^{n_{2}}\mathbf{Z}^{(1)}_{l_{1}l^{\prime}_{1}}x_{l_{1}l^{\prime}_{1}}+\sum_{l_{1}=1}^{n_{1}}\sum_{l^{\prime}_{1}=1}^{n_{2}}\sum_{l_{2}=l_{1}+1}^{n_{1}}\sum_{l^{\prime}_{2}=1}^{n_{2}}\mathbf{Z}^{(2)}_{l_{1}l^{\prime}_{1}l_{2}l^{\prime}_{2}}x_{l_{1}l^{\prime}_{1}}x_{l_{2}l^{\prime}_{2}}\\\
+\sum_{l_{1}=1}^{n_{1}}\sum_{l^{\prime}_{1}=1}^{n_{2}}\sum_{l_{2}=l_{1}+1}^{n_{1}}\sum_{l^{\prime}_{2}=1}^{n_{2}}\sum_{l_{3}=l_{2}+1}^{n_{1}}\sum_{l^{\prime}_{3}=1}^{n_{2}}\mathbf{Z}^{(3)}_{l_{1}l^{\prime}_{1}l_{2}l^{\prime}_{2}l_{3}l^{\prime}_{3}}x_{l_{1}l^{\prime}_{1}}x_{l_{2}l^{\prime}_{2}}x_{l_{3}l^{\prime}_{3}}+...\\\
+\sum_{l_{1}=1}^{n_{1}}\sum_{l^{\prime}_{1}=1}^{n_{2}}...\sum_{l_{n_{1}}=l_{n_{1}-1}+1}^{n_{1}}\sum_{l^{\prime}_{n_{1}}=1}^{n_{2}}\mathbf{Z}^{(n_{1})}_{l_{1}l^{\prime}_{1}\dots
l_{n_{1}}l^{\prime}_{n_{1}}}x_{l_{1}l^{\prime}_{1}}\dots
x_{l_{n_{1}}l^{\prime}_{n_{1}}}$ (20)

where $\mathcal{X}$ is defined:

$\mathcal{X}=\\{X\in\\{0,1\\}^{n_{1}\times n_{2}}:\forall
j,\sum_{i=1}^{n_{1}}x_{ij}\leq 1,\forall i\sum_{j=1}^{n_{2}}x_{ij}=1\\}$ (21)

###### Proof.

First, we will show EHGM converges, then it will be proven that the converged
solution is globally optimal.

The search begins with initializing queue $\mathbf{Q}_{1}=\mathbf{P}$. The
algorithm terminates with the exhaustion of $\mathbf{Q}_{1}$. Each set
$\mathbf{Q}_{m}\subset\mathbf{P}$ contains feasible k-assignments conditioned
on the assignment constraints and costs $\tilde{C},C^{*}$. Backtrack
(Algorithm 4) removes $\mathbf{x}_{m}$ from $\mathbf{Q}_{m-1}$ upon
enumeration of $\mathbf{Q}_{m}$. The recursion then falls back to selecting
from branch $m-1$, eventually exhausting $\mathbf{Q}_{m-1}$ just as in the
enumeration of $\mathbf{Q}_{m}$. This recursion continues until the first
branch $\mathbf{k}_{1}\in\mathbf{Q}_{1}$ is removed, signaling the exploration
of all assignments originating with the k-tuple $\mathbf{k}_{1}$. The
exploration is repeated for each $\mathbf{k}_{1}\in\mathbf{Q}_{1}$. Thus, all
possible assignments $X\in\mathcal{X}$ are explored via the branching scheme.

Assignments accrue a monotonically increasing cost $\tilde{C}$ to be compared
to $C^{*}$ with accompanying assignment $\mathbf{x}^{*}$ at each branch. A
complete assignment then drops the last k assignments from
$\tilde{\mathbf{x}}$, initializing the backwards recursion, emptying
$\mathbf{Q}_{m}$ until $\mathbf{k}^{(n_{k})}_{m-1}\in\mathbf{Q}_{m-1}$ is
exhausted. There are at most $|\mathbf{Q}_{m}|\leq n_{k}$ viable permutations
at branch m. Each possible branch is evaluated from
$\tilde{\mathbf{x}}=[\mathbf{k}_{1},\mathbf{k}_{2},\dots\mathbf{k}_{m-1}]$.
The $(m-1)^{st}$ branch $\mathbf{k}_{m-1}\in\mathbf{Q}_{m-1}$ is removed from
$\mathbf{Q}_{m-1}$ upon exhaustion of $\mathbf{Q}_{m}$:

$\displaystyle\tilde{\mathbf{x}}^{(m)}_{1}$
$\displaystyle=[\mathbf{k}^{(1)}_{1},\mathbf{k}^{(1)}_{2},\dots\mathbf{k}^{(1)}_{m-1},\mathbf{k}^{(1)}_{m}]$
$\displaystyle\tilde{\mathbf{x}}^{(m)}_{2}$
$\displaystyle=[\mathbf{k}^{(1)}_{1},\mathbf{k}^{(1)}_{2},\dots\mathbf{k}^{(1)}_{m-1},\mathbf{k}^{(2)}_{m}]$
$\displaystyle...$ $\displaystyle\tilde{\mathbf{x}}^{(m)}_{n_{k}}$
$\displaystyle=[\mathbf{k}^{(1)}_{1},\mathbf{k}^{(1)}_{2},\dots\mathbf{k}^{(1)}_{m-1},\mathbf{k}^{(n_{k})}_{m}]$

Each of the $n_{k}$ possible final branches from $\mathbf{k}^{(1)}_{m-1}$ is
explored, then $\mathbf{k}^{(1)}_{m-1}$ is removed from $\mathbf{Q}_{m-1}$.

The process follows for the $M^{th}$ branch, exhausting viable assignment sets
until $\mathbf{k}_{M-1}$ is removed. The recursion follows inductively back to
the exhaustion of $\mathbf{Q}_{1}$, signaling the end of the search. Thus, all
possible assignments $X\in\mathcal{X}$ are explored via the branching scheme.

The convergent and exhaustive algorithm will yield a globally optimal solution
$C^{*}=f(\mathbf{x}^{*})$ after exhausting $\mathbf{Q}_{1}$. As proven above
the additive decomposition of the cost structure (equation 7) is proven to be
satisfied by summing all selection and aggregation rule terms. Assume an
uninformed initialization $C^{*}=\infty$. Then the first pass will greedily
take the best permutation from the first set $\mathbf{Q}_{1}$:
$\mathbf{k}^{(1)}_{1}$, and the best from the second set given it does not
conflict with $\mathbf{k}^{(1)}_{1}$: $\mathbf{k}^{(2|1)}_{2}$. This process
will continue until the first complete assignment:
$\tilde{\mathbf{x}}=[\mathbf{k}^{(1)}_{1},\mathbf{k}^{(2|1)}_{2},\mathbf{k}^{(3|2,1)}_{3},\dots,\mathbf{k}^{(M|(M-1),\dots,1)}_{M}]$
with $\tilde{C}=f(\tilde{\mathbf{x}})$. The first Bracktrack removes
$\mathbf{k}^{(M|(M-1),\dots,1)}_{M}$, and the $M^{th}$ Visit call will exhaust
$\mathbf{Q}_{M}$. Subsequent Enqueue calls will limit only allow branches that
satisfy both the assignment constraints and the updated selection rule cost
(Algorithm 3). This follows that any k-tuple of assignments
$\mathbf{k}^{(j)}_{m}$ such that for
$\tilde{\mathbf{x}}=[\mathbf{k}^{(j_{1})}_{1},\mathbf{k}^{(j_{2}|j_{1})}_{2},\dots,\mathbf{k}^{(j_{d}|j_{(m-1)},\dots
j_{1})}_{m}]$:

$\tilde{C}+H_{m}(\tilde{\mathbf{x}}_{m-1},\mathbf{k}^{(j_{d}|j_{(m-1)},\dots
j_{1})}_{m})<C^{*}$

The additive decomposition of the objective paired with the assumed non-
negativity of the dissimilarity tensors $\mathbf{Z}^{(j)}$ results in each
branch monotonically increasing $\tilde{C}$:

$\tilde{C}+H_{m}(\tilde{\mathbf{x}}_{m-1},\mathbf{k}^{(j_{d}|j_{(m-1)},\dots
j_{1})}_{m})+I_{m}(\mathbf{k}^{(j_{1})}_{1},\mathbf{k}^{(j_{2}|j_{1})}_{2},\dots,\mathbf{k}^{(j_{d}|j_{(m-1)},\dots
j_{1})}_{m})\geq\tilde{C}$

The convergent search will thus eliminate all paths that are not globally
optimal. Incrementally updating the reserved solution $\mathbf{x}^{*}$ with
cost $C^{*}$ expedites convergence as each replacement is necessarily a better
solution. The resulting $\mathbf{x}^{*}$ and corresponding cost $C^{*}$ are
such that at no other full assignment $\tilde{\mathbf{x}}$ can replace
$\mathbf{x}^{*}$, by definition a globally optimal solution of f.

∎

### Model Fitting

Expert annotations are used to derive features such that the correct
assignment consistently achieves a minimal cost across the training set.
Features can be engineered and analyzed in context of point set matching just
as in traditional supervised learning tasks.

Features are expressed as attributes over hyperedge multiplicities
$d=1,2,\dots,n_{1}$. Hyperedge features $g^{(d)}_{s}$, $s=1,\dots,n_{d}$ are
given as input. Each feature $g_{s}^{(d)}$ assumes a Gaussian distribution,
and if $n_{d}\geq 2$ the features are modeled as a multivariate Gaussian
distribution. Measurements from the data are used to derive estimates of the
parameters of the Gaussian distributions: $\mathbf{mu}$ and $\mathbf{Sigma}$.
The most common application in heuristic approaches is to use the previous
frame’s feature values as the centers of the distributions. This standard
approach is effective for features that vary minimally, frame-to-frame.
However, certain angle measurements may vary greatly between frames. Mean
estimates across the training data can better account for macroscopic patterns
in features. The variances are then estimated from the feature values across
the training set.

The dissimilarity costs arise from the Mahalanbonis distance between a
hypothesized assignment’s feature measurements and the estimated template mean
values scaled by estimated covariance matrix. The dissimilarity tensors
$\mathbf{Z}^{(d)}$ are expressed as a function of the $n_{d}$ features of
hyperedge d. A partial assignment up to degree d:
$[(l_{1},\dots,l_{d})\mapsto(l^{\prime}_{1},\dots,l^{\prime}_{d})]$ invokes a
cost according to the $n_{d}$ features: $\sum_{s=1}^{n_{d}}g_{s}^{(d)}$. The
expected values: $\sum_{s=1}^{n_{d}}\bar{g}_{s}^{(d)}$ are calculated in
aggregate from training data for higher variance patterns:

$\bar{g}_{s}^{(d)}=\frac{\sum_{L=1}^{N}g_{s}^{(d)}(X_{L},\mathbf{X}_{L})}{N}$
(22)

where $X_{L}$ and $\mathbf{X}_{L}$ are the correct permutation and observed
point set, respectively, for sample $L$. The variance-covariance matrix uses
estimated means to estimate variances and covariances among feature
measurements in the annotated data:

$\hat{\sigma}^{(d)}_{a,b}=\sum_{L=1}^{N}(g_{a}^{(d)}(X_{L},\mathbf{X}_{L})-\bar{g}_{a}^{(d)})(g_{b}^{(d)}(X_{L},\mathbf{X}_{L})-\bar{g}_{b}^{(d)})$
(23)

$\hat{\mathbf{\Sigma}}^{(d)}_{g}=\left|\begin{bmatrix}\hat{\sigma}^{(d)}_{1,1}&\hat{\sigma}^{(d)}_{1,2}&...&\hat{\sigma}^{(d)}_{1,n_{d}}\\\
\hat{\sigma}^{(d)}_{2,1}&\hat{\sigma}^{(d)}_{2,2}&...&\hat{\sigma}^{(d)}_{2,n_{d}}\\\
...&...&...&...\\\
\hat{\sigma}^{(d)}_{n_{d},1}&\hat{\sigma}^{(d)}_{n_{d},2}&...&\hat{\sigma}^{(d)}_{n_{d},n_{d}}\end{bmatrix}\right|$

The selection rule tensor dissimilarity tensors $\mathbf{Z}^{(d)}\in
R^{\underbrace{n\times n,\dots,\times n}_{2d}}$ use both sets of estimates to
compute costs. The Mahalanobis distance is used to describe the scaled
distance between the observed attributed hyperedge to an estimated feature
description. Let
$\mathbf{g}^{(d)}=[g_{1}^{(d)},g_{2}^{(d)},\dots,g_{n_{d}}^{(d)}]^{\prime}$
and
$\bar{\mathbf{g}}^{(d)}=[\bar{g}_{1}^{(d)},\bar{g}_{2}^{(d)},\dots,\bar{g}_{n_{d}}^{(d)}]^{\prime}$

$\mathbf{Z}^{(d)}_{l_{1}l^{\prime}_{1}l_{2}l^{\prime}_{2}\dots
l_{d}l^{\prime}_{d}}=(\mathbf{g}^{(d)}-\bar{\mathbf{g}}^{(d)})^{\prime}(\hat{\mathbf{\Sigma}}^{(d)}_{g})^{-1}(\mathbf{g}^{(d)}-\bar{\mathbf{g}}^{(d)})$
(24)

The traditional approach uses the labeled coordinates in the prior frame to
build corresponding prior frame feature measurements. These prior frame
feature measurements then serve as the estimated center of the Gaussian
distribution. The covariance matrix estimation follows accordingly, in which
the variation in frame-to-frame differences is estimated from sample data.

### Posture Modeling

Embryonic C. elegans posture modeling used the aforementioned template
hypergraph for quantifying hypothesized seam cell identities throughout the
search process. The developing embryo elongates and as a result becomes more
coiled due to the constraining eggshell. As such, the template hypergraphs are
updated according to binned time intervals. Parameters are estimated from data
according to the point in development between first image and hatch.

Each image is captured at time $t$ with $n_{1}=20$ located nuclei centroids.
The coordinates can be stored as $\mathbf{X}\in R^{n\times 3}$ which
$\mathbf{X}_{i}=[x_{i},y_{i},z_{i}]$ representing the $i^{th}$ centroid in
$R^{3}$. The seam cells are ordered posterior to anterior: TL, TR, …, H0L,
H0R. Then let $1,3,5,\dots 19$ be in the indices of the left side, and
$2,4,6,\dots 20$ be the indices of the right side. The seam cells are paired
$(1,2)$ for the tail pair, then $(2,3),(4,5),(6,7),\dots(19,20)$ for the body
pairs. Let $\mathbf{L}=(\mathbf{l}_{1},\mathbf{l}_{2},\dots,\mathbf{l}_{10})$
denote the left side nuclei locations, and similarly $\mathbf{R}$ the right
side nuclei locations.

The sampled worm embryos develop at similar rates. However, the occurrence of
the first twitch, a point in development that triggers rapid physical changes,
varies slightly embryo-to-embryo. As a result, we apply a time normalization
in effort to compare feature measurements from images of the different
specimens. Each embryo’s time to first twitch is measured $s_{w}$, as well as
hatch time $h_{w}$ for sample embryo $w=1,2,\dots 16$. The time points for
each sample are indexed $k=1,2,\dots n_{w}$. Each volume’s imaging time
$t_{wk}$ is normalized to the $[0,1]$ range via Eq 25. Each normalized time
point $z_{wk}$ is scaled such that $z_{wk}=0$ represents first twitch, and
$z_{wk}=1$ hatching.

$z_{wk}=\frac{t_{wk}-s_{w}}{h_{w}-s_{w}}$ (25)

Features used in each of the three models: Sides, Pairs, and Posture are
described with accompanying plots of their distributions throughout both
development and within the embryo. All plots will feature the normalized time
of observation $z_{wk}$ on the horizontal axis. The vertical axis unit will
vary between angular measurements, distances (in $\mu m$), and ratios of each
which have no unit. Plots with multiple subplots in will measure features
calculated at using segments of the embryo posterior to anterior. The left-
most plot will depict calculations using the tail pair while the right-most
will depict the feature ending with the H0 pair at the anterior of the embryo.

#### Sides

The graphical model uses scaled distances between pairs of nuclei. The first
feature to analyze is the distance between paired cell nuclei (S1 Fig 1-A):

$\mathbf{PD}_{i}=\|\mathbf{L}_{i}-\mathbf{R}_{i}\|_{2}$ (26)

The distance between the left and right seam cell nuclei within a pair can be
interpreted as the width of the embryo measured at the sampled seam cell
locations. The first (left-most) subplot of S1 Fig 1-A illustrates distances
in microns between nuclei centroids of the tail pair for each observation. The
tail pair distance is used for the initial pair selection rule $H_{1}$ across
all models. The second set of distances form along the left and right sides of
the worm (S1 Fig 1-B). The lengths of chords between successive nuclei on each
side are calculated: $\|\mathbf{L}_{i+1}-\mathbf{L}_{i}\|_{2}$ and
$\|\mathbf{R}_{i+1}-\mathbf{R}_{i}\|_{2}$. Similar to the pair distances, side
length observations are highly variant.

Supp Fig. 1: Sides model features. A) Distances between nuclei of lateral
pairs. Notably, the tail pair distance (left-most panel) is constant
throughout imaging. The tail pair distance informs the initial pair selection
rule $H_{1}$. B) Chord lengths along left and right sides of the posture. Both
quadratic features show high variance.

#### Pairs

The Pairs model uses hyperedges connecting two or three pairs of nuclei (four
or six nuclei). The first four features measure pair-to-pair variation, while
the latter two features use triplets of pairs to measure angles formed by the
midpoints of the three pairs. The first pair-to-pair feature extends upon the
use of pair distances to better describe the coiled worm. The ratio of
sequential pair distances models the variation in width throughout the
assigned nuclei (Fig 2-A):

$PDR_{i}=\frac{PD_{i}}{PD_{i+1}}$ (27)

Each feature’s estimated mean is slightly greater than $1$, indicating that,
on average, the worm is widening from tail to head. Another easily interpreted
distance is the length of the chords connecting sequential pair midpoints.
This is a more robust measure of worm length as side lengths vary more based
upon the worm’s folding (S1 Fig 2-B):

$MD_{i}=\|\mathbf{M}_{i+1}-\mathbf{M}_{i}\|_{2}$ (28)

The length of the chords connecting sequential pair midpoints is a more robust
measure of worm length as side lengths vary more based upon the worm’s
folding. The cosine similarity is used to assess the degree to which
sequential sides are pointing in the same direction (S1 Fig 2-C):

$\phi_{i}=\frac{(\mathbf{R}_{i+1}-\mathbf{R}_{i})\cdot(\mathbf{L}_{i+1}-\mathbf{L}_{i})}{\|\mathbf{R}_{i+1}-\mathbf{R}_{i}\|_{2}\|\mathbf{L}_{i+1}-\mathbf{L}_{i}\|_{2}}\in[-1,1]$
(29)

The final two pair-to-pair Pairs features attempt to model two different types
of twist in the posture. The lateral and axial twists measures angles of
rotation from lateral and posterior views, respectively (S1 Fig 2-D).

$\displaystyle\mathbf{b}_{1}$
$\displaystyle=\frac{\mathbf{L}_{i+1}-\mathbf{L}_{i}}{\|\mathbf{L}_{i+1}-\mathbf{L}_{i}\|_{2}}$
(30) $\displaystyle\mathbf{b}_{2}$
$\displaystyle=\frac{\mathbf{L}_{i}-\mathbf{R}_{i}}{\|\mathbf{L}_{i}-\mathbf{R}_{i}\|_{2}}$
(31) $\displaystyle\mathbf{b}_{3}$
$\displaystyle=\frac{\mathbf{R}_{i}-\mathbf{R}_{i+1}}{\|\mathbf{R}_{i}-\mathbf{R}_{i+1}\|_{2}}$
(32) $\displaystyle\mathbf{b}_{4}$
$\displaystyle=\frac{\mathbf{R}_{i+1}-\mathbf{L}_{i+1}}{\|\mathbf{R}_{i+1}-\mathbf{L}_{i+1}\|_{2}}$
(33) $\displaystyle\mathbf{n}_{1}$
$\displaystyle=\mathbf{b}_{1}\times\mathbf{b}_{2}$ (34)
$\displaystyle\mathbf{n}_{2}$
$\displaystyle=\mathbf{b}_{2}\times\mathbf{b}_{3}$ (35)
$\displaystyle\mathbf{n}_{3}$
$\displaystyle=\mathbf{b}_{3}\times\mathbf{b}_{4}$ (36)
$\displaystyle\mathbf{c}_{1}$
$\displaystyle=<\mathbf{n}_{1}\times\mathbf{n}_{2},\mathbf{b}_{2}>$ (37)
$\displaystyle\mathbf{c}_{2}$ $\displaystyle=<\mathbf{n}_{1},\mathbf{n}_{2}>$
(38) $\displaystyle\psi_{i}$
$\displaystyle=\frac{1}{\pi}atan2(<\mathbf{n}_{1}\times\mathbf{n}_{2},\mathbf{b}_{2}>,<\mathbf{n}_{1},\mathbf{n}_{2}>))$
(39)

Axial twists present between a sequence of two pairs calculates the angle
obtained by projecting the chord linking pairs onto each other (S1 Fig 2-E):

$\displaystyle\mathbf{b}_{1}$
$\displaystyle=\frac{\mathbf{L}_{i}-\mathbf{L}_{i+1}}{\|\mathbf{L}_{i}-\mathbf{L}_{i+1}\|_{2}}$
(41) $\displaystyle\mathbf{b}_{2}$
$\displaystyle=\frac{\mathbf{R}_{i}-\mathbf{L}_{i}}{\|\mathbf{R}_{i}-\mathbf{L}_{i}\|_{2}}$
(42) $\displaystyle\mathbf{b}_{3}$
$\displaystyle=\frac{\mathbf{R}_{i+1}-\mathbf{R}_{i}}{\|\mathbf{R}_{i+1}-\mathbf{R}_{i}\|_{2}}$
(43) $\displaystyle\mathbf{b}_{4}$
$\displaystyle=\frac{\mathbf{L}_{i+1}-\mathbf{R}_{i+1}}{\|\mathbf{L}_{i+1}-\mathbf{R}_{i+1}\|_{2}}$
(44) $\displaystyle\mathbf{n}_{1}$
$\displaystyle=\mathbf{b}_{1}\times\mathbf{b}_{2}$ (45)
$\displaystyle\mathbf{n}_{2}$
$\displaystyle=\mathbf{b}_{2}\times\mathbf{b}_{3}$ (46)
$\displaystyle\mathbf{n}_{3}$
$\displaystyle=\mathbf{b}_{3}\times\mathbf{b}_{4}$ (47)
$\displaystyle\mathbf{c}_{1}$
$\displaystyle=<\mathbf{n}_{2}\times\mathbf{n}_{3},\mathbf{b}_{3}>$ (48)
$\displaystyle\mathbf{c}_{2}$ $\displaystyle=<\mathbf{n}_{2},\mathbf{n}_{3}>$
(49) $\displaystyle\tau_{i}$
$\displaystyle=\frac{1}{\pi}atan2(<\mathbf{n}_{2}\times\mathbf{n}_{3},\mathbf{b}_{3}>,<\mathbf{n}_{2},\mathbf{n}_{3}>))$
(50)

Angles along sides of the worm formed by triples of sequential nuclei
approximate bend in the worm along each side. These bend angles are highly
variant, especially frame-to-frame, in the same manner as side lengths in
Sides (S1 Fig 1-B). Angles formed by pair midpoints exacerbate the
computational burden as six nuclei are required, compare to three in a typical
angle calculation, but the midpoint based angles are less variant than angles
of each side (S1 Fig 2-F):

$\Theta_{i}=\frac{180}{\pi}\arccos{\frac{<\mathbf{M}_{i+1}-\mathbf{M}_{i},\mathbf{M}_{i+2}-\mathbf{M}_{i+1}>}{\|\mathbf{M}_{i+1}-\mathbf{M}_{i}\|_{2}\|\mathbf{M}_{i+2}-\mathbf{M}_{i+1}\|_{2}}}$
(52)

Each angle $\Theta\in[0,180]$ where $0$ would denote the worm perfectly folded
upon itself, and $180$ would define a flat worm. A second set of angles aims
to approximate the posterior to anterior bend in the worm. The angles
$\zeta_{i}$ are defines as the angles formed by fitted planes intersecting
between pair midpoints (S1 Fig 2-G):

$\zeta_{i}=\frac{180}{\pi}\frac{<(\mathbf{R}_{i+1}-\mathbf{L}_{i+1})\times(\mathbf{M}_{i+1}-\mathbf{M}_{i}),((\mathbf{R}_{i+1}-\mathbf{L}_{i+1})\times(\mathbf{M}_{i+1}-\mathbf{M}_{i+1}))>}{\|(\mathbf{R}_{i+1}-\mathbf{L}_{i+1})\times(\mathbf{M}_{i+1}-\mathbf{M}_{i})\|_{2}\|(\mathbf{R}_{i+1}-\mathbf{L}_{i+1})\times(\mathbf{M}_{i+1}-\mathbf{M}_{i+1})\|_{2}}$
(53)

Supp Fig. 2: Pairs model features. A) Ratios of pair distances (Eq 27). B)
Distance between successive pair midpoints (Eq 28). C) Cosine similarities
between successive left and right sides (Eq 29). D) Lateral axial twist angles
(Eq 30). E) Axial twist angles (Eq 41). F) Midpoint bend angles (Eq 52). G)
Planar intersection angles (Eq 53).

#### Posture

The Posture model is comprised of all Pairs features as well as the features
defined by the summations of each local feature measurement throughout the
hypothesized posture. Full posture features give insight into the changes in
the embryo’s shape throughout late-stage embryogenesis (S1 Fig 3). Worm length
follows an approximately logarithmic pattern.Total curvature follows a
negative exponential pattern. Earlier on the worm is fatter and cannot bend as
much. The worm elongates during development, allowing for sharper bends.

Supp Fig. 3: Posture model features include all Pairs features and posture-
wide versions of Pairs features. A) Summed ratios of pair distances (Eq 27).
B) Summed distances between successive pair midpoints (Eq 28). C) Summed
cosine similarities between successive left and right sides (Eq 29). D) Summed
lateral axial twist angles (Eq 30). E) Summed axial twist angles (Eq 41). F)
Summed midpoint bend angles (Eq 52). G) Summed planar intersection angles (Eq
53).