# Quantum Algorithm for Estimating Betti Numbers Using a Cohomology Approach

Nhat A. Nghiem  Department of Physics and Astronomy, State University of New
York at Stony Brook, Stony Brook, NY 11794-3800, USA  Xianfeng David Gu
Department of Computer Science, State University of New York at Stony Brook,
Stony Brook, NY 11794, USA  Department of Applied Mathematics & Statistics,
State University of New York at Stony Brook, Stony Brook, NY 11794, USA  Tzu-
Chieh Wei  Department of Physics and Astronomy, State University of New York
at Stony Brook, Stony Brook, NY 11794-3800, USA  C. N. Yang Institute for
Theoretical Physics, State University of New York at Stony Brook, Stony Brook,
NY 11794-3840, USA

###### Abstract

Topological data analysis has emerged as a powerful tool for analyzing large-
scale data. High-dimensional data form an abstract simplicial complex, and by
using tools from homology, topological features could be identified. Given a
simplex, an important feature is so-called Betti numbers. Calculating Betti
numbers classically is a daunting task due to the massive volume of data and
its possible high-dimension. While most known quantum algorithms to estimate
Betti numbers rely on homology, here we consider the ‘dual’ approach, which is
inspired by Hodge theory and de Rham cohomology, combined with recent advanced
techniques in quantum algorithms. Our cohomology method offers a relatively
simpler, yet more natural framework that requires exponentially less qubits,
in comparison with the known homology-based quantum algorithms. Furthermore,
our algorithm can calculate its $r$-th Betti number $\beta_{r}$ up to some
multiplicative error $\delta$ with running time
$\mathcal{O}\big{(}\log(c_{r})c_{r}^{2}/(c_{r}-\beta_{r})^{2}\delta^{2}\big{)}$,
where $c_{r}$ is the number of $r$-simplex. It thus works best when the $r$-th
Betti number is considerably smaller than the number of the $r$-simplex in the
given triangulated manifold.

##  I Introduction

Topology and geometry are lasting and vibrant areas in mathematics. Despite
being abstract, topology has laid the groundwork for many important tools that
have been widely applied in science and engineering wasserman2016topological ;
bubenik2015statistical ; gu2003genus ; gu2008computational ;
gu2020computational . Among them, topological data analysis (TDA) is gaining
much attention due to its utility in revealing important topological features
of datasets, which, in reality, can be sensitive to noise or sampling errors.
Persistent homology, which is built upon algebraic topology, is a powerful
method that can probe the topological structure of the underlying dataset.
Typically, a collection of high-dimensional data (i.e., vectors) forms an
abstract simplicial complex in which the connectivity between data points is
controlled by a quantity called length scale.

Figure 1: An (abstract) simplicial complex. Each point (red points), or
0-simplex, might be a vector in a very high dimension.

Given a simplicial complex, one is interested in its Betti numbers, which
reveal the number of connected components, loops, holes, etc., of a
topological structure. The Abelian property of homology allows the formulation
of problems into a linear-algebra framework, which can be more convenient to
work with practically. However, the massive volume of the dataset and its
exponential growth in dimensionality induce a very large computational cost.

Lloyd, Garnerone, and Zanard considered the problem of computing Betti numbers
in the quantum setting, and their proposed algorithm (the so-called LGZ
algorithm) claimed to yield an exponential speedup compared to classical
methods lloyd2016quantum . Their underlying approach essentially relies on
(simplicial) homology: given a simplicial complex $\Sigma$, the $k$-th Betti
number is the rank of the $k$-th homology group $H_{k}$. The rank of such a
group (or equivalently, the dimension of such space) is revealed by analyzing
the spectrum of the so-called boundary map $\partial$, which is a linear map
between chain spaces. In Ref. lloyd2016quantum , such analysis is then done by
applying quantum techniques, such as the quantum phase estimation
kitaev1995quantum  and sparse matrix simulation berry2007efficient , leading
to the LGZ quantum algorithm. Following Ref. lloyd2016quantum , a series of
works ubaru2021quantum ; mcardle2022streamlined ; gunn2019review ;
ameneyro2022quantum  have been substantially improving the running time of the
original one lloyd2016quantum . Most recently, a striking result from
schmidhuber2022complexity  clarified some assumptions and performance of the
LGZ-like algorithms lloyd2016quantum . In particular, Schmidhuber and Lloyd
schmidhuber2022complexity  showed that computing Betti numbers is #P-hard in
the general case and approximating them up to multiplicative error is NP-hard,
ruling out the generic efficient estimation of Betti numbers. The potentially
exponential advantage of the quantum algorithm can only be obtained if the
input is a specified complex instead of a list of vertices (which means that
we have to build the complex based on the length scale). Given a specified
complex $S$, the most efficient quantum method (combining ubaru2021quantum
and gunn2019review ) to estimate the $k$-th Betti number $\beta_{k}$ to some
multiplicative error $\delta$ is

$\displaystyle\mathcal{O}\Big{(}\frac{1}{\delta}n\kappa\sqrt{\frac{|c_{k}|}{\beta_{k}}}\Big{)},$

where $\kappa$ is the conditional number of the so-called Dirac operator
schmidhuber2022complexity , $n$ is the number of points (or vertices) in $S$,
and $|c_{k}|$ is the number of $k$-simplex at the given length scale. As
pointed out in schmidhuber2022complexity , the best classical algorithm can
compute the $k$-th Betti number at a complexity of
$\mathcal{O}\big{(}|c_{k}|^{3}\big{)}$. The exponential advantage is recovered
in the regime where $\beta_{k}\rightarrow|c_{k}|$ and $k\sim\mathcal{O}(n)$.

In this article, inspired by the exciting development of recent quantum
algorithms for TDA, we attempt to provide a ‘dual’ approach to
lloyd2016quantum , namely, via the cohomology approach. Another motivation for
our method is the discrete Hodge theory, which offers a powerful tool in
computational conformal geometry gu2008computational . Thus, it is interesting
to explore applying such a theory to TDA. We focus on the scenario where a
simplicial complex $S$ is a uniform triangulation of some manifold. By uniform
triangulation, we mean that the simplicial complex is built with uniform
simplices. For instance, a uniformly triangulated 2-manifold is formed by
equilateral triangles or 2-simplices properly gluing together. The same
construction holds for arbitrary dimensions, i.e., by gluing higher
dimensional objects. According to Schmidhuber and Lloyd
schmidhuber2022complexity , this is the case where an exponential speedup is
possible, as the description of a simplicial complex is given beforehand
instead of a set of vertices and pairwise distances.

As a brief summary, our ‘dual’ approach via cohomology is built upon De Rham
cohomology and Hodge theory, as they provide a direct link between the $k$-th
homology/cohomology group and a special group called the harmonic group.
Cohomology assigns a real number to each simplex, and therefore, we could
store all those numbers in a vector, which again could be stored using a
logarithmic number of (qu)bits. Other enabling elements for our algorithm are
the recent advances in quantum algorithmic techniques harrow2009quantum ;
lloyd2014quantum ; childs2017quantum ; nghiem2023improved ; low2017optimal ;
low2019hamiltonian ; gilyen2019quantum  that allow us to handle large matrix
operations, such as multiplication, inverting of matrices, and polynomial
transformation.

This article is organized as follows. First, in section II, we revise some
important quantum tools that would contribute as main recipes of our
subsequent quantum algorithm. In section III.1, we then give some background
of de Rham cohomology and Hodge theory in both smooth and discrete settings,
which contains the key insight into our algorithm. Next, in section III.2, we
describe the classical procedure for calculating the $k$-th Betti number in a
triangulated manifold. Subsequently, in section IV, we outline the
corresponding quantum algorithm for estimating the Betti number. Some analysis
and discussion follow subsequently, e.g., in Section VI, including the
generalization to different Betti numbers and higher dimensional triangulated
manifolds, e.g., in Section V. We conclude in Section VII with some summary
and remarks for future exploration.

##  II Some Preliminaries

In this section, we introduce the main quantum ingredients that are needed to
construct our subsequent algorithm. We formally recapitulate the main results
for brevity purposes and leave out the details.

######  Definition 1 (Block Encoding Unitary)

low2017optimal ; low2019hamiltonian ; gilyen2019quantum  Let $A$ be some
Hermitian matrix of size $N\times N$. Let a unitary $U$ have the following
form:

$\displaystyle U=\begin{pmatrix}A&\cdot\\\ \cdot&\cdot\\\ \end{pmatrix}.$

Then $U$ is said to be a block encoding of matrix $A$. Equivalently, we can
write:

$\displaystyle U=\ket{\bf{0}}\bra{\bf{0}}\otimes A+\cdots$

where $\ket{\bf{0}}$ denotes the first computational basis state in some
larger Hilbert space (of, e.g., multiple qubits). It is quite obvious from the
above definition that

$\displaystyle
A_{ij}=(\bra{\bf{0}}\otimes\bra{i})U(\ket{\bf{0}}\otimes\ket{j}),$ (1)

where $A_{ij}$ refers to the entry of $A$ at $i$-th row and $j$-th column.

The unitary block encoding of some matrix $A$ as described above allows us to
apply the tool of the so-called quantum signal processing low2017optimal ;
low2019hamiltonian ; gilyen2019quantum  to achieve arbitrary transformation of
$A$ with high efficiency. For instance, as shown in low2017optimal ;
low2019hamiltonian , the method yields an optimal algorithm for simulating a
Hamiltonian $\hat{H}$, i.e., implementing the operator $\exp(-i\hat{H}t)$ for
arbitrary $t$ given a black box access to entries of $\hat{H}$. We refer the
readers to these original works for details; instead, we recapitulate the
necessary recipe for our subsequent algorithm.

###### Lemma 1

low2017optimal ; gilyen2019quantum  Given unitary block encoding of some
matrix $A$ of dimension $n$, then it is possible to construct $\exp(-iAt)$ up
to accuracy $\epsilon$ in time

$\mathcal{O}\big{(}d_{A}(t+poly\log(1/\epsilon))\big{)},$

where $d_{A}$ is the required time complexity for preparing the unitary block
encoding of $A$.

Furthermore, the following result is also useful and will also be proved in
Appendix C.

######  Lemma 2 (Block Encoding of Product of Two Matrices gilyen2019quantum )

Given a unitary block encoding of two matrices $A_{1}$ and $A_{2}$, then there
exists an efficient procedure to construct a unitary block encoding of
$A_{1}A_{2}$.

######  Lemma 3 (Efficient Matrix Application nghiem2023improved )

