# Inverting the discrete curl operator: a novel graph algorithm to find a
vector potential of a given vector field

Silvano Pitassi<EMAIL_ADDRESS>Riccardo Ghiloni
<EMAIL_ADDRESS>Ruben Specogna<EMAIL_ADDRESS>University
of Udine, Polytechnic Department of Engineering and Architecture, EMCLab, via
delle scienze 206, 33100 Udine, Italy  University of Trento, Department of
Mathematics, via Sommarive 14, 38123 Povo-Trento, Italy

###### Abstract

We provide a novel framework to compute a discrete vector potential of a given
discrete vector field on arbitrary polyhedral meshes. The framework exploits
the concept of acyclic matching, a combinatorial tool at the core of discrete
Morse theory. We introduce the new concept of _complete acyclic matchings_ and
we show that they give the same end result of Gaussian elimination. Basically,
instead of doing costly row and column operations on a sparse matrix, we
compute equivalent cheap combinatorial operations that preserve the underlying
sparsity structure. Currently, the most efficient algorithms proposed in
literature to find discrete vector potentials make use of tree-cotree
techniques. We show that they compute a special type of complete acyclic
matchings. Moreover, we show that the problem of computing them is equivalent
to the problem of deciding whether a given mesh has a topological property
called _collapsibility_. This fact gives a topological characterization of
well-known termination problems of tree-cotree techniques. We propose a new
recursive algorithm to compute discrete vector potentials. It works directly
on basis elements of $1$\- and $2$-chains by performing elementary Gaussian
operations on them associated with acyclic matchings. However, the main
novelty is that it can be applied recursively. Indeed, the recursion process
allows us to sidetrack termination problems of the standard tree-cotree
techniques. We tested the algorithm on pathological triangulations with known
topological obstructions. In all tested problems we observe linear
computational complexity as a function of mesh size. Moreover, the algorithm
is purely graph-based so it is straightforward to implement and does not
require specialized external procedures. We believe that our framework could
offer new perspectives to sparse matrix computations.

###### keywords:

compatible discretizations, polyhedral meshes, vector potential, source fields

††journal: Journal of LaTeX Templates

##  1 Introduction

In this paper we present a novel framework to solve the discrete version of
the following potential problem: determine a vector field with specified curl,
i.e. a vector potential of a given vector field. To do so, an “inverse curl”
operator is required, necessarily including the specification of additional
conditions (a gauge condition in physical parlance) to uniquely define a
resulting solution. Let us now formally introduce the problem by first
considering the continuous case and then moving towards its discrete version,
where we define discrete counterparts of vector fields and differential
operators acting on them.

Let us consider a bounded domain $\Omega$ of $\mathbb{R}^{3}$, where a vector
field $\bm{J}$ is defined. Assume that $\Omega$ is topologically trivial,
i.e., it is homeomorphic to a closed $3$-ball (or, equivalently, to a cube).
We first determine necessary and sufficient conditions for assuring that
$\bm{J}$ is the curl of a vector field $\bm{H}$, i.e. $\bm{H}$ is a vector
potential of $\bm{J}$. The answer is well-known, being a classical problem in
vector analysis. A vector field $\bm{J}$ is the curl of a vector field if and
only if its divergence is zero and its flux is vanishing across all the (but
one) connected components of $\partial\Omega$. In our case, $\partial\Omega$
is connected, so the flux condition can be omitted since it is automatically
verified.

Let us consider a _mimetic discretization_ of our continuous problem. Mimetic
discretization methods like the Mimetic Finite Difference method (MFD) [1],
Discrete Geometric Approach (DGA) [2], Finite Integration Technique (FIT),
Discrete de Rham (DDR) methods [3] include the structure of exterior calculus,
thus retaining fundamental properties of the continuous theory.

We recall the main concepts of mimetic discretization in Section 2. We cover
the domain $\Omega$ of $\mathbb{R}^{3}$ with a polyhedral mesh $\mathcal{K}$,
namely a regular CW cell complex having cells being polyhedra. Continuous
vector fields and differential operators are replaced by their discrete
counterparts. We define a discrete vector field as a collection of _degrees of
freedom_ (DoFs). We introduce the discrete vector field $\bm{h}$, an array of
DoFs with values on each edge of $\mathcal{K}$; next, the discrete vector
field $\bm{i}$, an array of DoFs with values on each face of $\mathcal{K}$.
The central problem of this paper is to find a discrete vector potential
$\bm{h}$ such that

$\mathbb{C}\bm{h}=\bm{i},$ (1)

where $\mathbb{C}$ is the usual incidence matrix between the faces and the
edges of the mesh. In intimate analogy with the continuous case, a necessary
and sufficient condition to have a consistent discrete vector potential
$\bm{h}$ is that array $\bm{i}$ represents a _discrete solenoidal_ vector
field, i.e. it verifies $\mathbb{D}\bm{i}=\bm{0}$, where $\mathbb{D}$ is the
incidence matrix between cells and faces.

Our main motivation to solve problem Eq. 1 stems from the fact that this
algorithmic primitive is an enabling technology for solving many problems
arising in computational physics, from electromagnetism to elasticity and
fluid mechanics [4]. First, it can be used to solve the vector laplacian in
nearly linear time [5]. The idea is that, instead formulating the vector
laplacian by using a vector potential, the scalar potential [6] or the mixed-
hybrid [2] formulations could be used instead, which produce linear systems
that can be solved in nearly linear time by using algebraic multigrid methods.
Second, inverse discrete curl is at the root of efficient algorithms to
compute a cohomology basis and source fields for solving magnetostatics and
eddy current problems by mimetic or finite element methods [7, 8, 9, 10]. We
think $\bm{h}$ as a discrete magnetic field and $\bm{i}$ as a discrete
current; then Eq. 1 expresses the so-called discrete Ampère’s law. In the
electromagnetic literature, discrete fields $\bm{h}$ satisfying Eq. 1 are
often called _source fields_. A different application in computational
electromagnetics is to find a magnetic vector potential from a magnetic
induction field [11].

There are two analogous discrete potential problems: the problem of
determining a scalar potential with assigned gradient and a vector field with
assigned divergence. However, as we will see in our discourse, and as been
already pointed out in literature [12], these two problems are less
challenging than problem Eq. 1, since they are easily solved in linear worst-
case complexity using standard spanning tree constructions.

Many different algorithms have been proposed in literature to solve linear
system Eq. 1.

A naive solution to the problem would be to solve Eq. 1 by a linear system
solver. However, these techniques are not feasible from the practical point of
view since they show cubic worst-case computational complexity. Instead, the
most efficient methods to solve linear system Eq. 1 are based on the so-called
_tree-cotree decomposition_. Tree-cotree decomposition arises from graph
theory and consists in partitioning the edges of a graph into a _spanning
tree_ and its complement, referred to as the cotree. The basic idea, rooted in
the works [6, 7], goes as follows. First, the values of array $\bm{h}$
relative to a _spanning tree_ on the vertices and edges of $\mathcal{K}$ are
set to zero. Next, discrete Ampère’s law is iteratively enforced on each face
of the mesh. This principle is the core of the algorithms proposed in [13,
12].

Tree-cotree techniques are frequently claimed to be general in the literature.
However, proofs about their termination are not discussed at all [6, 7]. In
fact, such techniques are not guaranteed to converge. A careful analysis of
the termination properties can be found in [13], where it is shown that they
strongly depend on the choice of the spanning tree. In particular, there exist
situations involving a topologically trivial complex (for example, certain
meshes of a cube or a 3-ball) and a spanning tree on which such techniques do
not terminate.

To overcome termination issues of tree-cotree techniques, different approaches
have been proposed. The approach in [14] is based on the idea of symbolic
computations. Although being general, this approach is slower than previous
approaches. Moreover, it is difficult to implement in practice since it
requires specialized data structures to manage the symbolic computations; for
instance, it requires object oriented programming languages to implement the
symbolic computations. In [12], termination issues are solved using an
explicit formula based on a double integral computation. However, these double
integral computations require specialized algorithms that are time consuming.

The aim of this paper is to describe an efficient and easy to implement
algorithm to solve linear system Eq. 1. Our novel algorithm is based on
concepts of _discrete Morse theory_ [15]. This theory employs a construction
called _acyclic matching_ [16], which collects combinatorial operations
analogous to topological operations in the continuum. In Section 3 we review
basic concepts at the core of discrete Morse theory as well as our specialized
definitions.

Starting from acyclic matchings of discrete Morse theory, we derive three main
contributions.

Firstly, we provide a unified framework based on discrete Morse theory for the
solution of linear system Eq. 1 and we present in Section 4. We introduce the
novel concept of _complete acyclic matching_. The crucial fact of our
framework is the following: complete acyclic matchings give the same end
result of Gaussian elimination. What the whole procedure boils down to is
avoiding costly matrix algebra operations by performing equivalent cheap
combinatorial operations. By this procedure, the sparsity of linear system Eq.
1 is retained, whereas this would not be the case in the standard Gaussian
elimination. In fact, sparse systems become dense in intermediate steps due to
an inconvenient choice of the backward elimination algorithm.

Secondly, using the newly introduced framework, we show that tree-cotree
decomposition techniques are algorithms to compute specialized complete
acyclic Morse matchings. Although they use a different language to describe
the same actions, they provide the same end results. We show that the question
of finding this kind of specialized complete acyclic matchings is equivalent
to the following _collapsibility problem_ : decide whether a 3-dimensional
simplicial complex embedded in $\mathbb{R}^{3}$ with trivial topology and with
Lipschitz boundary is _collapsible_. This fact establish the topological
nature of termination problems of tree-cotree techniques and shows why they
arise in practical applications. Indeed, there exist examples of
triangulations of 3-balls that are not collapsible [17]. As a result, for
certain triangulations, avoiding termination problems by a careful choice of
the input spanning tree of tree-cotree techniques is an impossible task and
thus we need to resort to new approximation strategies.

Thirdly, we provide a new recursive algorithm to compute discrete vector
potentials we describe in Section 5. It consists of a new greedy heuristics to
construct acyclic matchings together with a recursive construction, which has
no analogous in classical discrete Morse theory. Indeed, we do not define a
new chain complex, the so-called _Morse complex_ , but instead we employ
acyclic matchings to express new basis elements during Gaussian elimination in
terms of the original basis elements of the vector space of 1-chains and the
vector space of 2-chains. Fundamentally, we do not tackle the problem of
computing specialized complete acyclic matchings like in tree-cotree
techniques, which, as discussed above, suffers from well-known topological
termination problems. Instead, we introduce a recursive approach whose basic
outline is as follows. We first try to find a complete acyclic matching on
linear system Eq. 1. If we do not succeed, we transform linear system Eq. 1
into a new smaller linear system, by considering a suitable subset of the
newly computed basis during Gaussian elimination. We show that if we can solve
this new smaller linear system, then we can get a solution of linear system in
Eq. 1. The crucial fact is that we can recursively apply the algorithm on the
new smaller linear system. More specifically, we try to find a complete
acyclic matching on the corresponding new smaller linear system. Again,
acyclic matchings provides a new basis during Gaussian elimination so that the
whole approach can be applied recursively. The fundamental computational
advantage of our recursive approach is that all the costly algebraic
operations are replaced by elementary cheap combinatorial operations on
significantly smaller instances of the original linear system Eq. 1.

We tested our recursive algorithm on challenging benchmark problems. These
include also pathological meshes with known topological obstructions. Our
algorithm exhibits linear computational complexity for all tested problems.
Moreover, the algorithm is also purely graph-theoretic, so straightforward to
implement and do not require any additional specialized procedure. We collect
simulation results in Section 6.

While obtaining a theoretical linear worst-case complexity bound is hard, our
algorithm solves all issues that are typically found in practice. Moreover,
the whole framework proves to be general and offers new perspectives to sparse
matrix computations. These observations are summarized in Section 7.

##  2 Notation

The domain of interest of this paper is a closed and bounded polyhedral domain
$\Omega$ of $\mathbb{R}^{3}$ with Lipschitz boundary. We assume that $\Omega$
has trivial topology, i.e., it is homeomorphic to a closed $3$-dimensional
ball or, equivalently, $\Omega$ is simply connected and its boundary
$\partial\Omega$ is connected ($\partial\Omega$ is homeomorphic to a
$2$-sphere indeed); see [18] (Section 6) and [19] (Section 3). We consider a
_polyhedral cell complex_ (or _polyhedral mesh_) subdivision $\mathcal{K}$ of
$\Omega$. Elements of $\mathcal{K}$ are called _cells_. A $k$-cell $\sigma$ is
a $k$-dimensional subset in $\mathbb{R}^{3}$ homeomorphic to a closed
$k$-dimensional ball. A $0$-cell is a point of $\mathbb{R}^{3}$. We equip each
$k$-cell with an inner orientation. We denote by $\dim{\sigma}$ the dimension
of the cell $\sigma\in\mathcal{K}$ and we write $\mathcal{K}_{k}$ the
subcollection of all $k$-cells in $\mathcal{K}$. We focus on the 3-dimensional
case, thus we have 3-cells (or volumes) in $\mathcal{K}_{3}$, 2-cells (or
faces) in $\mathcal{K}_{2}$, 1-cells (or edges) in $\mathcal{K}_{1}$ and
0-cells (or vertices) in $\mathcal{K}_{0}$. We denote by $c$ a generic volume,
by $f$ a face, by $e$ an edge and by $v$ a vertex. We denote by
$\mathfrak{v},\mathfrak{e},\mathfrak{f},\mathfrak{c}$ the cardinality of
$\mathcal{K}_{0},\mathcal{K}_{1},\mathcal{K}_{2},\mathcal{K}_{3}$,
respectively. A polyhedral cell complex $\mathcal{K}$ is _simplicial_ if all
its cells are simplicies and the boundary of each cell has the natural
simplicial decomposition, see [20]. If $\mathcal{K}$ is simplicial,
$\mathcal{K}$ is also called a triangulation of $\Omega$.

The mesh $\mathcal{K}$ has the structure of a _(regular) cell complex_ ,
namely, the following three conditions hold [20]. First, for each $k$-cell
$\sigma$ in $\mathcal{K}$ its _boundary_ $\partial\sigma$ is a union of
$(k-1)$-cells in $\mathcal{K}$ for $k\in\\{1,2,3\\}$ [20]. Second, given
distinct $k$-cells $\sigma,\tau$, their intersection $\sigma\cap\tau$ is
either empty or is a union of lower dimensional cells in $\mathcal{K}$. Third,
given a $l$-cell $\sigma$ and $k$-cell $\tau$ with $l\leq k$, $\sigma\neq\tau$
and $\sigma\cap\tau\neq\emptyset$, we have
$\sigma\cap\tau\subset\partial\tau$.

