remarkRemark exampleExample hypothesisHypothesis claimClaim problemProblem
Cholesky-like preconditioner for Hodge Laplacian A. Savostianov, F. Tudisco,
N. Guglielmi

# Cholesky-like Preconditioner for Hodge Laplacians via Heavy Collapsible
Subcomplex††thanks: Submitted to the editors on DATE.

Anton Savostianov Gran Sasso Science Institute, L’Aquila, Italy ()
<EMAIL_ADDRESS>Francesco Tudisco  Maxwell Institute & School of
Mathematics, University of Edinburgh, UK<EMAIL_ADDRESS>Nicola
Guglielmi Gran Sasso Science Institute, L’Aquila, Italy ()
<EMAIL_ADDRESS>
###### Abstract

Techniques based on $k$-th order Hodge Laplacian operators $L_{k}$ are widely
used to describe the topology as well as the governing dynamics of high-order
systems modeled as simplicial complexes. In all of them, it is required to
solve a number of least square problems with $L_{k}$ as coefficient matrix,
for example in order to compute some portions of the spectrum or integrate the
dynamical system. In this work, we introduce the notion of optimal collapsible
subcomplex and we present a fast combinatorial algorithm for the computation
of a sparse Cholesky-like preconditioner for $L_{k}$ that exploits the
topological structure of the simplicial complex. The performance of the
preconditioner is tested for conjugate gradient method for least square
problems (CGLS) on a variety of simplicial complexes with different dimensions
and edge densities. We show that, for sparse simplicial complexes, the new
preconditioner reduces significantly the condition number of $L_{k}$ and
performs better than the standard incomplete Cholesky factorization.

###### keywords:

simplicial complex, Hodge Laplacian, graph Laplacian, Cholesky preconditioner,
Gauss elimination, collapsible simplicial complex

65F08, 05C50, 57M15, 62R40

##  1 Introduction

Graph models of network systems are ubiquitous throughout the sciences.
However, a fast-growing line of work in modern network science has highlighted
the importance of higher-order models that take directly into account multi-
agent interactions rather than classical graph pairwise interactions.
Simplicial complexes are a popular higher-order generalization of graphs that
allows modeling non-diadic interactions by means of a sequence of nested nodal
simplicies. Given the nested structure of simplicial complexes, boundary
operators $B_{k}$, mapping each simplex to its border, and the associated
higher-order Laplacian operators $L_{k}$, [18], provide a fundamental
description of the topology of the system. For example, one may show the
convergence of $L_{k}$ in the thermodynamic limit to Hodge Laplacians on
manifolds, [4, 5]; the elements of the kernel of the higher-order Laplacian
operators correspond to $k$-dimensional holes in the complex (connected
components, one-dimensional holes, voids, etc.); the positive part of the
spectrum of $L_{k}$ can be used to measure the topological stability of the
system, [13]; the corresponding eigenvectors provide information about higher-
order topological features [8, 11] and can be used to describe the underlying
dynamics, e.g. in case of synchronization or higher-order diffusion, [9, 28].
Computing (part of) the spectrum of $L_{k}$ with standard iterative methods
heavily requires efficient solving a number of linear systems of the form
$L_{k}\boldsymbol{x}=\boldsymbol{f}$, see e.g. [7]. Moreover, the higher-order
system dynamics, for example in the case of random walks, PageRank, and
opinion spreading, [23, 24], are primarily defined in terms of the same type
of linear systems. While it was observed in [6] that for specifically
structured simplicial complexes one can efficiently solve the linear system
$L_{k}\boldsymbol{x}=\boldsymbol{f}$ directly, e.g. when the complex can be
embedded in a triangulation of a unit sphere in three dimensions, attempting
to extend that approach to arbitrary complex significantly increases the
complexity of the direct solver, [3].

In order to speed up the solution of linear systems with higher-order
Laplacians via iterative methods, in this work we propose a sparse and
computationally efficient Choelsky-like preconditioner for $L_{k}$ based on
the notion of weak collapsible complex. In the case of classical Laplacian
$L_{0}$, one can build a sparse stochastic Cholesky preconditioner based on
the underlying graph structure, [16]. By investigating the structure of the
Schur complements forming the Cholesky factorization method, we observe that
the same construction transfers to the higher-order setting only for
simplicial complexes with a specific weakly collapsible topological structure
[29]. As a result, we propose an algorithm that builds a fast lower-triangular
preconditioner based on the optimal collapsible subcomplex. While we present
the problem and the main ideas in the general setting, we will focus on the
case $k=1$ when designing the algorithm and testing its efficiency.

The rest of the paper is organized as follows: we introduce formal definitions
for simplicial complexes and $k$-Laplacian operators in Section 2; the
transition from the linear system $L_{k}\boldsymbol{x}=\boldsymbol{f}$ to the
least square problems on up-Laplacians $L_{k}^{\uparrow}$ is demonstrated in
Section 3, Theorem 3.1. Section 4 describes the Cholesky preconditioner and
the relation to the weak collapsibility of the simplicial complex, Section 5.
Given that, we introduce the heavy collapsible subcomplex preconditioner
HeCSin Section 6, based on the theoretically optimal subcomplex preconditioner
introduced in Lemmas 6.6 and 6.9. Finally, Section 7 provides numerical
results for HeCS in terms of condition number and performance of the
preconditioned conjugate gradient method, as compared to the unpreconditioned
system and the standard incomplete Cholesky factorization.

##  2  Simplicial complex and Homology groups

A simplicial complex $\mathcal{K}$ on the vertices $\\{v_{1},v_{2}\ldots
v_{n}\\}$ is a collection of simplices $\sigma$, sets of nodes with the
property that all the subsets of $\sigma$ are simplicies of $\mathcal{K}$ too.
We refer to a simplex made out of $k$ nodes $\sigma=[v_{i_{1}},\ldots
v_{i_{k+1}}]$ as being of order $k$, and write $\dim\sigma=k$; the set of all
the simplices of order $k$ in the complex $\mathcal{K}$ is denoted by
$\mathcal{V}_{k}(\mathcal{K})$. Thus, $\mathcal{V}_{0}(\mathcal{K})$ are the
vertices of $\mathcal{K}$, $\mathcal{V}_{1}(\mathcal{K})$ are edges between
pairs of vertices, $\mathcal{V}_{2}(\mathcal{K})$ triangles connecting three
vertices, and so on. We let $m_{k}=|\mathcal{V}_{k}(\mathcal{K})|$ denote the
cardinality of $\mathcal{V}_{k}(\mathcal{K})$.

Each set of simplices
$\mathcal{V}_{k}(\mathcal{K})=\left\\{\sigma_{1},\dots\sigma_{m_{k}}\right\\}$
induces a linear space of formal sums over the simplicies
$C_{k}(\mathcal{K})=\left\\{\sum_{i=1}^{m_{k}}\alpha_{i}\sigma_{i}\mid\alpha_{i}\in\mathds{R}\right\\}$
referred to as _chain space_ ; in particular, $C_{0}(\mathcal{K})$ is known as
the space of vertex states and $C_{1}(\mathcal{K})$ as the space of edge
flows. Simplices of different orders are related through the boundaries
operators $\partial_{k}$ mapping the simplex to its boundary; formally,
$\partial_{k}:C_{k}(\mathcal{K})\mapsto C_{k-1}(\mathcal{K})$ is defined
through the alternating sum:

$\partial_{k}[v_{1},v_{2},\ldots v_{k}]=\sum_{i=1}^{k}(-1)^{i-1}[v_{1},\ldots
v_{i-1},v_{i+1},\ldots v_{k}]$

By fixing an ordering for $\mathcal{V}_{k}(\mathcal{K})$ we can fix a
canonical basis for $C_{k}(\mathcal{K})$ and represent each boundary operator
as a matrix $B_{k}\in\mathrm{Mat}_{m_{k-1}\times m_{k}}$ with exactly $k$
nonzero entries per each column, being either $+1$ or $-1$. For these
matrices, the fundamental property of topology holds: _the boundary of the
boundary is zero_ , [18, Thm. 5.7]:

(1) $B_{k}B_{k+1}=0$

The matrix representation $B_{k}$ of the boundary operator $\partial_{k}$
requires fixing an ordering of the simplices in $\mathcal{V}_{k}(\mathcal{K})$
and $\mathcal{V}_{k-1}(\mathcal{K})$. As it will be particularly relevant for
the purpose of this work, we emphasize that we order triangles and edges as
follows: triangles in $\mathcal{V}_{2}(\mathcal{K})$ are oriented in such a
way that the first edge (in terms of the ordering of
$\mathcal{V}_{1}(\mathcal{K})$) in each triangle is positively acted upon by
$B_{2}$, i.e. the first non-zero entry in each column of $B_{2}$ is $+1$,
Figure 1.

Figure 1:  Example of the simplicial complex with ordering and orientation:
nodes from $\mathcal{V}_{0}(\mathcal{K})$ in orange, triangles from
$\mathcal{V}_{2}(\mathcal{K})$ in blue. Orientation of edges and triangles is
shown by arrows; the action of $B_{2}$ operator is given for both triangles.

The following definitions introduce the fundamental concepts of $k$-th
homology group and $k$-th order Laplacian. See [18] e.g. for more details.

######  Definition 2.1 (Homology group and higher-order Laplacian).

Since $\operatorname{im}B_{k+1}\subset\ker B_{k}$, the quotient space
$\mathcal{H}_{k}=\nicefrac{{\ker B_{k}}}{{\operatorname{im}B_{k+1}}}$, known
as $k$-th homology group, is correctly defined and the following isomorphisms
hold

$\mathcal{H}_{k}\cong\ker B_{k}\cap\ker
B_{k+1}^{\top}\cong\ker\left(B_{k}^{\top}B_{k}+B_{k+1}B_{k+1}^{\top}\right).$

The matrix $L_{k}=B_{k}^{\top}B_{k}+B_{k+1}B_{k+1}^{\top}$ is called the
$k$-th order _graph Laplacian_ ; the two terms
$L_{k}^{\downarrow}=B_{k}^{\top}B_{k}$ and
$L_{k}^{\uparrow}=B_{k+1}B_{k+1}^{\top}$ are referred to as the _down-
Laplacian_ and the _up-Laplacian_ , respectively.

The homology group $\mathcal{H}_{k}$ describes the $k$-th topology of the
simplicial complex $\mathcal{K}$: $\beta_{k}=\dim\mathcal{H}_{k}=\dim\ker
L_{k}$ coincides exactly with the number of $k$-dimensional holes in the
complex, known as the _ $k$-th Betti number_. In the case $k=0$, the operator
$L_{0}=L_{0}^{\uparrow}$ is exactly the classical graph Laplacian whose kernel
corresponds to the _connected components_ of the graph, while
$L_{0}^{\downarrow}=0$. For $k=1$ and $k=2$, the elements of $\ker L_{1}$ and
$\ker L_{2}$ describe the simplex 1-dimensional holes and voids respectively,
and are frequently used in the analysis of trajectory data, [23, 1].

Although more frequently found in their purely combinatorial form, the
definitions of simplicial complexes, homology groups, and higher-order
Laplacians admit a generalization to the weighted case. For the sake of
generality, in the rest of the work, we use the following notion of weighted
boundary operators (and thus weighted simplicial complexes), as considered in
e.g. [13].

######  Definition 2.2 (Weighted and normalised boundary matrices).

For _weight functions_
$w_{k}:\mathcal{V}_{k}(\mathcal{K})\mapsto\mathds{R}_{+}\cup\\{0\\}$, define
the diagonal weight matrix $W_{k}\in\mathrm{Mat}_{m_{k}\times m_{k}}$ as
$(W_{k})_{ii}=\sqrt{w_{k}(\sigma_{i})}$. Then the weighting scheme for the
boundary operators upholding the Hodge algebras (1) is given by:

(2) $B_{k}\mapsto W_{k-1}^{-1}B_{k}W_{k}$

Note that, with the weighting scheme Equation 2, the dimensionality of the
homology group is preserved, $\dim\ker L_{k}=\dim\ker\widehat{L}_{k}$, [13] as
well as the fundamental property of topology Equation 1. For the sake of
simplicity, in the following we will omit writing the weighting matrices
explicitly and will write $B_{k}$ to denote a generic possibly weighted
boundary operator $B_{k}=W_{k-1}^{-1}B_{k}W_{k}$.

##  3 Linear systems

Note that $k$-th order Laplacians $L_{k}$, $k>0$, are matrices with a large
number of zero entries since $m_{k}\ll m_{k-1}^{2}$ asymptotically (e.g. the
number of triangles $m_{2}$ is always bounded by $m_{1}^{3/2}$). In
particular, in analogy with the standard $k=0$ case, we say that the
simplicial complex $\mathcal{K}$ and the corresponding Laplacians $L_{k}$ as
sparse if $m_{k}=\mathcal{O}(m_{k-1}\log m_{k-1})$, that is the number of
simplices of higher order is comparable up to a constant times the number of
simplices of lower order times its logarithm. This is, for example, the case
for structured simplicies such as trees for the classical Laplacian, or
triangulations for $L_{1}$; similarly, one can show the existence of a sparser
simplicial complex approximating the original one with the
$m_{k}=\mathcal{O}(m_{k-1}\log m_{k-1})$ sparsity pattern, [21, 25].