Given a coherent oracle access to some s-sparse, Hermitian matrix $\hat{H}$ of
dimension $n\times n$, and a given $n\times 1$ state $\ket{b}$, then there is
a unitary $U_{H}$ that acts in the following way,

$\displaystyle
U_{H}\ket{0^{m}}\ket{b}=\ket{0^{m}}(\hat{H}/s)\ket{b}+\ket{\Phi_{\perp}},$

where $\ket{\Phi_{\perp}}$ is some unimportant state (not properly normalized)
that is orthogonal to $\ket{0^{m}}(H/s)\ket{b}$, i.e.,
$\ket{0^{m}}\bra{0^{m}}\otimes\mathbf{1}\ket{\Phi_{\perp}}=0$, and
$m=\log(n)+1$. The unitary $U_{H}$ runs in time

$\displaystyle\mathcal{O}\Big{(}\log(n),poly\big{(}\log(\frac{1}{\epsilon})\big{)}\Big{)},$

where $\epsilon$ is the error tolerance.

######  Lemma 4 (Matrix Inversion childs2017quantum )

Given an oracle access to some $s$-sparse, Hermitian matrix $\hat{H}$ of
dimension $n\times n$, and a given $n\times 1$ state $\ket{b}$. Then, there is
a unitary $U$ that acts as follows:

$\displaystyle
U_{H}\ket{0^{m}}\ket{b}=\ket{0^{m}}\hat{H}^{-1}\ket{b}+\ket{\Phi_{\perp}},$

where $\ket{\Phi_{\perp}}$ is some unimportant state (not properly normalized)
that is orthogonal to $\ket{0^{m}}(\hat{H}/s)\ket{b}$, i.e.,
$\ket{0^{m}}\bra{0^{m}}\otimes\mathbf{1}\ket{\Phi_{\perp}}=0$. The unitary
$U_{H}$ runs in time

$\mathcal{O}\Big{(}\kappa
s\log(n)poly\big{(}\log(\frac{1}{\epsilon})\big{)}\Big{)},$

where $\epsilon$ is the error tolerance, and $\kappa$ is the conditional
number of $\hat{H}$.

We further remark that while in harrow2009quantum , childs2017quantum  the
inverted matrix is assumed to be square and non-singular, some subsequent
works, such as gilyen2019quantum , were built upon an advanced technique
called unitary block encoding (or quantum signal processing low2017optimal ),
allowing the pseudo inverse if $\hat{H}$ is rectangular and/or has zero
eigenvalues.

##  III Classical Framework

This section introduces key results from differential geometry that underlie
classical and quantum algorithms.

###  III.1 de Rham Cohomology and Hodge Theory

De Rham cohomology is a very important tool that encompasses both algebraic
and differential topology and aims to probe the topological structure of
smooth manifolds. While in homology, the main objects are spaces (or groups)
of chains, in de Rham cohomology, the main objects are spaces of forms. The
homology group is formed via the equivalence classes of closed chains, whereas
in de Rham cohomology, the cohomology group is formed via the equivalence
classes of closed forms. Hodge’s theory is built on an important observation
that each cohomology class has a canonical representative, the so-called
harmonic forms. A standard result in Hodge’s theory is the Hodge decomposition
theorem:

######  Theorem 1 (Hodge Decomposition)

Suppose $M$ is an $n$-dimensional closed Riemannian manifold, then

$\displaystyle\Omega_{k}={\rm Img}(d^{k-1})\oplus{\rm Img}(\delta^{k+1})\oplus
H^{k}_{\Delta}(M),$ (2)

where $\Omega_{k}$ is the space of $k$-forms, Img denotes the image of a map,
$d^{k-1}$ is the exterior derivative map:
$\Omega_{k-1}(M)\rightarrow\Omega_{k}(M)$, $\delta^{k+1}$ is the
codifferential operator: $\Omega_{k+1}\rightarrow\Omega_{k}$, and
$H^{k}_{\Delta}$ is the space of harmonic $k$-form.

In other words, if $\omega\in\Omega_{k}$, it can be expressed and decomposed
as

$\displaystyle\omega=\delta\Omega+d\eta+h,$ (3)

where $\eta\in\Omega_{k-1}$, $\Omega\in\Omega_{k+1}$, anWe may abuse the
notation throughout this work a special property that it vanishes under the
action of $d$ and $\delta$, i.e., $dh=0$ and $\delta h=0$. Most importantly,
$h$ is unique for each cohomology class, i.e., if $\omega$ and
$\omega^{\prime}$ lie in the same cohomology class, then they have the same
$h$. Throughout this work, we may abuse the notation by writing $\delta$ and
$d$ without the superscript.

Generally, de Rham cohomology and Hodge theory work for the smooth setting.
But the extension to the discrete setting can be achieved by simply replacing
de Rham cohomology with simplicial cohomology by, e.g., identifying $k$-forms
with $k$-cochains, and so on. The discrete version has been developed and
applied extensively in real applications; see, e.g., Refs gu2008computational
; gu2003genus , as $d$ and $\delta$ operators can be represented as a linear
transformation on a set of corresponding simplices. We use the same notations
for both cases, and we treat groups and vector spaces in the same manner, as a
vector space is also an Abelian group, and we only work in the Abelian
(commutative) setting.

Now, we state two important results that provide a useful insight into our
algorithm, to be described below. As mentioned above, each cohomology class
has a unique representative, and, therefore, it directly implies that the
correspondence between two spaces (or groups), as stated in the two theorems
below.

###### Theorem 2

Given an $n$-dimensional closed Riemannian manifold. The $k$-th de Rham
cohomology group is isomorphic to the harmonic $k$-form group

$\displaystyle H^{k}_{dR}(M)\approx H^{k}_{\Delta}(M).$ (4)

Another standard and important result that we will employ is the following.

###### Theorem 3

The de Rham cohomology group is isomorphic to the simplicial cohomology group

$\displaystyle H^{k}_{dR}(M)\approx H^{k}(M).$ (5)

We remark that the above two theorems are standard results in the areas of
differential geometry and algebraic topology, which were rigorously
established in some standard texts, e.g., see bott1982differential .

###  III.2 Sketch of The Method for Calculating Betti Number

Theorems 2 and 3, plus the duality between the simplicial homology and
cohomology, show that all these groups are isomorphic, which reveals a
potential approach to calculate Betti numbers of a given simplicial complex
$\sum$. Given some $k$, the $k$-th Betti number is the rank of $k$-th homology
group, which is also the rank of the $k$-th cohomology group. If we regard
them as vector spaces, then the rank becomes the dimension. We remark that the
harmonic forms also form a vector space; therefore, the dimension of such
space could be revealed if we know the maximum number of linearly independent
vectors. The Hodge decomposition theorem 2 allows us to find the harmonic form
given an initial $k$-form $\omega$, as $h=\omega-d\eta-\delta\Omega$. In fact,
the Hodge decomposition theorem works with arbitrary forms, as we will discuss
in detail subsequently. We have the following algorithm 1 for finding the
$k$-th Betti number.

1:A set of randomized $M$ $k$-forms $\\{\omega_{i}\\}\in\Omega_{k}$. A
simplicial complex $\sum$ For each $i\in\\{M\\}$, do:

2:1\. Compute the coexact term $\delta\Omega_{i}$

3:2\. Compute the exact term $d\eta_{i}$

4:3\. Find the harmonic form $h_{i}=\omega_{i}-d\eta_{i}-\delta\Omega_{i}$

5:4\. Arrange $\\{h_{i}\\}$ as a matrix $\mathbb{H}$. Find the number
$\mathcal{K}$ of linearly independent harmonic forms among all $\\{h_{i}\\}$,
which is equivalently the rank of $\mathbb{H}$.

6:The $k$-th Betti number is $\mathcal{K}$.

Algorithm 1 Algorithm for Finding the $k$-th Betti Number

In the next two sections, we describe how to execute the above procedure in
the quantum setting and will elaborate on some details in the Appendix. For
convention, throughout the text, we use $\vec{x}$ to denote an arbitrary
vector, $\ket{\vec{x}}$ is the quantum state (or normalized vector)
corresponding to $\vec{x}$, and $|\vec{x}|$ denotes its length, i.e., the
$l_{2}$-norm. We remark that sometimes we make the following abuse of
notation: we typically use $\ket{\bf 0}$ to denote necessary ancillas required
to execute algorithms without specifying the actual number of qubits. That
being said, if two registers of states being $\ket{\bf 0}$ are provided, they
might not have the same dimension. The same convention holds for any
$\ket{Garbage}$ state since they are irrelevant to our quantum algorithm for
computing the rank.

##  IV Quantum Algorithm for Estimating Betti Numbers of Triangulated
2-Manifold

The classical algorithm 1, on which our quantum algorithm is based, can, in
principle, be used to find all the Betti numbers. We shall see how its quantum
version provides a speedup. To explain our quantum algorithm, we consider
specifically the first Betti number for simplicity and illustration. The next
section describes how to generalize to higher Betti numbers. The particular
simplicial complex that we consider here is the triangulation of a 2-dim
manifold, which means that it is composed of many 2-simplices glued together
(see Fig. 1). As we will show, the only difference with higher Betti numbers
is the specific entries or coefficients of linear equations. We remind readers
that our goal now is to construct a quantum version of Algorithm 1. To improve
the readability, we will construct each step one by one. We now begin with the
computation of the coexact term and then the exact term.

###  IV.1 Computing Coexact Terms $\delta\Omega$

Denote the number of 0-simplices, 1-simplices and 2-simplices as
$c_{0},c_{1}$, and $c_{2}$, respectively (and hene $c_{r}$ denote the number
of $r$-simplices). We could simply choose $m=\max(c_{0},c_{1},c_{2})$, for a
reason that would be clear later.

The first step in Algorithm 1 is to compute the coexact term $\delta\Omega$.
From the equation (3) we have

$\displaystyle d\omega(f)=d\delta\Omega(f),$ (6)

where $f$ is some 2-simplex. The above yields the following linear equation:

$\displaystyle A\cdot\vec{\Omega}=B,$ (7)