We can now define a new object, called a real $k$-chain. A $k$-chain $a$ of
$\mathcal{K}$ is a formal linear combination of $k$-cells
$a=\sum_{i=1}^{r}a_{i}\sigma_{i}$, where $\sigma_{i}$ are $k$-cells in
$\mathcal{K}$ and $a_{i}$ are real coefficients. The number $r$ denotes the
cardinality of the collection of $k$-cells in $\mathcal{K}$ and is any number
among $\mathfrak{v},\mathfrak{e},\mathfrak{f}$ or $\mathfrak{c}$. The set of
$k$-chains, equipped with the natural addition and scalar multiplication,
provides a real vector space. We denote it by $C_{k}(\mathcal{K})$. Note that
each $k$-cell is also a $k$-chain. If $\sigma$ is a $k$-cell, by $-\sigma$ we
denote the cell $\sigma$ but with opposite orientation. The set of all
$k$-cells form a basis for $C_{k}(\mathcal{K})$, which we call _canonical
basis_ for $C_{k}(\mathcal{K})$. We identify the boundary of each $k$-cell
with the linear combination of the $(k-1)$-cells in $\partial\sigma$ defined
by setting

$\partial\sigma\coloneqq\sum_{i=1}^{r}w_{i}\,\rho_{i},$ (2)

where $w_{i}$ is different from zero if and only if
$\rho_{i}\subset\partial\sigma$ and in this case, $w_{i}$ is equal to $+1$ if
$\rho_{i}$ has the orientation induced by that of $\sigma$ by using the right-
hand rule and $-1$ otherwise [21].

The real vector space of $k$-chains and the real vector space of
$(k-1)$-chains are connected by a linear map called _boundary operator_
$\partial_{k}:C_{k}(\mathcal{K})\to C_{k-1}(\mathcal{K})$. We define the
boundary operator by linearity on the space of chains by setting

$\partial_{k}a\coloneqq\sum_{i=1}^{r}a_{i}\,\partial\sigma_{i},$ (3)

where $a=\sum_{i=1}^{r}a_{i}\sigma_{i}$ as above. Note that Eq. 3 is well-
defined since $\mathcal{K}$ is a cell complex. Since $\mathcal{K}$ is a cell
complex, it can verified that $\partial_{k-1}\circ\partial_{k}=0$ for
$k\in\\{1,2,3\\}$; see, for example [20].

Let us now consider the concept of a _$k$ -cochain_. A $k$-cochain $b$ acts on
a $k$-chain to produce a real number and therefore $k$-cochains are elements
of the dual space of $C_{k}(\mathcal{K})$. We define the vector space of
$k$-cochains $C^{k}(\mathcal{K})$ to be the dual space of linear functionals
$b:C_{k}(\mathcal{K})\to\mathbb{R}$. We denote the value of a $k$-chain $a$
under a $k$-cochain $b$ as $\langle{b}\,,{a}\rangle\coloneqq b(a)$. Let us
consider the canonical basis $\\{\sigma_{i}\in C_{k}(\mathcal{K})\mid
i=1,\dots,r\\}$ of the vector space of $k$-chains $C_{k}(\mathcal{K})$. From
basic linear algebra, there exist unique linear functionals $\\{\sigma^{i}\in
C^{k}(\mathcal{K})\mid i=1,\dots,r\\}$ such that

$\langle{\sigma^{i}}\,,{\sigma_{j}}\rangle=\delta_{ij},$ (4)

where $\delta_{ij}$ is the Kronecker delta. The set $\\{\sigma^{i}\in
C^{k}(\mathcal{K})\mid i=1,\dots,r\\}$ defined by Eq. 4 form a basis for the
vector space of $k$-cochains $C^{k}(\mathcal{K})$, which is called _canonical
dual basis_. We have established a one-to-one correspondence between chains
and cochains. This chain-cochain natural duality yields the real linear
isomorphism $\phi_{k}:C_{k}(\mathcal{K})\to C^{k}(\mathcal{K})$ sending each
$\sigma_{i}$ to $\sigma^{i}$. We will write a generic $k$-cochain $b\in
C^{k}(\mathcal{K})$ as a sum $b=\sum_{i=1}^{r}b_{i}\sigma^{i}$ with real
coefficients $b_{i}$.

For $k$-cochains, in intimate analogy with chains, we can define a _coboundary
operator_ $\delta^{k}:C^{k}(\mathcal{K})\to C^{k+1}(\mathcal{K})$ as the _dual
of the boundary operator_ , i.e., it is defined by requiring that, for every
$b\in C^{k}(\mathcal{K})$ and $a\in C_{k+1}(\mathcal{K})$, the following
identity holds

$\langle{\delta^{k}b}\,,{a}\rangle=\langle{b}\,,{\partial_{k+1}a}\rangle.$ (5)

In mimetic methods the coboundary operator $\delta^{k}$ acts as a discrete
counterpart of the continuous differential operators [22, 21]. Specifically,
$\delta^{0}$ acts as the discrete gradient, $\delta^{1}$ as the discrete curl
and $\delta^{2}$ as the discrete divergence.

A straightforward calculation using Eq. 5 shows that
$\delta^{k}\circ\delta^{k-1}=0$ for $k\in\\{1,2\\}$. These relations mimic the
structure of continuous differential operators [22, 21]. In particular,
discrete differential operators form a chain complex

$C^{*}=\cdots\xrightarrow{\delta^{k-1}}C^{k}(\mathcal{K})\xrightarrow{\delta^{k}}C^{k+1}(\mathcal{K})\xrightarrow{\delta^{k+1}}\cdots,$
(6)

where $C^{k}(\mathcal{K})=0$ if $k<0$ or $k>3$. Since the domain $\Omega$ is
topologically trivial the sequence is _exact_ for $k\neq 0$, i.e. it satisfies
$\mathrm{im}(\delta^{k-1})=\mathrm{ker}(\delta^{k})$ for $k\neq 0$.

In the case of $k$-chains, there is a natural choice of a basis given by the
canonical basis. Using the isomorphism $\phi_{k}:C_{k}(\mathcal{K})\to
C^{k}(\mathcal{K})$, we have also fixed a canonical dual basis for
$C^{k}(\mathcal{K})$. Since the coboundary operator $\delta^{k}$ is a linear
map between $C^{k}(\mathcal{K})$ and $C^{k+1}(\mathcal{K})$, it can be
represented, using the fixed bases of $C^{k}(\mathcal{K})$ and
$C^{k+1}(\mathcal{K})$, as a matrix. Thus, to represent the coboundary
operator as a matrix, we must always explicitly state which bases are chosen
and, in fact, we will soon see the benefits of changing the bases.

Let us now consider an arbitrary basis for the vector space of $k$-chains
$C_{k}(\mathcal{K})$. We denote it by
$\mathcal{B}_{k}=\\{\dots,\xi_{i},\dots\\}$. Using the isomorphism
$\phi_{k}:C_{k}(\mathcal{K})\to C^{k}(\mathcal{K})$,
$\phi_{k}(\mathcal{B}_{k})=\\{\dots,\xi^{i}=\phi_{k}(\xi_{i}),\dots\\}$ is a
basis for $C^{k}(\mathcal{K})$. In what follows, we take this process of
choosing a basis for $C^{k}(\mathcal{K})$ for granted. When this is done, we
say that we have chosen a basis $\mathcal{B}=\bigcup_{k}\mathcal{B}_{k}$ for
the entire chain complex $C^{*}$, i.e. a basis for each $C^{k}(\mathcal{K})$.
We write $(C^{*},\mathcal{B})$ to denote a chain complex with a basis. We
define the _canonical basis_ $\widehat{\mathcal{B}}$ for $C^{*}$ to be the
basis $\widehat{\mathcal{B}}=\bigcup_{k}\widehat{\mathcal{B}}_{k}$ where each
$\widehat{\mathcal{B}}_{k}$ is the canonical basis for $C_{k}(\mathcal{K})$.

Having chosen bases in $C_{k}(\mathcal{K})$ and hence in $C^{k}(\mathcal{K})$,
we denote by $\mathbb{D}_{k}$ the matrix associated with $\delta^{k}$ for
$k\in\\{0,1,2\\}$. We introduce a more common notation from mimetic methods
$\mathbb{G}\coloneqq\mathbb{D}_{0}$, $\mathbb{C}\coloneqq\mathbb{D}_{1}$ and
$\mathbb{D}\coloneqq\mathbb{D}_{2}$. We can also define $\mathbb{D}_{3}$ as
the null operator from $C_{3}(\mathcal{K})$ to $0$.

We represent $k$-chains and $k$-cochains by vectors of size $r$ that contain
the real numbers with respect to the ordered bases. A $k$-chain
$a=\sum_{i=1}^{r}a_{i}\xi_{i}$ in a basis $\\{\dots,\xi_{i},\dots\\}$ is
represented by the column vector $\mathbb{R}^{r}\ni\bm{a}=(a_{1}\cdots
a_{r})^{T}$ and a $k$-cochain $b=\sum_{i=1}^{r}b_{i}\xi^{i}$ in the basis
$\\{\dots,\xi^{i}=\phi_{k}(\xi_{i}),\dots\\}$ is represented by the column
vector $\mathbb{R}^{r}\ni\bm{b}=(b_{1},\cdots b_{r})^{T}$.

As explained in the introduction, the aim of this paper is to devise an
efficient algorithm for the solution of the following problem: find array
$\bm{h}\in\mathbb{R}^{\mathfrak{e}}$ such that

$\mathbb{C}\bm{h}=\bm{i},$ (7)

where $\bm{i}\in\mathbb{R}^{\mathfrak{f}}$ satisfy $\mathbb{D}\bm{i}=\bm{0}$.

The rank of the matrix $\mathbb{C}$ is not maximal, thus linear system Eq. 7
has an infinite number of solutions. In fact, if $\bm{h}$ is a solution of Eq.
7, then array $\bm{h}^{\prime}=\bm{h}+\mathbb{G}\bm{\psi}$ is also a solution
of Eq. 7 since we have

$\mathbb{C}\bm{h}^{\prime}=\mathbb{C}(\bm{h}+\mathbb{G}\bm{\psi})=\mathbb{C}\bm{h}=\bm{i},$
(8)

where $\bm{\psi}\in\mathbb{R}^{\mathfrak{v}}$ and we have used the chain
complex property $\mathbb{C}\mathbb{G}=0$ of Eq. 6.

In what follows, the concept of a partition of a given index set will play an
important role. Let $I$ be a finite index set. A partition of $I$ is a family
of disjoint subsets $\\{I_{1},\dots,I_{p}\\}$ of $I$ such that
$\bigcup_{l=1}^{p}I_{l}=I$. The _subvector of $\bm{v}=(v_{i})_{i\in
I}\in\mathbb{R}^{I}$ induced by $I_{l}$_ is

${{\bm{v}}|_{I_{l}}}\coloneqq(v_{i})_{i\in I_{l}},$ (9)

for $l\in\\{1,\dots,p\\}$. A representation of the vector $\bm{v}$ as a block
vector is given by

$\bm{v}=({{\bm{v}}|_{I_{l}}})_{l\in\\{1,\dots,p\\}}.$ (10)

Let us consider a product of index sets $I$ and $J$. We will need a
corresponding notion of Eq. 9 for a matrix whose entries are indexed by
elements in $I\times J$. Let us consider partitions of $I$ and $J$ as
$\\{I_{1},\dots,I_{p}\\}$ and $\\{J_{1},\dots,J_{q}\\}$, respectively. We have
a corresponding partition of $I\times J$ as a family of disjoint subsets
$\\{O_{1},\dots,O_{n}\\}$ such that $O_{i}=I_{l}\times J_{m}$ for some
$l\in\\{1,\dots,p\\}$, $m\in\\{1,\dots,q\\}$ and $I\times
J=\bigcup_{i=1}^{n}O_{i}$. The _submatrix of
$\mathbb{A}=(\mathbb{A}_{i,j})_{i\in I,j\in J}\in\mathbb{R}^{I\times J}$
induced by $I_{l}\times J_{m}$_ is

${{\mathbb{A}}|_{I_{l}\times J_{m}}}\coloneqq(\mathbb{A}_{i,j})_{i\in
I_{l},j\in J_{m}}\in\mathbb{R}^{I_{l}\times J_{m}}.$ (11)

A representation of the matrix $\mathbb{A}$ as a block matrix is given by

$\mathbb{A}=({{\mathbb{A}}|_{I_{l}\times
J_{m}}})_{l\in\\{1,\dots,p\\},m\in\\{1,\dots,q\\}}.$ (12)

##  3 Discrete Morse Theory

The underlying principle of our construction follows an ad hoc reformulation
of Forman’s Discrete Morse theory [15] given by Kozlov [16], where the basic
tool is a combinatorial object called acyclic matching. Several special cases
of our construction have already appeared in literature. We present a
formulation of discrete Morse theory adapted to our purposes, along with
smaller, more illustrative instances, which will provide insights on the
structure of our algorithm.

###  3.1 Informal introduction to discrete Morse theory

The first concept is that of _elementary collapse_. One may view discrete
Morse theory as a generalization of the theory of simplicial collapses. The
concept of collapse, originated in Whitehead’s work [23], provides a
combinatorial operation that is analogous to the continuous operation called
deformation retraction, i.e., the operation of continuously shrinking a
topological space to a subset. More specifically, let $(\sigma,\tau)$ be a
pair of cells such that $\sigma\subset\tau$ and $\dim\sigma=\dim\tau-1$. For
this pair, to induce an elementary collapse, we require $\tau$ to be a cell of
maximal dimension in $\mathcal{K}$ and the only one cell of $\mathcal{K}$
containing $\sigma$; we refer to this as saying that the pair $(\sigma,\tau)$
is _free_ in $\mathcal{K}$. Equivalently, we also say that $\sigma$ is free in
$\mathcal{K}$; see Fig. 1(a). We say that $\mathcal{K}$ _collapses to_
$\mathcal{L}$ if one could get from $\mathcal{K}$ to $\mathcal{L}$ in a finite
sequence of elementary collapses. If $\mathcal{K}$ is equivalent to a single
vertex, then we say that $\mathcal{K}$ is _collapsible_ ; in this case there
exists a sequence of elementary collapses leaving a single vertex.