The spectrum of the $k$-th order Laplacian defines the $k$-th topology of the
complex $\mathcal{K}$. In most applications, one is actually interested in
just a few nonzero eigenvalues, possibly in the leftmost or the rightmost part
of the spectrum [12, 13, 22]. For large and sparse simplicial complexes, the
methods of choice for the computation of these eigenspaces are iterative
eigensolvers that require solving a (possibly large) number of linear systems
(least-squares problems) of the form $L_{k}\boldsymbol{x}=\boldsymbol{f}$, for
suitable $\boldsymbol{x},\boldsymbol{f}\bot\ker L_{k}$.

The following elementary result shows that solving
$L_{k}\boldsymbol{x}=\boldsymbol{f}$ can be brought down to the solution of
$L_{k}^{\uparrow}\boldsymbol{x}=\boldsymbol{f}$. We provide the proof for
completeness.

######  Theorem 3.1 (Joint $k$-Laplacian solver).

The linear system $L_{k}\boldsymbol{x}=\boldsymbol{f}$ can be reduced to a
sequence of consecutive least square problems for isolated up-Laplacians.
Precisely, $\boldsymbol{x}$ is a solution of system (3),

(3)
$L_{k}\boldsymbol{x}=\boldsymbol{f}\qquad\text{s.t.}\qquad\boldsymbol{x},\boldsymbol{f}\perp\ker
L_{k}$

if and only if it can be written as
$\boldsymbol{x}=B_{k}^{\top}\boldsymbol{u}+\boldsymbol{x}_{2}$, where:

$\\!\widehat{\boldsymbol{u}}=\underset{\boldsymbol{z}}{\operatorname{arg\,min}}\left\|L_{k-1}^{\uparrow}\boldsymbol{z}-B_{k}\boldsymbol{f}_{1}\right\|,\quad\\!\\!\boldsymbol{u}=\underset{\boldsymbol{z}}{\operatorname{arg\,min}}\left\|L_{k-1}^{\uparrow}\boldsymbol{z}-\widehat{\boldsymbol{u}}\right\|,\quad\\!\\!\boldsymbol{x}_{2}=\underset{\boldsymbol{y}}{\operatorname{arg\,min}}\left\|L_{k}^{\uparrow}\boldsymbol{y}-\boldsymbol{f}_{2}\right\|$

and $\boldsymbol{f}=\boldsymbol{f}_{1}+\boldsymbol{f}_{2}$ with
$\boldsymbol{f}_{1}=B_{k}^{\top}\boldsymbol{z}_{1}$,
$\boldsymbol{z}_{1}=\underset{\boldsymbol{z}}{\operatorname{arg\,min}}\left\|L_{k-1}^{\uparrow}\boldsymbol{z}-B_{k}\boldsymbol{f}\right\|$.

######  Proof 3.2.

For the $k$-th order Hodge Laplacian $L_{k}$, it holds (see [18] e.g.):

(4) $\mathds{R}^{m_{k-1}}=\hbox
to0.0pt{$\displaystyle\overbrace{\phantom{\operatorname{im}B_{k}^{\top}\oplus\ker
L_{k}}}^{\ker
B_{k+1}^{\top}}$\hss}\operatorname{im}B_{k}^{\top}\oplus\underbrace{\ker
L_{k}\oplus\operatorname{im}B_{k+1}}_{\ker B_{k}}\,.$

Thus,
$\boldsymbol{x}=\boldsymbol{x}_{1}+\boldsymbol{x}_{2}=B_{k}^{\top}\boldsymbol{u}+B_{k+1}\boldsymbol{v}$
for some $\boldsymbol{u}$, $\boldsymbol{v}$. Then, the system
$L_{k}^{\uparrow}\boldsymbol{x}=\boldsymbol{f}$ is equivalent to:

$\begin{cases}L_{k}^{\downarrow}B_{k}^{\top}\boldsymbol{u}=B_{k}^{\top}B_{k}B_{k}^{\top}\boldsymbol{u}=\boldsymbol{f}_{1}\\\
L_{k}^{\uparrow}B_{k+1}\boldsymbol{v}=B_{k+1}B_{k+1}^{\top}B_{k+1}\boldsymbol{v}=\boldsymbol{f}_{2}\\\
\boldsymbol{f}_{1}+\boldsymbol{f}_{2}=\boldsymbol{f}\end{cases}$

where $\boldsymbol{f}_{1}$ and $\boldsymbol{f}_{2}$ are the Hodge
decomposition of the right hand side $\boldsymbol{f}$, with
$\boldsymbol{f}_{1}\in\operatorname{im}B_{k}^{\top}$ and
$\boldsymbol{f}_{2}\in\operatorname{im}B_{k+1}$. Then
$\boldsymbol{f}=B_{k}^{\top}\boldsymbol{z}_{1}+B_{k+1}\boldsymbol{z}_{2}$ and,
after multiplication by $B_{k}$,
$B_{k}B_{k}^{\top}\boldsymbol{z}_{1}=L_{k-1}^{\uparrow}\boldsymbol{z}_{1}=B_{k}\boldsymbol{f}\iff\min_{\boldsymbol{z}_{1}}\left\|L_{k-1}^{\uparrow}\boldsymbol{z}_{1}-B_{k}\boldsymbol{f}\right\|$
and $\boldsymbol{f}_{2}=\boldsymbol{f}-B_{k}^{\top}\boldsymbol{z}_{1}$.

Finally, we note that equation
$B_{k}^{\top}B_{k}B_{k}^{\top}\boldsymbol{u}=\boldsymbol{f}_{1}\iff(L_{k-1}^{\uparrow})^{2}\boldsymbol{u}=B_{k}\boldsymbol{f}_{1}$
can be solved by two consecutive LS-problems:

$\min_{\widehat{\boldsymbol{u}}}\left\|L_{k-1}^{\uparrow}\widehat{\boldsymbol{u}}-B_{k}\boldsymbol{f}_{1}\right\|,\qquad\min_{{\boldsymbol{u}}}\left\|L_{k-1}^{\uparrow}{\boldsymbol{u}}-\widehat{\boldsymbol{u}}\right\|$

which corresponds to the solution of the down part of the original system.

Theorem 3.1 reduces the problem of solving a linear system (least-squares
problem) for the entire Laplacian matrix $L_{k}$ to the problem of solving a
system for the up-Laplacian $L_{k}^{\uparrow}$. Thus, in the next sections, we
propose an efficient combinatorial Cholesky-like preconditioner for
$L_{k}^{\uparrow}$.

##  4 Efficient solver for the up-Laplacian

When the up-Laplacian least-square problem
$\min_{\boldsymbol{x}}\|L_{k}^{\uparrow}\boldsymbol{x}-\boldsymbol{f}\|$ is
solved with standard fast iterative methods such as CGLS, [2, 14], the
convergence rate and overall behavior of the solver highly depend on the
condition number of $L_{k}^{\uparrow}$, defined as

$\kappa(L_{k}^{\uparrow})=\frac{\sigma_{\max}(L_{k}^{\uparrow})}{\sigma_{\min}(L_{k}^{\uparrow})}\,,$

where $\sigma_{\max}$ and $\sigma_{\min}$ are the largest and smallest
positive singular values of $L_{k}^{\uparrow}$.

In order to reduce $\kappa(L_{k}^{\uparrow})$, we move from
$\min_{\boldsymbol{x}}\|L_{k}^{\uparrow}\boldsymbol{x}-\boldsymbol{f}\|$ to
the symmetrically preconditioned system

$\min_{\boldsymbol{x}}\|\left(C^{+}L_{k}^{\uparrow}(C^{+})^{\top}\right)\left(C^{\top}\boldsymbol{x}\right)-C^{+}\boldsymbol{f}\|$

where $C^{+}$ denotes the Moore-Penrose pseudoinverse of $C$, and $C$ is
chosen so that (a) the transition between unconditioned and preconditioned
systems is bijective, (b) the matrix $C^{+}L_{k}^{\uparrow}(C^{+})^{\top}$ is
better conditioned than the initial $L_{k}^{\uparrow}$, and (c) the pseudo-
inverse of the preconditioner $C$ can be efficiently implemented. In
particular, a cheap matvec of the pseudo-inverse $C^{+}\boldsymbol{f}$ is
guaranteed in the case of $C$ being a lower-triangular matrix. This is the
case, for example, for the _Cholesky_ and the _incomplete Cholesky_
preconditioners [10, 15, 20].

Following this line of approach, a variety of techniques have been developed
for preconditioning the Laplacian matrix $L_{0}^{\uparrow}=L_{0}$. In
particular, one of the most prominent approaches is based on the celebrated
spectral sparsification result by Spielman and Srivastava, [25]. This result
has been recently extended to the higher-order Laplacians, [21]. However,
while the step from sparsification to preconditioner is well understood in the
$k=0$ case, when $k>0$, efficiently computing an effective preconditioner from
a sparsified simplicial complex is not trivial.

In the next subsection, we review how a sparse approximate Cholesky
preconditioner is computed from the spectral sparsifier in the graph case. We
then highlight the difficulties for $k>0$ and will devote the rest of the
paper to design of an efficient approach for $k=1$.

###  4.1 Schur Complements and Cholesky preconditioner for $k=0$

We start with a useful observation about the structure of the Gaussian
elimination process to form the Cholesky factor. Let $A$ be a real symmetric
positive definite matrix. Then its _Schur complements_ $S_{i}$ obtained in the
process of Cholesky factorization via Gaussian elimination can be defined
recursively as follows: Let $\delta_{i}$ be the $i$-th canonical vector,
$(\delta_{i})_{i}=1$ and $(\delta_{i})_{j}=0$, $i\neq j$, then set $S_{0}=A$
and for $i=1,2,3,\dots$

(5) $\displaystyle
S_{i}=S_{i-1}-\frac{1}{\delta_{i}^{\top}S_{i-1}\delta_{i}}S_{i-1}\delta_{i}\delta_{i}^{\top}S_{i-1}^{\top}$
$\displaystyle\boldsymbol{c}_{i}=\frac{1}{\sqrt{\delta_{i}^{\top}S_{i-1}\delta_{i}}}S_{i-1}\delta_{i}$

With these definitions, the Cholesky factor $C$ such that $A=CC^{\top}$ is
formed by the columns $\boldsymbol{c}_{i}$, namely

$C=[\boldsymbol{c}_{1}\,\boldsymbol{c}_{2}\,\ldots\boldsymbol{c}_{m}].$

Using induction, one can directly verify that $C$ defined above is indeed
lower triangular. Moreover, since
$S_{i}-S_{i-1}=\boldsymbol{c}_{i}\boldsymbol{c}_{i}^{\top}$ and $S_{m}=0$, we
have
$A=S_{0}=S_{0}-S_{m}=\sum(S_{i}-S_{i-1})=\sum\boldsymbol{c}_{i}\boldsymbol{c}_{i}^{\top}=CC^{\top}$.

The case of classical Laplacian $L_{0}$ is uniquely suited for Schur
complements: indeed, in the case $k=0$, each Schur complement $S_{i}$ has the
form $S_{i}=H_{i}+K_{i}$ where both $H_{i}$ and $K_{i}$ are graph Laplacians
on the same set of vertices. By this key observation, one can construct a
stochastic Cholesky preconditioner for $L_{0}$: each $i$th step of Gauss
elimination substitutes the dense Laplacian $K_{i}$ with a sparser Laplacian
on a subgraph around the vertex $i$, [16], to produce fast and sparse
approximate Cholesky multiplier. This property unfortunately holds only for
$L_{0}$ and does not carry over to $L_{k}$ for $k>0$, as we demonstrate in the
next subsection.

###  4.2 The structure of the Schur complements $S_{i}$ for $k=1$

One of the keys to the decomposition $S_{i}=H_{i}+K_{i}$ for the Schur
complements of $L_{0}$ is the decomposition $L_{0}=\sum_{e}L_{0}(e)$, where
the sum goes over all the edges of the graph and $L_{0}(e)$ is the rank-1
graph Laplacian corresponding to the edge $e$,
$L_{0}(e)=w_{1}(e)\boldsymbol{e}\boldsymbol{e}^{\top}$ with $\boldsymbol{e}$
being the corresponding incidence column from the operator $W_{0}^{-1}B_{1}$.
One can directly show an analogous decomposition into rank-1 terms for any up-
Laplacian:

######  Lemma 4.1 (Rank-1 decomposition of $L_{k}^{\uparrow}$).

For a simplicial complex $\mathcal{K}$, the following decomposition into
rank-1 up-Laplacian terms holds for any $L_{k}^{\uparrow}$:

$L_{k}^{\uparrow}=\sum_{t\in\mathcal{V}_{k}(\mathcal{K})}L_{k}^{\uparrow}(t)=\sum_{t\in\mathcal{V}_{k}(\mathcal{K})}w_{k}(t)\boldsymbol{e}_{t}\boldsymbol{e}_{t}^{\top}$