where $A$ is a sparse matrix $\in R^{c_{2}\times c_{2}}$ which encodes the
linear action of $d\delta$, $\vec{\Omega}$ is a column vector containing all
$\\{\Omega(f_{i})\\}$, and $B$ encodes the action of $d\omega$. Both
$\vec{\Omega}$ and $B$ are $\in R^{c_{2}\times 1}$. We want to start from a
1-form vector $\vec{\omega}$ (this is a vector of size $c_{1}\times 1$ that
contains the values of 1-form $\omega$ on all 1-simplices), then $B$ can be
written in a simpler form as $B=C\cdot\vec{\omega}$ where $C\in R^{c_{2}\times
c_{1}}$, which encodes explicitly the operation $d$ with a known expression
and is also sparse, and $\vec{\omega}\in R^{c_{1}\times 1}$. Therefore, the
above equation becomes

$\displaystyle A\cdot\vec{\Omega}=C\cdot\vec{\omega}.$ (8)

In order to find $\vec{\Omega}$, in principle, we could use the quantum linear
solver harrow2009quantum . However, both $A$ and $C$ are not Hermitian, and do
not share the same dimension, and therefore, in order to apply the method, we
need to modify the system as follows harrow2009quantum ; wiebe2012quantum ,

$\displaystyle\begin{bmatrix}{\bf{0}}_{c_{2}\times c_{2}}&A\\\
A^{T}&{\bf{0}}_{c_{2}\times c_{2}}\\\
\end{bmatrix}\cdot\begin{bmatrix}{\bf{0}}_{c_{2}\times 1}\\\ \vec{\Omega}\\\
\end{bmatrix}=\begin{bmatrix}{\bf{0}}_{c_{2}\times c_{2}}&C\\\
C^{T}&{\bf{0}}_{c_{1}\times c_{1}}\\\
\end{bmatrix}\cdot\begin{bmatrix}{\bf{0}}_{c_{2}\times 1}\\\ \vec{w}\\\
\end{bmatrix}.$ (9)

By doing so we obtain a new matrix that is square and Hermitian. We denote the
equation for the new system as

$\displaystyle
A^{\prime}\cdot\vec{\Omega^{\prime}}=C^{\prime}\cdot\vec{\omega^{\prime}},$
(10)

where $A^{\prime},C^{\prime}\in R^{c\times c}$ with $c=c_{1}+c_{2}$, and
$\vec{\Omega^{\prime}}$ contains the solution $\vec{\Omega}$ in its last
$c_{2}$ entries, and zeros in the remaining. On the right-hand side,
$\vec{\omega^{\prime}}$ contains the 1-form $\omega$ in its last $c_{1}$
entries. (We note that if one desires the dimensions $c_{1}$ and $c_{2}$ to be
powers of 2, one can always pad 0’s in the corresponding matrices and vectors
above.) Then the solution has the explicit form:

$\displaystyle\vec{\Omega^{\prime}}=A^{\prime-1}\cdot
C^{\prime}\cdot\vec{\omega^{\prime}},$ (11)

Given access to entries of $A$ and $C$ (and hence, $A^{\prime}$ and
$C^{\prime}$), we can then solve for $\vec{\Omega^{\prime}}$ using the quantum
linear solver harrow2009quantum ; childs2017quantum . We remark that in this
case, aside from inverting a matrix, we also need to apply $C^{\prime}$ to
$\vec{\omega^{\prime}}$. Such application is in fact a modification of the
inverting eigenvalue technique originally introduced in harrow2009quantum  and
was also encountered in wiebe2012quantum . More recently, a highly efficient
technique for matrix multiplication was outlined in nghiem2023improved , which
is built upon the Chebyshev polynomial approach in childs2017quantum . Here,
we simply integrate the two steps together. We remark that we need to find the
coexact term $\delta\Omega$, which can be computed by another round of matrix
application, $\vec{\delta\Omega}=R_{1}\cdot\vec{\Omega}$, where $R_{1}$ is a
sparse matrix of size $c_{1}\times c_{2}$ that is the discrete version of
$\delta$. In order to apply the same techniques, we made the same modification
as above and obtain:

$\displaystyle\vec{\delta\Omega^{\prime}}$
$\displaystyle=R_{1}^{\prime}\cdot\vec{\Omega^{\prime}}$ (12)
$\displaystyle=R_{1}^{\prime}\cdot A^{\prime-1}\cdot
C^{\prime}\cdot\vec{\omega^{\prime}},$ (13)

where $R_{1}^{\prime}$ again is the isometry embedding of $R_{1}$, in a
similar form as in Eq. (11). Note that $\vec{\delta\Omega^{\prime}}$ now
contains the coexact values $\vec{\delta\Omega}$ in its first $c_{1}$ entries.
Since $w^{\prime}$ could be arbitrary, we can safely work with its quantum
state $\ket{w^{\prime}}$ as the length does not matter. Suppose we are
equipped with $\ket{w^{\prime}}$, consecutive application of matrices and
matrix inversion (see Eqn 13), which we denote as unitary $U_{\delta\Omega}$,
yields the following operation:

$\displaystyle U_{\delta\Omega}\ket{w^{\prime}}\ket{\bf
0}=\ket{\psi}=(|\vec{\delta\Omega^{\prime}}|/s_{\Omega})\cdot\ket{\bf
0}\ket{\vec{\delta\Omega^{\prime}}}+\ket{{\rm garbage}_{1}},$ (14)

where $\ket{{\rm garbage}_{1}}$ is some unimportant state that is orthogonal
to $\ket{\vec{\delta\Omega^{\prime}}}\ket{\bf 0}$; and
$s_{\Omega}=s_{R_{1}^{\prime}}s_{C^{\prime}}$ which is the product of the
sparsity of $R_{1}^{\prime}$ and $C^{\prime}$, as the application of Lemma 3
contains the sparsity of the matrix, and we remark that as we will show in
appendix A, the sparsity is known. We remind that the notation $\ket{\bf 0}$
is an abbreviation of extra registers that is required in the application of
childs2017quantum ; nghiem2023improved .

###  IV.2 Computing exact term $d\eta$

Now we proceed to the second step in Algorithm 1, which is computing the exact
term $d\eta$. This step is very similar to what we had for $\delta\Omega$.
From equation (3), we have:

$\displaystyle\delta\omega(v)=\delta d\eta(v),$ (15)

where $v$ is some 0-simplex. The above yields the following linear equation:

$\displaystyle K\cdot\vec{\eta}=D\cdot\vec{\omega},$ (16)

where $K$ is a matrix that encodes the linear action of $\delta d$. Proceeding
similarly as in the $\delta\Omega$ case, we embed the above system to a square
$c\times c$ system (note the dimension) and denote the resultant equation in
the enlarged system as

$\displaystyle\vec{d\eta^{\prime}}=R_{2}^{\prime}\cdot K^{\prime-1}\cdot
D^{\prime}\cdot\vec{\omega^{\prime}},$ (17)

where $\vec{d\eta^{\prime}}$ contains the exact term $d\eta$ in its first
$c_{1}$ entries. Application of matrix multiplication and matrix inversion in
a similar manner to what we did in the previous part results in a unitary
$U_{d\eta}$ that, given $\ket{w^{\prime}}$, yields the following:

$\displaystyle U_{d\eta}\ket{w^{\prime}}\ket{\bf
0}=\ket{\phi}=(|\vec{d\eta^{\prime}}|/s_{\eta})\cdot\ket{\bf
0}\ket{\vec{d\eta^{\prime}}}+\ket{{\rm garbage}_{2}}.$ (18)

where again, we note that $s_{\eta}$ is known.

###  IV.3 Preparing 1-form state $\ket{w^{\prime}}$

The above two subsections show that the first step in both is to prepare the
1-form state $\ket{w^{\prime}}$, which is the state of dimension $c$ that
contains the real 1-form $w$ in the last $c_{1}$ entries. We note that there
are $M$ different initial 1-form $w_{j}$, so we need to be able to efficiently
prepare many initial states $\ket{w^{\prime}_{j}}$’s.

According to the input in Algorithm 1, we note that for any $j$, entries of
$\ket{w^{\prime}_{j}}$ are randomly chosen. Therefore, we can pick a short-
depth random unitary $U_{w}$ of dimension $c\times c$. For each index state
$\ket{j}$, we have that $U_{w}\ket{j}$ is the $j$-th column of $U_{w}$, which
is also randomized. In order to filter out those top entries, i.e., making
these first $c-c_{1}$ entries to be zero, and keeping the last $c_{1}$ entries
non-zero, we can multiply the vector $U_{w}\ket{j}$ with a matrix $A$ that is
zero everywhere except those last $c_{1}$ entries on the diagonal being 1. In
other words, we have $A_{ii}=0$ for $i<c-c_{1}$ and $A_{ii}=1$ for any $i\geq
c-c_{1}$. The matrix $A$ is apparently row/column computable, which is
equivalent to being coherently accessbile. Therefore, we can apply Lemma 3 to
obtain $AU_{w}\ket{j}$, which is our desired 1-form state
$\ket{w^{\prime}_{j}}$. We have the following result:

######  Lemma 5 (Preparation of $\ket{w^{\prime}}$)

Given an input state $\ket{j}$ plus ancilla initialized in $\ket{\bf 0}$, then
it is highly efficient to construct the unitary $U_{W}$ that performs the
following:

$\displaystyle U_{W}\ket{\bf 0}\ket{j}=\ket{\bf
0}\ket{w^{\prime}_{j}}+\ket{\rm Garbage},$ (19)

The running time of $U_{W}$ is

$\mathcal{O}\Big{(}\log(c),poly(\log(1/\epsilon))\Big{)},$

where $\ket{\rm Garbage}$ is orthogonal to $\ket{\bf 0}\ket{w^{\prime}_{j}}$
and $\epsilon$ is the error tolerance.

Important Remark: In the following section, for clarity, we make a subtle
convention that we focus on the subspace spanned by $\ket{\bf 0}$ in the first
register, which means that any further stated operations on
$\ket{w^{\prime}_{j}}$ (plus possible ancilla) are actually being controlled
by $\ket{\bf 0}$ in the Equation 5 of Lemma 5. For convenience, we therefore
suppress writing out this register unless necessary.

###  IV.4 Computing Harmonic Form $h$

Now we proceed to describe our algorithm to compute $h$. We first note that
while $M$ can be arbitrary large, we choose $M=c_{1}$ in our algorithm, as it
is sufficient for our purpose.

Suppose we begin with $\ket{w^{\prime}_{j}}\ket{\bf 0}$, where $\ket{\bf 0}$
is an (additional) ancilla system (we recall the convention regarding the
working subspace that we set before). Basically we append extra $\ket{\bf 0}$
to the right of Equation 5 and consider only the subspace spanned by $\ket{\bf
0}$ that appears on the left of $\ket{w^{\prime}_{j}}$ in Eqn 5. As we said
previously, this can be done by executing any unitary on the first register of
$\ket{w^{\prime}_{j}}\ket{\bf 0}$, controlled by the second register being
$\ket{\bf 0}$ on the left). We then append an extra register initialized as:

$\frac{1}{2}\big{(}\ket{00}-\ket{01}-\ket{10}+\ket{11}\big{)},$

which can be prepared by two Hadamard gates acting on $\ket{1}\ket{1}$. Our
whole system is then:

$\displaystyle\frac{1}{2}\Big{(}\ket{00}-\ket{01}-\ket{10}+\ket{11}\Big{)}\ket{w^{\prime}_{j}}\ket{\bf
0}.$ (20)

We note that the quantum 1-form $\ket{w^{\prime}_{j}}$ has dimension $c$, and
the actual 1-form $w_{j}$ is stored in the last $c_{1}$ entries, and we need
to relocate them to the top $c_{1}$ entries so that we can do the subtraction
properly (note that $h=w-\delta\Omega-d\eta$). A simple way to achieve this
goal is to perform the basis permutation, as we only change the corresponding
basis. However, it is unknown (at least to us) whether there is any universal
and efficient way to realize arbitrary permutation. Instead, we propose a
solution that relies on matrix multiplication. Let $A$ be some symmetric
matrix such that $A_{ij}=1$ if $i\leq c_{1},j>c-c_{1}$ (note that we need to
assume $c\geq 2c_{1}$, but this is possible as we only need to enlarge the
system and set the extra coefficients to be zero), and 0 otherwise. Then, a
simple multiplication of $A$ with $\ket{w^{\prime}_{j}}$ will exactly relocate
those last $c_{1}$ entries. Since $A$ is easily computable, applying Lemma 3
allows us to implement such quantum multiplication. We have the following
lemma:

###### Lemma 6

Given a quantum state $\ket{x^{\prime}}$ of dimension $c$ that is non-zero in
the last $c_{1}$ rows, and zero otherwise. Let $\vec{x}$ be the vector with
the same entries as $\ket{x^{\prime}}$ but instead have those non-zero entries
in the first $c_{1}$ rows (with the same order of entries). Then, the
following operation can be achieved efficiently:

$\displaystyle U_{C}\ket{x^{\prime}}\ket{\bf 0}=|\vec{x}|\ket{x}\ket{\bf
0}+\ket{Garbage},$ (21)

where $\ket{Garbage}$ is some state that is orthogonal to $\ket{x}\ket{\bf
0}$.

Now we specifically look at the part $\ket{00}\ket{w^{\prime}_{j}}\ket{\bf
0}$. Using the above Lemma, we can apply the $U_{C}$, being controlled by the
first two qubits being $\ket{00}$, to the state $\ket{w^{\prime}_{j}}$ to
produce the shifted version of the vector $\ket{w^{\prime}_{j}}$, for which we
denote simply as $\ket{w_{j}}$. (Without confusion, we will use the language
that $U_{c}$ is applied to $\ket{w^{\prime}_{j}}$ controlled by $\ket{00}$.)
To be more specific, the state $\ket{00}\ket{w^{\prime}_{j}}$ is transformed
to:

$\displaystyle\ket{\phi_{1}}=\ket{00}(|\vec{w_{j}}|\ket{w_{j}}\ket{\bf
0}+\ket{Garbage}).$ (22)

Next, we consider the part $\ket{01}\ket{w^{\prime}_{j}}\ket{\bf 0}$. We aim
to apply $U_{\delta\Omega}$ to $\ket{w^{\prime}_{j}}\ket{\bf 0}$ being
controlled by the first two qubits being $\ket{01}$. We then obtain the
following state:

$\displaystyle\ket{\phi_{2}}=\ket{01}\Big{(}|\vec{\delta\Omega^{\prime}_{j}}|/s_{\Omega}\cdot\ket{\vec{\delta\Omega_{j}^{\prime}}}\ket{\bf
0}+\ket{{\rm garbage}}\Big{)}.$ (23)

Similarly, we consider the part $\ket{10}\ket{w_{j}^{\prime}}\ket{\bf 0}$. We
apply $U_{d\eta}$ being controlled by the first two qubits being $\ket{10}$ to
obtain:

$\displaystyle\ket{\phi_{3}}=\ket{10}\Big{(}|\vec{d\eta^{\prime}_{j}}|/s_{\eta}\cdot\ket{\vec{d\eta_{j}^{\prime}}}\ket{\bf
0}+\ket{{\rm garbage}}\Big{)}.$ (24)

For the last part $\ket{11}\ket{w^{\prime}_{j}}\ket{\bf 0}$, we simply flip
the first bit in the register $\ket{\bf 0}$ controlled by the register
$\ket{11}$. In other words, we transform $\ket{11}\ket{w^{\prime}_{j}}\ket{\bf
0}$ into

$\ket{11}\ket{w^{\prime}_{j}}\ket{1\bf\tilde{0}},$

where ${\bf\tilde{0}}$ denotes one less 0 from ${\bf 0}$.

So far we have transformed the equation 20 into:

$\displaystyle\ket{\phi}=\frac{1}{2}\Big{(}\ket{\phi_{1}}-\ket{\phi_{2}}-\ket{\phi_{3}}+\ket{11}\ket{w^{\prime}_{j}}\ket{1\bf\tilde{0}}\Big{)}.$
(25)

Denote the whole unitary process that transforms the initial state $\ket{\bf
0}\ket{j}$ (which actually began from Lemma 5) to the above state as
$U_{\phi}$, which is $|00\rangle\langle 00|\otimes U_{C}\otimes
I+|01\rangle\langle 01|\otimes U_{d\eta}\otimes I+|10\rangle\langle 10|\otimes
U_{d\eta}\otimes I+|11\rangle\langle 11|\otimes I\otimes X_{1}$. As a summary,
by explicitly including the previously suppressed register and all the
ancillas, we have:

$\displaystyle U_{\phi}\ket{\bf 0}\ket{00}\ket{j}\ket{\bf 0}=\ket{\bf
0}\ket{\phi}+\ket{\rm Garbage}.$ (26)

The reason why we have $\ket{\phi}$ entangled with $\ket{\bf 0}$ on the r.h.s.
is due to the convention that we have made after Lemma 5, as everything was
done being (additionally) controlled by $\ket{\bf 0}$ in Eqn. 5.

Now we use a different procedure with the state $\ket{\bf
0}\ket{00}\ket{i}\ket{\bf 0}$. Let $s_{m}$ be some integer value that is
greater than both $s_{\Omega}+s_{\eta}$, e.g., $s_{m}=s_{\Omega}+s_{\eta}+1$.
We use controlled rotation gates to transform the register $\ket{00}$ to

$\frac{1}{s_{m}}\ket{00}+\frac{s_{\Omega}}{s_{m}}\ket{01}+\frac{s_{\eta}}{s_{m}}\ket{10}+G\ket{11},$

where $G$ in the above state refers to the required normalization factor,
i.e.,

$G=\sqrt{1-\frac{1}{s_{m}^{2}}-\frac{s_{\Omega}^{2}}{s_{m}^{2}}-\frac{s_{\eta}^{2}}{s_{m}^{2}}}.$

We again note that as the linear matrices that implement $d$ and $\delta$ are
known, their sparsities are also known. The above state can be obtained from
$\ket{00}$ by the following procedure. We first use rotation gate to transform
$\ket{00}$ to:

$\ket{0}\Big{(}\sqrt{1/s_{m}^{2}+(s_{\eta}/s_{m})^{2}}\ket{0}+\sqrt{(s_{\Omega}/s_{m})^{2}+G^{2}}\ket{1}\Big{)}.$

We then perform two controlled rotation gates. The first one is controlled by
$\ket{0}$, and the second is controlled by $\ket{1}$ in the second register.
More specifically, we transform:

$\displaystyle\sqrt{1/s_{m}^{2}+(s_{\eta}/s_{m})^{2}}\ket{00}\rightarrow\sqrt{1/s_{m}^{2}+(s_{\eta}/s_{m})^{2}}\Big{(}\frac{1/s_{m}}{\sqrt{1/s_{m}^{2}+(s_{\eta}/s_{m})^{2}}}\ket{00}+\frac{s_{\eta}/s_{m}}{\sqrt{1/s_{m}^{2}+(s_{\eta}/s_{m})^{2}}}\ket{10}\Big{)},$
(27)

and

$\displaystyle\sqrt{(s_{\Omega}/s_{m})^{2}+G^{2}}\ket{01}\rightarrow\sqrt{(s_{\Omega}/s_{m})^{2}+G^{2}}\Big{(}\frac{s_{\Omega}/s_{m}}{\sqrt{(s_{\Omega}/s_{m})^{2}+G^{2}}}\ket{01}+\frac{G}{\sqrt{(s_{\Omega}/s_{m})^{2}+G^{2}}}\ket{11}\Big{)}.$
(28)

We thus obtain our desired two-qubit state and the other part. Denote for
convenience the whole process as $U$, i.e.,

$\displaystyle U\ket{\bf 0}\ket{00}\ket{i}\ket{\bf 0}=\ket{\bf
0}\Big{(}\frac{1}{s_{m}}\ket{00}+\frac{s_{\Omega}}{s_{m}}\ket{01}+\frac{s_{\eta}}{s_{m}}\ket{10}+G\ket{11}\Big{)}\ket{i}\ket{\bf
0}\ket{\bf 0}\equiv\ket{\Phi}.$ (29)

We then make the following crucial observation:

$\displaystyle\bra{\bf 0}\bra{i}U^{\dagger}U_{\phi}\ket{\bf
0}\ket{j}=\braket{\Phi,\phi}=\frac{1}{s_{m}}\Big{(}\bra{i}\big{(}\vec{w_{j}}-|\vec{\delta\Omega^{\prime}_{j}}|\cdot\ket{\vec{\delta\Omega_{j}^{\prime}}}-|\vec{d\eta^{\prime}_{j}}|\cdot\ket{\vec{d\eta_{j}^{\prime}}}\big{)}\Big{)},$
(30)

which allows us to compute the components in the above combination of vectors,

$\vec{w_{j}}-|\vec{\delta\Omega^{\prime}_{j}}|\cdot\ket{\vec{\delta\Omega_{j}^{\prime}}}-|\vec{d\eta^{\prime}_{j}}|\cdot\ket{\vec{d\eta_{j}^{\prime}}},$