Dropping the uniqueness condition on $\tau$, we obtain what we refer to as an
_internal collapse_ , see Fig. 1(b).

Geometrically, in both cases, we obtain a collapse of the pair $(\sigma,\tau)$
by contracting the whole cell $\tau$ onto $\partial\tau\setminus\sigma$.

Figure 1: (a) Elementary collapse of free pair $(\sigma,\tau)$. (b) Internal
collapse of pair $(\sigma,\tau)$; the resulting cell complex is not more
simplicial.

In intimate analogy with elementary collapses, we may combine many internal
collapses to form a sequence of internal collapses, again without affecting
the homotopy type.

We thus have a family of pairs
$\\{(\sigma_{1},\tau_{1}),\dots,(\sigma_{n},\tau_{n})\\}$ to be collapsed, in
this order. One may view the set of all such pairs as a _matching_ on
$\mathcal{K}$. Accordingly, we refer to cells contained in some pair as
_matched_ and other cells as _unmatched_ or _critical_ (with respect to the
matching).

Let $\mathcal{K}^{\,(i)}$ be the resulting cell complex after the first $i$
collapses. For the pairs to form a sequence of elementary collapses, we
require that each new pair $(\sigma_{i},\tau_{i})$ is free in
$\mathcal{K}^{\,(i-1)}$. For generic collapses we apply the same requirement,
except that we restrict our attention to the family of matched cells.
Specifically, we do not require $(\sigma_{i},\tau_{i})$ to be free in
$\mathcal{K}^{\,(i-1)}$, but $\tau_{i}$ must be the only matched cell of
$\mathcal{K}^{\,(i-1)}$ containing $\sigma_{i}$. Equivalently, for each $i$,
we should have that $\sigma_{i}$ is not contained in
$\tau_{i+1},\dots,\tau_{n}$ for $i\in\\{1,\dots,n-1\\}$. We refer to a
matching on $\mathcal{K}$ admitting an ordering with this property as
_acyclic_. We formalize all this concepts in Section 3.2.

The main theorem of discrete Morse theory states that an acyclic matching
induces a homotopy equivalence between $\mathcal{K}$ and the so-called Morse
complex, a cell complex formed by critical cells only [16] (Theorem 11.13
(b)).

###  3.2 Acyclic matchings

We start our exposition by examining acyclic matchings from a purely
combinatorial point of view without any reference to topology. Indeed, our
interest is in using discrete Morse theory to develop a fast algorithm for the
solution of linear system Eq. 7 to be applied to cell complexes arising from
experimental or numerical meshes of real case problems. We give a specific
version of combinatorial discrete Morse theory by Kozlov [16] that is adapted
to our purposes.

For any $\sigma,\tau\in\widehat{\mathcal{B}}_{k}$, define
$\langle{\sigma}\,,{\tau}\rangle$ to be $1$ if $\sigma=\tau$ and $0$
otherwise. Extend, by linearity, $\langle{\cdot}\,,{\cdot}\rangle$ to a scalar
product on $C_{k}(\mathcal{K})$. Note that we can identify the scalar product
$\langle{\cdot}\,,{\cdot}\rangle$ with the duality product between chains and
cochains in Section 2 via the isomorphism $\phi_{k}:C_{k}(\mathcal{K})\to
C^{k}(\mathcal{K})$ in Eq. 4, i.e.
$\langle{\sigma}\,,{\tau}\rangle=\langle{\phi_{k}(\sigma)}\,,{\tau}\rangle=\phi_{k}(\sigma)(\tau)$.

Let us consider the chain complex $(C^{*},\mathcal{B})$ with basis
$\mathcal{B}$. We define a relation $\prec$ on $\mathcal{B}$ as follows. Given
distinct basis elements $\sigma\in\mathcal{B}_{k}$ and
$\tau\in\mathcal{B}_{k+1}$,

$\sigma\prec\tau\iff\langle{\sigma}\,,{\partial_{k+1}\tau}\rangle\neq 0.$ (13)

If $\sigma\prec\tau$, then we say that $\sigma$ and $\tau$ are _incident_.

We introduce the _boundary set_ and _coboundary set_ of
$\sigma\in\mathcal{B}_{k}$ as

$\mathrm{bd}_{\mathcal{B}}{(\sigma)}\coloneqq\\{\,\rho\in\mathcal{B}_{k-1}\mid\rho\prec\sigma\\},$
(14)

and

$\mathrm{cobd}_{\mathcal{B}}{(\sigma)}\coloneqq\\{\,\rho\in\mathcal{B}_{k+1}\mid\sigma\prec\rho\\},$
(15)

respectively.

Let $\sigma\in\mathcal{B}_{k}$. If the cardinality of
$\mathrm{cobd}_{\mathcal{B}}{(\sigma)}$ is one, then we say that $\sigma$ is
_free_. In this case, there exists a unique basis element $\tau$ such that
$\sigma\prec\tau$, and we also say that the pair $(\sigma,\tau)$ is free. If
the cardinality of $\mathrm{cobd}_{\mathcal{B}}{(\sigma)}$ is greater than
one, then we say that $\sigma$ in _internal_. In this case, if
$\tau\in\mathrm{cobd}_{\mathcal{B}}{(\sigma)}$, then we also say that the pair
$(\sigma,\tau)$ is internal.

######  Definition 3.1 (Matching, acyclic matching).

A _matching_ $\mathcal{M}$ on $\mathcal{B}$ is a family of pairs
$\\{(\sigma,\tau)\\}$ with $\sigma,\tau\in\mathcal{B}$ such that:

  1. 1.

$(\sigma,\tau)\in\mathcal{M}$ implies $\sigma\prec\tau$.

  2. 2.

each $\sigma\in\mathcal{B}$ is the first component of at most one pair
$(\sigma,\tau)$ in $\mathcal{M}$.

A matching $\mathcal{M}$ is called _acyclic_ if there does not exist a cycle

$\tau_{1}\succ\sigma_{1}\prec\tau_{2}\succ\cdots\prec\tau_{h}\succ\sigma_{h}\prec\tau_{1},$
(16)

with $h\geq 2$, $(\sigma_{i},\tau_{i})\in\mathcal{M}$ for all
$i\in\\{1,\dots,h\\}$ and all $\tau_{i}\in\mathcal{B}$ being distinct.

A matching $\mathcal{M}_{k}$ of $k$-chains on $\mathcal{B}$ is a matching such
that if $(\sigma,\tau)\in\mathcal{M}_{k}$ then $\sigma\in\mathcal{B}_{k}$.

The following result is a reformulation of Theorem 11.2 in [16] by Kozlov. It
describes the crucial combinatorial property that characterizes acyclic
matchings. Its proof can be obtained by a suitable adaption of the mentioned
Theorem 11.2, see pages 181-182 of [16].

######  Theorem 1.

A matching $\mathcal{M}$ on $\mathcal{B}$ is acyclic if and only if there
exists a total order of pairs of $\mathcal{M}$ as
$\\{(\sigma_{1},\tau_{1}),\dots,(\sigma_{n},\tau_{n})\\}$ such that, for every
$i\in\\{1,\dots,n-1\\}$, $\sigma_{i}$ is not incident to any
$\tau_{i+1},\dots,\tau_{n}$.

In what follows, we will write an acyclic matching $\mathcal{M}$ as a sequence
$\\{(\sigma_{1},\tau_{1}),\dots,(\sigma_{n},\tau_{n})\\}$, where it is
understood that the total order is chosen according to Theorem 1.

Let $\mathcal{M}$ be a matching on $\mathcal{B}$. We say that a basis element
in $\mathcal{B}$ is _matched_ (with respect to $\mathcal{M}$) if it is
contained in some pair in $\mathcal{M}$ (both as first or second component of
the pair), otherwise, we say that it is _unmatched_ or _critical_ (with
respect to $\mathcal{M}$).

We denote by $\mathcal{U}_{k}$ the set of all $\tau\in\mathcal{B}_{k}$ such
that $\tau$ is matched with some $\sigma\in\mathcal{B}_{k-1}$. Similarly, we
denote by $\mathcal{D}_{k}$ the set of all $\sigma\in\mathcal{B}_{k}$ such
that $\sigma$ is matched with some $\tau\in\mathcal{B}_{k+1}$. Given
$\mathcal{D}_{k}$ and $\mathcal{U}_{k}$, there is a corresponding set of
_critical_ $k$-chains

$\mathcal{C}_{k}\coloneqq\mathcal{B}_{k}\setminus(\mathcal{D}_{k}\cup\mathcal{U}_{k}).$
(17)

Finally, we set $\mathcal{U}\coloneqq\bigcup_{k}\mathcal{U}_{k}$,
$\mathcal{D}\coloneqq\bigcup_{k}\mathcal{D}_{k}$ and
$\mathcal{C}\coloneqq\bigcup_{k}\mathcal{C}_{k}$. It easy to see that the sets
$\mathcal{U},\mathcal{D},\mathcal{C}$ provide a partition of $\mathcal{B}$

$\mathcal{B}=\mathcal{U}\cup\mathcal{D}\cup\mathcal{C}.$ (18)

Our definition of a matching is related to the presentation of the
combinatorial Morse theory of Forman [15] and in particular the more recent
formulation given by Kozlov [16]. In earlier presentations, elements in
$\mathcal{U}$ and $\mathcal{D}$ are not explicitly introduced since what is
important is only the bijective pairing between their elements. Instead, in
our setting they will play a fundamental role since we use discrete Morse
theory from a purely combinatorial point of view and elements in $\mathcal{U}$
and $\mathcal{D}$ will be used to select suitable submatrices. The set of
critical elements is present also in classical discrete Morse theory. The set
of critical elements will play a fundamental role in Algorithm 3, where
critical elements become the new input for subsequent iterations. In Algorithm
3, it is also essential to be able to express new basis elements as a function
of the previous ones. Hence, we need to keep track of basis structure at each
iteration.

An important difference with classical discrete Morse theory is that in Forman
[15] a new chain complex, the so-called Morse complex, is constructed. Our
version of discrete Morse theory operates directly on basis elements by
performing elementary operations on it. This is needed to describe how
incidence matrices transform in the new bases due to our combinatorial
operations.

###  3.3 Basis transformations associated with an acyclic matching

We begin by considering simple examples to develop geometric intuition behind
the general definitions. When simplifying a cell complex, the effect of a
collapse is that of changing the structure of the basis $\mathcal{B}$, by
performing elementary operations on it. There are three elementary operations
to obtain a new basis from a previous one. If
$\mathcal{B}_{k}=\\{\dots,\xi_{i},\xi_{j},\dots\\}$ is a basis of
$C_{k}(\mathcal{K})$, then a new basis may be obtained by

  1. 1.

Exchanging elements $\xi_{i}$ and $\xi_{j}$.

  2. 2.

Multiplying $\xi_{i}$ by $-1$.

  3. 3.

Replacing $\xi_{j}$ by $\xi_{j}+q\xi_{i}$ with $q\in\mathbb{R}$.

Let us consider a prototype example of a simplicial complex $\mathcal{K}$ in
Fig. 2. The set of all 1-chains is generated by the canonical basis
$\widehat{\mathcal{B}}_{1}=\\{e_{1},e_{2},e_{3},e_{4},e_{5}\\}$ and the set of
all 2-chains by the canonical basis
$\widehat{\mathcal{B}}_{2}=\\{f_{1},f_{2}\\}$. In Fig. 2 edges
$e_{1},e_{2},e_{3},e_{4}$ are free since each of them is incident to exactly
one face in $\mathcal{K}$. Instead, edge $e_{5}$ is internal since is the
common edge of $f_{1}$ and $f_{2}$. Depending on whether a collapse is
elementary or internal we have corresponding elementary operations on
$\mathcal{B}$.

Figure 2: (a) The simplicial complex $\mathcal{K}$. (b) Elementary collapse of
the free pair $(e_{1},f_{2})$. (c) Internal collapse of the pair
$(e_{5},f_{2})$; note that the geometric realization of the resulting cell
complex is not more simplicial.

######  Example 2 (Elementary collapse).

Let us consider the elementary collapse of the free pair $(e_{1},f_{2})$. The
obtained cell complex in Fig. 2(b) is generated by the set of critical basis
elements $\mathcal{C}$. In fact, we have
$\mathcal{C}_{1}=\\{e_{2},e_{3},e_{4},e_{5}\\}$ and
$\mathcal{C}_{2}=\\{f_{1}\\}$. We get a new basis $\mathcal{B}^{\prime}$ of
$\mathcal{K}$ as $\mathcal{B}^{\prime}_{1}=\mathcal{D}_{1}\cup\mathcal{C}_{1}$
and $\mathcal{B}^{\prime}_{2}=\mathcal{U}_{2}\cup\mathcal{C}_{2}$, where
$\mathcal{D}_{1}=\\{e_{1}\\}$ and $\mathcal{U}_{2}=\\{f_{2}\\}$. We see that,
in the case of an elementary collapse, we get a new basis
$\mathcal{B}^{\prime}$ of $\mathcal{K}$ by performing elementary operations of
type 1.

######  Example 3 (Interior collapse).

Let us now consider the collapse of the pair $(e_{5},f_{2})$. Contrary to the
previous case, edge $e_{5}$ is not free, so we cannot consider an elementary
collapse of the pair $(e_{5},f_{2})$. However, we can collapse $(e_{5},f_{2})$
as an internal collapse. The obtained cell complex in Fig. 2(c) is not
generated by critical basis elements in $\mathcal{C}$ as in the previous case.
In fact, we have $\mathcal{C}_{1}=\\{e_{1},e_{2},e_{3},e_{4}\\}$ and
$\mathcal{C}_{2}=\\{f_{1}\\}$. Instead, it is generated by a new set of basis
elements obtained from $\mathcal{C}$ by adding a linear combination of other
basis elements. We consider the linear transformation $\mathcal{C}_{2}\ni
f_{1}\mapsto f_{1}+f_{2}$. The new set of critical basis elements is
$\mathcal{C}_{1}=\\{e_{1},e_{2},e_{3},e_{4}\\}$ and
$\mathcal{C}_{2}=\\{f_{1}+f_{2}\\}$. We get a different basis
$\mathcal{B}^{\prime}$ of the cell complex $\mathcal{K}$ as
$\mathcal{B}^{\prime}_{1}=\mathcal{D}_{1}\cup\mathcal{C}_{1}$ and
$\mathcal{B}^{\prime}_{2}=\mathcal{U}_{2}\cup\mathcal{C}_{2}$, where
$\mathcal{D}_{1}=\\{e_{5}\\}$ and $\mathcal{U}_{2}=\\{f_{2}\\}$. We see that,
in the case of an internal collapse, we get a new basis $\mathcal{B}^{\prime}$
of $\mathcal{K}$ by performing elementary operations of type 1, 2 and 3.