where, for each $t\in\mathcal{V}_{k}(\mathcal{K})$,
$L_{k}^{\uparrow}(t)=w_{k}(t)\boldsymbol{e}_{t}\boldsymbol{e}_{t}^{\top}$ is
the rank-1 up-Laplacian for each triangle $t$ such that (a)
$w_{k}(t)=[W_{k}^{2}]_{tt}$ is the weight of the simplex $t$, and (b)
$\boldsymbol{e}_{t}$ is the vector obtained as the action of the boundary
operator $W_{k}^{-1}B_{k}$ on $t$, i.e. the $t$-th column of the matrix
$W_{k}^{-1}B_{k}$.

Using Lemma 4.1, we obtain the following characterization for the first Schur
complement $S_{1}$ for $L_{1}^{\uparrow}$:

######  Lemma 4.2 (First Schur complement $S_{1}$ for $L_{1}^{\uparrow}$).

For the up-Laplacian $L_{1}^{\uparrow}$, the following derivation of the first
Schur complement holds:

$S_{1}=\sum_{t\mid 1\notin
t}w_{2}(t)\boldsymbol{e}_{t}\boldsymbol{e}_{t}^{\top}+\frac{1}{2\Omega_{1}}\sum_{\begin{subarray}{c}t_{1}\mid
1\in t_{1}\\\ t_{2}\mid 1\in
t_{2}\end{subarray}}w_{2}(t_{1})w_{2}(t_{2})\big{[}\boldsymbol{e}_{t_{1}}-\boldsymbol{e}_{t_{2}}\big{]}\big{[}\boldsymbol{e}_{t_{1}}-\boldsymbol{e}_{t_{2}}\big{]}^{\top}$

where $\Omega_{1}$ is the total weight of all triangles adjacent to the edge
$1$, $\Omega_{1}=\sum_{t\mid 1\in t}w_{2}(t)$.

######  Proof 4.3.

Following LABEL:eq:gauss, one needs to compute the constant
$\delta_{1}^{\top}S_{0}\delta_{1}$ and the rank-1 matrix
$S_{0}\delta_{1}\delta_{1}^{\top}S_{0}$ where $S_{0}=L_{1}^{\uparrow}$. Note
that
$\boldsymbol{e}_{t}^{\top}\delta_{1}=\frac{1}{\sqrt{w_{1}(1)}}\mathds{1}_{1\in
t}$, the indicator of the triangle $t$ having the edge $1$, since $1$ is
necessarily the first (in the chosen order) edge in the triangle; hence,
$\sum_{t}w_{2}(t)\boldsymbol{e}_{t}^{\top}\delta_{1}=\frac{1}{\sqrt{w_{1}(1)}}\sum_{t\mid
1\in t}w_{2}(t)$:

$\displaystyle\delta_{1}^{\top}S_{0}\delta_{1}$
$\displaystyle=\delta_{1}^{\top}\left(\sum_{t}w_{2}(t)\boldsymbol{e}_{t}\boldsymbol{e}_{t}^{T}\right)\delta_{1}=\sum_{t}w_{2}(t)\left(\delta_{1}^{\top}\boldsymbol{e}_{t}\right)^{2}=\frac{1}{w_{1}(1)}\sum_{t\mid
1\in t}w_{2}(t)=\frac{\Omega_{\\{1\\}\mid\varnothing}}{w_{1}(1)}$
$\displaystyle S_{0}\delta_{1}\delta_{1}^{\top}S_{0}$
$\displaystyle=\sum_{t_{1},t_{2}}w_{2}(t_{1})w_{2}(t_{2})\boldsymbol{e}_{t_{1}}\boldsymbol{e}_{t_{1}}^{\top}\delta_{1}\delta_{1}^{\top}\boldsymbol{e}_{t_{2}}\boldsymbol{e}_{t_{2}}^{\top}=\frac{1}{w_{1}(1)}\sum_{\begin{subarray}{c}t_{1}\mid
1\in t_{1}\\\ t_{2}\mid 1\in
t_{2}\end{subarray}}w_{2}(t_{1})w_{2}(t_{2})\boldsymbol{e}_{t_{1}}\boldsymbol{e}_{t_{2}}^{\top}\,.$

By symmetry,

$\sum_{\begin{subarray}{c}t_{1}\mid 1\in t_{1}\\\ t_{2}\mid 1\in
t_{2}\end{subarray}}w(t_{1})w(t_{2})\boldsymbol{e}_{t_{1}}\boldsymbol{e}_{t_{2}}^{\top}=\frac{1}{2}\sum_{\begin{subarray}{c}t_{1}\mid
1\in t_{1}\\\ t_{2}\mid 1\in
t_{2}\end{subarray}}w(t_{1})w(t_{2})\left(\boldsymbol{e}_{t_{1}}\boldsymbol{e}_{t_{2}}^{\top}+\boldsymbol{e}_{t_{2}}\boldsymbol{e}_{t_{1}}^{\top}\right),$

one can note by a straightforward arithmetic transformation that

$\boldsymbol{e}_{t_{1}}\boldsymbol{e}_{t_{2}}^{\top}+\boldsymbol{e}_{t_{2}}\boldsymbol{e}_{t_{1}}^{\top}=-\big{[}\boldsymbol{e}_{t_{1}}-\boldsymbol{e}_{t_{2}}]\big{[}\boldsymbol{e}_{t_{1}}-\boldsymbol{e}_{t_{2}}]^{\top}+\boldsymbol{e}_{t_{1}}\boldsymbol{e}_{t_{1}}^{\top}+\boldsymbol{e}_{t_{2}}\boldsymbol{e}_{t_{2}}^{\top}$

Finally, $\sum_{\begin{subarray}{c}t_{1}\mid 1\in t_{1}\\\ t_{2}\mid 1\in
t_{2}\end{subarray}}w(t_{1})w(t_{2})\boldsymbol{e}_{t_{1}}\boldsymbol{e}_{t_{1}}^{\top}=\Omega_{1}\sum_{t\mid
1\in t}w(t)\boldsymbol{e}_{t}\boldsymbol{e}_{t}^{\top}$.

The structure of the first Schur complement $S_{1}$ is reminiscent of the
classical graph Laplacian case as $S_{1}$ is formed by the sum of two terms:

$H_{1}=\sum_{t\mid 1\notin
t}w_{2}(t)\boldsymbol{e}_{t}\boldsymbol{e}_{t}^{\top},\qquad
K_{1}=\frac{1}{2\Omega_{1}}\sum_{\begin{subarray}{c}t_{1}\mid 1\in t_{1}\\\
t_{2}\mid 1\in
t_{2}\end{subarray}}w_{2}(t_{1})w_{2}(t_{2})\big{[}\boldsymbol{e}_{t_{1}}-\boldsymbol{e}_{t_{2}}\big{]}\big{[}\boldsymbol{e}_{t_{1}}-\boldsymbol{e}_{t_{2}}\big{]}^{\top}.$

As in the $k=0$ case, the term $H_{1}$, which corresponds to the portion of
the original $L_{1}^{\uparrow}$ not adjacent to the edge being eliminated
(edge $1$), is an up-Laplacian. The remaining term of the original
$L_{1}^{\uparrow}$ which consists of up-Laplacians adjacent to edge $1$, is
transformed here into the matrix $K_{1}$, which we refer to as the _cyclic
term_. Unfortunately, unlike the $0$-Laplacian case, it is easy to realize
that the cyclic term $K_{1}$ is generally not an up-Laplacian. We show this
with a simple illustrative example below.

######  Example 4.4.

Assume the simplicial complex formed by two triangles,
$\mathcal{V}_{0}(\mathcal{K})=\\{1,2,3,4\\}$,
$\mathcal{V}_{1}(\mathcal{K})=\\{12,13,14,23,24\\}$ and
$\mathcal{V}_{2}(\mathcal{K})=\\{123,124\\}$, adjacent by the first edge $12$
with $W_{1}=W_{2}=I$. Then, vectors
$\boldsymbol{e}_{t_{1}}=\begin{pmatrix}1&-1&0&1&0\end{pmatrix}^{\top}$,
$\boldsymbol{e}_{t_{2}}=\begin{pmatrix}1&0&-1&0&1\end{pmatrix}^{\top}$ and
$K_{1}=\begin{pmatrix}0&-1&1&1&-1\end{pmatrix}\cdot\begin{pmatrix}0&-1&1&1&-1\end{pmatrix}^{\top}$.
Note that $K_{1}$ is denser than $L_{1}^{\uparrow}$ and has lost the
structural balance of $L_{1}^{\uparrow}$:
$\operatorname{diag}K_{1}=K_{1}\boldsymbol{1}=\boldsymbol{0}$ where
$\operatorname{diag}L_{1}^{\uparrow}=L_{1}^{\uparrow}\boldsymbol{1}=\left(\Omega_{i}\right)_{i=1}^{m_{2}}$.

Further computation of Schur complements $S_{2},S_{3},S_{4}\dots$ unavoidably
introduces a fastly growing number of similar cyclic terms. Thus, to build an
efficient preconditioner for $L_{1}^{\uparrow}$, one needs to deal with the
cyclic terms arising on each step $S_{i}$.

To this end, note that $K_{1}$ is zero if and only if there is a unique
triangle $t\in\mathcal{V}_{2}(\mathcal{K})$ adjacent to the edge $1$.
Maintaining this property for all Schur complements is equivalent to the
concept of weak collapsibility which we introduce next.

##  5 Collapsibility of a Simplicial Complex

In this section we borrow the terminology from [29] to introduce the concept
of collapsibility and then weak collapsibility of a simplicial complex. The
simplex $\sigma\in\mathcal{K}$ is _free_ if it is a face of exactly one
simplex $\tau=\tau(\sigma)\in\mathcal{K}$ of higher order (maximal face). The
_collapse_ $\mathcal{K}\backslash\\{\sigma\\}$ of $\mathcal{K}$ at a free
simplex $\sigma$ is the operation of reducing $\mathcal{K}$ to
$\mathcal{K}^{\prime}$, where $\mathcal{K}^{\prime}=\mathcal{K}-\sigma-\tau$;
namely, this is the operation of removing a simplex $\tau$ having an
accessible (not included in another simplex) face $\sigma$. A sequence of
collapses done at the simplicies $\Sigma=\\{\sigma_{1},\sigma_{2},\ldots\\}$
is called a _collapsing sequence_ ; formally:

######  Definition 5.1 (Collapsing sequence).

Let $\mathcal{K}$ be a simplicial complex. We say that
$\Sigma=\\{\sigma_{1},\sigma_{2},\ldots\\}$ is a _collapsing sequence_ for
$\mathcal{K}$ if $\sigma_{1}$ is free in $\mathcal{K}=\mathcal{K}^{(1)}$ and
each $\sigma_{i}$, $i>1$, is free at
$\mathcal{K}^{(i)}=\mathcal{K}^{(i-1)}\backslash\\{\sigma_{i}\\}$. The
resulting complex $\mathcal{L}$ obtained collapsing $\mathcal{K}$ at $\Sigma$
is denoted by $\mathcal{L}=\mathcal{K}\backslash\Sigma$.

Note that, by definition, every collapsing sequence $\Sigma$ has a
corresponding sequence
$\mathds{T}=\\{\tau(\sigma_{1}),\tau(\sigma_{2}),\ldots\\}$ of maximal faces
being collapsed at every step. The notion of _collapsible simplicial complex_
is defined in [29] as follows

######  Definition 5.2 (Collapsible simplicial complex).

The simplicial complex $\mathcal{K}$ is collapsible if there exists a
collapsing sequence $\Sigma$ such that $\mathcal{K}$ collapses to a single
vertex, i.e. $\mathcal{K}\backslash\Sigma=\\{v\\}$ for some
$v\in\mathcal{V}_{0}(\mathcal{K})$.

While least square problems with collapsible simplicial complexes can be
solved directly in an efficient way, [6], collapsibility is a strong
requirement for a simplicial complex. In fact, determining whether the complex
is collapsible is in general _NP-complete_ , [27], even though it an be almost
linear for a specific set of families of $\mathcal{K}$, [6]. Moreover,
simplicial complexes are rarely collapsible, as we discuss in the following.

Next, we recall the concept of a $d$-Core, [26]:

######  Definition 5.3 ($d$-Core).

A $d$-Core is a subcomplex of $\mathcal{K}$ such that every simplex of
dimension $d-1$ belongs to at least $2$ $d$-simplices.

So, for example, a $2$-Core of a $2$-skeleton $\mathcal{K}$, is a subcomplex
of the original complex $\mathcal{K}$ such that every edge from
$\mathcal{V}_{1}(\mathcal{K})$ belongs to at least $2$ triangles from
$\mathcal{V}_{2}(\mathcal{K})$. Finally, we say that $\mathcal{K}$ is
$d$-collapsible if it can be collapsed only by collapses at simplices
$\sigma_{i}$ of order smaller than $d$, i.e. $\dim\sigma_{i}\leq d-1$. Then we
have the following criterium:

######  Lemma 5.4 ([19]).

$\mathcal{K}$ is $d$-collapsible if and only if it does not contain any
$d$-core.