which is exactly the harmonic form $h_{j}$ associated with the initial 1-form
$w_{j}$. Therefore, the inner product $\braket{\Phi,\phi}$ is exactly the
$i$-th entry of the $j$-th harmonic form $h_{j}$. By definition, the unitary
$U^{\dagger}U_{\phi}$ is exactly the unitary block encoding of $\mathbb{H}$
(where $\mathbb{H}$ was defined in algorithm 1), scaled down by a factor
$s_{m}$. We note that the scaling-down factor would not affect the subsequent
algorithm, as the kernel of a scaled-down matrix is the same. In fact, in the
Appendix, we will show that $s_{m}$ is actually very small, of order unity.
Therefore, we omit the factor $s_{m}$ in subsequent discussion.

Given the unitary encoding of $\mathbb{H}$, it is trivial to obtain the
unitary encoding of $\mathbb{H}^{\dagger}$, as it is just a matter of
transposition (note that $\mathbb{H}$ is real). Therefore, it is
straightforward to apply Lemma 2 to obtain the unitary block encoding of
$\mathbb{H}^{\dagger}\mathbb{H}$. Since all the form $w_{j}$’s (and hence,
their harmonic forms) are real, it is safe to use $\mathbb{H}^{T}\mathbb{H}$
simply. We will use $\mathbb{H}^{T}\mathbb{H}$ instead of $\mathbb{H}$ because
$\mathbb{H}$ is not generally symmetric, while $\mathbb{H}^{T}\mathbb{H}$ is
apparently symmetric. Lemma 1 then allows us to efficiently simulate
$\exp(-i\mathbb{H}^{T}\mathbb{H}t)$. The next goal is to estimate the
dimension of $Ker(\mathbb{H}^{T}\mathbb{H})$, where $Ker$ refers to the kernel
space. This is similar to what was done in lloyd2016quantum , where the
authors also employed the simulation of the boundary map combined with quantum
phase estimation method kitaev1995quantum  to extract the kernel of the
corresponding boundary operator, which then reveals the Betti numbers.
Basically, we prepare the following mixed state:

$\rho=\frac{1}{c_{1}}\sum_{j=1}^{c_{1}}\ket{j}\bra{j},$

by first preparing $\sum_{j=1}^{c_{1}}\ket{j}/\sqrt{c_{1}}$ and using CNOT
gates to copy the bitstring to another ancilla system initialized in $\ket{\bf
0}$. Tracing out either system yields $\rho$. We then run the quantum phase
estimation algorithm with $\exp(-i\mathbb{H}^{T}\mathbb{H}t)$ and the input
state $\rho$. Note that kernel space corresponds to the eigenvectors with zero
eigenvalues; therefore, for those eigenvectors, the quantum phase estimation
algorithm will output the zero-bit string on the extra register that holds the
outcome of phase estimation. Since the mixture $\rho$ is uniform, the
probability of measuring zeros on the extra register in the quantum phase
estimation circuit is:

$\displaystyle p_{0}=\frac{{\rm dim}Ker(\mathbb{H}^{T}\mathbb{H})}{c_{1}}.$
(31)

Therefore, we can estimate dim $Ker(\mathbb{H}^{T}\mathbb{H})$ by repeating
the algorithm to extract the measurement outcomes. In order to estimate
$p_{0}$ to additive accuracy $\epsilon$, we need to repeat the process
$\mathcal{O}(1/\epsilon^{2})$ times. However, ${\rm
dim}Ker(\mathbb{H}^{T}\mathbb{H})$ is an integer, and we require the ability
to estimate it to some multiplicative error $\delta_{d}$. One can choose

$\epsilon=\delta_{d}\frac{{\rm dim}Ker(\mathbb{H}^{T}\mathbb{H})}{c_{1}}$

to achieve such a goal. As the last step in our algorithm, the following
result shows that dim $Ker(\mathbb{H}^{T}\mathbb{H})$ is exactly dim
$Ker(\mathbb{H})$.

###### Lemma 7

Let $X=A^{T}A$ be some matrix of size $d\times d$ (A is not necessary to be
square). Then Ker($X$) = Ker($A$).

Proof: ($\rightarrow$) It is obvious that Ker($A$) $\subset$ Ker($X$), since
if $Ax=0$ (i.e., $x\in Ker(A)$) then $A^{T}Ax=0$.
($\leftarrow$) Let $x\in Ker(X)$, hence $A^{T}Ax=0$. We consider the inner
product

$x^{T}A^{T}Ax=x^{T}(A^{T}Ax)=0.$

But $x^{T}A^{T}=(Ax)^{T}$, therefore, the above equality is equivalent to

$(Ax)^{T}Ax=0.$

The above product is also the definition of a norm of a real vector $Ax$,
which is, by definition, equal to 0 if and only if $Ax=0$. It implies that
$x\in Ker(A)$. Our proof is thus completed.

The result sheds some light on how the Betti numbers can be calculated. We
remark that dim $Ker(\mathbb{H})$ is not the Betti number itself but the rank
of $\mathbb{H}$ which is $c_{1}-{\rm dim}Ker(\mathbb{H})$ is equal to
$\beta_{1}$ – the first Betti number. Therefore, the ability to find ${\rm
dim}Ker(\mathbb{H}^{T}\mathbb{H})$ (which is $c_{1}-\beta_{1})$ via the
quantum phase estimation method outlined earlier yields the first Betti
number.

As a brief summary, our algorithm mainly employs the matrix inversion and
matrix multiplication technique to transform the initial 1-form
$\ket{w^{\prime}}$ to its harmonic representation. Quantum signal processing
is then employed, combined with the quantum phase estimation plus repeating
measurements to extract the Betti number. The main computational expenses come
from matrix multiplication and inversion, as well as the number of
measurements required to estimate Betti number to high accuracy. For the
matrix multiplication/ inversion step, the matrices that appeared to have a
dimension $c\times c$; therefore, these steps take $\mathcal{O}(\log(c))$
times (we ignore the sparsity dependent because it is very small, as we will
see later). We emphasize a quite important point that, per Lemma 4, its
running time also involves the conditional number, which can be large in
general. This issue has been resolved in clader2013preconditioned , where the
authors provided a generalized quantum inversion method that allows a high
conditional number matrix to be efficiently inverted. While the method in
clader2013preconditioned  induces a larger scaling on the sparsity, we note
that the sparsity in our case is very small, as we will explicitly show in
Appendix A. Therefore, we can safely integrate the result
clader2013preconditioned  into our work without substantial further scaling.
Further recall that $c=c_{1}+c_{2}\leq 2\max(c_{0},c_{1},c_{2})$, where
$c_{0,1,2}$’s are the number of 0-simplices, 1-simplices and 2-simplices in
the given complex, respectively. It is straightforward to see that
$\mathcal{O}(\log(c))=\mathcal{O}(\log(n))$ where $n$ is the number of points.
We state our main result in the following theorem.

######  Theorem 4 (Estimating 1st Betti number)

Let $\sum$ be a simplicial complex with $n$ points, corresponding to the
triangulation of a 2-manifold. The 1st Betti number $\beta_{1}$ can be
estimated to multiplicative accuracy $\delta_{d}$ in time

$\displaystyle\mathcal{O}\Big{(}\frac{\log(n)}{\delta_{d}^{2}}\cdot\frac{c^{2}_{r}}{(c_{r}-\beta_{r})^{2}}\Big{)}.$

##  V Generalization to Different Betti Numbers and Higher Dimensional
Manifolds

Given a triangulated 2-manifold, aside from the first, there are also the
zeroth and the second Betti numbers $\beta_{0},\beta_{2}$. The zeroth Betti
number is always one here, as we already assume that the given manifold is a
connected graph that is triangulated. We have presented a quantum algorithm
for the first Betti number. In order to compute $\beta_{2}$, we first
randomize a different 2-form $\Omega$. We then deform them to the harmonic
form by using the Hodge decomposition (2). The difference is that for a
2-manifold, there is no higher than 2-simplex; therefore, there will be no co-
exact term. The decomposition is thus

$\displaystyle\Omega=d\omega+h,$ (32)

where $\omega$ is a 1-form, and $h$ is the corresponding harmonic 2-form of
$\Omega$. We then follow the same procedure as we did for calculating
$\beta_{1}$.

The next important and somewhat subtle point is about higher dimensional
cases. As we have emphasized from the beginning, the underlying groundwork of
our method is the discrete Hodge theory. In order to apply this, for example,
using Eqn. 2, we need to be able to formulate the operator $d\delta$ for all
forms (see Eqn. 7). This formula depends on the dimension of the given
manifold. The main reason is that, while the operator $d$ is easily formulated
in terms of a matrix given the simplex, the formulation of the operation
$\delta$ relies on the Poincaré duality, which uses the concept of dual cell.

Figure 2: Dual cell complexes in 2 and 3 dimensions. Left figure: In a
triangulated 2-manifold, the dual to 2-simplex $[v_{i},v_{j},v_{k}]$ is a
point o. The dual to an edge $[v_{i},v_{j}]$ is an edge $o\bar{o}$. Right
figure: In a triangulated 3-manifold, the dual to 3-simplex
$[v_{i},v_{j},v_{k},v_{m}]$ is a point o. The dual to a face (2-simplex)
$[v_{i},v_{j},v_{k}]$ is an edge $o\bar{o}$.

For a triangulated 2-manifold, a dual cell to an edge (1-simplex) is an edge.
However, for a triangulated 3-manifold, a dual cell to an edge is no longer an
edge but a face (2-simplex). It means that even for the same first Betti
number $\beta_{1}$, the matrix coefficients (see Eqn. 7) are different.
Fortunately, the computation procedure is the same, as it still relies on
deforming a given $k$-form to its harmonic part. We can proceed with the
computation of arbitrary Betti numbers and in arbitrary dimensional manifold
in exactly the same way as we did in the above algorithm, provided the matrix
form of $\delta$ (and hence $d\delta$) is given. Despite that there is no
exact expression for $\delta$, its matrix form can be efficiently calculated
given the triangulation of the manifold (see Appendix B). Therefore, given a
triangulated $m$-dimensional manifold $\sum$ having $n$ points, its $r$-th
Betti number $\beta_{r}$ can be estimated (with a multiplicative error
$\delta_{e}$) with time complexity

$\displaystyle\mathcal{O}\Big{(}\frac{\log(c_{r})}{\delta_{e}^{2}}\cdot\frac{c_{r}^{2}}{(c_{r}-\beta_{r})^{2}}\Big{)},$