Let $\mathcal{M}_{k}=\\{(\sigma_{1},\tau_{1}),\dots,(\sigma_{n},\tau_{n})\\}$
be an acyclic matching of $k$-chains on $\mathcal{B}$. We now give a recursive
definition of the _change of basis associated with $\mathcal{M}_{k}$_ and we
denote it by $\mathcal{B}\cdot\mathcal{M}_{k}$. To start with, given the
matched pair $(\sigma,\tau)$, we define the basis
$\mathcal{B}\cdot(\sigma,\tau)$ as follows. $\mathcal{B}\cdot(\sigma,\tau)$ is
obtained from $\mathcal{B}$ by performing two actions. First, we consider the
partition of $\mathcal{B}=\mathcal{U}\cup\mathcal{D}\cup\mathcal{C}$, in
particular, we have

$\mathcal{B}_{k}=\mathcal{D}_{k}\cup\mathcal{C}_{k},$ (19)
$\mathcal{B}_{k+1}=\mathcal{U}_{k+1}\cup\mathcal{C}_{k+1}.$ (20)

Second, the pair $(\sigma,\tau)$ acts on the set of critical elements
$\mathcal{C}$ as follows

$\mathcal{C}_{k}\ni\sigma^{\prime}\mapsto\sigma^{\prime},$ (21)
$\mathcal{C}_{k+1}\ni\tau^{\prime}\mapsto\tau^{\prime}-\frac{\langle{\sigma}\,,{\partial_{k+1}\tau^{\prime}}\rangle}{\langle{\sigma}\,,{\partial_{k+1}\tau}\rangle}\tau.$
(22)

We see that $\mathcal{B}\cdot(\sigma,\tau)$ is again a basis of $C^{*}$.
Indeed, it is obtained from $\mathcal{B}$ by adding linear combinations of
other basis elements. We define $\mathcal{B}\cdot\mathcal{M}_{k}$ recursively
by the rule

$\displaystyle\mathcal{B}\cdot\\{(\sigma_{1},\tau_{1}),\dots,(\sigma_{n},\tau_{n})\\}\coloneqq\big{(}\mathcal{B}\cdot\\{(\sigma_{1},\tau_{1}),\dots,(\sigma_{n-1},\tau_{n-1})\\}\big{)}\cdot(\sigma_{n},\tau_{n}),$
(23)

with $n\geq 2$. We see that Eq. 23 is well-defined. Indeed, if
$\mathcal{B}\cdot\\{(\sigma_{1},\tau_{1}),\dots,(\sigma_{n-1},\tau_{n-1})\\}$
is basis, then
$\mathcal{B}\cdot\\{(\sigma_{1},\tau_{1}),\dots,(\sigma_{n},\tau_{n})\\}$ is
obtained by adding linear combinations of basis elements, hence it is a basis.
Moreover, $\mathcal{M}_{k}$ is an acyclic matching on
$\mathcal{B}\cdot\\{(\sigma_{1},\tau_{1}),\dots,(\sigma_{n-1},\tau_{n-1})\\}$,
since transformations Eq. 21, Eq. 22 leaves matched elements in $\mathcal{B}$
invariant.

Note that, when $(\sigma,\tau)$ is a free pair, the transformation Eq. 22
leaves elements in $\mathcal{C}_{k+1}$ invariant. Hence, no algebraic
operations are required and the new basis is simply a permutation of the
previous one. Instead, when $(\sigma,\tau)$ is internal, we have to apply at
least one transformation Eq. 22.

##  4 Acyclic matchings and Gaussian elimination

We shall now head towards an algorithm to reduce the matrix $\mathbb{C}$ to a
row echelon form by means of elementary operations on the basis $\mathcal{B}$.
The whole procedure boils down to a standard train of thought used in basic
linear algebra. When the coboundary operator $\delta^{1}$ is given as a finite
matrix $\mathbb{C}$, the bases and orders are already determined. However, we
can get any other bases by applying elementary row and column operations on
the matrix $\mathbb{C}$. Our algorithm will produce a change of bases in such
a way that the matrix $\mathbb{C}$ in the new bases has an invertible upper
triangular submatrix (i.e. a matrix with non-zeros only in its upper triangle
and main diagonal) and thus it can be transformed in row echelon form. Hence,
we can fast solve the system by processing the unknown variables in reverse
order, a standard process known as back substitution.

###  4.1 Acyclic matchings and Gaussian elimination

The crucial observation is the following general novel result.

######  Lemma 1.

Denote by $\mathbb{D}_{k}$ one among the matrices $\mathbb{G}$, $\mathbb{C}$
or $\mathbb{D}$ for $k$ equal to $0$,$1$ or $2$, respectively. Let
$\mathcal{M}_{k}$ be an acyclic matching of $k$-chains on $\mathcal{B}$. Then,
${{\mathbb{D}_{k}}|_{\mathcal{U}_{k+1}\times\mathcal{D}_{k}}}$, the submatrix
of $\mathbb{D}_{k}$ induced by $\mathcal{U}_{k+1}\times\mathcal{D}_{k}$, is
upper triangular and invertible.

###### Proof.

Since the matching
$\mathcal{M}_{k}=\\{(\sigma_{1},\tau_{1}),\dots,(\sigma_{n},\tau_{n})\\}$ is
acyclic, for every $i\in\\{1,\dots,n-1\\}$, the basis element $\sigma_{i}$ is
not incident to any basis element $\tau_{i+1},\dots,\tau_{n}$. Thus, all the
non-zero entries in each column are above the diagonal as the rows and columns
are arranged in the total order induced by the matching. ∎

By Lemma 1, for every acyclic matching there is a corresponding upper
triangular submatrix
${{\mathbb{D}_{k}}|_{\mathcal{U}_{k+1}\times\mathcal{D}_{k}}}$ of matched rows
and columns of $\mathbb{D}_{k}$. However, Lemma 1 is decisive only when the
number of matched pairs in $\mathcal{M}_{k}$ is equal to the rank of
$\mathbb{D}_{k}$. In fact, in this case,
${{\mathbb{D}_{k}}|_{\mathcal{U}_{k+1}\times\mathcal{D}_{k}}}$ is an
invertible submatrix of $\mathbb{D}_{k}$ of order equal to the rank
$\mathbb{D}_{k}$, hence, we can write a solution of Eq. 7 after setting some
free variables to zero. For this reason, we introduce the new concept of
complete acyclic matching.

######  Definition 4.1 (Complete acyclic matching).

We say that a matching $\mathcal{M}_{k}$ of $k$-chains on $\mathcal{B}$ is
_complete_ if the number of matched pairs in $\mathcal{M}_{k}$ is equal to the
rank of $\mathbb{D}_{k}$.

Now, we illustrate how the upper triangular submatrix
${{\mathbb{D}_{k}}|_{\mathcal{U}_{k+1}\times\mathcal{D}_{k}}}$ induced by a
complete acyclic matching allows us to operationally obtain a solution of Eq.
7 using back substitution.

Let us consider a complete acyclic matching $\mathcal{M}_{1}$ of 1-chains. By
applying Lemma 1, the action of $\mathcal{M}_{1}$ on Eq. 7 is equivalent to
Gaussian elimination. It is thus sufficient to invert submatrix
${{\mathbb{C}}|_{\mathcal{U}_{2}\times\mathcal{D}_{1}}}$ after setting free
variables in $\mathcal{C}_{1}$ to zero. We write a discrete potential $\bm{h}$
solution of Eq. 7 as

$\displaystyle{{\bm{h}}|_{\mathcal{D}_{1}}}={{\mathbb{C}}|_{\mathcal{U}_{2}\times\mathcal{D}_{1}}}^{-1}{{\bm{i}}|_{\mathcal{U}_{2}}},$
(24) $\displaystyle{{\bm{h}}|_{\mathcal{C}_{1}}}=\bm{0}.$ (25)

Since ${{\mathbb{C}}|_{\mathcal{U}_{2}\times\mathcal{D}_{1}}}$ is upper
triangular, we can evaluate
${{\mathbb{C}}|_{\mathcal{U}_{2}\times\mathcal{D}_{1}}}^{-1}{{\bm{i}}|_{\mathcal{U}_{2}}}$
by back substitution in linear time [24].

The aim of the back substitution is to determine the coefficient values of
$\bm{h}$. Let $\mathcal{M}_{1}=\\{(e_{1},f_{1}),\dots,(e_{n},f_{n})\\}$. The
process of back substitution is so-called because one determines the
coefficient values backwards, by first computing $\bm{h}_{e_{n}}$, then
substituting back into the previous equation to solve for $\bm{h}_{e_{n-1}}$
and repeating through $\bm{h}_{e_{n-1}}$. A naive combination of these
coefficients in one step for each $\bm{h}_{e_{i}}$ is very time consuming
since large intermediate expressions are generated. This is avoided by
combining coefficients pairwise as in the following standard back substitution
algorithm.

Algorithm 1 Back substitution process

1:procedure Backsubstitution($\mathcal{M}_{1},\mathcal{B}$)

2: for $i$ ranging from $n$ down to $1$ do

3:
$\bm{h}_{e_{i}}=(\mathbb{C}_{f_{i},e_{i}})^{-1}(\bm{i}_{f_{i}}-\sum_{e\in\mathrm{bd}_{\mathcal{B}}{(f_{i})}}\mathbb{C}_{f_{i},e}\bm{h}_{e})$

4: end for

5: return $\bm{h}$;

6:end procedure

In the above discussion, although we focused on matrix $\mathbb{C}$, we have
actually detailed the proof of the following theorem.

######  Theorem 4.

Given a complete acyclic matching $\mathcal{M}_{k}$ of $k$-chains on
$\mathcal{B}$, for $k\in\\{0,1,2\\}$ we can find in linear time a discrete
potential $\bm{v}\in C^{k}(\mathcal{K})$ of $\bm{w}\in C^{k+1}(\mathcal{K})$
with $\mathbb{D}_{k+1}\bm{w}=\bm{0}$, namely a solution of
$\mathbb{D}_{k}\bm{v}=\bm{w}$.

###  4.2 On the problem of constructing complete acyclic matchings: the case
of tree-cotree techniques

We now turn to the theoretical issue of constructing complete acyclic
matchings $\mathcal{M}_{k}$ of $k$-chains on the canonical basis
$\widehat{\mathcal{B}}$.

Case $k=0$. We need to construct a complete acyclic matching $\mathcal{M}_{0}$
of 0-chains on $\widehat{\mathcal{B}}$. By definition, the number of matched
elements in $\mathcal{M}_{0}$ has to be equal to the rank of $\mathbb{G}$,
which is $\mathfrak{v}-1$. Let us consider a _spanning tree_ $T$ on
$\mathcal{K}$. It can be constructed in worst-case linear time using standard
graph algorithms, for instance _breadth-first search_ (BFS) algorithm [25].
There is a standard reasoning to define an acyclic matching corresponding to a
spanning tree $T$ [26]. It is obtained by mimicking a spanning tree traversal
process. The construction goes as follows. Since $T$ is a tree, there exists
at least a _leaf_ in $T$, namely a vertex $v\in\mathcal{K}_{0}$ with only one
incident edge in $T$. Pick a leaf $v$ and pair it with the unique edge $e$
containing it. Then, add the pair $(v,e)$ to $\mathcal{M}_{0}$ and remove
$v,e$ from $T$. Since $v$ is a leaf, the obtained graph is again a tree. By
iteratively repeating this process, we define a complete acyclic matching
$\mathcal{M}_{0}$. The matching $\mathcal{M}_{0}$ is well-defined: if a vertex
$v$ is matched during the process, it is removed from $T$ and thus cannot
appear in any other pair; it is acyclic since every tree is. The matching
$\mathcal{M}_{0}$ is complete, since $\mathcal{K}$ is connected, $T$ is
spanning, namely all vertices of $\mathcal{K}$ are in $T$ and they are
eventually added to $\mathcal{M}_{0}$ in the above process, except the last
one.

Case $k=2$. Before considering the case $k=1$, we show the similarity of the
construction for $k=2$ with the case $k=0$. We need to construct a complete
acyclic matching $\mathcal{M}_{2}$ of 2-chains on basis
$\widehat{\mathcal{B}}$. By definition, the number of matched elements in
$\mathcal{M}_{2}$ has to equal to the rank of $\mathbb{D}$, which is
$\mathfrak{c}$. Note that the cell complex $\mathcal{K}$ defines a manifold
with boundary. It implies that if $c,c^{\prime}$ are distinct volumes and $f$
is a face such that $f\subset\partial c\cap\partial c^{\prime}$, then
$c,c^{\prime}$ are the only volumes that contain $f$. Thanks to the manifold
condition, it is well-defined the so-called _complete dual graph_ of
$\mathcal{K}$. The dual graph
$\tilde{G}=(\tilde{V}\cup\tilde{V}_{\infty},\tilde{E}\cup\tilde{E}_{\infty})$
of $\mathcal{K}$ is a graph with set of vertices $\tilde{V}$ given by volumes
of $\mathcal{K}$ and $\\{c,c^{\prime}\\}\in\tilde{E}$ if $\dim(c\cap
c^{\prime})=2$. There is an additional vertex
$\tilde{v}_{\infty}\in\tilde{V}_{\infty}$ and there are additional edges
$\tilde{E}_{\infty}$ in $\tilde{G}$. Each edge in $\tilde{E}_{\infty}$
corresponds to a volume $c$ whose boundary $\partial c$ contains a boundary
face $f$; thus, if $f$ is a boundary face and $c$ is the unique volume
incident to it, then $\\{\tilde{v}_{\infty},c\\}\in\tilde{E}_{\infty}$. Let us
consider a spanning tree $\tilde{T}$ on $\tilde{G}$. By repeating the same
construction detailed for $k=0$ on $\tilde{T}$, we obtain a finite sequence
$S$ of pairs of the form $(c,\tilde{e})$, where $c$ is a vertex in
$\tilde{V}\cup\tilde{V}_{\infty}$ (i.e. a volume of $\mathcal{K}$ or
$\tilde{v}_{\infty}$) and $\tilde{e}=\\{c,c^{\prime}\\}$ is an edge in
$\tilde{E}\cup\tilde{E}_{\infty}$. For each pair $(c,\tilde{e})$, we choose a
face $f$ of $\mathcal{K}$ such that $f\subset\partial c\cap\partial
c^{\prime}$ if $\tilde{v}_{\infty}\notin\tilde{e}=\\{c,c^{\prime}\\}$,
$f\subset\partial c$ if $c^{\prime}=\tilde{v}_{\infty}$ and $f\subset\partial
c^{\prime}$ if $c=\tilde{v}_{\infty}$. Replace in $S$ each pair
$(c,\tilde{e})$ with either $(f,c)$ if $c\neq\tilde{v}_{\infty}$ or
$(f,c^{\prime})$ if $c=\tilde{v}_{\infty}$. In this way, we obtain a complete
acyclic matching $\mathcal{M}_{2}$ of 2-chains of $\mathcal{K}$. Thus, all
vertices of the dual graph are eventually added to $\mathcal{M}_{2}$, except
the vertex $\tilde{v}_{\infty}$, and all volumes in $\mathcal{K}$ are matched.