The $d$-Core is the generalization of the cycle for the case of
$1$-collapsibility of a classical graph, and finding a $d$-Core inside a
complex $\mathcal{K}$ is neither trivial nor computationally cheap. Note that
a $d$-Core is in general dense, due to its definition, and does not have a
prescribed structure. We illustrate simple exemplary cores in the case of
$d=2$ in Figure 2, hinting at the combinatorial many possible configurations
for a general $d$-Core, for $d\geq 2$.

Figure 2:  $2$-Core, examples: all $3$-cliques in graphs are included in
corresponding $\mathcal{V}_{2}(\mathcal{K})$.

Due to the appearance of $d$-Core being a local phenomenon, an arbitrary
simplicial complex $\mathcal{K}$ tends to contain $2$-Cores as long as
$\mathcal{K}$ is denser than the trivially collapsible case.

While collapsibility is a strong requirement, in the next section, we show
that a weaker condition is enough to efficiently design a preconditioner for
any “sparse enough” simplicial complex in the $2$-skeleton case. Whilst the
extension to $k>2$ is not hard, a careful reconsideration of necessary
modifications in the proof of Theorem 5.6 may be required.

###  5.1 Weak collapsibility

Let the complex $\mathcal{K}$ be restricted up to its $2$-skeleton,
$\mathcal{K}=\mathcal{V}_{0}(\mathcal{K})\cup\mathcal{V}_{1}(\mathcal{K})\cup\mathcal{V}_{2}(\mathcal{K})$,
and assume $\mathcal{K}$ is collapsible. Then a collapsing sequence $\Sigma$
necessarily involves collapses at simplices $\sigma_{i}$ of different orders:
at edges (eliminating _edges_ and _triangles_) and at vertices (eliminating
_vertices_ and _edges_). One can show that for a given collapsing sequence
$\Sigma$ there is a reordering $\tilde{\Sigma}$ such that
$\dim\tilde{\sigma_{i}}$ in the reordered sequence are non-increasing, [6,
Lemma 2.5]. Namely, if such a complex is collapsible, then there is a
collapsible sequence $\Sigma=\\{\Sigma_{1},\Sigma_{0}\\}$ where $\Sigma_{1}$
contains all the collapses at edges first and $\Sigma_{0}$ is composed of
collapses at vertices. Note that the partial collapse
$\mathcal{K}\backslash\Sigma_{1}=\mathcal{L}$ eliminates all the triangles in
the complex, $\mathcal{V}_{2}(\mathcal{L})=\varnothing$; otherwise, the whole
sequence $\Sigma$ is not collapsing $\mathcal{K}$ to a single vertex. Since
$\mathcal{V}_{2}(\mathcal{L})=\varnothing$, the associated up-Laplacian
$L_{1}^{\uparrow}(\mathcal{L})=0$.

######  Definition 5.5 (Weakly collapsible complex).

A simplicial complex $\mathcal{K}$ restricted to its $2$-skeleton is called
_weakly collapsible_ , if there exists a collapsing sequence $\Sigma_{1}$ such
that the simplicial complex $\mathcal{L}=\mathcal{K}\backslash\Sigma_{1}$ has
no simplices of order $2$, i.e. $\mathcal{V}_{2}(\mathcal{L})=\varnothing$ and
$L_{1}^{\uparrow}(\mathcal{L})=0$.

Note that while a collapsible complex is necessarily weakly collapsible, the
opposite does not hold. Consider the following example in Figure 3: the
initial complex is weakly collapsible either by a collapse at $[3,4]$ or at
$[2,4]$. After this, the only available collapse is at the vertex $[4]$
leaving the uncollapsible $3$-vertex structure.

Figure 3: Example of weakly collapsible but not collapsible simplicial complex

For a given simplicial complex $\mathcal{K}$, one can use the following greedy
algorithm in order to find a collapsing sequence $\Sigma$ and test for
collapsibility: at each iteration perform any of possible collapses; in the
absence of free edges, the complex should be considered not collapsible. This
greedy procedure is illustrated in Algorithm 1.

Specifically, let $\Delta_{\sigma}$ be a set of triangles of $\mathcal{K}$
containing the edge $\sigma$, $\Delta_{\sigma}=\\{t\mid
t\in\mathcal{V}_{2}(\mathcal{K})\text{ and }\sigma\in t\\}$. Then the edge
$\sigma$ is free and has only one adjacent triangle $\tau=\tau(\sigma)$ iff
$|\Delta_{\sigma}|=1$; let $F=\\{\sigma\mid|\Delta_{e}|=1\\}$ is a set of all
free edges. Note that $|\Delta_{e}|\leq m_{0}-2=\mathcal{O}(m_{0})$.

Algorithm 1 GREEDY_COLLAPSE($\mathcal{K}$): greedy algorithm for the weak
collapsibility

1:initial set of free edges $F$, adjacency sets
$\\{\Delta_{\sigma_{i}}\\}_{i=1}^{m_{1}}$

2:$\Sigma=[\;],\;\mathds{T}=[\;]$ $\triangleright$ initialize the collapsing
sequence

3:while $F\neq\varnothing$ and  $\mathcal{V}_{2}(\mathcal{K})\neq\varnothing$
do

4: $\sigma\leftarrow\texttt{pop}(F)$, $\;\tau\leftarrow\tau(\sigma)$
$\triangleright$ pick a free edge $\sigma$

5: $\mathcal{K}\leftarrow\mathcal{K}\backslash\\{\sigma\\}$,
$\;\Sigma\leftarrow[\,\Sigma\;\;\sigma\,],\;\mathds{T}\leftarrow[\mathds{T}\;\tau]$
$\triangleright$ $\tau$ is a triangle being collapsed;
$\tau=[\sigma,\sigma_{1},\sigma_{2}]$

6: $\Delta_{\sigma_{1}}\leftarrow\Delta_{\sigma_{1}}\backslash\tau$,
$\;\Delta_{\sigma_{2}}\leftarrow\Delta_{\sigma_{2}}\backslash\tau$
$\triangleright$ remove $\tau$ from adjacency lists

7: $F\leftarrow F\cup\\{\sigma_{i}\,|\,i=1,2\text{ and
}|\Delta_{\sigma_{i}}|=1\\}$ $\triangleright$ update $F$ if any of
$\sigma_{1}$ or $\sigma_{2}$ has become free

8:end while

9:return $\mathcal{K},\,\Sigma,\,\mathds{T}$

Algorithm 1 recursively picks up a free edge $\sigma$ from the set $F$,
performs a collapse $\mathcal{K}\leftarrow\mathcal{K}\backslash\\{\sigma\\}$
and updates $F$ for the collapsed subcomplex. Then, the greedy approach may
fail to find the collapsing sequence only if it gets stuck on the collapsible
complex, so the order of collapses matters. We demonstrate the validity of the
greedy Algorithm 1 in the following theorem:

######  Theorem 5.6.

Algorithm 1 finds a weakly collapsible sequence of $2$-skeleton simplicies in
polynomial time. Additionally, it finds a collapsing sequence if and only if
the simplicial complex is collapsible.

######  Proof 5.7.

Clearly, Algorithm 1 runs polynomially with respect to the number of simplexes
in $\mathcal{K}$, so its consistency automatically yields polynomiality.

The failure of the greedy algorithm would indicate the existence of a weakly
collapsible complex $\mathcal{K}$ such that the greedy algorithm gets stuck at
a $2$-Core, which is avoidable for another possible order of collapses. Among
all such complexes, let $\mathcal{K}$ be any minimal one with respect to the
number of triangles $m_{2}$. Then there exiss a free edge
$\sigma\in\mathcal{V}_{1}(\mathcal{K})$ such that
$\mathcal{K}\backslash\\{\sigma\\}$ is _collapsible_ and another free edge
$\sigma^{\prime}\in\mathcal{V}_{2}(\mathcal{K})$ such that
$\mathcal{K}\backslash\\{\sigma^{\prime}\\}$ is _not collapsible_.

Note that if $\mathcal{K}$ is minimal then any pair of free edges $\sigma_{1}$
and $\sigma_{2}$ belong to the same triangle:
$\tau(\sigma_{1})=\tau(\sigma_{2})$.

Indeed, for any $\tau(\sigma_{1})\neq\tau(\sigma_{2})$,
$\mathcal{K}\backslash\\{\sigma_{1},\sigma_{2}\\}=\mathcal{K}\backslash\\{\sigma_{2},\sigma_{1}\\}$.
Let $\tau(\sigma_{1})\neq\tau(\sigma_{2})$ for at least one pair of
$\sigma_{1}$ and $\sigma_{2}$; in our assumption, either (1) both
$\mathcal{K}\backslash\\{\sigma_{1}\\}$ and
$\mathcal{K}\backslash\\{\sigma_{2}\\}$, (2) only
$\mathcal{K}\backslash\\{\sigma_{1}\\}$ or (3) none are collapsible.

In the first case, either $\mathcal{K}\backslash\\{\sigma_{1}\\}$ or
$\mathcal{K}\backslash\\{\sigma_{2}\\}$ is a smaller example of the complex
satisfying the assumption above (since a bad edge $\sigma^{\prime}$ can not be
either $\sigma_{1}$ or $\sigma_{2}$ and belongs to collapsed complexes),
hence, violating the minimality.

If only $\mathcal{K}\backslash\\{\sigma_{1}\\}$ is collapsible, then
$\mathcal{K}\backslash\\{\sigma_{2},\sigma_{1}\\}$ is not collapsible (since
$\mathcal{K}\backslash\\{\sigma_{2}\\}$ is not collapsible and $\sigma_{1}$ is
free in $\mathcal{K}$ and in $\mathcal{K}\backslash\\{\sigma_{2}\\}$ assuming
$\tau(\sigma_{1})\neq\tau(\sigma_{2})$); hence,
$\mathcal{K}\backslash\\{\sigma_{1},\sigma_{2}\\}$ is not collapsible (since
$\mathcal{K}\backslash\\{\sigma_{1},\sigma_{2}\\}=\mathcal{K}\backslash\\{\sigma_{2},\sigma_{1}\\}$
as we stated above), so $\mathcal{K}\backslash\\{\sigma_{1}\\}$ is a smaller
example of a complex satisfying the assumption.

Finally, if both $\mathcal{K}\backslash\\{\sigma_{1}\\}$ and
$\mathcal{K}\backslash\\{\sigma_{2}\\}$ are collapsible, then for known
$\sigma^{\prime}$ such that $\mathcal{K}\backslash\\{\sigma^{\prime}\\}$ is
not collapsible, $\tau(\sigma^{\prime})\neq\tau(\sigma_{1})$ or
$\tau(\sigma^{\prime})\neq\tau(\sigma_{2})$, which can be treated as the
previous point.

As a result, for $\sigma$ ( $\mathcal{K}\backslash\\{\sigma\\}$ is
collapsible) and for $\sigma^{\prime}$ (
$\mathcal{K}\backslash\\{\sigma^{\prime}\\}$ is not collapsible ) it holds
that
$\tau(\sigma)=\tau(\sigma^{\prime})\Rightarrow\sigma\cap\sigma^{\prime}=\\{v\\}$,
so after collapses $\mathcal{K}\\{\sigma\\}$ and
$\mathcal{K}\backslash\\{\sigma^{\prime}\\}$ we arrive at two identical
simplicial complexes besides the hanging edge ($\sigma^{\prime}$ or $\sigma$)
irrelevant for the weak collapsibility. Since a simplicial complex can not be
simultaneously collapsible and not collapsible, the question of weak
collapsibility can always be resolved by the greedy algorithm which has
polynomial complexity.

### Computational cost of the greedy algorithm

The complexity of Algorithm 1 rests upon the precomputed
$\sigma\mapsto\Delta_{\sigma}$ structure that de-facto coincides with the
boundary operator $B_{2}$ (assuming $B_{2}$ is stored as a sparse matrix, the
adjacency structure describes its non-zero entries). Similarly, the initial
$F$ set can be computed alongside the construction of $B_{2}$ matrix. Another
concession is needed for the complexity of the removal of elements from
$\Delta_{\sigma_{i}}$ and $F$, which may vary from $\mathcal{O}(1)$ on average
up to guaranteed $\log(|\Delta_{\sigma_{i}}|)$. As a result, given a pre-
existing $B_{2}$ operator, Algorithm 1 runs linearly, $\mathcal{O}(m_{1})$, or
almost linearly depending on the realisation, $\mathcal{O}(m_{1}\log m_{1})$.

##  6  Preconditioning through the Subsampling of the $2$-Core

In Section 4 we have shown that the efficient computation of the Cholesky
multiplier for $L_{1}^{\uparrow}$ is complicated by arising cyclic terms. At
the same time, the absence of the cyclic terms in the Schur complements is a
characteristic of weakly collapsible complexes. Using this observation, in
this section, we develop a Cholesky-like preconditioning scheme based on an
efficient Cholesky multiplier for weakly collapsible complexes.

First, we demonstrate that a weakly collapsible simplicial complex
$\mathcal{K}$ immediately yields an exact Cholesky decomposition for its up-
Laplacian:

######  Lemma 6.1.

Assume the 2-skeleton simplicial complex $\mathcal{K}$ is weakly collapsible
through the collapsing sequence $\Sigma$. Let $\mathds{T}$ be the
corresponding sequence of maximal faces and let $P_{\Sigma}$ and
$P_{\mathds{T}}$ be the permutation matrices of the two sequences, i.e. such
that $\left[P_{\Sigma}\right]_{ij}=1\iff j=\sigma_{i}$, and similarly for
$P_{\mathds{T}}$. Then $C=P_{\Sigma}B_{2}P_{\mathds{T}}$ is an exact Cholesky
mulitplier for $P_{\Sigma}L_{1}^{\uparrow}(\mathcal{K})P_{\Sigma}^{\top}$,
i.e. $P_{\Sigma}L_{1}^{\uparrow}(\mathcal{K})P_{\Sigma}^{\top}=CC^{\top}$.

######  Proof 6.2.

Note that the sequences $\Sigma$ and $\mathds{T}$ (and the multiplication by
the corresponding permutation matrices) impose only the reordering of
$\mathcal{V}_{1}(\mathcal{K})$ and $\mathcal{V}_{2}(\mathcal{K})$,
respectively; after such reordering the $i$-th edge collapses the
$i$-triangle. Hence, the first $(i-1)$ entries of the $i$-th column of the
matrix $B_{2}$ (
$\left[B_{2}\right]_{:,i}=\sqrt{w(t_{i})}\boldsymbol{e}_{t_{i}}$ ) are $0$,
otherwise one of the previous edges is not free. As a result, $C$ is lower-
triangular and by a direct computation one has
$CC^{\top}=P_{\Sigma}L_{1}^{\uparrow}(\mathcal{K})P_{\Sigma}^{\top}$.

An arbitrary simplicial complex $\mathcal{K}$ is generally not weakly
collapsible (see Figure 2). Specifically, weak collapsibility is a property of
sparse simplicial complexes with the sparsity being measured by the number of
triangles $m_{2}$ (in the weakly collapsible case $m_{2}<m_{1}$ since each
collapse at the edge eliminates exactly one triangle); hence, the removal of
triangles from $\mathcal{V}_{2}(\mathcal{K})$ can potentially destroy the
2-Core structure inside $\mathcal{K}$ and make the complex weakly collapsible.

With this observation, in order to find a cheap and effective preconditioner
for $L_{1}^{\uparrow}(\mathcal{K})$, one may search for a weakly collapsible
subcomplex $\mathcal{L}\subseteq\mathcal{K}$ and use its exact Cholesky
multiplier $C$ as a preconditioner for $\mathcal{K}$.

Specifically, such subcomplex $\mathcal{L}$ should satisfy the following
properties:

  1. (1)

it has the same set of nodes and edges,
$\mathcal{V}_{0}(\mathcal{L})=\mathcal{V}_{0}(\mathcal{K})$ and
$\mathcal{V}_{1}(\mathcal{L})=\mathcal{V}_{1}(\mathcal{K})$;

  2. (2)

triangles in $\mathcal{L}$ are subsampled,
$\mathcal{V}_{2}(\mathcal{L})\subseteq\mathcal{V}_{2}(\mathcal{K})$;

  3. (3)

$\mathcal{L}$ is weakly collapsible through some collapsing sequence $\Sigma$
and corresponding sequence of maximal faces $\mathds{T}$;

  4. (4)

$\mathcal{L}$ has the same 1-homology as $\mathcal{K}$, that is $\ker
L_{1}(\mathcal{K})=\ker L_{1}(\mathcal{L})$;

  5. (5)

the Cholesky multiplier $C=P_{\Sigma}B_{2}(\mathcal{L})P_{\mathds{T}}$
improves the condition number of $L_{1}^{\uparrow}(\mathcal{K})$, namely
$\kappa_{+}(C^{+}P_{\Sigma}L_{1}^{\uparrow}(\mathcal{K})P_{\Sigma}^{\top}(C^{+})^{\top})\ll\kappa_{+}(L_{1}^{\uparrow}(\mathcal{K}))$.

Let us comment on the conditions above. Conditions (1) and (2) are
automatically met when a subcomplex $\mathcal{L}$ is obtained from
$\mathcal{K}$ through the elimination of triangles. Condition (3) is a
structural requirement on $\mathcal{L}$ and can be guaranteed by the design of
the subcomplex using the proposed Algorithm 2. Condition (4) guarantees that
the preconditioning strategy is bijective, as we show in the next Lemma 6.3.
Finally, condition (5) asks for a better condition number and is checked
numerically in Section 7. However, whilst one can not guarantee improvement in
preconditioning quality, we can provide an explicit formula for the condition
number
$\kappa_{+}(C^{+}P_{\Sigma}L_{1}^{\uparrow}(\mathcal{K})P_{\Sigma}^{\top}(C^{+})^{\top})$
assuming an arbitrary subcomplex forms the preconditioner $C$. We provide such
a formula in the next Section 6.1 and discuss how to use it to construct a
preconditioner via heavy collapsible subcomplex.

######  Lemma 6.3 (On the conservation of the 1-homology red and condition
(4)).

For any $\mathcal{L}$ subcomplex of $\mathcal{K}$, the following statements
about 1-homology hold:

  1. (i)

if $\mathcal{L}$ satisfies conditions (1) and (2) then it can only extend the
kernel, i.e. we have $\ker L_{1}(\mathcal{K})\subseteq\ker
L_{1}(\mathcal{L})$;

  2. (ii)

if additionally $\mathcal{L}$ satisfies condition (4), then the
preconditioning scheme is bijective, in the sense that the solution to the
original least square problem
$\min_{\boldsymbol{x}}\|L_{1}^{\uparrow}(\mathcal{K})\boldsymbol{x}-\boldsymbol{f}\|$
and the one preconditioned with the Cholesky factor of $\mathcal{L}$ coincide.

######  Proof 6.4.

For (i) it is sufficient to note that the elimination of the triangle
$t\in\mathcal{V}_{2}(\mathcal{K})$ lifts the restriction
$\boldsymbol{e}_{t}^{\top}\boldsymbol{x}=0$ for $\boldsymbol{x}\in\ker
L_{1}(\mathcal{K})$; hence, if $\boldsymbol{x}\in\ker L_{1}(\mathcal{K})$,
then $\boldsymbol{x}\in\ker L_{1}(\mathcal{L})$. For (ii), note that bijection
between the systems $L_{k}^{\uparrow}\boldsymbol{x}=\boldsymbol{f}$ and
$\left(C^{+}L_{k}^{\uparrow}(C^{+})^{\top}\right)C^{\top}\boldsymbol{x}=C^{+}\boldsymbol{f}$
is guaranteed by $\ker C^{\top}=\ker L_{k}^{\uparrow}=\ker B_{k+1}^{\top}$
(assuming $\boldsymbol{x}\perp\ker L_{k}^{\uparrow}$). Then, by the spectral
inheritance principle, [13, Thm. 2.7], $\ker
L_{k}^{\uparrow}(\mathcal{X})=\ker L_{k}(\mathcal{X})\oplus
B_{k}^{\top}\cdot\operatorname{im}L_{k-1}^{\uparrow}$. The second part,
$B_{k}^{\top}\cdot\operatorname{im}L_{k-1}^{\uparrow}$, consists of the action
of $B_{k}^{\top}$ on non-zero related eigenvectors of $L_{k-1}^{\uparrow}$ and
is not dependent on $\mathcal{V}_{k+1}(\mathcal{K})$ (triangles, in case
$k=1$), hence remains preserved in the subcomplex $\mathcal{L}$. Since by
condition (4) $\ker L_{1}(\mathcal{K})=\ker L_{1}(\mathcal{L})$, the same
statement holds for up-Laplacians, $\ker L_{1}^{\uparrow}(\mathcal{K})=\ker
L_{1}^{\uparrow}(\mathcal{L})$. Since $C$ is an exact Cholesky multiplier for
$L_{1}^{\uparrow}(\mathcal{L})$, $\ker L_{1}^{\uparrow}(\mathcal{L})=\ker
C^{\top}$ and $\ker L_{1}^{\uparrow}(\mathcal{K})=\ker B_{2}^{\top}$ yielding
$\ker C^{\top}=\ker B_{2}^{\top}$ and bijectivity.

###  6.1 Preconditioning quality by the Subcomplex

Note that subcomplex $\mathcal{L}$ is fully defined by the subset $\mathds{T}$
of subsampled triangles, $\mathds{T}\subset\mathcal{V}_{2}(\mathcal{K})$, so
$\mathcal{L}=\mathcal{V}_{0}(\mathcal{K})\cup\mathcal{V}_{1}(\mathcal{K})\cup\mathds{T}$).
We introduce the following matrix notation corresponding to the subsampling:

######  Definition 6.5 (Subsampling matrix).

Let $\mathds{T}$ be a subset of triangles,
$\mathds{T}\subset\mathcal{V}_{2}(\mathcal{K})$, then $\Pi$ is a subsampling
matrix if

  * •

$\Pi$ is diagonal, $\Pi\in\mathds{R}^{m_{2}\times m_{2}}$;

  * •

$(\Pi)_{ii}=1\iff i\in\mathds{T}$; otherwise, $(\Pi)_{ii}=0$.

Assuming the subset $\mathds{T}$ of subsampled triangles (or, equivalently,
the subsampling matrix $\Pi$) is given, one needs only the triangle weight
matrix $\widehat{W}$ in order to obtain $L_{1}^{\uparrow}(\mathcal{L})$ and
the corresponding Cholesky multiplier $C$.

Generally speaking, (squared) weights $\widehat{W}_{2}^{2}$ of sampled
triangles $\mathds{T}$ may differ from the original weights $W_{2}^{2}$. Let
$\widehat{W}_{2}^{2}=W_{2}^{2}+\Delta W_{2}$, where $\Delta W_{2}$ is still
diagonal, but entries are not necessarily positive. Then one can formulate the
question of the optimal weight redistribution as the optimization problem:

$\min_{\Delta
W_{2}}\left\|L_{1}^{\uparrow}(\mathcal{L})-L_{1}^{\uparrow}(\mathcal{K})\right\|=\min_{\Delta
W_{2}}\left\|B_{2}\left[\Pi(W_{2}^{2}+\Delta W_{2})\Pi-
W_{2}^{2}\right]B_{2}^{\top}\right\|$

Here and only here, since we manipulate weights, we use the unweighted $B_{2}$
matrix so one can have explicit access to the weight matrix $W_{2}$.

######  Lemma 6.6 (Optimal weight choice for the subcomplex).

Let $\mathcal{L}$ be subcomplex of $\mathcal{K}$ with fixed corresponding
subsampling matrix $\Pi$. Then the optimal weight perturbation for the
subsampled triangles is

$\Delta W_{2}\equiv 0,$

so the best choice of weights of subcomplex is $\widehat{W}_{2}=W_{2}\Pi$.

######  Proof 6.7.

Let $\Delta W_{2}=\Delta W_{2}(t)$ where $t$ is a time parametrization; we can
compute the gradient $\nabla_{\Delta
W_{2}}\sigma_{1}\left(L_{1}^{\uparrow}(\mathcal{L})-L_{1}^{\uparrow}(\mathcal{K})\right)$
through the derivative
$\frac{d}{dt}\sigma_{1}\left(L_{1}^{\uparrow}(\mathcal{L})-L_{1}^{\uparrow}(\mathcal{K})\right)$:

$\displaystyle\dot{\sigma}_{1}$
$\displaystyle=\boldsymbol{x}^{\top}B_{2}\Pi\dot{\Delta W_{2}}\Pi
B_{2}^{\top}\boldsymbol{x}=\left\langle B_{2}\Pi\dot{\Delta W_{2}}\Pi
B_{2}^{\top},\,\boldsymbol{x}\boldsymbol{x}^{\top}\right\rangle=\mathrm{Tr}\left(B_{2}\Pi\dot{\Delta
W_{2}}\Pi B_{2}^{\top}\boldsymbol{x}\boldsymbol{x}^{\top}\right)=$
$\displaystyle=\left\langle\Pi
B_{2}^{\top}\boldsymbol{x}\boldsymbol{x}^{\top}B_{2}\Pi,\,\dot{\Delta
W_{2}}\right\rangle=\left\langle\nabla_{\Delta W_{2}}\sigma_{1},\dot{\Delta
W_{2}}\right\rangle$

By projecting onto the diagonal structure of the weight perturbation,

$\nabla_{\operatorname{diag}\Delta
W_{2}}\,\sigma_{1}=\operatorname{diag}\left(\Pi
B_{2}^{\top}\boldsymbol{x}\boldsymbol{x}^{\top}B_{2}\Pi\right).$