where $c_{r}$ is the number of $r$-simplices in $\sum$. In the low $r$ limit
($r\ll n$), then $\log(c_{k})\in\mathcal{O}(\log(n))$; whereas for the limit
$r\sim n$, then $\log(c_{r})\sim r\log(n)$. Our algorithm achieves its best
performance when $\beta_{r}\ll c_{r}$, which is opposite to that of the LGZ
lloyd2016quantum ; schmidhuber2022complexity .

##  VI Analysis and Discussion

To make a comparison, we remind some of the results regarding the calculation
of the first Betti number. As mentioned in schmidhuber2022complexity , given a
specified simplicial complex $\sum$ with $n$ points, the best classical
algorithm to compute the first Betti number $\beta_{1}$ takes time
$\mathcal{O}(c_{1}^{3})$, where $c_{1}$ is the number of 1-simplices in
$\sum$. As also mentioned in schmidhuber2022complexity , the best running time
of the improved-LGZ algorithm to estimate Betti number to multiplicative error
$\delta$ is

$\displaystyle\mathcal{O}\Big{(}\frac{1}{\delta}n\kappa\sqrt{\frac{c_{1}}{\beta_{1}}}\Big{)}.$
(33)

The advantage of the improved-LGZ algorithm compared to the classical
algorithm is obtained in the regime $\beta_{1}\rightarrow c_{1}$ (and
generally, $\beta_{k}\rightarrow c_{k}$), as they emphasize in their work
schmidhuber2022complexity . More specifically, in such a regime for a
triangulated 2-manifold as in our specific case, the classical running time is
$\mathcal{O}(n^{3})$, since for a triangulation of $n$ points, the number of
edges (1-simplices) is bounded by $\mathcal{O}(n)$. The running time of the
improved LGZ is (we ignore the $\kappa$ factor) $\mathcal{O}(n/\delta)$.
Hence, their quantum algorithm yields a cubic speedup.

In the same regime, our algorithm turns out to be slower than both above
approaches, as $c_{1}^{2}$ will have the order $\sim\mathcal{O}(n^{4})$. On
the other hand, in the regime where $\beta_{1}$ is small compared to $c_{1}$,
our algorithm will have a running time

$\mathcal{O}\Big{(}\frac{\log(n)}{\delta^{2}}\Big{)},$

which is an exponential speedup compared to both the classical algorithm and
that of the quantum algorithms (LGZ and its improved version) lloyd2016quantum
; schmidhuber2022complexity . Such an interesting opposite performance is a
somewhat unexpected outcome from the duality between cohomology and homology.

Another subtle point that we would like to discuss is the assumption that we
made regarding the shape of the simplicial complex. In our work, we only deal
with a uniformly triangulated manifold, i.e., all the composing simplices are
similar in size, and the manifold is built by properly gluing them together.
In some previous works, such as lloyd2016quantum ; mcardle2022streamlined ;
schmidhuber2022complexity ; ubaru2021quantum , the setting is seemingly more
general where the shape could be arbitrary, as two points only get connected
if their distance is smaller than a known threshold. One may wonder if our
assumption would severely limit the practicality of the outlined algorithm.
The answer is no, as the topology of the underlying manifold only depends on
the connectivity but not the actual distance between arbitrary two points. As
an example, let us consider the 2-dimensional case. One can imagine that,
given a triangle with arbitrary angles/lengths, it is deformable or
topologically equivalent to an equilateral triangle. Therefore, the
topological space formed by the union, or by gluing different triangles
together, is topologically equivalent to the union of equilateral triangles,
which form the uniformly triangulated manifold. Since they are topologically
equivalent, their underlying topological properties are the same.
Consequently, performing computation on the uniformly triangulated manifold is
more convenient, as we discuss further in Appendix B, where we provide an
explicit formula for the codifferential operator. Finally, we remark that in
general, Hodge theory (see Theorem 2) works with closed manifold (in both
smooth and discrete setting). If our initial configuration is an open
triangulated manifold, e.g, with a boundary, then a simple trick is to double
cover the simplicial complex while preserving the symmetry. In this way, we
can recover the closeness condition and proceed with our algorithm.

##  VII Conclusion

Our work has provided a ‘dual’ approach to lloyd2016quantum  and is built upon
the (discrete) Hodge theory and de Rham cohomology. There are a few major
reasons underlying the advantage of our approach compared to lloyd2016quantum
. In lloyd2016quantum , the authors basically quantized the homology approach,
associating the chain group to a vector space and finding singular values/
vectors of the boundary map $\partial$. The key step in lloyd2016quantum  is
the identification of the chain group as a vector space, whereby each simplex
is represented by a basis state, which means that the vector space needs to be
at least as large as the number of simplices contained in $\sum$. By doing
this way, the resource scales as a polynomial of $n$. In reality, the high
value of $n$ is usually desired (large-scale analysis), which means that
poly($n$) could be very high, and hence the algorithm induces a very high
computational cost. The cohomology approach we have adopted here fits nicely
in such a large-scale setting. We recall that in (simplicial) cohomology, a
cochain is a map $C\rightarrow R$ where $C$ is some chain group (space).
Equivalently, we can imagine that each simplex is associated with a real
number, and, therefore, we can have a more efficient way of storing our data,
as we only need $\sim\log(M)$ qubits to store an $M$-dimensional vector, which
further reduces the resource needed for processing. Another major point is
that in lloyd2016quantum , it is required to generate the proper simplicial
complex state, which contributes substantially to the computational cost and,
in certain cases, has a high failing probability (for details, see
lloyd2016quantum ). Here, in the cohomology approach, the initial 1-form state
can, in fact, be chosen arbitrarily, which is more convenient. We further make
an important remark that usually, the first and second Betti numbers are
important enough for us to infer the underlying structure of a dataset, which
implies that for those low Betti numbers, our algorithm actually has running
time $\sim\mathcal{O}(\log(n))$. It is worth noting that our dual approach
does not perform well in the regime of high Betti numbers, for which the LGZ
algorithm works very well. Whether there is a mixed approach that can give
rise to better performance in the intermediate regime is left for future
exploration.

###### Acknowledgements.

This work was supported in part by the US Department of Energy, Office of
Science, National Quantum Information Science Research Centers, Co-design
Center for Quantum Advantage (C2QA) under contract number DE-SC0012704
(T.-C.W.), and by the National Science Foundation under Grants No.
CMMI-1762287 (X.D.G.) and No. FAIN-2115095 (X.D.G.), as well as NIH
3R01LM012434-05S1 (X.D.G.) and NIH 1R21EB029733-01A1 (X.D.G.). We also
acknowledge the support from a Seed Grant from Stony Brook University’s Office
of the Vice President for Research.

## References

  * (1) Larry Wasserman.  Topological data analysis.  arXiv preprint arXiv:1609.08227, 2016.
  * (2) Peter Bubenik et al.  Statistical topological data analysis using persistence landscapes.  J. Mach. Learn. Res., 16(1):77–102, 2015.
  * (3) Xianfeng Gu, Yalin Wang, Tony F Chan, Paul M Thompson, and Shing-Tung Yau.  Genus zero surface conformal mapping and its application to brain surface mapping.  IEEE Transaction on Medical Imaging (TMI), 23(8):949–958, 2004\.
  * (4) Xianfeng David Gu and Shing-Tung Yau.  Computational conformal geometry, volume 1.  International Press Somerville, MA, 2008.
  * (5) Xianfeng Gu, Feng Luo, and Shing Tung Yau.  Computational conformal geometry behind modern technologies.  Notices of the American Mathematical Society, 67(10):1509–1525, 2020.
  * (6) Seth Lloyd, Silvano Garnerone, and Paolo Zanardi.  Quantum algorithms for topological and geometric analysis of data.  Nature communications, 7(1):1–7, 2016.
  * (7) A Yu Kitaev.  Quantum measurements and the abelian stabilizer problem.  arXiv preprint quant-ph/9511026, 1995.
  * (8) Dominic W Berry, Graeme Ahokas, Richard Cleve, and Barry C Sanders.  Efficient quantum algorithms for simulating sparse hamiltonians.  Communications in Mathematical Physics, 270(2):359–371, 2007.
  * (9) Shashanka Ubaru, Ismail Yunus Akhalwaya, Mark S Squillante, Kenneth L Clarkson, and Lior Horesh.  Quantum topological data analysis with linear depth and exponential speedup.  arXiv preprint arXiv:2108.02811, 2021.
  * (10) Sam McArdle, András Gilyén, and Mario Berta.  A streamlined quantum algorithm for topological data analysis with exponentially fewer qubits.  arXiv preprint arXiv:2209.12887, 2022.
  * (11) Sam Gunn and Niels Kornerup.  Review of a quantum algorithm for betti numbers.  arXiv preprint arXiv:1906.07673, 2019.
  * (12) Bernardo Ameneyro, Vasileios Maroulas, and George Siopsis.  Quantum persistent homology.  arXiv preprint arXiv:2202.12965, 2022.
  * (13) Alexander Schmidhuber and Seth Lloyd.  Complexity-theoretic limitations on quantum algorithms for topological data analysis.  arXiv preprint arXiv:2209.14286, 2022.
  * (14) Aram W Harrow, Avinatan Hassidim, and Seth Lloyd.  Quantum algorithm for linear systems of equations.  Physical review letters, 103(15):150502, 2009.
  * (15) Seth Lloyd, Masoud Mohseni, and Patrick Rebentrost.  Quantum principal component analysis.  Nature Physics, 10(9):631–633, 2014.
  * (16) Andrew M Childs, Robin Kothari, and Rolando D Somma.  Quantum algorithm for systems of linear equations with exponentially improved dependence on precision.  SIAM Journal on Computing, 46(6):1920–1950, 2017.
  * (17) Nhat A Nghiem and Tzu-Chieh Wei.  An improved method for quantum matrix multiplication.  Quantum Information Processing, 22(8):299, 2023.
  * (18) Guang Hao Low and Isaac L Chuang.  Optimal hamiltonian simulation by quantum signal processing.  Physical review letters, 118(1):010501, 2017.
  * (19) Guang Hao Low and Isaac L Chuang.  Hamiltonian simulation by qubitization.  Quantum, 3:163, 2019.
  * (20) András Gilyén, Yuan Su, Guang Hao Low, and Nathan Wiebe.  Quantum singular value transformation and beyond: exponential improvements for quantum matrix arithmetics.  In Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing, pages 193–204, 2019.
  * (21) Raoul Bott, Loring W Tu, et al.  Differential forms in algebraic topology, volume 82.  Springer, 1982.
  * (22) Nathan Wiebe, Daniel Braun, and Seth Lloyd.  Quantum algorithm for data fitting.  Physical review letters, 109(5):050505, 2012.
  * (23) B David Clader, Bryan C Jacobs, and Chad R Sprouse.  Preconditioned quantum linear system algorithm.  Physical review letters, 110(25):250504, 2013.
  * (24) Mikio Nakahara.  Geometry, topology and physics.  CRC press, 2018.