By combining Theorem 4 with the above results, we state the following theorem
which solves the discrete potential problem for $k\in\\{0,2\\}$.

######  Theorem 5.

For $k\in\\{0,2\\}$, we can find in linear time a discrete potential
$\bm{v}\in C^{k}(\mathcal{K})$ of $\bm{w}\in C^{k+1}(\mathcal{K})$ with
$\mathbb{D}_{k+1}\bm{w}=\bm{0}$, namely a solution of
$\mathbb{D}_{k}\bm{v}=\bm{w}$.

Case $k=1$. With no surprise, it turns out to be the most challenging case. To
begin with, we show the following result.

######  Theorem 6.

Let $\mathcal{K}$ be a 3-dimensional topologically trivial simplicial complex
embedded in $\mathbb{R}^{3}$ with Lipschitz boundary as in Section 2. Then,
$\mathcal{K}$ is collapsible if and only if there exists a complete acyclic
matching $\mathcal{M}_{1}$ of 1-chains on the canonical basis
$\widehat{\mathcal{B}}$.

###### Proof.

Let us consider a sequence of elementary collapses $\mathcal{M}$ of
$\mathcal{K}$ leading from $\mathcal{K}$ to a vertex. We define an acyclic
matching $\mathcal{M}_{1}$ of 1-chains (on $\widehat{\mathcal{B}}$) by
selecting, according to the total order of $\mathcal{M}$, all pairs made by
matched edges and faces in $\mathcal{M}$. Note that $\mathcal{M}_{1}$ is
necessarily acyclic since $\mathcal{M}$ is. Next, since $\mathcal{K}$
collapses to a vertex, all faces in $\mathcal{K}$ are matched in
$\mathcal{M}$. In particular, all volumes in $\mathcal{K}$ are matched in
$\mathcal{M}$ and there are $\mathfrak{c}$ corresponding matched faces. Thus,
there are $\mathfrak{f}-\mathfrak{c}$ matched faces in $\mathcal{M}_{1}$. It
follows that $\mathcal{M}_{1}$ is complete since the rank of $\mathbb{C}$ is
exactly $\mathfrak{f}-\mathfrak{c}$.

To prove the converse result, let us consider a complete acyclic matching of
1-chains $\mathcal{M}_{1}$.

We need the following construction. First let us consider the complete dual
graph
$\tilde{G}=(\tilde{V}\cup\tilde{V}_{\infty},\tilde{E}\cup\tilde{E}_{\infty})$
of $\mathcal{K}$ as defined above for the case $k=2$. Consider the subgraph
$\tilde{G}_{\mathcal{M}_{1}}=(\tilde{V}\cup\tilde{V}_{\infty},\tilde{E}_{\mathcal{M}_{1}})$
of $\tilde{G}$ which has the same vertices of $\tilde{G}$ and
$\tilde{E}_{\mathcal{M}_{1}}$ includes an edge
$\\{c,c^{\prime}\\}\in\tilde{E}$ of $\tilde{G}$ if the unique face $f$ of
$\mathcal{K}$ such that $f\subset c\cap c^{\prime}$ is critical with respect
to $\mathcal{M}_{1}$. We now show that $\tilde{G}_{\mathcal{M}_{1}}$ has
$\mathfrak{c}$ edges (i.e, the number of vertices of $\tilde{G}$ minus one)
and is connected. First, since $\mathcal{M}_{1}$ is a complete acyclic
matching the number of critical faces, which is equal to the number of edges
in $\tilde{E}_{\mathcal{M}_{1}}$, is
$\mathfrak{f}-(\mathfrak{f}-\mathfrak{c})=\mathfrak{c}$. Second, assume that,
for the sake of contradiction, $\tilde{G}_{\mathcal{M}_{1}}$ is disconnected.
Let $\tilde{V}^{*}$ be the set of nodes in a fixed connected component of
$\tilde{G}_{\mathcal{M}_{1}}$. Let
$\tilde{E}_{\tilde{V}^{*}}\subset\tilde{E}\cup\tilde{E}_{\infty}$ be the
subset of _cut edges_ , namely edges of $\tilde{G}$ with one vertex in
$\tilde{V}^{*}$ and one vertex in its complement
$(\tilde{V}\cup\tilde{V}_{\infty})\setminus\tilde{V}^{*}$. Note that since
$\tilde{G}$ is connected, $\tilde{E}_{\tilde{V}^{*}}$ is not empty. Moreover,
by definition of $\tilde{G}_{\mathcal{M}_{1}}$, to each
$\tilde{e}\in\tilde{E}_{V^{*}}$ corresponds a matched face $f$ of
$\mathcal{K}$ with respect to $\mathcal{M}_{1}$ forming the pair
$(e,\tilde{e})=(e,f)\in\mathcal{M}_{1}$ for a unique edge $e$ of
$\mathcal{K}$.

Now we construct a cycle of the form $\tilde{e}_{1}\succ
e_{1}\prec\tilde{e}_{2}\succ\dots\prec\tilde{e}_{h}\succ
e_{h}\prec\tilde{e}_{1}$ with $h\geq 2$,
$(e_{i},\tilde{e}_{i})\in\mathcal{M}_{1}$ for all $i\in\\{1,\dots,h\\}$ and
all $\tilde{e}_{i}\in\tilde{E}_{V^{*}}$ being distinct. Start with a cut edge
$\tilde{e}_{1}$ and consider the unique edge $e_{1}$ of $\mathcal{K}$ forming
the matched pair $(e_{1},\tilde{e}_{1})\in\mathcal{M}_{1}$. Then, $e_{1}$ is
contained in at least one other cut edge $\tilde{e}_{2}\in\tilde{E}_{V^{*}}$,
otherwise $\tilde{e}_{1}$ cannot be a cut edge. By iteratively repeating this
process, we obtain a sequence $S$ of pairs
$(e_{1},\tilde{e}_{1}),\dots,(e_{i},\tilde{e}_{i})$ such that
$\tilde{e}_{1}\succ e_{1}\prec\tilde{e}_{2}\succ\dots\prec\tilde{e}_{i-1}\succ
e_{i}\prec\tilde{e}_{i}$. Since we have finite graphs, eventually we will run
out of those cut edges which do not appear as the second component of a pair
in $S$. Thus, at some step $i$ of the above process, we must get a cut edge
that is the second component of a pair in $S$ and we get a cycle of the form
Eq. 16, which is a contradiction since $\mathcal{M}_{1}$ is acyclic. We have
proved that $\tilde{G}_{\mathcal{M}_{1}}$ is a connected subgraph of
$\tilde{G}$ with $\mathfrak{c}$ edges and thus it is a spanning tree of
$\tilde{G}$. We define a complete acyclic matching $\mathcal{M}_{2}$ of
2-chains by using the spanning tree construction detailed above for $k=2$ on
$\tilde{G}_{\mathcal{M}_{1}}$.

Now, $\mathcal{M}_{2}$ and $\mathcal{M}_{1}$ define a sequence of elementary
collapses leading from $\mathcal{K}$ to a single vertex as follows. Since
every volume of $\mathcal{K}$ is matched with respect to $\mathcal{M}_{2}$,
the set of critical cells of $\mathcal{K}$ with respect to $\mathcal{M}_{2}$
form a 2-dimensional subcomplex $\mathcal{K}^{\,(1)}$ of $\mathcal{K}$. By
applying Theorem 11.13 (a) in [16], there exists a sequence of elementary
collapses leading from $\mathcal{K}$ to $\mathcal{K}^{\,(1)}$. Next, observe
that, by construction, the set of matched faces in $\mathcal{K}$ with respect
to $\mathcal{M}_{2}$ is exactly the set of critical faces of $\mathcal{K}$
with respect to $\mathcal{M}_{1}$. Thus, the set of critical cells of
$\mathcal{K}^{\,(1)}$ with respect to $\mathcal{M}_{1}$ form a 1-dimensional
subcomplex $\mathcal{K}^{\,(2)}$ of $\mathcal{K}^{\,(1)}$. By applying Theorem
11.13 (a) in [16] there exists a sequence of elementary collapses leading from
$\mathcal{K}^{\,(1)}$ to $\mathcal{K}^{\,(2)}$. $\mathcal{K}^{\,(2)}$ is a
topologically trivial 1-dimensional subcomplex, thus it is a spanning tree on
$\mathcal{K}$. Thus $\mathcal{K}$ collapses to a single vertex. ∎

Now we focus on the _Spanning Tree Technique_ (STT) algorithm introduced in
[13] and implicitly used in many works [7, 8, 12, 6]. The algorithm, given as
input a spanning tree $T$ on $\mathcal{K}$, computes a discrete vector
potential $\bm{h}$ as follows. To start with, STT sets the value
$\bm{h}_{e}=0$ if $e\in T$, otherwise the value $\bm{h}_{e}$ at this stage is
unknown. Next, all faces $f$ of $\mathcal{K}$ are loaded into a list $L$. The
main loop of STT works until there are no more faces in $L$. In each
iteration, we randomly search for a face $f$ with two boundary edges
$e_{1},e_{2}\subset\partial f$ such that the values
$\bm{h}_{e_{1}},\bm{h}_{e_{2}}$ are known. Then, the value of $\bm{h}_{e}$ of
the remaining boundary edge $e\subset\partial f$ is determined by

$\bm{h}_{e}=(\mathbb{C}_{f,e})^{-1}(\bm{i}_{f}-\mathbb{C}_{f,e_{1}}\bm{h}_{e_{1}}-\mathbb{C}_{f,e_{2}}\bm{h}_{e_{2}})$
(26)

and face $f$ is removed from the list $L$. In the case when $L$ is non-empty
and there is no available face $f$ satisfying the above property, then STT
does not terminate since it stalls in a infinite loop. As shown in [13], STT
termination depends on the choice of the input spanning tree $T$.

We now show that the STT algorithm boils down to a procedure to construct
complete acyclic matchings of 1-chains on the canonical basis
$\widehat{\mathcal{B}}$. This result completes the picture of the equivalence
between tree-cotree techniques and the problem of finding complete acyclic
matching of $k$-chains. The idea behind this observation has it root on the
fact that tree-cotree techniques describe the same actions of acyclic
matchings although using a different language.

To formally state the next theorem we introduce the following notations. Given
a spanning tree $T$ on $\mathcal{K}$, we say that STT _terminates_ (with input
$T$) if it does not stall in an infinite loop. We say that STT _uses the pair
$(e,f)$_ if STT determines the value of edge $e\subset\partial f$ via Eq. 26
during its main loop execution. Finally, we say that STT _uses edge $e$_ if
STT uses the pair $(e,f)$ for some face $f\in L$.

######  Theorem 7.

Let $\mathcal{K}$ be a simplicial complex as in Theorem 6. There exists a
spanning tree on $\mathcal{K}$ for which STT terminates if and only if there
exists a complete acyclic matching of 1-chains on the canonical basis
$\widehat{\mathcal{B}}$.

###### Proof.

Let $T$ be a spanning tree such that STT terminates. We now define an acyclic
matching $\mathcal{M}_{T}$ of 1-chains on $\widehat{\mathcal{B}}$ such that
every _cotree edge_ of $T$, namely an edge $e\in\mathcal{K}_{1}$ which is not
an edge of $T$, is matched in $\mathcal{M}_{T}$. $\mathcal{M}_{T}$ is
constructed during STT execution as follows. We initialize
$\mathcal{M}_{1}=\emptyset$. Next, if STT uses a pair $(e,f)$ during its
execution then we add $(e,f)$ to $\mathcal{M}_{1}$. In this way we get a total
order of pairs in $\mathcal{M}_{1}$ where a pair
$(e_{i},f_{i})\in\mathcal{M}_{1}$ comes before than a pair
$(e_{j},f_{j})\in\mathcal{M}_{1}$ in this total order with $i<j$ if STT uses
$(e_{i},f_{i})$ before than $(e_{j},f_{j})$. The set $\mathcal{M}_{T}$ of all
such pairs is a complete acyclic matching. First, we see that
$\mathcal{M}_{T}$ is a matching since if STT uses the pair $(e,f)$, then the
value of $\bm{h}_{e}$ set by STT can never be reassigned. Second, the matching
$\mathcal{M}_{T}$ is acyclic since if STT uses the pair $(e_{i},f_{i})$, then
$e_{i}$ is not contained in the boundary of any matched face $f_{j}$ in
$\mathcal{M}_{1}$ for $1\leq j<i$. Third, since STT terminates, every cotree
edge $e$ of $T$ is matched in $\mathcal{M}_{1}$. Using Euler’s formula, the
fact that $\mathcal{K}$ is topologically trivial and that $T$ is a spanning
tree on $\mathcal{K}$ we get
$\mathfrak{e}-(\mathfrak{v}-1)=\mathfrak{f}-\mathfrak{c}$. Thus
$\mathcal{M}_{1}$ is complete.