Note that $\operatorname{diag}\left(\Pi
B_{2}^{\top}\boldsymbol{x}\boldsymbol{x}^{\top}B_{2}\Pi\right)_{ii}=|\Pi
B_{2}^{\top}\boldsymbol{x}|^{2}_{i}$; then the stationary point is
characterized by $\Pi B_{2}^{\top}\boldsymbol{x}=0\iff\boldsymbol{x}\in\ker
L_{1}^{\uparrow}$. The latter is impossible since $\boldsymbol{x}$ is the
eigenvector corresponding to the largest eigenvalue; hence, since
$\Pi(W_{2}^{2}+\Delta W_{2})\Pi\neq W_{2}^{2}$, the optimal perturbation is
$\Delta W_{2}\equiv 0$.

We have established the optimal choice of the weights provided the subsampling
matrix $\Pi$. As a result, assuming $B_{2}$ is weighted, optimal
$L_{1}^{\uparrow}(\mathcal{L})=B_{2}\Pi B_{2}^{\top}$. Now we proceed to
characterize the quality of the preconditioning in terms of the matrix $\Pi$
in Theorem 6.9, starting with a necessary technical relation between
$\operatorname{im}B_{2}^{\top}$ and $\ker\Pi$.

######  Remark 6.8.

Knowing the optimal weight for sampled triangles from Lemma 6.6, one needs to
preserve the kernel of subsampled Laplacian

$\ker\left(B_{2}\Pi B_{2}^{\top}\right)=\ker\left(B_{2}B_{2}^{\top}\right)$

to form a correct preconditioner $C$. Since we have $\Pi=\Pi^{2}$ and $\ker
L_{1}^{\uparrow}=\ker B_{2}^{\top}$, then $\ker\left(B_{2}\Pi
B_{2}^{\top}\right)=\ker\left(\Pi B_{2}^{\top}\right)$. Additionally, $\ker
B_{2}^{\top}\subseteq\ker\left(\Pi B_{2}^{\top}\right)$, so
$\ker\left(B_{2}\Pi
B_{2}^{\top}\right)\neq\ker\left(B_{2}B_{2}^{\top}\right)\iff$ there exists
$\boldsymbol{y}\in\operatorname{im}B_{2}^{\top}$ such that
$B_{2}^{\top}\boldsymbol{y}\neq 0$ and $B_{2}^{\top}\boldsymbol{y}\in\ker\Pi$.
Then in order to preserve the kernel, one needs
$\operatorname{im}B_{2}^{\top}\cap\ker\Pi=\\{0\\}$.

######  Theorem 6.9 (Conditioning by the Subcomplex).

Let $\mathcal{L}$ be a weakly collapsible subcomplex of $\mathcal{K}$ defined
by the subsampling matrix $\Pi$ and let $C$ be a Cholesky multiplier of
$L_{1}^{\uparrow}(\mathcal{L})$ defined as in Lemma 6.1. Then the conditioning
of the symmetrically preconditioned $L_{1}^{\uparrow}$ is given by:

$\kappa_{+}\left(C^{+}P_{\Sigma}L_{1}^{\uparrow}P_{\Sigma}^{\top}\left(C^{+}\right)^{\top}\right)=\left(\kappa_{+}\left(\left(S_{1}V_{1}^{\top}\Pi\right)^{+}S_{1}\right)\right)^{2}=\left(\kappa_{+}(\Pi
V_{1})\right)^{2},$

where $V_{1}$ forms the orthonormal basis on $\operatorname{im}B_{2}^{\top}$.

######  Proof 6.10.

By Lemma 6.6, $W_{2}(\mathcal{L})=\Pi W_{2}$; then let us consider the lower-
triangular preconditioner $C=P_{\Sigma}B_{2}\Pi P_{\mathds{T}}$ for
$P_{\Sigma}L_{1}^{\uparrow}P_{\Sigma}^{\top}$; then the preconditioned matrix
is given by:

$\displaystyle
C^{+}\left(P_{\Sigma}L_{1}^{\uparrow}P_{\Sigma}^{\top}\right)(C^{+})^{\top}$
$\displaystyle=\left(P_{\Sigma}B_{2}\Pi
P_{\mathds{T}}\right)^{+}\left(P_{\Sigma}L_{1}^{\uparrow}P_{\Sigma}^{\top}\right)\left(\left(P_{\Sigma}B_{2}\Pi
P_{\mathds{T}}\right)^{+}\right)^{\top}=$
$\displaystyle=P_{\mathds{T}}^{\top}\left(B_{2}\Pi\right)^{+}L_{1}^{\uparrow}\left(\left(B_{2}\Pi\right)^{+}\right)^{\top}P_{\mathds{T}}$

Note that $P_{\mathds{T}}$ is unitary, so
$\kappa_{+}(P_{\mathds{T}}XP_{\mathds{T}}^{\top})=\kappa_{+}(X)$, hence the
principle matrix is
$\left(B_{2}\Pi\right)^{+}L_{1}^{\uparrow}\left(\left(B_{2}\Pi\right)^{+}\right)^{\top}=\left(B_{2}\Pi\right)^{+}B_{2}B_{2}^{\top}\left(\left(B_{2}\Pi\right)^{+}\right)^{\top}$.
Since $\kappa_{+}(X^{\top}X)=\kappa_{+}^{2}(X)$, then in fact one needs to
consider $\kappa_{+}\left(\left(B_{2}\Pi\right)^{+}(B_{2})\right)$. Let us
consider the SVD-decomposition for $B_{2}=USV^{\top}$; more precisely,

$B_{2}=USV^{\top}=\begin{pmatrix}U_{1}&U_{2}\end{pmatrix}\begin{pmatrix}S_{1}&0\\\
0&0\end{pmatrix}\begin{pmatrix}V_{1}^{\top}\\\
V_{2}^{\top}\end{pmatrix}=U_{1}S_{1}V_{1}^{\top}$

where $S_{1}$ is a diagonal invertible matrix. Note that $U$ and $U_{1}$ have
orthonormal columns and $S_{1}$ is diagonal and invertible, so

$(B_{2}\Pi)^{+}B_{2}=\left(SV^{\top}\Pi\right)^{+}SV^{\top}=\left(S_{1}V_{1}^{\top}\Pi\right)^{+}S_{1}V_{1}^{\top}$

By the definition of the condition number $\kappa_{+}$, one needs to compute
$\sigma^{+}_{\min}$ and $\sigma^{+}_{\max}$ where:

$\sigma^{+}_{\min}=\underset{\boldsymbol{x}\perp\ker\left(\left(S_{1}V_{1}^{\top}\Pi\right)^{+}S_{1}V_{1}^{\top}\right)}{\min}\frac{\left\|\left(S_{1}V_{1}^{\top}\Pi\right)^{+}S_{1}V_{1}^{\top}\boldsymbol{x}\right\|}{\|\boldsymbol{x}\|}$

and analogously for $\sigma^{+}_{\max}$. Here we prove the result for
$\sigma^{+}_{\min}$. The proof transfers unchanged to $\sigma^{+}_{\max}$ by
changing $\min$ to $\max$.

Note that
$\operatorname{im}B_{2}^{\top}=\operatorname{im}V_{1}=\operatorname{im}V_{1}S_{1}$,
so by Remark 6.8, $\ker\Pi\cap\operatorname{im}V_{1}S_{1}=\\{0\\}$, hence
$\ker\Pi V_{1}S_{1}=\ker V_{1}S_{1}$. Since $\ker
V_{1}S_{1}\cap\operatorname{im}S_{1}V_{1}^{\top}=\\{0\\}$, one gets $\ker\Pi
V_{1}S_{1}\cap\operatorname{im}S_{1}V_{1}^{\top}=\\{0\\}$. By the properties
of the pseudo-inverse we have that $\ker\Pi
V_{1}S_{1}=\ker\left(S_{1}V_{1}^{\top}\Pi\right)^{\top}=\ker\left(S_{1}V_{1}^{\top}\Pi\right)^{+}$;
as a result,
$\ker\left(\left(S_{1}V_{1}^{\top}\Pi\right)^{+}S_{1}V_{1}^{\top}\right)=\ker
S_{1}V_{1}^{\top}$. Since $S_{1}$ is invertible,
$\ker\left(\left(S_{1}V_{1}^{\top}\Pi\right)^{+}S_{1}V_{1}^{\top}\right)=\ker
V_{1}^{\top}$.

For $\boldsymbol{x}\in\ker
V_{1}^{\top}\Rightarrow\boldsymbol{x}\in\operatorname{im}V_{1}$, so
$\boldsymbol{x}=V_{1}\boldsymbol{y}$. Since $V_{1}^{\top}V_{1}=I$,
$\|\boldsymbol{x}\|=\|V_{1}\boldsymbol{y}\|$ and:

$\sigma^{+}_{\min}=\underset{\boldsymbol{y}}{\min}\frac{\left\|\left(S_{1}V_{1}^{\top}\Pi\right)^{+}S_{1}\boldsymbol{y}\right\|}{\|\boldsymbol{y}\|}\underset{\boldsymbol{z}=S_{1}\boldsymbol{y}}{=}\underset{\boldsymbol{z}}{\min}\frac{\left\|\left(S_{1}V_{1}^{\top}\Pi\right)^{+}\boldsymbol{z}\right\|}{\left\|S_{1}^{-1}\boldsymbol{z}\right\|}$

Note that
$\boldsymbol{v}=\left(S_{1}V_{1}^{\top}\Pi\right)^{+}\boldsymbol{z}\iff\begin{cases}S_{1}V_{1}^{\top}\Pi\boldsymbol{v}=\boldsymbol{z}\\\
\boldsymbol{v}\perp\ker S_{1}V_{1}^{\top}\Pi\end{cases}$ and $\ker
S_{1}V_{1}^{\top}\Pi=\ker V_{1}^{\top}\Pi$, so:

$\sigma^{+}_{\min}=\underset{\boldsymbol{v}\perp\ker
V_{1}^{\top}\Pi}{\min}\frac{\|\boldsymbol{v}\|}{\|V_{1}^{\top}\Pi\boldsymbol{v}\|}$

Hence
$\kappa_{+}\left(C^{+}P_{\Sigma}L_{1}^{\uparrow}P_{\Sigma}^{\top}(C^{+})^{\top}\right)=\kappa_{+}^{2}(V_{1}^{\top}\Pi)=\kappa_{+}^{2}(\Pi
V_{1})$.

###  6.2 Algorithm: Preconditioner via Heavy Collapsible Subcomplex

Following Theorem 6.9, note that $C$ is a perfect preconditioner for
$P_{\Sigma}L_{1}^{\uparrow}P_{\Sigma}^{\top}$ if $\Pi=I$, since
$\kappa_{+}(V_{1})=1$ and we would compute the exact Cholesky decomposition.
However, this is computationally prohibitive. Thus, it is natural to try to
find a sparser $\Pi$ so that $\Pi V_{1}$ is as close to $V_{1}$ as possible.
Multiplication by $\Pi$ cancels rows in $V_{1}$ corresponding to the
eliminated triangles; at the same time, rows in $V_{1}$ are scaled by the
weights of the triangles since
$\mathrm{span}V_{1}=W_{2}\operatorname{im}B_{2}^{\top}$, so we expect to $\Pi$
eliminating triangles with lowest weights to be a good choice.

Based on this observation and Theorem 6.9, we provide an algorithm for
preconditioning $L_{1}^{\uparrow}(\mathcal{K})$ which aims to eliminate
triangles of the lowest weight thus constructing what we call a heavy weakly
collapsible subcomplex $\mathcal{L}$ with largest possible total weight of
triangles. The exact Cholesky multiplier $C$ of
$L_{1}^{\uparrow}(\mathcal{L})$ is cheap to compute and is used as a
preconditioner for $L_{1}^{\uparrow}(\mathcal{K})$.

The proposed Algorithm 2 works as follows: start with an empty subcomplex
$\mathcal{L}$; then, at each step try to extend $\mathcal{L}$ with the
heaviest unconsidered triangle $t$: $\mathcal{L}\to\mathcal{L}\cup\\{t\\}$ –
here the extension includes the addition of the triangle $t$ with all its
vertices and edges to the complex $\mathcal{L}$. If the extension
$\mathcal{L}\cup\\{t\\}$ is weakly collapsible, it is accepted as the new
$\mathcal{L}$, otherwise $t$ is rejected; in either case the triangle $t$ is
removed from the list of unconsidered triangles, i.e. $t$ is not considered
for a second time.

Algorithm 2 HEAVY_SUBCOMPLEX$(\mathcal{K},W_{2})$: construction a heavy
collapsible subcomplex

1:the original complex $\mathcal{K}$, weight matrix $W_{2}$

2:$\mathcal{L}\leftarrow\varnothing,\,\mathds{T}\leftarrow\varnothing$
$\triangleright$ initial empty subcomplex

3:while there is unprocessed triangle in $\mathcal{V}_{2}(\mathcal{K})$ do

4: $t\leftarrow\texttt{nextHeaviestTriangle}(\mathcal{K},W_{2})$
$\triangleright$ e.g. iterate through $\mathcal{V}_{2}(\mathcal{K})$ sorted by
weight

5: if $\mathcal{L}\cup\\{t\\}$ is weakly collapsible then $\triangleright$ run
GREEDY_COLLAPSE$(\mathcal{L}\cup\\{t\\})$ (Algorithm 1)