##  Appendix A Further elaboration on linear equation

Some necessary background on topology and geometry can be found in any
standard textbook, such as [24]. In this section, we discuss in detail how to
solve for $\Omega$ and $\eta$ that appear in the Hodge decomposition (3),
which constitutes the linear equation similar to Eq. (7) that we will solve
with a quantum means. The Hodge decomposition holds for arbitrary dimensions,
but since we deal with a triangulated 2-manifold and present how to compute
the first Betti number $\beta_{1}$, we proceed in this section with the case
of 1-forms. Further generalization to higher dimensions will be left for the
next section, which allows for the adoption of our quantum procedure.

Let $\omega$ be some 1-form. Recall that the Hodge decomposition allows
$\omega$ to be written as

$\displaystyle\omega=\delta\Omega+d\eta+h,$ (34)

where $\Omega$ is a 2-form, $\eta$ is a 0-form and $h$ is called the harmonic
form, $d$ is the exterior derivative, which maps from the space of 0-form to
1-form, and $\delta$ is the codifferential operator, which maps a 2-form to a
1-form. We note that $h$ vanishes under $d$ and $\delta$: $dh=0$ and $\delta
h=0$. These two operators have a special property:

$d^{2}\equiv d\circ d=0,$

and

$\delta^{2}\equiv\delta\circ\delta=0.$

Applying the derivative map $d$ to both sides of Eqn. 34, and using the above
properties of both the derivative map and the harmonic form, we obtain:

$\displaystyle d\omega=d\delta\Omega.$ (35)

We have mentioned that in the discrete context, an arbitrary $k$-form is
approximated by a $k$-cochain. Let $[v_{1},v_{2},v_{3}]$ be a 2-simplex, then
we have:

$\displaystyle d\omega[v_{1},v_{2},v_{3}]=d\delta\Omega[v_{1},v_{2},v_{3}].$
(36)

Figure 3:  A part of triangulated 2-manifold. Each face (2-simplex) is labeled
with red color.

From the above equation, we have

$d\omega[v_{1},v_{2},v_{3}]=\omega[v_{1},v_{2}]+\omega[v_{2},v_{3}]+\omega[v_{3},v_{1}].$

If $\omega$ is known, then we can compute $d\omega[v_{1},v_{2},v_{3}]$.
Likewise:

$d\delta\Omega[v_{1},v_{2},v_{3}]=\delta\Omega[v_{1},v_{2}]+\delta\Omega[v_{2},v_{3}]+\delta\Omega[v_{3},v_{1}].$

A standard procedure using the Poincaré duality to formulate $\delta$ is done
in the standard textbook, such as [4]. Here we quote the results. From Fig. 3,
$f_{1}$ and $f_{2}$ are two faces (2-simplices) adjacent to $[v_{1},v_{2}]$.
We then have the following formula for $\delta\Omega[v_{1},v_{2}]$:

$\displaystyle\delta\Omega[v_{1},v_{2}]=\frac{1}{w_{12}}\Big{(}\frac{\Omega(f_{1})}{|f_{1}|}-\frac{\Omega(f_{2})}{|f_{2}|}\Big{)},$
(37)

where $|f_{1,2}|$ is the area of the corresponding faces, or 2-simplices, and
$w_{12}$ is the cotangent edge weight (see in Fig. 3 for the illustration),

$w_{12}=\frac{1}{2}\big{(}\cot(\angle v_{1}v_{3}v_{2})+\cot(\angle
v_{1}v_{4}v_{2})\big{)}.$

The above formula for the cotangent edge weight needs to be slightly modified
if the edge $[v_{1},v_{2}]$ lies on the boundary (for example, edge
$[v_{1},v_{4}]$ on Fig. 3), which means that there is only one face adjacent
to it. In such a case, then we only include the cotangent of one angle, which
is opposite to $[v_{1},v_{4}]$ ($\angle v_{1}v_{2}v_{4}$).

Applying formula 37, we then obtain a full formula for
$d\delta\Omega[v_{1},v_{2},v_{3}]$, which is:

$\displaystyle
d\delta\Omega[v_{1},v_{2},v_{3}]=\frac{1}{w_{12}}\Big{(}\frac{\Omega(f_{1})}{|f_{1}|}-\frac{\Omega(f_{2})}{|f_{2}|}\Big{)}+\frac{1}{w_{23}}\Big{(}\frac{\Omega(f_{3})}{|f_{3}|}-\frac{\Omega(f_{2})}{|f_{2}|}\Big{)}+\frac{1}{w_{13}}\Big{(}\frac{\Omega(f_{4})}{|f_{4}|}-\frac{\Omega(f_{2})}{|f_{2}|}\Big{)}.$
(38)

We recall that $d\omega[v_{1},v_{2},v_{3}]=d\delta\Omega[v_{1},v_{2},v_{3}]$,
and as $\omega$ is known, we obtain an algebraic equation for $\Omega$. We
note that the Eqn. 38 contains the first order of $\Omega$ (acting on certain
2-simplices), and, therefore, repeating equation 36 eventually yields a linear
equation (Eqn. 7), and the solution is unique [4]. We note that, the
coefficients of matrix $A$ in Eqn. 7 in general will include terms of the
following form

$\frac{1}{w\cdot|f|},$

where $w$ is a cotangent edge weight, and $|f|$ is the area of some 2-simplex.
In general, these quantities can be computed provided the pairwise distance
between vertices that form the corresponding simplices. In our problem, we are
working with a uniform triangulated manifold, which means that all the
2-simplices or triangles are equilateral. Therefore, the cotangent edge weight
is simply:

$w=\frac{1}{2}(\frac{1}{\sqrt{3}}+\frac{1}{\sqrt{3}})=\frac{1}{\sqrt{3}}$

for all cotangent edge weights (including the internal ones) and the value of
$|f|$ is

$|f|=\frac{\sqrt{3}a^{2}}{2},$

where $a$ is the distance between arbitrary two connected points, or
0-simplices. The value $a$ depends on the given dataset and, in principle, it
can be arbitrary, as the length does not affect the topology of a given
configuration.

Now, we elaborate the computation of 0-form $\eta$ in Eqn. 34. By applying
$\delta$ to both sides of Eqn. 34, we have:

$\delta\omega=\delta d\eta.$

Figure 4: Illustration of dual cell of 0-simplex $v_{0}$. The whole region
inside the red-colored boundaries is the dual cell of $v_{0}$, with area
denoted as $\tilde{|v_{0}|}$. This area could be computed given the pairwise
distance between $v_{0}$ and its adjacent vertices. In our case where the
given triangulation is uniform and the distance between two vertices is $a$,
it is easy to calculate the area $|\tilde{v_{0}}|$, which is
$pa^{2}/4\sqrt{3}$, where $p$ is the number of adjacent points to $v_{0}$.

As $\delta\omega$ is a 0-form, using the same method of the Poincaré duality
[4], the formula for $\delta\omega$ acting on some 0-simplex or a point
$v_{0}$ is:

$\displaystyle\delta\omega[v_{0}]=-\frac{1}{\tilde{|v_{0}|}}\sum_{v_{i}\sim
v_{0}}w_{0j}\cdot\omega[v_{0},v_{j}],$ (39)

where the symbol $\sim$ denotes the adjacency (i.e., vertex $v_{j}$ is
connected to $v_{0}$ via an edge), $w_{0,j}$ is the cotangent edge weight
(defined above), and $|\tilde{v_{0}}|$ is the area of the dual cell complex of
$v_{0}$ (we will elaborate further the dual cell complex in the subsequent
section).

The formulation of $\delta d\eta$ is the same, as $d\eta$ is also a 1-form.
More explicitly, we have:

$\delta d\eta[v_{0}]=-\frac{1}{\tilde{|v_{0}|}}\sum_{v_{i}\sim
v_{0}}w_{0j}\cdot d\eta[v_{0},v_{j}]=-\frac{1}{\tilde{|v_{0}|}}\sum_{v_{i}\sim
v_{0}}w_{0j}\cdot(\eta[v_{j}]-\eta[v_{0}]).$

Since $\delta\omega=\delta d\eta$, the area term $\tilde{|v_{0}|}$ cancels
out, fortunately. By applying the above formulas for all vertices, eventually,
we obtain a linear equation for $\eta$ with known coefficients $w$’s, i.e.,
from the cotangent edge weights.

##  Appendix B General Version of Discrete Hodge Theory

As we have emphasized, going to higher dimensional manifolds induces more
complication, as the rule for $d\delta$ and $\delta d$ that we described in
the previous section will no longer hold exactly the same form. The main
reason is due to the use of the Poincaré duality. For a different dimension,
the dual cell to each simplex is generally different. Hence, it affects the
way that we formulate the codifferential operator $\delta$ (and its
composition with $d$). We remark that the underlying theory of our algorithm
relies on results from non-trivial geometry and topology, which means that we
do not expect to cover everything here, and we refer the readers to the
literature, such as [24], and particularly, more relevant to our perspective
from computational algorithms, Ref. [4]. In this section, we aim to exemplify
the subtle points regarding the generalization of our algorithm to higher
dimensions in more details than what we described in Sec. V.

###  B.1 Poincare Dual Cell Complex

We first review concepts of Poincare dual cell complex. The key point is that,
given a triangulated manifold of some dimension $n$, the dual cell of some
$k$-simplex is not necessarily a proper simplex.

Figure 5: Dual cell complex in 2-d case

As we can see from Fig. 5, the order of the given simplex (0,1,or 2-simplex)
and of its dual cell match the dimension of given manifold, which is 2. In the
2d case, the dual cell to a 2-simplex, denoted $f$ (rightmost ones) is a point
$\bar{f}$ (of zero order), which is the center of the circle that passes all
three points of the original 2-simplex $f$. In the middle two figures, the
1-simplex $e$ has the dual cell $\bar{e}$. The way $\bar{e}$ is constructed is
that, first, we construct the dual cell of the two triangles, or 2-simplices,
that are adjacent to $e$, which are two points. Then we connect them, hence
forming a line $\bar{e}$. Likewise for the leftmost figures, the dual cell of
a given point or 0-simplex $v$ is formed by connecting the dual cell of all
triangles (2-simplex) that $v$ belongs to, and including the interior. Hence,
according the bottom left figure, the dual cell of $v$, $\bar{v}$ is the
region inside the red lines, including the boundary. Similar construction
holds for a triangulated manifold of arbitrary dimension.