Conversely, let $\mathcal{M}_{1}$ be a complete acyclic matching of 1-chains
on $\widehat{\mathcal{B}}$. We need the following construction, which mimics
the one described in Theorem 6. Consider the subgraph
$G_{\mathcal{M}_{1}}=(V,E_{\mathcal{M}_{1}})$ of $\mathcal{K}$ with set of
vertices given by vertices of $\mathcal{K}$ and $E_{\mathcal{M}_{1}}$ includes
an edge $e\in\mathcal{K}_{1}$ if $e$ is critical with respect to
$\mathcal{M}_{1}$. We now show that $G_{\mathcal{M}_{1}}$ has $\mathfrak{v}-1$
edges and is connected. First, since $\mathcal{M}_{1}$ is a complete acyclic
matching the number of critical edges, which is equal to the number of edges
in $E_{\mathcal{M}_{1}}$, is
$\mathfrak{e}-(\mathfrak{f}-\mathfrak{c})=\mathfrak{v}-1$ where we have used
Euler’s formula $\mathfrak{v}-\mathfrak{e}+\mathfrak{f}-\mathfrak{c}=1$ and
the fact that $\mathcal{K}$ is topologically trivial. Second, assume that, for
the sake of contradiction, $G_{\mathcal{M}_{1}}$ is disconnected. Let $V^{*}$
be the set of nodes in a fixed connected component of $G_{\mathcal{M}_{1}}$.
Let $E_{V^{*}}\subset\mathcal{K}_{1}$ be the subset of _cut edges_ , namely
edges of $G$ with one vertex in $V^{*}$ and one vertex in its complement
$V\setminus V^{*}$. Note that since $\mathcal{K}$ is connected, $E_{V^{*}}$ is
not empty. Moreover, by definition of $G_{\mathcal{M}_{1}}$, to each $e\in
E_{V^{*}}$ corresponds a matched edge $e$ of $\mathcal{K}$ with respect to
$\mathcal{M}_{1}$ forming the pair $(e,f)\in\mathcal{M}_{1}$ for a unique face
$f$ of $\mathcal{K}$.

Now we construct a cycle of the form $f_{1}\succ e_{1}\prec
f_{2}\succ\dots\prec f_{h}\succ e_{h}\prec f_{1}$ with $h\geq 2$,
$(e_{i},f_{i})\in\mathcal{M}_{1}$ for all $i\in\\{1,\dots,h\\}$ and all
$e_{i}\in E_{V^{*}}$ being distinct. Start with a cut edge $e_{1}$ and
consider the unique face $f_{1}$ of $\mathcal{K}$ forming the matched pair
$(e_{1},f_{1})\in\mathcal{M}_{1}$. Then, $f_{1}$ contains at least one other
cut edge $e_{2}\in E_{V^{*}}$, otherwise $e_{1}$ cannot be a cut edge. By
iteratively repeating this process, we obtain a sequence $S$ of pairs
$(e_{1},f_{1}),\dots,(e_{i},f_{i})$ such that $f_{1}\succ e_{1}\prec
f_{2}\succ\dots\prec f_{i-1}\succ e_{i}\prec f_{i}$. Since we have finite
graphs, eventually we will run out of those cut edges which do not appear as
the first component of a pair in $S$. Thus, at some step $i$ of the above
process, we must get a cut edge that is the first component of a pair in $S$
and we get a cycle of the form Eq. 16, which is a contradiction since
$\mathcal{M}_{1}$ is acyclic. We have proved that $G_{\mathcal{M}_{1}}$ is a
connected subgraph of $\mathcal{K}$ with $\mathfrak{v}-1$ edges and thus it is
a spanning tree on $\mathcal{K}$.

Now we prove that STT with input $G_{\mathcal{M}_{1}}$ terminates. Note that,
by definition of $G_{\mathcal{M}_{1}}$, if $(e,f)\in\mathcal{M}_{1}$ then $e$
is a _cotree edge_ of $G_{\mathcal{M}_{1}}$, namely $e$ is not an edge of
$G_{\mathcal{M}_{1}}$. Denote by $C$ the set of all cotree edges of
$G_{\mathcal{M}_{1}}$. Moreover, since $\mathcal{M}_{1}$ is complete,
reasoning as above using Euler’s formula, we get that every cotree edge of
$G_{\mathcal{M}_{1}}$ is matched with respect to $\mathcal{M}_{1}$. Since
$\mathcal{M}_{1}=\\{(e_{1},f_{1}),\dots,(e_{n},f_{n})\\}$ is acyclic, we can
order pairs in $\mathcal{M}_{1}$ in such a way that, for every
$i\in\\{1,\dots,n-1\\}$, $f_{i}$ is not incident to any $e_{i+1},\dots,e_{n}$.
Let $E_{\text{STT}}$ be the set of cotree edges used by STT during its
execution. Note that $E_{\text{STT}}$ is not empty since STT can use at least
the pair $(e_{1},f_{1})$, thanks to the total order chosen on
$\mathcal{M}_{1}$. Suppose that STT does not terminate. This means that,
during its execution, it does not use any cotree edge and $C\setminus
E_{\text{STT}}\neq\emptyset$. We will show that this is impossible. Let $i$ be
the minimum integer such that the cotree edge $e_{i}$ belongs to $C\setminus
E_{\text{STT}}$. Since pairs in $\mathcal{M}_{1}$ are ordered as described
above, $i>1$ and $f_{i}$ can be only incident to cotree edges $e_{j}$ with
$j\leq i$. Moreover, by definition of $i$, there is no $j$ with $j<i$ such
that the cotree edge $e_{j}$ belongs to $C\setminus E_{\text{STT}}$. Hence,
STT should have at least used the pair $(e_{i},f_{i})$, i.e. $e_{i}$ belongs
to $E_{\text{STT}}$ as well. This gives the desired contradiction and
completes the proof. ∎

The proof of Theorem 7 shows that, if the STT algorithm terminates for a given
spanning tree $T$, starting from $T$ we can construct a complete acylclic
matching $\mathcal{M}_{T}$ of 1-chains on $\widehat{\mathcal{B}}$. But also
the other way around, namely, if we have a complete acyclic matching
$\mathcal{M}_{1}$ of 1-chains on $\widehat{\mathcal{B}}$, starting from
$\mathcal{M}_{1}$ we can construct a spanning tree $G_{\mathcal{M}_{1}}$ for
which STT terminates.

By combining Theorem 6 and Theorem 7 we now state the following result which
gives a topological characterization of termination problems of tree-cotree
techniques.

######  Theorem 8.

Let $\mathcal{K}$ be a simplicial complex as in Theorem 6. Then, there exists
a spanning tree on $\mathcal{K}$ for which STT terminates if and only if
$\mathcal{K}$ is collapsible.

There are known examples of triangulations of 3-balls which are not
collapsible [17]. Hence, using Theorem 8, there are triangulations
$\mathcal{K}$ of 3-balls such that, for every possible spanning tree of
$\mathcal{K}$ given as input, STT does not terminate.

We remark that it is NP-complete to decide whether a given 3-dimensional
simplicial complex (embedded or not) is collapsible [27]. However, to the
authors knowledge, the related question for the case of 3-dimensional
simplicial complexes embedded in $\mathbb{R}^{3}$ and with Lipschitz boundary
is still an open problem. Yet, numerical evidence shows that this problem is
very difficult in general although good heuristics exists [28].

These results are the cause of the well-known termination problems of tree-
cotree techniques [13]; see also the discussion in Section 6.

###  4.3 The case where we cannot find a complete acyclic matching

As shown by Theorem 6 and the discussion at the end of Section 4.2, for cell
complex $\mathcal{K}$ given as input, we cannot find in general a complete
acyclic matching $\mathcal{M}_{1}$ of 1-chains on the canonical basis
$\widehat{\mathcal{B}}$. Consequently, Theorem 4 cannot be applied in general
for $\mathcal{B}=\widehat{\mathcal{B}}$.

Let us now assume that acyclic matching $\mathcal{M}_{1}$ is not complete. We
shall now show that, after solving another linear system, we can still get a
discrete vector potential solution $\bm{h}$ solution of Eq. 7 by exploiting
back substitution.

As shown in Section 3.3, an acyclic matching acts on a basis $\mathcal{B}$ by
performing elementary operations on it. Let us consider the new basis
$\mathcal{B}^{\prime}=\mathcal{B}\cdot\mathcal{M}_{1}=\mathcal{U}\cup\mathcal{D}\cup\mathcal{C}$.
There is a corresponding block partition of linear system Eq. 7 as

$\begin{pmatrix}{{\mathbb{C}}|_{\mathcal{U}_{2}\times\mathcal{D}_{1}}}&{{\mathbb{C}}|_{\mathcal{U}_{2}\times\mathcal{C}_{1}}}\\\
{{\mathbb{C}}|_{\mathcal{C}_{2}\times\mathcal{D}_{1}}}&{{\mathbb{C}}|_{\mathcal{C}_{2}\times\mathcal{C}_{1}}}\end{pmatrix}\begin{pmatrix}{{\bm{h}}|_{\mathcal{D}_{1}}}\\\
{{\bm{h}}|_{\mathcal{C}_{1}}}\end{pmatrix}=\begin{pmatrix}{{\bm{i}}|_{\mathcal{U}_{2}}}\\\
{{\bm{i}}|_{\mathcal{C}_{2}}}\end{pmatrix}.$ (27)

The crucial fact turns out to be that
${{\mathbb{C}}|_{\mathcal{C}_{2}\times\mathcal{D}_{1}}}$ is a zero matrix.
Thus, we can determine a discrete vector potential
$\bm{h}=({{\bm{h}}|_{\mathcal{D}_{1}}},{{\bm{h}}|_{\mathcal{C}_{1}}})^{T}$
solution of Eq. 7 by solving, in order,

${{\mathbb{C}}|_{\mathcal{C}_{2}\times\mathcal{C}_{1}}}{{\bm{h}}|_{\mathcal{C}_{1}}}={{\bm{i}}|_{\mathcal{C}_{2}}},$
(28)
${{\mathbb{C}}|_{\mathcal{U}_{2}\times\mathcal{D}_{1}}}{{\bm{h}}|_{\mathcal{D}_{1}}}={{\bm{i}}|_{\mathcal{U}_{2}}}-{{\mathbb{C}}|_{\mathcal{U}_{2}\times\mathcal{C}_{1}}}{{\bm{h}}|_{\mathcal{C}_{1}}},$
(29)

where Eq. 29 is solved by exploiting back substitution as in Section 4.1.

Let us prove that ${{\mathbb{C}}|_{\mathcal{C}_{2}\times\mathcal{D}_{1}}}$ is
a zero matrix. Let
$\mathcal{M}_{1}=\\{(\sigma_{1},\tau_{1})\dots,(\sigma_{n},\tau_{n})\\}$. Let
us consider bases $\\{\dots,\xi^{i}=\phi_{k}(\xi_{i}),\dots\\}$ and
$\\{,\dots,\pi_{i},\dots,\\}$ of $C^{1}(\mathcal{K})$ and
$C_{2}(\mathcal{K})$, respectively. We see that the $(i,j)$-entry of
$\mathbb{C}$ is
$\mathbb{C}_{i,j}=\langle{\delta^{1}(\xi^{j})}\,,{\pi_{i}}\rangle=\langle{\xi^{j}}\,,{\partial_{2}\pi_{i}}\rangle=\langle{\xi_{j}}\,,{\partial_{2}\pi_{i}}\rangle$,
where we have used Eq. 5 and the isomorphism $\phi_{k}:C_{k}(\mathcal{K})\to
C^{k}(\mathcal{K})$ in Eq. 4. Thus, to prove that
${{\mathbb{C}}|_{\mathcal{C}_{2}\times\mathcal{D}_{1}}}$ is a zero matrix we
have to show that every $\sigma_{i}\in\mathcal{D}_{1}$ with
$i\in\\{1,\dots,n\\}$ is not incident to any $\rho\in\mathcal{C}_{2}$.

We proceed by induction on $n$.

Suppose that $n=1$. Necessarily, $\tau_{1}\in\mathcal{U}_{2}$ is the only
basis element in
$\mathcal{B}\cdot\mathcal{M}_{1}=\mathcal{U}\cup\mathcal{D}\cup\mathcal{C}$
incident on $\sigma_{1}$. In fact, if $\tau\in\mathcal{B}$ is incident on
$\sigma_{1}$, then its image under the transformation Eq. 22 is not incident
on $\sigma_{1}$. It follows that $\sigma_{1}$ is not incident to any
$\rho\in\mathcal{C}_{2}$.

We now assume the statement true for $n-1$ and we prove it for $n$. Let us
consider the basis
$\mathcal{B}\cdot\mathcal{M}_{1}=\mathcal{U}\cup\mathcal{D}\cup\mathcal{C}$.
Proceeding as above, we see that $\sigma_{n}$ is not incident to any
$\rho\in\mathcal{C}_{2}$. To conclude, it is sufficient to show that also each
$\sigma_{i}$ with $i\in\\{1,\dots,n-1\\}$ is not incident to any
$\rho\in\mathcal{C}_{2}$. Basis elements in $\mathcal{B}\cdot\mathcal{M}_{1}$
are obtained from that of
$\mathcal{B}\cdot\\{(\sigma_{1},\tau_{1}),\dots,(\sigma_{n-1},\tau_{n-1})\\}$
by applying transformation Eq. 22 with $(\sigma,\tau)=(\sigma_{n},\tau_{n})$.
Using the induction hypothesis we see that the first term in Eq. 22 (i.e.,
$\tau^{\prime}\in C_{2}(\mathcal{K})$) is not incident on any $\sigma_{i}$
with $i\in\\{1,\dots,n-1\\}$. Moreover, each $\sigma_{i}$ with
$i\in\\{1,\dots,n-1\\}$ is not incident to $\tau_{n}$ as pairs in
$\mathcal{M}_{1}$ are ordered as in Theorem 1. This completes the proof.

We now state the following theorem which combines the results of this section
with those of Section 4.1. It can be thought as a specific algebraic version
of Forman’s discrete Morse complex construction [15]. We think that our
presentation and terminology shed light on the linear algebra behind the more
abstract discrete Morse theory constructions.

######  Theorem 9.

Denote by $\mathbb{D}_{k}$ one among the matrices $\mathbb{G}$, $\mathbb{C}$
or $\mathbb{D}$ for $k$ equal to $0$,$1$ or $2$, respectively. Let
$\mathcal{M}_{k}$ be an acyclic matching of $k$-chains on $\mathcal{B}$. Then,
there is a corresponding block partition of $\mathbb{D}_{k}$ as

$\mathbb{D}_{k}=\begin{pmatrix}{{\mathbb{D}_{k}}|_{\mathcal{U}_{k+1}\times\mathcal{D}_{k}}}&{{\mathbb{D}_{k}}|_{\mathcal{U}_{k+1}\times\mathcal{C}_{k}}}\\\
0&{{\mathbb{D}_{k}}|_{\mathcal{C}_{k+1}\times\mathcal{C}_{k}}}\end{pmatrix},$
(30)