6:
$\mathcal{L}\leftarrow\mathcal{L}\cup\\{t\\},\,\mathds{T}\leftarrow[\mathds{T}\;t]$
$\triangleright$ extend $\mathcal{L}$ by $t$

7: end if

8:end while

9:return $\mathcal{L},\,\mathds{T},\,\Sigma$ $\triangleright$ here $\Sigma$ is
the collapsing sequence of $\mathcal{L}$

######  Remark 6.11.

We show next that a subcomplex $\mathcal{L}$ sampled with Algorithm 2
satisfies properties (1)–(4) above: indeed,
$\mathcal{V}_{0}(\mathcal{K})=\mathcal{V}_{0}(\mathcal{L})$,
$\mathcal{V}_{1}(\mathcal{K})=\mathcal{V}_{1}(\mathcal{L})$ and $\mathcal{L}$
is weakly collapsible by construction. It is less trivial to show that the
subsampling $\mathcal{L}$ does not increase the dimensionality of $1$-homology
(see Lemma 6.3).

Assuming the opposite, the subcomplex $\mathcal{L}$ cannot have any additional
$1$-dimensional holes in the form smallest-by-inclusion cycles of more than
$3$ edges: since this cycle is not present in $\mathcal{K}$, it is “covered”
by at least one triangle $t$ which necessarily has a free edge, so
$\mathcal{L}$ can be extended by $t$ and remain weakly collapsible.
Alternatively, if the only additional hole corresponds to the triangle $t$ not
present in $\mathcal{L}$; then, reminiscent of the proof for Theorem 5.6, let
us consider the minimal by inclusion simplicial complex $\mathcal{K}$ for
which it happens. Then the only free edges in $\mathcal{L}$ are the edges of
$t$, otherwise $\mathcal{K}$ is not minimal. At the same time, in such setups
$t$ is not registered as a hole since it is an outer boundary of the complex
$\mathcal{L}$, e.g. consider the exclusion of exactly one triangle in the
tetrahedron case, Figure 2, which proves that $\mathcal{L}$ cannot extend the
1-homology of $\mathcal{K}$.

The complexity of Algorithm 2 is $\mathcal{O}(m_{1}m_{2})$ at worst which
could be considered comparatively slow: the algorithm passes through every
triangle in $\mathcal{V}_{2}(\mathcal{K})$ and performs collapsibility check
via Algorithm 1 on $\mathcal{L}$ which never has more than $m_{1}$ triangles
since it is weakly collapsible. Note that Algorithm 2 and Theorem 6.9 do not
depend on $\mathcal{K}$ being a 2-Core; moreover, the collapsible part of a
generic $\mathcal{K}$ is necessarily included in the subcomplex $\mathcal{L}$
produced by Algorithm 2. Hence a prior pass of GREEDY_COLLAPSE$(\mathcal{K})$
reduces the complex to a smaller 2-Core $\mathcal{K}^{\prime}$ with faster
HEAVY_SUBCOMPLEX$(\mathcal{K}^{\prime},W_{2})$ since
$\mathcal{V}_{1}(\mathcal{K}^{\prime})\subset\mathcal{V}_{1}(\mathcal{K})$ and
$\mathcal{V}_{2}(\mathcal{K}^{\prime})\subset\mathcal{V}_{2}(\mathcal{K})$.

We summarise the whole procedure for computing the preconditioner next:

  * •

reduce a generic simplicial complex $\mathcal{K}$ to a 2-Core
$\mathcal{K}^{\prime}$ through the collapsing sequence $\Sigma_{1}$ and the
corresponding sequence of triangles $\mathds{T}_{1}$ through the greedy
Algorithm 1;

  * •

form a heavy weakly connected subcomplex $\mathcal{L}$ from
$\mathcal{K}^{\prime}$ with the collapsing sequence $\Sigma_{2}$ and the
corresponding sequence of triangles $\mathds{T}_{2}$ using Algorithm 2;

  * •

form the preconditioner $C$ by permuting and subsampling $B_{2}$ using the
subset of triangles $\mathds{T}=\mathds{T}_{1}\cup\mathds{T}_{2}$ (that
determines the subsampling matrix $\Pi$) and the associated collapsing
sequence $\Sigma=\Sigma_{1}\cup\Sigma_{2}$, via $C=P_{\Sigma}B_{2}\Pi
P_{\mathds{T}}$.

original complex $\mathcal{K}$2-Core
$\mathcal{K}^{\prime}$$\mathcal{K}^{\prime}\subset\mathcal{K}$weakly
collapsible heavy subcomplex
$\mathcal{L}$HEAVY_SUBCOMPLEX$(\mathcal{K}^{\prime},W_{2})$ Algorithm
2GREEDY_COLLAPSE($\mathcal{K}$) Algorithm 1 projection matrix
$\Pi$preconditioner $C$

Figure 4:  The scheme of the simplicial complex transformation: from the
original $\mathcal{K}$ to the heavy weakly collapsible subcomplex
$\mathcal{L}$.

We refer to the preconditioner built in this way (see also Figure 4) as a
_heavy collapsible subcomplex_ preconditioner (HeCS-preconditioner).

##  7 Numerical experiments

We present here a number of numerical experiments to validate the performance
of the proposed preconditioning strategy. All the experiments are run using
julia on Apple M1 CPU and can be reproduced with the code available at
https://github.com/COMPiLELab/HeCS.

###  7.1 Conjugate Gradient Least Square Method

The preconditioning performance for the least square problem
$\min_{\boldsymbol{x}\perp\ker
L_{1}^{\uparrow}}\|L_{1}^{\uparrow}\boldsymbol{x}-\boldsymbol{f}\|$ is
measured on the _conjugate gradient least square method_ (CGLS), [2, 14]. The
method requires computing matvec operations for the matrix $L_{1}^{\uparrow}$
and its preconditioned alternatives; CGLS converges as
$(\sqrt{\kappa_{+}(A)}-1)/(\sqrt{\kappa_{+}(A)}+1)$ and we run it until the
infinity norm of the residual
$\boldsymbol{r}_{i}=L_{1}^{\uparrow}\boldsymbol{x}_{i}-\boldsymbol{f}$ falls
below a given threshold, i.e. $\|\boldsymbol{r}_{i}\|_{\infty}\leq\epsilon$.

###  7.2 Shifted incomplete Cholesky preconditioner

$L_{1}^{\uparrow}$ is a singular matrix. Assuming $U$ is an orthogonal basis
of $\ker L_{1}^{\uparrow}$, one can move to $L_{1}^{\uparrow}\to
L_{1}^{\uparrow}+\alpha UU^{\top}$ which can be preconditioned by non-singular
methods. Specifically, we use
$C_{\alpha}=\texttt{ichol}(L_{1}^{\uparrow}+\alpha UU^{\top})$ to compare with
the HeCS preconditioner, Figure 4.

Calculating such a shift requires efficiently computing $\ker
L_{1}^{\uparrow}$, which in principle has a complexity comparable to the
original system. On the other hand, in our setting, given the spectral
inheritance principle (see e.g. [13, Thm. 2.7]), an orthogonal basis $U$ can
be formed directly using the vectors $B_{1}^{\top}\boldsymbol{x}$,
$\boldsymbol{x}\in(\bf 1)^{\perp}$, when $\mathcal{K}$ has trivial 0- and
1-homologies (i.e. it is formed by one connected component, $\ker
L_{0}=\mathrm{span}\\{\bf 1\\}$, and has no 1-dimensional holes, $\ker
L_{1}=0$).

Note that the HeCS preconditioner instead works without requiring any
triviality of the topology of the complex.

###  7.3 Problem setting: Enriched triangulation as a simplicial complex

To illustrate the behaviour of the preconditioned system
$C^{+}P_{\Sigma}L_{1}^{\uparrow}P_{\Sigma}^{\top}(C^{+})^{\top}$, we consider
a sparse simplicial complex $\mathcal{K}$, i.e. we assume
$m_{2}=\mathcal{O}(m_{1}\ln m_{1})$. Note that the developed routine can be
applied in denser cases, although one can expect a certain loss of efficiency.
To generate problem settings within this range, $\mathcal{K}$ is synthesized
as an enriched triangulation of $N$ points on the unit square with a
prescribed edge sparsity pattern $\nu$ as follows:

  1. (1)

$\mathcal{V}_{0}(\mathcal{K})$ is formed by the corners of the unit square and
$(N-4)$ points sampled uniformly at random from $U\left([0,1]^{2}\right)$;

  2. (2)

the Delaunay triangulation of $\mathcal{V}_{0}(\mathcal{K})$ is computed; all
edges and $3$-cliques of the produced graph are included in
$\mathcal{V}_{1}(\mathcal{K})$ and $\mathcal{V}_{2}(\mathcal{K})$
respectively;

  3. (3)

$d\geq 0$ edges (excluding the outer boundary) are chosen at random and
eliminated from $\mathcal{V}_{1}(\mathcal{K})$; triangles adjacent to the
chosen edges are eliminated from $\mathcal{V}_{2}(\mathcal{K})$. As result,
produced complex $\mathcal{K}$ has a non-trivial $1$-homology;

  4. (4)

the sparsity pattern is defined as $\nu=m_{2}/q(m_{1})$ where
$q(m_{1})=\mathcal{O}(m_{1}\ln m_{1})$ is the highest density of triangles for
the sparse case; additional edges on $\mathcal{V}_{0}(\mathcal{K})$ are added
to $\mathcal{V}_{1}(\mathcal{K})$ alongside with new appearing $3$-cliques to
reach the target $m_{2}/q(m_{1})$ value. The initial sparsity pattern of the
triangulation is denoted by $\nu_{\Delta}$.

###  7.4 Heavy subcomplex and triangle weight profile

Algorithm 2, aims to build a heavy weakly collapsible subcomplex $\mathcal{L}$
such that the total weight of triangles in $\mathcal{L}$ is close to the total
weight of triangles in the original complex $\mathcal{K}$. At the same time,
the number of triangles in $\mathcal{L}$ is limited,
$m_{2}(\mathcal{L})<m_{1}(\mathcal{L})=m_{1}(\mathcal{K})$, due to the weak
collapsibility, while the number of triangles in $\mathcal{K}$ can go up to
$q(m_{1}(\mathcal{K}))$. Hence, the quality of the preconditioner is
determined by the triangle weight distribution $w(\cdot)$ on
$\mathcal{V}_{2}(\mathcal{K})$: namely, if $w_{2}(t)$ are distributed
uniformly and independently, the quality of the preconditioning falls rapidly
after $\nu>\nu_{\Delta}$ and, at the same time, the original matrix
$L_{1}^{\uparrow}$ remains well-conditioned in such configurations. Instead,
we observe that unbalanced weight distributions lead to ill-conditioned
$L_{1}^{\uparrow}$ and, thus, here we consider the following two situations:

  * •

the weights of triangles are random variables that are _independent_ of each
other and distributed as heavy-tailed Cauchy distributions or bi-modal
Gaussian distribution $\mathcal{N}(1,\sigma_{1})+\mathcal{N}(1/3,\sigma_{2})$.
This way we generate a sufficiently large number of heavy triangles and a
cluster of reducible triangles with small $w_{2}(t)$;

  * •

the weights of triangles are _topologically dependent_ , that is $w_{2}(t)$ is
a function of the weights of the neighboring triangles. A way to implement
this dependence is to set $w_{2}(t)=f(w_{1}(e_{1}),w_{1}(e_{2}),w_{1}(e_{3}))$
for the triangle $t=(e_{1},e_{2},e_{3})$. Two well-known common choices of $f$
are the min-rule $w_{2}(t)=\min(w_{1}(e_{1}),w_{1}(e_{2}),w_{1}(e_{3}))$, [13,
17], and the product rule $w_{2}(t)=w_{1}(e_{1})w_{1}(e_{2})w_{1}(e_{3})$, [4,
5]. In this way, the edge weight profile $w_{1}(\cdot)$ on
$\mathcal{V}_{1}(\mathcal{K})$ is transformed to an unbalanced distribution
$w_{2}(t)$.

###  7.5 Timings

One needs to separately discuss the time cost of the computation of the
preconditioning operator (as a tuplet $C$ and the permutation $P_{\Sigma}$)
and of the matvec computation for the preconditioned operator
$C^{+}P_{\Sigma}L_{1}^{\uparrow}P_{\Sigma}^{\top}(C^{+})^{\top}$.

Note that matvec of $L_{1}^{\uparrow}$ has the complexity of the number of
non-zero elements, $\mathcal{O}(m_{2})$; the HeCS preconditioner $C$ has $\leq
3m_{1}$ non-zero elments and lower triangular structure, so matvecs of either
$C^{+}$ and $(C^{+})^{\top}$, as well as the permutation matrix $P_{\Sigma}$,
have complexity $\mathcal{O}(m_{1})$. Hence, the complexity of each
preconditioned CGLS iteration is $\mathcal{O}(m_{2}+m_{1})$, as opposed to the
original $\mathcal{O}(m_{2})$; asymptotically one expects $m_{2}\gg m_{1}$, so
the preconditioning scheme is efficient. On the other hand, the shifted ichol
preconditioner $C_{\alpha}$ loses the sparsity due to the shift; as a result,
the application of $C_{\alpha}^{+}$ costs $\mathcal{O}(m_{0}m_{1}+m_{2})$
since $\mathrm{rank}\,U=\mathcal{O}(m_{0})$ in all considered scenarios.