###  B.2 Codifferential Operator

Let $M$ be a triangulated $n$-dimensional manifold. Now we examine the
codifferential operator $\delta:\Omega^{k}\rightarrow\Omega^{k-1}$, where
$\Omega^{k}$ is the space of $k$-form $(k<n)$. Let $\sigma$ be some $(k-1)$
simplex and $\bar{\sigma}$ be its corresponding Poincaré dual-cell complex.
Additionally, let $\partial\bar{\sigma}=\sum_{i}\gamma_{i}$ be the boundary of
$\bar{\sigma}$, where in general, $\gamma_{i}$ is some simplex (the
orientation is included). Let $\Delta$ be some $k$-form. We recall that, as we
have mentioned in the main text, in discrete setting, a form is identical as a
cochain. We have thus the following general formula:

$\displaystyle\delta\Delta(\sigma)=(-1)^{nk+n+1}(-1)^{(n-k+1)(k-1)}\frac{|\sigma|}{\bar{|\sigma|}}\sum_{i}\frac{|\gamma_{i}|}{|\bar{\gamma_{i}}|}\Delta(\bar{\gamma_{i}}).$
(40)

The above formula is somewhat complicated, but it is straightforward, since in
the previous section, we have provided a systematic way to construct a dual
cell complex. Given an $n$-dim triangulated manifold, we can, in principle,
construct the corresponding dual cell complex and can then carry out the
formula for $\delta$ algorithmically. Given distances between points on the
original manifold $M$, in principle, we can compute corresponding length,
area, volume elements of simplices of $M$ and their dual cell complex. For the
uniform triangulation, the calculation is much easier, as we have computed
some quantities like $w$ and $|f|$ in the previous section. Therefore, all the
linear systems (such as in Eqn. (7)) are entry-wise accessible.

We will exemplify the formula 40 further for both 2-dim and 3-dim cases by
direct calculation. In the previous section A, we already encounter the
formulation. We will try to reproduce the same formulas for $\delta\Omega$ and
$\delta w$ in 2-d case.

2-dimensional Triangulated Manifold

Figure 6: Dual cell complex of a vertex and an edge.

Let $\Omega$ be 2-form (or 2-cochain) and $\omega$ be 1-form (or 1-cochain).
First we try to apply (40) to compute $\delta\Omega(e)$ on the right of Fig.
6. According to the figure, the dual cell of the triangle
$f_{1}=[v_{0},v_{1},v_{2}]$ is denoted as $\tilde{o}$, and the dual cell of
triangles $f_{2}=[v_{0},v_{1},v_{3}]$ is $o$. The dual cell of $e$, or
1-simplex $[v_{0},v_{1}]$, is $\bar{e}=o\bar{o}$. Applying Eqn. 40, we have

$\delta\Omega(e)=\frac{|e|}{|\bar{e}|}\Big{(}\frac{|o|}{|\bar{o}|}\Omega(\bar{o})-\frac{|\tilde{o}|}{|\bar{\tilde{o}}|}\Omega(\bar{\tilde{o}})\Big{)}.$

Since both $o$ and $\tilde{o}$ are points, $|o|$ and $|\tilde{o}|$ are equal
to 1. Furthermore, the dual of $o$ and $\tilde{o}$ are
$[v_{0},v_{1},v_{3}],[v_{0},v_{1},v_{2}]\equiv f_{1},f_{2}$, respectively, as
we defined in the previous section. A result from elementary geometry [4]
shows that

$\frac{|\bar{e}|}{|e|}=\frac{1}{2}\big{(}\cot(\angle
v_{1}v_{3}v_{0})+\cot(\angle v_{1}v_{2}v_{0})\big{)}=w_{12},$

which is defined as the cotangent edge weight. Combining everything, we
obtain:

$\displaystyle\delta\Omega(e)=\frac{1}{w_{12}}\Big{(}\frac{1}{f_{1}}\Omega(f_{1})-\frac{1}{f_{2}}\Omega(f_{2})\Big{)},$
(41)

which is exactly what we got in Eqn. 37. Similar calculation can be done for
$\delta\omega$, which leads to a similar form as in Eqn. 39.

3-dimensional Triangulated Manifold
In the 3d case, the situation is a bit more complicated because the dual of
each simplex is no longer the same as in the 2d case. However, the formulation
is still straightforward, as we can systematically construct the dual cell and
use the same formula 40. We work out the following example; see figure 7.

Figure 7: Example in 3-D

We denote the 2-simplex $\sigma=[v_{i},v_{j},v_{k}]$ and its dual cell
$\bar{\sigma}=o\bar{o}$; we further denote 3-simplex
$[v_{n},v_{i},v_{j},v_{k}],[v_{m},v_{i},v_{j},v_{k}]$ as $\alpha,\beta$ . Let
$\Omega$ be some 3-form. The boundary of $\bar{\sigma}$ is:

$\partial\bar{\sigma}=\bar{o}-o$

The action of $\delta\Omega(\sigma)$ is:

$\displaystyle\delta\Omega(\sigma)=\frac{|\sigma|}{|o\bar{o}|}\Big{(}\frac{1}{|\alpha|}\Omega(\alpha)-\frac{1}{|\beta|}\Omega(\beta)\Big{)}.$
(42)

Since the 3-simplex is a tetrahedron, its volume is
$V=|\alpha|=|\beta|=a^{3}/6\sqrt{2}$, and the length of $o\bar{o}$ is
$a/\sqrt{6}$; the area of $\sigma$ is $|\sigma|=\sqrt{3}a^{2}/4$.

##  Appendix C Product of Block Encoding

Here we explicitly show that, given a unitary block encoding for matrices
$A_{1}$ and $A_{2}$, respectively, one can construct the unitary block
encoding of $A_{1}A_{2}$, or $A_{2}A_{1}$ with a few modest steps. While it
has been established in the original work [20], we present a solution here for
completeness.

Let $U_{1},U_{2}$ be block encoding of $A_{1},A_{2}$, i.e, we have:

$\displaystyle U_{1,2}=\begin{pmatrix}A_{1,2}&\cdot\\\ \cdot&\cdot\\\
\end{pmatrix}.$

Equivalently, we can write $U$ (we drop the subscript for now, as it should
not cause any confusion):

$\displaystyle U=\ket{\bf{0}}\bra{\bf{0}}\otimes A+\cdots$

We now observe the following property:

$\displaystyle U\ket{\bf 0}\ket{\phi}=\ket{\bf
0}A\ket{\phi}+\sum_{j\neq\textbf{0}}\ket{j}\ket{\phi_{j}},$ (43)

where $\ket{\phi}$ and all $\ket{\phi_{j}}$ (the latter being essentially
garbage states) have the same dimension as matrix $A$. For a reason that will
be clear later on, we borrow an extra qubit initialized in $\ket{0}$ and
rewrite the above equation as:

$\displaystyle\mathbb{I}\otimes U\ket{\bf 0}\ket{\phi}=\ket{0}\ket{\bf
0}A\ket{\phi}+\ket{0}\sum_{j\neq\textbf{0}}\ket{j}\ket{\phi_{j}}.$ (44)

Now we use $X$ gate to flip the ancilla qubit to obtain the state:

$\displaystyle\ket{1}\ket{\bf
0}A\ket{\phi}+\ket{1}\sum_{j\neq\textbf{0}}\ket{j}\ket{\phi_{j}}.$

Denote the whole above unitary process, which maps $\ket{0}\ket{\bf
0}\ket{\phi}$ to $\ket{1}\ket{\bf
0}A\ket{\phi}+\ket{1}\sum_{j\neq\textbf{0}}\ket{j}\ket{\phi_{j}}$ as
$\mathcal{U}$. For matrix $A_{1}$, given a computational basis state
$\ket{i}$, we have:

$\displaystyle\mathcal{U}_{1}\ket{0}\ket{\bf 0}\ket{i}=\ket{1}\ket{\bf
0}A_{1}\ket{i}+\ket{1}\sum_{j\neq\textbf{0}}\ket{j}\ket{\phi_{j}}.$ (45)

For matrix $A_{2}$, we would use the conjugate transpose of the original
$U_{2}$ (without the extra CNOT step that flips the ancilla) to obtain:

$\displaystyle U_{2}^{\dagger}\ket{0}\ket{\bf 0}\ket{k}=\ket{0}\ket{\bf
0}A_{2}^{\dagger}\ket{k}+\ket{0}\sum_{j\neq\textbf{0}}\ket{j}\ket{\phi_{j}}.$
(46)

Now for the above state, we use the second register (associated with $\ket{\bf
0}$) as controlling qubits to flip the first qubit $\ket{0}$ to $\ket{1}$,
conditioned on it being $\ket{\bf 0}$. In other words, from the above state,
we obtain:

$\displaystyle\ket{1}\ket{\bf
0}A_{2}^{\dagger}\ket{k}+\ket{0}\sum_{j\neq\textbf{0}}\ket{j}\ket{\phi_{j}}.$
(47)

We denote the process containing of $U_{2}$ and the controlled-X step above as
$\mathcal{U}_{2}$, i.e.,

$\displaystyle\mathcal{U}_{2}\ket{0}\ket{\bf 0}\ket{k}=\ket{1}\ket{\bf
0}A_{2}^{\dagger}\ket{k}+\ket{0}\sum_{j\neq\textbf{0}}\ket{j}\ket{\phi_{j}}.$
(48)

It is straightforward to observe that, thanks to the orthogonality of
computational basis states, the inner product:

$\displaystyle\bra{0}\bra{\bf
0}\bra{k}\mathcal{U}_{2}^{\dagger}\mathcal{U}_{1}\ket{0}\ket{\textbf{0},i}=\bra{k}A_{2}A_{1}\ket{i}\equiv(A_{2}A_{1})_{ki},$
(49)

which is exactly the definition of unitary block encoding. Therefore, we have
successfully constructed the block encoding of $A_{2}A_{1}$. For the reverse
order $A_{1}A_{2}$, the procedure is exactly the same, as we just need to
reverse the role of $\mathcal{U}_{1,2}$ from the above procedure.