where ${{\mathbb{D}_{k}}|_{\mathcal{U}_{k+1}\times\mathcal{D}_{k}}}$ is upper
triangular and invertible.

##  5 Algorithm description

The goal of this section is a recursive algorithm that reduces matrix
$\mathbb{C}$ into a row echelon form by means of elementary operations on the
basis $\mathcal{B}=\bigcup_{k}\mathcal{B}_{k}$ and we present in Section 5.2.
We first describe in Section 5.1 our novel greedy procedure to construct
acyclic matchings.

###  5.1 Greedy approach to construct acyclic matchings

To minimize computational effort of change of basis in Eq. 28, we have to
carefully choose how to construct the acyclic matching $\mathcal{M}_{1}$. It
is visible from Section 3.3 that constructing acyclic matchings by internal
collapses can get more complicated than by elemetary collapses since changes
of basis of type 2 and 3 are involved. In this case, we have to express the
new basis elements as a linear combination of the previous ones and hence they
cannot be removed from the data structure after each collapse.

The above discussion motivates the concept of _degree_ of a basis element
$\sigma\in\mathcal{B}$. We define $\mathrm{deg}(\sigma)$ as the cardinality of
the coboundary of $\sigma$, i.e. the set
$\mathrm{cobd}_{\mathcal{B}}{(\sigma)}$ defined in Eq. 15. Note that if
$\mathrm{deg}(\sigma)$ is 1 then $\sigma$ is free. In our approach, also the
case $\mathrm{deg}(\sigma)=2$ will play a fundamental role. We define
$\sigma\in\mathcal{B}$ to be _flat_ if $\mathrm{deg}(\sigma)=2$.

To reduce the amount of computation during Gaussian elimination, it is wise to
first search for basis elements with lowest degree. Indeed, if
$\mathrm{deg}(\sigma)$ is 1, then $\sigma$ is free and the new basis is a
selection of the previous one and no algebraic operations are needed. In this
case, it should be possible to avoid all the matrix algebra computation and
efficiently organize the basis so that matrices are in triangular form as in
Lemma 1.

Keeping in mind the above heuristics, we search for basis elements having
smallest degree. Specifically, we do not strictly choose basis elements with
minimum degree but instead we proceed in a sequential manner. First, we search
for all free basis elements until exhaustion. Next, we search for all flat
basis elements until exhaustion. We have pursued this method because been
motivated by its practical implementation and performance on test problems,
rather than by following the best theoretical greedy approach; see the
discussion in Section 6.

We construct acyclic matchings using a standard elementary collapse greedy
procedure, where we search for collapsing sequences of free basis elements in
a monotone-like fashion [28]. That is, we proceed in sequential order with
respect to the dimension of the basis elements by first collapsing 2-chains
and then 1-chains.

Algorithm 2 Random strategy to construct acyclic matchings

1:procedure ConstructAcyclicMatchings($\mathcal{B}$)

2: for $k$ ranging from $2$ down to $1$ do

3: while there exists a free pair $(\sigma,\tau)$ in $\mathcal{B}_{k}$ do

4: elementary collapse $(\sigma,\tau)$

5: insert $(\sigma,\tau)$ to $\mathcal{M}_{k}$

6: end while

7: while there exists a flat pair $(\sigma,\tau)$ in $\mathcal{B}_{k}$ do

8: internal collapse of $(\sigma,\tau)$

9: insert $(\sigma,\tau)$ to $\mathcal{M}_{k}$

10: end while

11: end for

12: return $\mathcal{M}_{2}\cup\mathcal{M}_{1}$

13:end procedure

It is clear that Algorithm 2 always terminates. Moreover, the obtained
matching is acyclic since at each iteration collapses are performed. Note that
Algorithm 2 requires no backtracking since new free or flat pairs can only
appear after each new collapse. Thus, the worst-case complexity is linear
using a suitable algorithm implementation that employs a list data structure.

As proved in Section 4.2, we can always find a complete acyclic matching of
2-chains by using a standard spanning tree construction. The proof of this
result implies that the order in which we collapse 2-chains in Algorithm 2 is
not important. Therefore, there always exists a complete acyclic matching of
2-chains and we can get one by performing collapses in a random fashion.

###  5.2 Recursive algorithm

We now present a recursive construction of acyclic matchings.

Normally one aims at finding a complete acyclic matching, namely a matching
that reaches the needed number of matched pairs so that we can apply Theorem
4. However, as proved in Theorem 6, determining a complete acyclic matching
$\mathcal{M}_{1}$ of 1-chains is a hard algorithmic problem and we do not
tackle it. In Algorithm 2 we employed a greedy strategy that randomly selects
an acyclic matching but it can easily happen that the obtained acyclic
matching $\mathcal{M}_{1}$ is not complete.

Let us consider the case where $\mathcal{M}_{1}$ is not complete. After
applying basis transformations associated with $\mathcal{M}_{1}$ to the
current basis $\mathcal{B}$, we get a new basis
$\mathcal{B}^{\prime}=\mathcal{B}\cdot\mathcal{M}_{1}=\mathcal{U}\cup\mathcal{D}\cup\mathcal{C}$.
In particular, the set of critical basis elements $\mathcal{C}$ in
$\mathcal{B}^{\prime}$ forms a linearly independent subset of
$\mathcal{B}^{\prime}$. As pointed out in Section 4.3, the focus now shifts to
the set of critical basis elements $\mathcal{C}$.

The novel idea is that we may call the whole routine recursively, where the
output set of critical basis elements obtained from the previous iteration
becomes the input basis for the next iteration. Fundamentally, we think linear
system Eq. 28 as an instance of the original linear system Eq. 7, although
considering a subset of the previous basis. This operation is well-defined
since we have the partition
$\mathcal{B}^{\prime}=\mathcal{U}\cup\mathcal{D}\cup\mathcal{C}$ and the set
$\mathcal{C}$ of critical basis elements form a linearly independent subset of
$\mathcal{B}$.

If at a certain recursion stage we get a complete acyclic matching
$\mathcal{M}_{1}$ of 1-chains, then we can find a discrete vector potential
$\bm{h}$ solution of Eq. 7 by recursively applying the reasoning described in
Section 4.3. Otherwise, we recursively apply the same routine on the obtained
set of critical basis elements $\mathcal{C}$.

Algorithm 3 Construction of a discrete vector potential $\bm{h}$ solution of
Eq. 7

1:procedure ConstructDiscreteVectorPotential($\mathcal{B}$) $\triangleright$
at the beginning $\mathcal{B}$ is the canonical basis $\widehat{\mathcal{B}}$

2: construct an acyclic matching $\mathcal{M}_{2}\cup\mathcal{M}_{1}$ on
$\mathcal{B}$ using Algorithm 2

3: if $\mathcal{M}_{1}$ is complete then

4: set ${{\bm{h}}|_{\mathcal{C}_{1}}}=\bm{0}$

5: determine ${{\bm{h}}|_{\mathcal{D}_{1}}}$ by back substitution using
Algorithm 1

6: return $({{\bm{h}}|_{\mathcal{D}_{1}}},{{\bm{h}}|_{\mathcal{C}_{1}}})^{T}$

7: else

8: compute new basis $\mathcal{B}^{\prime}$ from $\mathcal{B}$ by change of
basis associated with $\mathcal{M}_{2}\cup\mathcal{M}_{1}$

9: if $\mathcal{M}_{1}$ is empty then

10: determine ${{\bm{h}}|_{\mathcal{C}_{1}}}$ in Eq. 28 using a sparse linear
system solver

11: return ${{\bm{h}}|_{\mathcal{C}_{1}}}$

12: else

13: recursively call this routine with input the set of critical elements
$\mathcal{C}$ of
$\mathcal{B}^{\prime}=\mathcal{U}\cup\mathcal{D}\cup\mathcal{C}$

14: set ${{\bm{h}}|_{\mathcal{C}_{1}}}$ to the value returned by the recursive
process

15: determine ${{\bm{h}}|_{\mathcal{D}_{1}}}$ by back substitution using
Algorithm 1

16: return $({{\bm{h}}|_{\mathcal{D}_{1}}},{{\bm{h}}|_{\mathcal{C}_{1}}})^{T}$

17: end if

18: end if

19:end procedure

For a given input basis $\mathcal{B}$, Algorithm 3 produces a sequence of
bases
$\mathcal{B}=\mathcal{B}^{\,(0)},\mathcal{B}^{\,(1)},\dots,\mathcal{B}^{\,(n)}$,
where index $i$ keeps track of the recursion depth. Each
$\mathcal{B}^{\,(i+1)}$ is constructed from $\mathcal{B}^{\,(i)}$ by recursion
as follows. We have the partition of
$\mathcal{B}^{\prime\,(i)}=\mathcal{B}^{\,(i)}\cdot(\mathcal{M}_{2}^{\,(i)}\cup\mathcal{M}_{1}^{\,(i)})=\mathcal{U}^{\,(i)}\cup\mathcal{D}^{\,(i)}\cup\mathcal{C}^{\,(i)}$,
associated with acyclic matching
$\mathcal{M}_{2}^{\,(i)}\cup\mathcal{M}_{1}^{\,(i)}$ for
$i\in\\{0,\dots,n\\}$. Note that from the discussion at the end of Section 5.1
it follows that $\mathcal{M}_{2}^{\,(i)}$ is not empty only for $i=0$. We set
$\mathcal{B}^{\,(i+1)}=\mathcal{C}^{\,(i)}$.

It is clear that the cardinality of each basis $\mathcal{B}^{\,(i)}$ decreases
as long as there are some collapses to be made by Algorithm 2. If no free or
flat basis elements are available in Algorithm 2, then the recursion is
stopped and a sparse linear system solver is used. So the algorithm always
terminates.

At the end of the recursion, we can determine a discrete vector potential
$\bm{h}$ solution of Eq. 7 by recursively applying the reasoning described in
Section 4.3. We write the general form a discrete vector potential $\bm{h}$
solution of Eq. 7 as

$\bm{h}=({{\bm{h}}|_{\mathcal{D}_{1}^{\,(0)}}},\dots,{{\bm{h}}|_{\mathcal{D}_{1}^{\,(n)}}},{{\bm{h}}|_{\mathcal{C}_{1}^{\,(n)}}})^{T},$
(31)

where each ${{\bm{h}}|_{\mathcal{D}_{1}^{\,(i)}}}$ is the subvector of
$\bm{h}$ induced by $\mathcal{D}_{1}^{\,(i)}$, see Fig. 3. We start by
possibly determining ${{\bm{h}}|_{\mathcal{C}_{1}^{\,(n)}}}$ using a linear
system solver. Next, we determine each
${{\bm{h}}|_{\mathcal{D}_{1}^{\,(0)}}},\dots,{{\bm{h}}|_{\mathcal{D}_{1}^{\,(n)}}}$
by back substitution starting from ${{\bm{h}}|_{\mathcal{D}_{1}^{\,(n)}}}$
down to ${{\bm{h}}|_{\mathcal{D}_{1}^{\,(1)}}}$. This proves the correctness
of Algorithm 3.

######  Theorem 10.

Given as input the canonical basis $\widehat{\mathcal{B}}$, Algorithm 3
returns a discrete vector potential $\bm{h}$ solution of Eq. 7.

Figure 3: General block structure of matrix $\mathbb{C}$ produced by Algorithm
3.

Concerning the computational complexity of Algorithm 3, if at a certain
recursion stage no free or flat basis elements are available in Algorithm 2, a
linear system solver is employed. Thus, the worst-case complexity is cubical
with respect to the size of the input mesh. Yet, the average complexity of
Algorithm 3 in all tested problems has been linear. This is because, in
practice, only one recursively call of Algorithm 3 is needed. Indeed, the
first recursive call of Algorithm 3 always finds a complete acyclic matching
$\mathcal{M}_{1}$ with respect to new basis $\mathcal{B}^{\,(1)}$. In other
words, there is no need in practice to solve a linear system of the form Eq.
28 with matrix
${{\mathbb{C}}|_{\mathcal{C}_{2}^{\,(1)}\times\mathcal{C}_{1}^{\,(1)}}}$ but
instead we can recursively apply back substitution to determine
${{\bm{h}}|_{\mathcal{D}_{1}^{\,(1)}}},{{\bm{h}}|_{\mathcal{D}_{1}^{\,(0)}}}$,
in this order, after setting free variables in $\mathcal{C}_{1}^{\,(1)}$ to
zero.

##  6 Numerical results

In this section we illustrate the performance of Algorithm 3. We consider
different sets of test problems. In the first set, we focus on simple
triangulations that appear in practical boundary value problems. Next, we
present more complicated benchmark triangulations.

For all triangulations that appear in practical boundary value problems, our
greedy procedure in Algorithm 2 always finds a complete acyclic matching
$\mathcal{M}_{1}^{\,(0)}$. Only for the more complicated benchmark problems it
is necessary to exploit the novel recursive procedure in Algorithm 3. However,
to compute a discrete vector potential it is enough, in all tested problems,
only one recursively call of Algorithm 3.

The algorithm has been implemented in C++. All the numerical computations have
been performed in a Intel Core i7-3720QM, with a processor at 2.60 GHz in a
laptop with 16 GB of RAM.

###  6.1 Triangulations coming from real case boundary value problems

We consider a triangulation coming from a computational electromagnetics
application. As an example, the computation of the source magnetic field for
the TEAM problem 7 has been addressed [29].

Table 1 contains information on the number of cells of triangulations of
different sizes together with the time (in milliseconds) required to compute
the discrete vector potential $\bm{h}$ using Algorithm 3. It is worth noticing
that in a triangulation with about 2 million tetrahedra our procedure computes
a discrete vector potential under a second. We run our Algorithm 2 with
different triangulations of the same metal plate and on each example Algorithm
2 finds a complete acyclic matching $\mathcal{M}_{1}^{\,(0)}$. Thus, there is
no recursive call of Algorithm 3. For the considered examples in Table 1 we
can clearly see the linear behaviour of the computational time with respect to
the size of the triangulations.

Table 1: Running times of Algorithm 3 for the modified TEAM benchmark example for triangulations of decreasing size. Name | Tetrahedra $(\mathfrak{c})$ | Faces $(\mathfrak{f})$ | Edges $(\mathfrak{e})$ | Vertices $(\mathfrak{v})$ | Time [ms]
---|---|---|---|---|---
Mesh 1 | 1,851,493 | 3,871,379 | 2,419,350 | 399,465 | 992
Mesh 2 | 1411688 | 2847256 | 1683787 | 248220 | 756
Mesh 3 | 529,664 | 1,065,104 | 626,566 | 91,127 | 284
Mesh 4 | 186264 | 378588 | 226584 | 34261 | 101