In Figure 5(a) we compare the performance of the two preconditioners for the
enriched triangulation on $m_{0}=32$ vertices and a varying number of edges
$m_{1}$: the cost of one CGLS iteration for HeCS preconditioner is higher than
the original system but asymptotically approaches the matvec cost of the
unpreconditioned system, whilst the ichol-preconditioner $C_{\alpha}$ is an
order larger.

(a) Single iteration timing: the average time of matvec computation for the
original system (blue), shifted ichol (orange) and HeCS preconditioner
(green).

(b) Computation time for the heavy subcomplex preconditioner in case of
enriched triangulations on $m_{0}$ vertices

Figure 5: Timings of HeCS-perconditioner

Additionally, we demonstrate the time complexity for the computation of the
HeCS preconditioner for enriched triangulation on $m_{0}=14,20,26,32$ vertices
and varying total edge number $m_{1}$, Figure 5(b) (here $m_{\Delta}$ denotes
the number of edges in the initial triangulation); as $\mathcal{K}$ becomes a
denser simplicial complex, HEAVY_COMPLEX($\mathcal{K}$) follows the expected
complexity $\mathcal{O}(m_{1}\cdot m_{2})$. In comparison with the complexity
of Cholesky decomposition which is $\mathcal{O}(m_{1}^{3})$, in the sparse
case $m_{2}=\mathcal{O}(m_{1}\ln m_{1})$ and the overall cost of HeCS
computation is always upper-bounded by $\mathcal{O}(m_{1}^{2}\ln(m_{1}))$.

###  7.6 Performance of the preconditioner

We demonstrate the quality of HeCS preconditioner for enriched triangulations
on $m_{0}$ vertices with $d=2$ initially eliminated edges and for varying
total number of edges $m_{1}$ such that the initial sparsity
$\faktor{m_{2}}{q(m_{1})}$ is increased until the induced number of triangles
$m_{2}$ reaches $q(m_{1})=9C^{2}m_{1}\ln(4m_{1})$ with $C=\frac{1}{2}$. This
quantity is chosen in accordance to [21, 25], where it is shown to be the
highest density for the sparse case. For each pair of parameters
$\left(m_{0},\faktor{m_{2}}{q(m_{1})}\right)$, $N=25$ simplicial complexes are
generated; triangle weight profile $w_{2}(t)$ is given by the following two
scenarios:

  1. (1)

indendent triangle weights with bi-modal imbalance, where
$w_{2}(t)\sim\mathcal{N}(1,1/3)$ for the original triangulation
$t\in\mathcal{K}_{\Delta}$ and $w_{2}(t)\sim\mathcal{N}(1/2,1/6)$ otherwise,
Figure 6;

  2. (2)

dependent triangle weights through the min-rule:
$w_{2}(t)=\min\left\\{\right.w_{1}(e_{1}),\allowbreak w_{2}(e_{2}),\allowbreak
w_{3}(e_{3})\left.\right\\}$ for $t=(e_{1},e_{2},e_{3})$ and edge weights are
folded normal variables, $w_{1}(e_{i})\sim|\mathcal{N}|[0,1]$, Figure 7.

For each weight profile we measure the condition number
($\kappa_{+}(C^{+}P_{\Sigma}L_{1}^{\uparrow}P_{\Sigma}^{\top}(C^{+})^{\top})$
vs $\kappa_{+}(L_{1}^{\uparrow})$), Figures 6 and 7, left, and the
corresponding number of CGLS iteration, Figures 6 and 7, right. In the case of
the min-rule, we provide a high-performance test for matrices up to $10^{5}$
in size.

In the case of the independent triangle weights, Figure 6, HeCS
preconditioning shows gains in $\kappa_{+}$ for the first, sparser, part of
the simplicial complexes; conversely, for the min-rule profile induced by the
folded normal edges’ weights, Figure 7, developed method outperforms the
original system for all tested sparsity patterns $\faktor{m_{2}}{q(m_{1})}$.
Noticeably, HeCS preconditioning performs better in terms of the actual CGLS
iterations, Figures 6 and 4, right, than in terms of $\kappa_{+}$, and, hence,
significantly speeds up the iterative solver for $L_{1}^{\uparrow}$.

Figure 6: Preconditioning quality for enriched triangulations with a varying
number of vertices $m_{0}=16,25,50,100$ and sparsity patterns
$\faktor{m_{2}}{q(m_{1})}$ and independent bi-modal weight profile: condition
numbers $\kappa_{+}$ on the left and the number of CGLS iterations on the
right. Average results among $25$ generations are shown in solid (HeCS) and in
dash (original system); colored areas around the solid line show the
dispersion among the generated complexes.  Figure 7: Preconditioning quality
for enriched triangulations with a varying number of vertices
$m_{0}=25,100,400,1600$ and sparsity patterns $\faktor{m_{2}}{q(m_{1})}$ and
dependent min-rule weight profile with folded normal edge weights: condition
numbers $\kappa_{+}$ on the left and the number of CGLS iterations on the
right. Average results among $25$ generations are shown in solid (HeCS) and in
dash (original system); colored areas around the solid line show the
dispersion among the generated complexes.

Finally, we demonstrate the comparative performance with the shifted
incomplete Cholesky preconditioner, $C_{\alpha}$, Figure 8; here we are forced
to guarantee trivial $0$\- and $1$-homologies, so no edges are eliminated in
the triangulation, $d=0$, and we check the kernel of $L_{1}$ for triviality
after the generation of $\mathcal{K}$. Similarly to the previous results,
preconditioning with the shifted ichol $C_{\alpha}$ is more efficient than
HeCS preconditioning for the “densest part” of the considered simplicial
complexes which means that our developed method still performs better in case
of the sparser $\mathcal{K}$. Moreover, the applicability of the shifted ichol
is limited to the cases of trivial homologies which is not the case for HeCS
preconditioning.

Figure 8:  Comparison of the preconditioning quality between HeCS(solid),
shifted ichol (semi-transparent) and original system (dashed) for the enriched
triangulation on $m_{0}=25$ vertices and varying sparsity patterns $\nu$ and
dependent min-rule weight profile with uniform edge weights: condition numbers
$\kappa_{+}$ on the left and the number of CGLS iterations on the right.
Average results among $25$ generations are shown in solid (HeCS and ichol) and
in dash (original system); colored areas around the solid line show the
dispersion among the generated complexes.

### Acknowledgments

N.G. acknowledges that his research was supported by funds from the Italian
MUR (Ministero dell’Università e della Ricerca) within the PRIN 2022 Project
“Advanced numerical methods for time dependent parametric partial differential
equations with applications” and the Pro3 joint project entitled “Calcolo
scientifico per le scienze naturali, sociali e applicazioni: sviluppo
metodologico e tecnologico”. N.G. and F.T. acknowledge support from MUR-PRO3
grant STANDS and PRIN-PNRR grant FIN4GEO. F.T. also acknowledges support from
INdAM – GNCS Project, number CUP_E53C22001930001. The authors are members of
INdAM - GNCS which they thank for support.

## References

  * [1] A. R. Benson, D. F. Gleich, and J. Leskovec, Higher-order organization of complex networks, Science, 353 (2016), pp. 163–166.
  * [2] Å. Björck, T. Elfving, and Z. Strakos, Stability of conjugate gradient and lanczos methods for linear least squares problems, SIAM Journal on Matrix Analysis and Applications, 19 (1998), pp. 720–736.
  * [3] M. Black, W. Maxwell, A. Nayyeri, and E. Winkelman, Computational topology in a collapsing universe: Laplacians, homology, cohomology, in Proceedings of the 2022 Annual ACM-SIAM Symposium on Discrete Algorithms (SODA), SIAM, 2022, pp. 226–251.
  * [4] Y.-C. Chen and M. Meila, The decomposition of the higher-order homology embedding constructed from the $k$-laplacian, Advances in Neural Information Processing Systems, 34 (2021), pp. 15695–15709.
  * [5] Y.-C. Chen, M. Meila, and I. G. Kevrekidis, Helmholtzian eigenmap: Topological feature discovery & edge flow learning from point cloud data, CoRR, abs/2103.07626 (2021).
  * [6] M. B. Cohen, B. T. Fasy, G. L. Miller, A. Nayyeri, R. Peng, and N. Walkington, Solving 1-Laplacians in nearly linear time: Collapsing and expanding a topological ball, in Proceedings of the Twenty-Fifth Annual ACM-SIAM Symposium on Discrete Algorithms, SIAM, 2014, pp. 204–216.
  * [7] J. W. Demmel, Iterative Methods for Eigenvalue Problems, ch. Applied Numerical Linear Algebra, pp. 361–387.
  * [8] S. Ebli and G. Spreemann, A notion of harmonic clustering in simplicial complexes, in 2019 18th IEEE International Conference On Machine Learning And Applications (ICMLA), IEEE, 2019, pp. 1083–1090.
  * [9] L. V. Gambuzza, F. Di Patti, L. Gallo, S. Lepri, M. Romance, R. Criado, M. Frasca, V. Latora, and S. Boccaletti, Stability of synchronization in simplicial complexes, Nature communications, 12 (2021), p. 1255.
  * [10] G. H. Golub and C. F. Van Loan, Matrix computations, JHU press, 2013\.
  * [11] V. P. Grande and M. T. Schaub, Disentangling the spectral properties of the Hodge Laplacian: Not all small eigenvalues are equal, arXiv:2311.14427, (2023).
  * [12] V. P. Grande and M. T. Schaub, Topological point cloud clustering, arXiv:2303.16716, (2023).
  * [13] N. Guglielmi, A. Savostianov, and F. Tudisco, Quantifying the Structural Stability of Simplicial Homology, Journal of Scientific Computing, 97 (2023), p. 2.
  * [14] M. R. Hestenes, E. Stiefel, et al., Methods of conjugate gradients for solving linear systems, Journal of research of the National Bureau of Standards, 49 (1952), pp. 409–436.
  * [15] N. J. Higham, Reliable Numerical Computation, Oxford University Press, 1990, ch. Analysis of the Cholesky decomposition of a semi-definite matrix.
  * [16] R. Kyng and S. Sachdeva, Approximate gaussian elimination for Laplacians: fast, sparse, and simple, in 2016 IEEE 57th Annual Symposium on Foundations of Computer Science (FOCS), IEEE, 2016, pp. 573–582.
  * [17] H. Lee, M. K. Chung, H. Kang, H. Choi, S. Ha, Y. Huh, E. Kim, and D. S. Lee, Coidentification of group-level hole structures in brain networks via Hodge Laplacian, in Medical Image Computing and Computer Assisted Intervention–MICCAI 2019: 22nd International Conference, Shenzhen, China, October 13–17, 2019, Proceedings, Part IV 22, Springer, 2019, pp. 674–682.
  * [18] L.-H. Lim, Hodge laplacians on graphs, Siam Review, 62 (2020), pp. 685–715.
  * [19] D. Lofano and A. Newman, The worst way to collapse a simplex, Israel Journal of Mathematics, 244 (2021), pp. 625–647.
  * [20] T. A. Manteuffel, An incomplete factorization technique for positive definite linear systems, Mathematics of computation, 34 (1980), pp. 473–497.
  * [21] B. Osting, S. Palande, and B. Wang, Spectral sparsification of simplicial complexes for clustering and label propagation, Journal of computational geometry, 11 (2022).
  * [22] E. Ribando-Gros, R. Wang, J. Chen, Y. Tong, and G.-W. Wei, Combinatorial and Hodge Laplacians: Similarity and difference, 2023.
  * [23] M. T. Schaub, A. R. Benson, P. Horn, G. Lippner, and A. Jadbabaie, Random walks on simplicial complexes and the normalized Hodge 1-Laplacian, SIAM Review, 62 (2020), pp. 353–391.
  * [24] M. T. Schaub, J.-B. Seby, F. Frantzen, T. M. Roddenberry, Y. Zhu, and S. Segarra, Signal processing on simplicial complexes, in Higher-Order Systems, Springer, 2022, pp. 301–328.
  * [25] D. A. Spielman and N. Srivastava, Graph sparsification by effective resistances, vol. 40, SIAM, 2011, pp. 1913–1926.
  * [26] M. Tancer, d-collapsibility is np-complete for d greater or equal to 4, Chicago Journal OF Theoretical Computer Science, 3 (2010), pp. 1–28.
  * [27] M. Tancer, Recognition of collapsible complexes is NP-complete, Discrete & Computational Geometry, 55 (2016), pp. 21–38.
  * [28] J. J. Torres and G. Bianconi, Simplicial complexes: higher-order spectral dimension and dynamics, Journal of Physics: Complexity, 1 (2020), p. 015002.
  * [29] J. H. C. Whitehead, Simplicial spaces, nuclei and m-groups, Proceedings of the London mathematical society, 2 (1939), pp. 243–327.