###  6.2 Bing’s House

A Bing’s House is now considered [30]. The simplicial complex, homeomorphic to
a 3-dimensional ball, can be obtained by replacing every surface in the Bing’s
House by a thick wall made of 3-cells. At the end of this procedure we obtain
the polyhedron in Fig. 4. Although we can informally identify two “chambers”,
it can be demonstrated that the Bing’s House is homeomorphic to the three-
dimensional ball.

As in the previous set of tests, Table 2 contains information about the number
of cells of the considered triangulations together with the computational time
required to compute a discrete vector potential. We have found that in almost
all runs of Algorithm 2 we get a complete acyclic matching
$\mathcal{M}_{1}^{\,(0)}$. Only in a few cases we need to resort to a
recursive call of Algorithm 3. However, to compute a discrete vector potential
it is enough, in all these cases, only one recursive call of Algorithm 3,
given that in the first recursion we always find a complete acyclic matching
$\mathcal{M}_{1}^{\,(1)}$.

To measure the complexity of the first recursive call, we consider the
cardinality of the basis $\mathcal{B}_{2}^{\,(1)}$, namely the output basis of
2-chains becoming the input for the first recursive call of Algorithm 3. We
have found that the cardinality of $\mathcal{B}_{2}^{\,(1)}$ is always less
than 10 on thousands of algorithm runs with different choices of the acyclic
matching $\mathcal{M}_{2}^{\,(0)}$ of 2-chains. Accordingly, as reported in
Table 2, we observe no influence of the recursive call in Algorithm 3 on the
linear behaviour of the running times with respect to the size of the
triangulations.

Table 2: Running times of Algorithm 3 for various triangulations of the thick Bing’s House. Name | Tetrahedra $(\mathfrak{c})$ | Faces $(\mathfrak{f})$ | Edges $(\mathfrak{e})$ | Vertices $(\mathfrak{v})$ | Time [ms]
---|---|---|---|---|---
Bing 1 | 800,020 | 1,600,537 | 937,631 | 137,115 | 429
Bing 2 | 87,221 | 175,317 | 102,212 | 14,117 | 47

###  6.3 Knot-theoretic obstructions

We consider 3-balls of $\mathbb{R}^{3}$ which admit non-collapsible
triangulations.

Obstructions coming from short knots have been considered first in the works
[31, 32].

Figure 4: Two views of the considered 3-dimensional thickening of a Bing’s
house with two rooms.

In [31], Bing proved, using knot theory, that some triangulations of the
3-ball are not collapsible. Bing’s construction works as follows. One starts
with a triangulated 3-ball $\Omega$ and introduces a “knotted spanning arc” in
its 1-skeleton. A knotted spanning arc is an arc as in Fig. 5. We dig a knot-
shaped tubular hole inside $\Omega$ starting from the top and we stop digging
one step before the tunnel go through the bottom of $\Omega$. In this way we
obtain a 3-ball $\Omega^{\prime}$ containing a knot having all its edges on
the boundary of $\Omega^{\prime}$, except for a single interior edge.

Figure 5: A knotted spanning arc in a 3-ball $\Omega$ at the core of Furch’s
construction.

If the knot is sufficiently complicated (like a double, or a triple trefoil),
Bing’s ball cannot be collapsible [31, 17]. In contrast, if the knot is simple
enough (like a single trefoil), then Bing’s ball may be collapsible. The
construction also appears in the 1924 work [33] of Furch and for the present
discussion we refer to it as _Furch’s knotted ball_ ; see [34] (Section 3.1)
for an historical account.

Firstly, we consider the simplest case of Furch’s knotted ball with only one
trefoil knot. Table 3 summarizes the geometrical information of this
triangulation named Furch 1. In each run of Algorithm 3, we do not find a
complete acyclic matching $\mathcal{M}_{1}^{\,(0)}$ of 1-chains. Thus,
Algorithm 3 is recursively called. However, only one recursive call is needed
in all the considered runs of Algorithm 3, given that in the first recursion
we always find a complete acyclic matching $\mathcal{M}_{1}^{\,(1)}$.

To measure the complexity of the first recursive call, we consider, as done in
Section 6.2, the cardinality of the basis $\mathcal{B}_{2}^{\,(1)}$. We have
found that the cardinality of $\mathcal{B}_{2}^{\,(1)}$ is always less than 30
on thousands of algorithm runs with different choices of the acyclic matching
$\mathcal{M}_{2}^{\,(0)}$ of 2-chains.

As the last benchmark problem, we consider a more complicated obstruction. We
dig one hundred trefoil knots in a parallel-like fashion starting from the top
of $\Omega$. Table 3 reports again the results for this triangulation named
Furch 100. We found that the cardinality of $\mathcal{B}_{2}^{\,(1)}$ is
always less than 360 on thousands of algorithm runs with different choices of
the acyclic matching $\mathcal{M}_{2}^{\,(0)}$ of 2-chains. Also in this case,
only one recursive call is needed in all considered runs of Algorithm 3, given
that in the first recursion we always find a complete acyclic matching
$\mathcal{M}_{1}^{\,(1)}$.

We observe no influence of the recursive call in Algorithm 3 on the linear
behaviour of the running times with respect to the size of the triangulations.
This is because of the small cardinality of $\mathcal{B}_{2}^{\,(1)}$.

This example show the effectiveness and generality of our method. In fact,
using the approach in [12], we have computed more than one hundred double
integral evaluations. Similarly, using the approach in [14], we have
constructed and then solved a linear system of more than one hundred equations
having as unknowns the symbolic variables employed in the approach. Therefore,
Furch 100 is an explicit example of a triangulation on which the approaches
[12, 14] perform poorly compared to Algorithm 3. To have a provably good
method which is reliable in practice, Algorithm 3 is expressly needed.

An important point is that in all tested problems we do not see a dependence
between the cardinality of the basis $\mathcal{B}_{2}^{\,(1)}$ and the number
of cells of the input triangulations. Thus, this quantity can be used as an
indicator of how easy it is to find a complete acyclic matching
$\mathcal{M}_{1}^{\,(0)}$ on a given input triangulation, namely, it
quantifies the “topological complexity” of the triangulation.

Table 3: Running times of Algorithm 3 for the Furch’s knotted balls. Name | Tetrahedra $(\mathfrak{c})$ | Faces $(\mathfrak{f})$ | Edges $(\mathfrak{e})$ | Vertices $(\mathfrak{v})$ | Time [ms]
---|---|---|---|---|---
Furch 100 | 243062 | 506619 | 311547 | 47991 | 129
Furch 1 | 31189 | 63830 | 38445 | 5805 | 17

##  7 Conclusions

The novel algorithm presented in this paper was proved to be general,
straightforward to implement and outperforms competing state-of-art algorithms
in the class of admissible meshes while showing linear average complexity with
respect to the input mesh size. By applying discrete Morse theory, we have
shown that for the important class of simplicial triangulations we achieve
linear computational complexity for all considered test problems. These
include, besides real case triangulations having simple topological
properties, also pathological triangulations. A challenging test case made of
one hundred trefoil knots have been considered and yet the proposed algorithm
succeeds in computing a discrete vector potential. Thus, we expect that our
algorithm works for every practical mesh that one encounters in practical
applications.

Worst-case complexity analysis can be misleading in the analysis of our
Algorithm 3. Worst-case complexity analysis of our algorithm assumes that at
certain iteration no new matched pairs are possible. In this case we need to
employ a sparse linear system solver, which leads to a cubical worst-case
complexity with respect to mesh size parameters. Yet, the average complexity
is linear in all tested problems. Our point of view is that the reason we see
linear computational complexity is because our examples restrict to
3-dimensional cell complexes decomposing bounded domains of $\mathbb{R}^{3}$
with sufficiently regular boundaries. We have observed that our recursive
strategy based on algebraic discrete Morse theory is efficient at solving
potential topological obstructions that can appear in 3-dimensional cases.
However, to prove that linear worst-case times are guaranteed is still an open
problem and will be subject to future work.

## References

  * Lipnikov et al. [2014] K. Lipnikov, G. Manzini, M. Shashkov,  Mimetic finite difference method,  Journal of Computational Physics 257 (2014) 1163–1227. Physics-compatible numerical methods.
  * Pitassi et al. [2021] S. Pitassi, R. Ghiloni, F. Trevisan, R. Specogna,  The role of the dual grid in low-order compatible numerical schemes on general meshes,  Journal of Computational Physics 436 (2021) 110285.
  * Di Pietro and Droniou [2021] D. A. Di Pietro, J. Droniou,  An arbitrary-order method for magnetostatics on polyhedral meshes based on a discrete de rham sequence,  Journal of Computational Physics 429 (2021) 109991.
  * Tran-Cong [1990] T. Tran-Cong,  On the potential of a solenoidal vector field,  Journal of Mathematical Analysis and Applications 151 (1990) 557–580.
  * Cohen et al. [2014] M. B. Cohen, B. T. Fasy, G. L. Miller, A. Nayyeri, R. Peng, N. J. Walkington,  Solving 1-laplacians in nearly linear time: Collapsing and expanding a topological ball,  in: SODA, 2014.
  * Webb and Forghani [1989] J. P. Webb, B. Forghani,  A single scalar potential method for 3d magnetostatics using edge elements,  International Magnetics Conference (1989) JD1–JD1.
  * Le Menach et al. [1998] Y. Le Menach, S. Clenet, F. Piriou,  Determination and utilization of the source field in 3d magnetostatic problems,  IEEE Transactions on Magnetics 34 (1998) 2509–2512.
  * Rodríguez et al. [2013] A. A. Rodríguez, E. Bertolazzi, R. Ghiloni, A. Valli,  Construction of a finite element basis of the first de Rham cohomology group and numerical solution of 3d magnetostatic problems,  SIAM Journal on Numerical Analysis 51 (2013) 2380–2402.
  * Dłotko and Specogna [2013] P. Dłotko, R. Specogna,  Physics inspired algorithms for (co)homology computations of three-dimensional combinatorial manifolds with boundary,  Computer Physics Communications 184 (2013) 2257–2266.
  * Alonso Rodríguez et al. [2015] A. Alonso Rodríguez, E. Bertolazzi, R. Ghiloni, A. Valli,  Finite element simulation of eddy current problems using magnetic scalar potentials,  Journal of Computational Physics 294 (2015) 503–523.
  * Silberman et al. [2019] Z. J. Silberman, T. R. Adams, J. A. Faber, Z. B. Etienne, I. Ruchlin,  Numerical generation of vector potentials from specified magnetic fields,  Journal of Computational Physics 379 (2019) 421–437.
  * Rodríguez and Valli [2015] A. A. Rodríguez, A. Valli,  Finite element potentials,  Applied Numerical Mathematics 95 (2015) 2–14.
  * Dlotko and Specogna [2010] P. Dlotko, R. Specogna,  Critical analysis of the spanning tree techniques,  SIAM J. Numer. Anal. 48 (2010) 1601–1624.
  * Dlotko and Specogna [2011] P. Dlotko, R. Specogna,  Efficient generalized source field computation for h-oriented magnetostatic formulations,  European Physical Journal-applied Physics 53 (2011) 20801.
  * Forman [1998] R. Forman,  Morse theory for cell complexes,  Advances in Mathematics 134 (1998) 90–145.
  * Kozlov [2008] D. Kozlov,  Combinatorial algebraic topology,  in: Algorithms and computation in mathematics, 2008\.
  * Benedetti and Lutz [2013] B. Benedetti, F. Lutz,  Knots in collapsible and non-collapsible balls,  Electron. J. Comb. 20 (2013) P31.
  * Cantarella et al. [2002] J. Cantarella, D. DeTurck, H. Gluck,  Vector calculus and the topology of domains in 3-space,  The American Mathematical Monthly 109 (2002) 409 – 442.
  * Benedetti et al. [2010] R. Benedetti, R. Frigerio, R. Ghiloni,  The topology of Helmholtz domains,  arXiv: Geometric Topology (2010).
  * Christiansen [2008] S. Christiansen,  A construction of spaces of compatible differential forms on cellular complexes,  Mathematical Models and Methods in Applied Sciences 18 (2008) 739–757.
  * Tonti [2013] E. Tonti,  The mathematical structure of classical and relativistic physics: A general classification diagram,  2013\.
  * Lipnikov et al. [2014] K. Lipnikov, G. Manzini, M. Shashkov,  Mimetic finite difference method,  J. Comput. Phys. 257 (2014) 1163–1227.
  * Whitehead [1939] J. Whitehead,  Simplicial spaces, nuclei and m‐groups,  Proceedings of The London Mathematical Society (1939) 243–327.
  * Strang [1993] G. Strang,  Introduction to linear algebra,  1993\.
  * Cormen et al. [2009] T. Cormen, C. Leiserson, R. Rivest, C. Stein,  Introduction to algorithms, third edition,  2009\.
  * Lewiner et al. [2003] T. Lewiner, H. Lopes, G. Tavares,  Optimal discrete morse functions for 2-manifolds,  Comput. Geom. 26 (2003) 221–233.
  * Tancer [2016] M. Tancer,  Recognition of collapsible complexes is np-complete,  Discrete & Computational Geometry 55 (2016) 21–38.
  * Benedetti and Lutz [2014] B. Benedetti, F. H. Lutz,  Random discrete Morse theory and a new library of triangulations,  Experimental Mathematics 23 (2014) 66 – 94.
  * Fujiwara and Nakata [1990] K. Fujiwara, T. Nakata,  Results for benchmark problem 7 (asymmetrical conductor with a hole),  Compel-the International Journal for Computation and Mathematics in Electrical and Electronic Engineering 9 (1990) 137–154.
  * Cohen [1973] M. M. Cohen,  A course in simple-homotopy theory,  1973\.
  * Bing [1964] R. Bing,  Some aspects of the topology of 3-manifolds related to the poincare conjecture (1964).
  * Goodrick [1968] R. Goodrick,  Non-simplicially collapsible triangulations of In,  1968\.
  * Furch [1924] R. Furch,  Zur grundlegung der kombinatorischen topologie,  Abh.Math.Semin.Univ.Hambg. 3 (1924) 69–88.
  * Ziegler [1998] G. M. Ziegler,  Shelling polyhedral 3-balls and 4-polytopes,  Discrete & Computational Geometry 19 (1998) 159–174.