# Asymptotic Homology of Brownian Motion on a Riemannian Manifold

Alberto Verjovsky and Ricardo F. Vila-Freyer  Instituto de Matemáticas
Universidad Nacional Autónoma de México
Apdo. Postal 273, Admon. de correos #3, C.P. 62251 Cuernavaca, Morelos,
México<EMAIL_ADDRESS>CIMAT
Apdo. Postal 402, Guanajuato, Gto., C. P. 36000, México<EMAIL_ADDRESS>
(Date: August 30, 2024)

###### Abstract.

We prove, using the celebrated result by Spitzer about winding of planar
Brownian motion, and the existence of harmonic morphisms
$f:M\to{\mathbb{S}}^{1}$ representing cohomology classes in
$\text{H}^{1}(M,{\mathbb{Z}})$, that there is a stochastic process
$H_{t}:{\mathcal{C}}(M)\to{\text{Hom}(\text{H}^{1}(M;{\mathbb{R}}),{\mathbb{R}})}\simeq{\text{H}_{1}(M;{\mathbb{R}})}$
($t\in[0,\infty)$), where ${\mathcal{C}}(M)=\\{\alpha:[0,\infty)\to
M:\alpha\,\,\text{is continuous}\\}$, which has a multivariate Cauchy
distribution i.e., such that for each nontrivial cohomology class
$[\omega]\in{\text{H}^{1}(M;{\mathbb{R}}),{\mathbb{R}})}$, represented by a
closed 1-form $\omega$, in the de Rham cohomology, the process
$A^{\omega}_{t}:{\mathcal{C}}(M)\to{\mathbb{R}}\,$ ($t\in[0,\infty)$) with
$A^{\omega}_{t}(B)=H_{t}(B)([\omega]),\,B\in{\mathcal{C}}(M)$ converges in
distribution, with respect to Wiener measure on ${\mathcal{C}}(M)$, to a
Cauchy’s distribution, with parameter 1. The process describes the
“homological winding” of the Brownian paths in $M$, thus it can be regarded as
a generalization of Spitzer result. The last section discusses the asymptotic
behavior of holonomy along Brownian paths.

###### Key words and phrases:

Brownian motion, asymptotic cycles, asymptotic linking number, asymptotic
homology

###### 2010 Mathematics Subject Classification:

58J65, 60J65, 58J65

###### Contents

  1. 1 Introduction
  2. 2 Brownian motion on a Riemannian manifold (a short review).
  3. 3 The asymptotic behavior of homology associated with Brownian motion
  4. 4 Lifting to Holonomy

##  1\. Introduction

The purpose of this paper is to describe the asymptotic probability laws of
homological cycles defined by random paths for any compact Riemannian
manifold. There is previous work on probability theory in the study of the
winding number of planar Brownian motion, which motivated our point of view,
due to P. Lévy, F. Spitzer, R. Durrett, and Pitman and Yor; among others. (See
respectively, [17], [25], [5], and [23]).

We use the idea of Schwartzman [24], who studies long segments of orbits,
$O(T,x)=\\{f_{t}(x)\,:\,t\in[0,T]\\}$, of a flow $f_{t}:M\to{M}$, which
preserves the measure $\mu$ on a manifold $M$ . After joining the endpoints
$x,\,f_{T}(x)$ of $O(T,x)$ with a “short curve” (for instance a geodesic arc
realizing the distance between the endpoints), one obtains an oriented cycle
$C_{T}(x)$ that represents an element in
$\text{H}_{1}(M;{\mathbb{R}})=\text{H}_{1}(M;{\mathbb{Z}})\otimes_{\mathbb{Z}}{\mathbb{R}}$.
He proves that the limit $\underset{t\to\infty}{\lim}\frac{C_{T}(x)}{T}:=C(x)$
exists for $\mu$-almost all $x$ . He calls the 1 -dimensional cycle $C(x)$ the
asymptotic cycle through $x$. The map
$x\overset{{\mathbb{A}}}{\mapsto}{C(x)}\in\text{H}_{1}(M;{\mathbb{R}})$ is a
measurable map which is invariant under the flow:
${\mathbb{A}}(x)={\mathbb{A}}(f_{t}(x))$ (i.e., $C(x)=C(f_{t}(x)))$. The
integral $A:=\int_{M}{\mathbb{A}}(x)d\mu\in\text{H}_{1}(M;{\mathbb{R}})$ is
called the asymptotic cycle of the flow with respect to $\mu$. Thus we may
think of a measure-preserving flow as a “diffuse cycle”.

If the flow is a differentiable flow on a differentiable manifold we may also
thing of the pair $(\\{{f_{t}}\\}_{t\in{\mathbb{R}}},\mu)$ as a de Rham
current or foliated cycle [4].

If $X$ is the vector field tangent to the flow, then the map
$A_{dR}:\text{H}^{1}(M;{\mathbb{R}})\to{\mathbb{R}}$ given by
$A_{dR}([\omega])=\int_{M}\,\omega(X)d\mu$, where $\omega$ is a closed 1-form,
is a linear form on cohomology which represents the asymptotic cycle if $\mu$
es ergodic.

We aim to define an asymptotic homological probability law associated with
Brownian motion on a Riemannian manifold by analogy with the case of
Schwartzman described in the previous paragraph. The construction is based on
an approximation to Brownian motion by piecewise smooth curves, following Itô
[12]: giving an interpretation in terms of a Stochastic Integral. (Of course,
there are other possible ways to do this. For example, Pinsky [21] constructs
Brownian motion by Isotropic Transport Process and this construction also
gives a nice way to approximate Brownian motion by a piecewise smooth path).
To understand better the interpretation of the asymptotic cycle associated
with Brownian motion in the first homology group with real coefficients, we
review the Eilenberg-Bruschlinsky-de Rham homology theory as in [24].
Fundamental to motivate our construction has been the result by Spitzer about
the winding number of planar Brownian motion mentioned above. In fact, by
Poincaré duality, for de Rham cohomology and with the Eilenberg-Bruschlinsky
interpretation of cohomology and homology theory, as homotopy classes of maps
to the circle and maps to the circle modulo those that can be lifted to the
real, respectively, we can understand the asymptotic cycle as a measure of the
winding of a random path around the circle by composing the path with a
harmonic map to the circle which is not homotopic to a constant. This together
with some results of harmonic morphisms allows us to prove that the asymptotic
cycle associated with Brownian motion is a random function
$H:{\mathcal{C}}(M)\to\text{H}_{1}(M;{\mathbb{R}})$ such that
$\varphi\circ{H}$ has a Cauchy distribution with respect to Wiener measure on
${\mathcal{C}}(M)$, for any linear functional
$\varphi:\text{H}_{1}(M;{\mathbb{R}})\to{\mathbb{R}}$ (Definition 1 and
Theorem 1 below). We remark that all the constructions depend on the choice of
a Riemannian metric on the manifold, since the construction of Brownian motion
depends on the latter.

For more geometrical properties of Brownian motion on a Riemannian manifold,
and as an introduction to the topic, we refer to the beautiful survey by
Kendall [15] on stochastic differential geometry, and of course to the already
classical book by McKean [18].

Since the paper is written by geometers, and will probably be read by both
probabilists and geometers, we find it useful to review a few concepts.
Finally, we want to thank L. Gorostiza and J. M. Gambaudo for helpful
suggestions and references, and also we want to acknowledge conversations with
R. Iturriaga and B. Deroin who pointed out what the asymptotic cycle for the
Poincaré disc should be, and the results by J. J. Prat, [22] and F. Ledrappier
[16] which helped the authors clarify the correct averaging.

##  2\. Brownian motion on a Riemannian manifold (a short review).

The existence of Brownian motion is equivalent to the existence of the Wiener
measure on the space of continuous paths. We review the properties of the
Wiener measure. We first describe the Wiener measure on ${\mathbb{R}}^{N}$. We
will denote the Wiener measure associated to the Brownian motion starting at
$0\in{\mathbb{R}}^{N}$ by $W_{0}$. It is a probability measure on

${\mathcal{C}}_{0}({\mathbb{R}}^{N})=\\{\alpha:[0,\infty)\longrightarrow{\mathbb{R}}^{N},\,\,\alpha(0)=0,\,\,\alpha\,\,\text{continuous}\\},$

constructed in the following way.

Let $\Delta$ be the Laplace operator on ${\mathbb{R}}^{N}$, then the heat
operator $L=\Delta-\partial/{\partial t}$. It has a fundamental solution
$h(x,y,t)=h_{t}(x,y)$ with the following properties:

  1. (1)

$h:({\mathbb{R}}^{N}\times{\mathbb{R}}^{N}-\text{Diagonal})\times{\mathbb{R}}^{(\geq
0)}\mapsto{\mathbb{R}}^{(\geq 0)}$ is smooth.

  2. (2)

$h(x,y,t)=h(y,x,t).$

  3. (3)

$L_{x}h=0=L_{y}h$ where $L_{x}$ is $L$ operating on the first variable of $L$.

  4. (4)

The Chapman-Kolmogoroff identity

$h_{t+s}(x,y)=\int_{{\mathbb{R}}^{N}}h_{s}(x,z)h_{t}(z,y)dz.$

  5. (5)

$\int_{{\mathbb{R}}^{N}}h_{t}(x,y)dy=1$ for all $x\in{\mathbb{R}}^{N}$ and
$t>0.$

The Borel $\sigma$-algebra is generated by the sets
$\rho_{\tau}^{-1}(B)\subset{\mathcal{C}}_{0}({\mathbb{R}}^{N})$ where
$\tau=(t_{1},\dots,t_{n})$ with $0<t_{1}<\dots<t_{n}$, and
$\rho_{\tau}:{\mathcal{C}}_{0}({\mathbb{R}}^{N})\mapsto{\mathbb{R}}^{N}\times\dots\times{\mathbb{R}}^{N}$
(n-copies of ${\mathbb{R}}^{N}$) is the evaluation map
$\rho_{\tau}(\alpha)=(\alpha(t_{1}),\dots\alpha(t_{n}))$, and $B$ is a Borel
subset of ${\mathbb{R}}^{N}\times\dots\times{\mathbb{R}}^{N}$ ($n$-times,
where $n$ is an arbitrarily changing integer). We have

$W_{0}(\rho_{\tau}^{-1}(B))=\int_{B}h_{t_{1}}(0,x_{1})h_{t_{2}-t_{1}}(x_{1},x_{2})\cdots
h_{t_{n}-t_{n-1}}(x_{n-1},x_{n})dx_{1}dx_{2}\cdots dx_{n}$

It is known that $W_{0}$ determines a probability Borel measure on
${\mathcal{C}}_{0}({\mathbb{R}}^{N})$ i.e., with
$W_{0}({\mathcal{C}}_{0}({\mathbb{R}}^{N}))=1$ (a construction essentially due
to Wiener, see [7] or [8] and the references therein)

If $(M,g)$ is a compact Riemannian manifold, we can construct Brownian motion
on $M$ by “rolling” M along a Brownian path in its tangent space. This
construction is formalized below. We remark that the existence of Brownian
motion is equivalent to the existence of the Wiener measure on the space of
continuous paths with a fixed starting point
${\mathcal{C}}_{p}=\\{\alpha:[0,\infty)\mapsto
M,\,\,\alpha(0)=p,\,\,\alpha\,\,continuous\\}$.

Assume that we have a compact Riemannian manifold $(M,g)$, and we have
Brownian motion $B_{p}(t)$ starting at a point $p$. Following Itô [12, 13],
this is constructed by solving a stochastic differential equation.

We assume as given Brownian motion $\beta(t)$ on ${\mathbb{R}}^{n}$ starting
at $0$, identify the tangent space $T_{p}(M)$ isometrically with
${\mathbb{R}}^{n}$ using the Riemannian metric $g_{p}$ on $T_{p}(M)$ and the
usual Euclidean metric on ${\mathbb{R}}^{n}$. Then we obtain the Brownian
motion on the manifold $M$ starting at a point $p$ by the development map at
the point. Formally, one solves a stochastic differential equation as follows
(see [12, 13]). In local coordinates on $M$, we denote by $g_{ij}$ the
components of the Riemannian tensor, $g^{ik}g_{kj}=\delta_{j}^{i}$ and
$\Gamma_{ij}^{k}$ the Christoffel symbols of the Riemannian connection. Let

$\Delta=\frac{1}{2}g^{ij}\nabla_{i}\nabla_{j}$
$=\frac{1}{2}g^{ij}\frac{\partial^{2}}{\partial x^{i}\partial
x^{j}}-\frac{1}{2}g^{ij}\Gamma_{ij}^{k}\frac{\partial}{\partial x^{k}}$

be the Laplace-Beltrami operator.

In the same coordinates, we define $\sigma_{k}^{i}$ and $m^{k}$ by

$\sum_{k}\sigma_{k}^{i}\sigma_{k}^{j}=g^{ij}$
$m^{k}=-\frac{1}{2}g^{ij}\Gamma_{ij}^{k}$

and solve the stochastic differential equation:

$dB^{i}(t)=\sigma_{k}^{i}(B(t))d\beta^{k}(t)+m^{i}(B(t))dt$ (1)

with $B^{i}(0)=p$. We obtain the paths of Brownian motion starting at $p\in
M$. (See the nice intuitive explanation in [26] or the more formal one in
[15]).

######  Definition 1.

Let ${\mathcal{C}}(M)=\\{\alpha:[0,\infty)\longrightarrow
M,\,\,\,\alpha\,\,continuous\\}.$ Let $\Pi:{\mathcal{C}}(M)\to{M}$ be defined
by $\Pi(\alpha)=\alpha(0)$. Then, as shown by [7], ${\mathcal{C}}(M)$ is a
Fréchet manifold and $\Pi$ is a differentiable locally trivial fibration with
fiber ${\mathcal{C}}_{p}=\Pi^{-1}(p)$, for each $p\in{M}$.

Let us define the measure $W$ on ${\mathcal{C}}(M)$ as the measure whose
restriction to a local trivialization
$U\times{\mathcal{C}}_{p}\simeq\Pi^{-1}(U)$ is the volume measure of $M$
restricted to the open set $U\subset{M}$ times $W_{p}$. More precisely, if
${\mathcal{U}}\in{\mathcal{C}}(M)$ is a Borel subset, we have

$W({\mathcal{U}})=\int\limits_{M}W_{p}[{\mathcal{U}}\cap{\mathcal{C}}_{p}(M)]dv_{g}(p),$

where $dv_{g}$ denotes the volume measure on $M$ induced by the Riemannian
metric $g$.

######  Definition 2.

The measure $W$ is the Wiener measure on ${\mathcal{C}}(M)$.

If we define the time shift semigroup
$F_{s}:{\mathcal{C}}(M)\longrightarrow{\mathcal{C}}(M),\,\,s\geq 0,$ as
$F_{s}(\alpha)(t)=\alpha(s+t)$, then this semigroup preserves the measure $W$.
Therefore the triple $({\mathcal{C}}(M),W,F_{s})$ is a measure-preserving non-
reversible dynamical system.

Since almost all paths with respect to Wiener measure are not differentiable
and we will need to integrate closed 1-forms over paths, in this paper we will
follow Itô’s construction in [12], Theorem 1. We review here the way he
approximates the Brownian path by a piecewise smooth curve. Given a Brownian
path $B(t)$ on $M$, consider a partition
${\mathcal{P}}=\\{0=t_{0},t_{1},t_{2},\ldots,t_{n}=T\\}$ of the interval
$[0,T]$. Join each pair of points $B(t_{i})$, $B(t_{i+1})$ with the unique
geodesic path joining these two points, obtaining a piecewise smooth curve
which we denote by $B(t,{\mathcal{P}})$. Itô shows that as the mesh of the
partition tends to zero, the piecewise smooth curve obtained tends to the
given Brownian path. We will use several times below this approximation.

We will also consider from now on, Brownian motion on the interval
$[0,\infty)$, even if in the constructions below, we only use for each path
$B(t)$, the trajectory up to a finite time $T$

######  Remark 1.

Since we will need an approximation of the Brownian curve by piecewise smooth
curves, it is also instructive to describe Pinsky’s method [21] to obtain a
random path on a manifold using the so-called isotropic transport process. The
idea is as follows: we give a partition
${\mathcal{P}}=\\{0=t_{0},\,t_{1},\dots,t_{n}=T\\}$ of the interval $[0,T]$
($n$ arbitrary). With initial point $\gamma(0)$, choose successively an
arbitrary (random) direction in the unit sphere of each $T_{\gamma(t_{i})}M$,
with the uniform distribution in the unit sphere defined using the Riemannian
metric in the corresponding tangent space. Then consider the geodesic path
$\gamma:[t_{i},t_{i+1}]\rightarrow M$ starting at the base point
$\gamma(t_{i})$ in the direction of the unit tangent vector chosen. Pinsky
proves that when the manifold $M$ is (i) compact or (ii) it has bounded
negative curvature, the union of all these geodesics will define a piecewise
smooth curve that converges to a random curve as the size of the partition
${\mathcal{P}}$ tends to zero at an appropriate rate. This is another way to
construct Brownian motion on $M$.

##  3\. The asymptotic behavior of homology associated with Brownian motion

In this section, we use the existence of Brownian motion on a Riemannian
manifold to construct a process
$H_{t}:{\mathcal{C}}(M)\to{\text{Hom}(\text{H}^{1}(M;{\mathbb{R}}),{\mathbb{R}})}\simeq{\text{H}_{1}(M;{\mathbb{R}})}$
($t\in[0,\infty)$) associated to it. From now on we will assume that our
Manifold will be Stochastically Complete, namely that Brownian motion exists
for all nonnegative time $t\geq 0$, more precisely, with probability 1,
Brownian motion starting at a point $X\in M$ does not explode (E. Hsu, page
107 [11]).

We will assume we have chosen a priori, for every pair of points $q_{1}$ and
$q_{2}$ of $M$, a collection of smooth short curves
$\\{\gamma_{q_{1},q_{2}}:[0,1]\to{M}\\}$ of uniformly bounded length joining
the points. Such a collection of curves exists and depends in a measurable way
on $q_{1}$ and $q_{2}$. (See Arnold and Khesin [1], page 155). We give a more
precise definition of a collection of short curves in section 4 below). Our
construction will be independent of the choice of this collection of short
curves.

######  Definition 3.

The Brownian cycle map is the measurable map:

$[\,\,\cdot\,\,]:{\mathcal{C}}(M)\times{\mathbb{R}}^{+}\rightarrow\text{H}_{1}(M,{\mathbb{R}}),$
(2)

defined by the formula $(B,T)\mapsto[\Gamma(B,T)]$, where $\Gamma(B,T)$ is the
closed curve obtained by joining the endpoints of the restriction of a
Brownian trajectory to the interval $[0,T],\,$ $B(t):[0,T]\rightarrow M$ with
the curve $\gamma_{{}_{{B(T)},{B(0)}}}$ of our collection of short curves. By
concatenation we may assume that $\Gamma(B,T)$ is parametrized by a function
$F_{(B,T)}:[0,T+1]\to{M}$ as follows: $F_{(B,T)}(t)=B(t)$ if $t\in[0,T]$ and
$F_{(B,T)}(t)=\gamma_{{}_{{B(T)},{B(0)}}}(t-T)$ if $t\in[T,T+1]$. This defines
a singular cycle $[\Gamma(B,T)]$ in the homology, group
$\text{H}_{1}(M;{\mathbb{R}})=\text{H}_{1}(M;{\mathbb{Z}})\otimes{\mathbb{R}}$.
This cycle is the value of our map $[\,\,\cdot\,\,]$.

To see the meaning of the Brownian cycle map, we use the approximation of the
random path and Poincaré duality as follows. Given a partition ${\mathcal{P}}$
of the interval $[0,T]$, we denote by $B(T,{\mathcal{P}})$ the piecewise
smooth curve as constructed above that approximates $B(t)$. By joining the
endpoints in the same way with $\gamma_{B(T),B(0)}$, we get a cycle that
depends on the random path $B(t)$, the length of the interval $T$, and the
partition ${\mathcal{P}}$. We denote it by $\Gamma^{{\mathcal{P}}}(B,T)$. This
is a piecewise smooth cycle, over which we can integrate a closed 1-form $w$:

$\int\limits_{\Gamma^{{\mathcal{P}}}(B,T)}w\,.$

Then the following limit:

$\int\limits_{\Gamma(B,T)}w=\lim_{max|\Delta t|\to
0}\int\limits_{\Gamma^{{\mathcal{P}}}(B,T)}w\,,$

well-defined. By Stokes’ theorem, it only depends on the cohomology class
$[\omega]\in{\text{H}^{1}(M;{\mathbb{R}})}$ represented by $\omega$. Poincaré
duality permits us to interpret the random cycle $[\Gamma(B,T)]$ as the linear
map $A_{(B,T)}:\text{H}^{1}(M;{\mathbb{R}})\to{\mathbb{R}}$ given by:

$A_{(B,T)}([w])=\underset{\Gamma(B,T)}{\int}w\in{\mathbb{R}}.$ (3)

######  Remark 2.

$A_{(B,T)}([w])$ in formula 3 does not depend on the parametrization of
$\Gamma(B,T)$.

Now consider the time shift map $F_{s}$ acting in the image of our map
$[\,\,\cdot\,\,]$, which is just a translation by time $s$ on our Brownian
path. For $s\geq 0$, it acts as follows. We will first describe it for the
piecewise smooth approximation. Given the path $B(t)$, consider the curve
given by the composition $F_{s}(B(t))=B(s+t)$. Now, construct in the same way,
the piecewise smooth approximation to this new random path (which starts in
$B(s)$), and the cycle $(F_{s})_{*}[\Gamma^{{\mathcal{P}}}(B,T)]$ as the cycle
obtained with the piecewise approximation of the curve $F_{s}\cdot
B:[s,s+T]\mapsto M$, closed with a curve $\gamma_{B(T+s),B(s)}$ from our fixed
collection of short curves. Of course, the map behaves well after taking
limits, giving a time translation on the cycle constructed above,
$(F_{s})_{*}[\Gamma(B,T)]$.

######  Lemma 1.

For a fixed closed 1-form $\omega$ and $B\in{\mathcal{C}}(M)$, the map
$T\mapsto{A_{(B,T)}([\omega])}$ is measurable, and subadditive with respect to
the semigroup action $B\mapsto(F_{s})(B),\,\,s\geq 0,\,B\in{\mathcal{C}}(M)$.

###### Proof.

Suppose we denote by ${\mathcal{P}}$ a partition of the interval $[0,T+S]$
such that both $S$ and $T$ are on the partition. Since we will need to
consider the cycle $\Gamma^{{\mathcal{P}}}(B,T+S)$ in two pieces, first going
to time $S$, then from $S$ to time $T+S$. Let us denote by
$\Gamma^{{\mathcal{P}}}(B,S)$, and consider the curve $B$ under the shift
$F_{S}$ by $S$ units, denoted by $B_{F_{S}}$. Then it defines a cycle
$\Gamma^{{\mathcal{P}}}(B_{F_{S}},T)$. (To simplify notation we are using the
same notation ${\mathcal{P}}$ for the part of the partition on each of the
corresponding subintervals). We have:

$\int\limits_{\Gamma^{{\mathcal{P}}}(B,S+T)}w=\int\limits_{\Gamma^{{\mathcal{P}}}(B,S)}w\,+\int\limits_{\Gamma^{{\mathcal{P}}}(F_{S}(B),T)}w\,-\int\limits_{\gamma_{B(S),B(0)}}w.$

Since the curves $\gamma_{p,q}$ are uniformly bounded,
$\left\lvert\int\limits_{\gamma_{{}_{B(S),B(0)}}}w\right\rvert\leq{c},\,\forall\,\gamma_{{}_{B(S),B(0)}}$,
where $c>0$ is a constant, depending only on the 1-form $w$. Consequently,

$\int_{\Gamma^{{\mathcal{P}}}(B,S+T)}w\leq\left[\int_{\Gamma^{{\mathcal{P}}}(B,S)}w\right]+\left[\int_{\Gamma^{{\mathcal{P}}}(F_{S}(B),T)}w\right]+c$

Now take limits as the mesh of the partition ${\mathcal{P}}$ tends to zero, we
obtain.

$A_{(B,T+S)}([\omega])=\left[\int_{\Gamma(B,S+T)}w\right]\leq\left[\int_{\Gamma(B,S)}w\right]+\left[\int_{\Gamma(F_{S}(B),T)}w\right]+c,$

Therefore we have the following subadditive property:

$A_{(B,T)}([\omega])\leq A_{\Gamma(B,T)}w+A_{\Gamma(F_{S}(B),T)}w+c$ (4)

Measurability is a direct consequence of the properties of Brownian motion. ∎

The subadditive ergodic theorem ([28], Theorem 10.1, page 231), applied to the
function $T\mapsto{A_{(B,T)}([\omega])}$ and the dynamical system
$F_{s}:{\mathcal{C}}(M)\to{\mathcal{C}}(M)$, implies the following corollary:

######  Corollary 1.

For each fixed homology class $[\omega]$ represented by the closed 1-form
$\omega$ the limit

${\mathcal{A}}(B,[\omega]):=\lim_{T\to\infty}\frac{A_{(B,T)}([\omega])}{T}$
(5)

is equal to 0 for almost every $B\in{\mathcal{C}}(M)$, with respect to
Wiener’s measure.

######  Remark 3.

In the following, at the end of Theorem 1, we will show that for every closed
one form $\omega$ and almost every $B\in{\mathcal{C}}(M)$, with respect to the
Wiener measure, the limit in equation 5 is zero.

However, if instead of dividing by $T$ in equation 5 we divide by
$\log(T^{\prime})$ (where $T^{\prime}$ is an appropriate random change of time
of Brownian motion on $M$) we will have a nontrivial random function. In fact,
we will show the following:

######  Theorem 1.

Let us consider the stochastic process

$A_{t}:{\mathcal{C}}(M)\times{\text{H}^{1}(M,{\mathbb{R}})}\to{\mathbb{R}},\,t\geq
0,$

given by the rule $A_{t}(B,\alpha)=A_{(B,t)}([w])$ where $\omega$ is a closed
1-form representing the class $\alpha\in\text{H}^{1}(M,{\mathbb{R}})$,
$B\in{\mathcal{C}}(M)$ and $A_{(B,t)}([w])$ is given by formula 3. Then, for
each fixed $B\in{\mathcal{C}}(M)$ and $t\geq 0$, the map
$h_{(t,B)}:\text{H}^{1}(M,{\mathbb{R}})\to{\mathbb{R}}\,,$
$h_{(t,B)}(\omega)\overset{def}{=}A_{t}(B,\omega)$ is a linear map and so, by
duality, is an element of $\text{H}_{1}(M;{\mathbb{R}})$. This produces for
each $t\in[0,\infty)$ a measurable map
$H_{t}:{\mathcal{C}}(M)\to{\text{H}_{1}(M;{\mathbb{R}})},\,B\overset{H_{t}}{\longmapsto}{h_{(t,B)}}$.

  1. (1)

Therefore, the process $\\{{H_{t}}\\}_{t\geq 0}$ is a multivariate process on
${\mathcal{C}}(M)$, with values in
${\text{H}_{1}(M;{\mathbb{R}})\simeq{\mathbb{R}}^{\beta_{1}}}$ (where
$\beta_{1}$ is the first Betti number of $M$). For each fixed homology class
$\alpha\in{\text{H}_{1}(M;{\mathbb{R}})},\,\alpha\neq 0$, represented by the
closed 1-form $\omega\neq 0$, and $t\geq 0$ one has the random variable
$H_{t}(\omega):{\mathcal{C}}(M)\to{\mathbb{R}}$,
$H_{t}(\omega)(B)=A_{(B,t)}([\omega])$. Then, after a measurable change of
time in the Brownian motion, with the new clock with time $t^{\prime}$ and
every $[\omega]\in{\text{H}^{1}(M;{\mathbb{R}})}$, the following holds:

$\lim_{t^{\prime}\to\infty}\,W\left(\frac{H_{t^{\prime}}(\omega)}{\log(t^{\prime})}\leq{y}\right)=\lim_{t^{\prime}\to\infty}\,W\left(\frac{A_{(B,t^{\prime})}([\omega])}{\log(t^{\prime})}\leq{y}\right)=\frac{1}{\pi}\int^{y}_{-\infty}\frac{dx}{1+x^{2}},$

where $W$ is Wiener measure on ${\mathcal{C}}(M)$ (definition 2).

  2. (2)

Therefore, the process $\\{{H_{t}}\\}_{t\geq 0}$ is a multivariate process on
${\mathcal{C}}(M)$, which converges in distribution to a multivariate Cauchy’s
distribution with parameter 1.

Since $\log{T^{\prime}}=o(T^{\prime})$, as $T^{\prime}\to\infty$ it follows
that for any $[\omega]$ the limit in equation 5 is 0 for almost all
$B\in{\mathcal{C}}(M)$, as stated in Corollary 1.

The proof of Theorem 1 consists essentially in reducing the theorem to the
particular case when $M={\mathbb{S}}^{1}$, because in this case, the theorem
is essentially Spitzer’s theorem [25]. Therefore, before the proof of the
theorem, we will state the theorem of Spitzer in a form that is adapted to our
proof.

Let ${\mathbb{S}}^{1}=\\{z\in{\mathbb{C}}\,:\,|z|=1\\}$ be the unit circle in
the complex plane. Let
${\mathcal{C}}({\mathbb{S}}^{1})=\\{\gamma:[0,\infty)\to{\mathbb{S}}^{1}\,:\,\,\text{such
that }\,\gamma\,\text{is continuous}\\}$. Let $P_{\bullet}$ be the standard
Wiener measure on ${\mathcal{C}}({\mathbb{S}}^{1})$ (i.e., paths starting at
any point on the circle 2).

Let $\beta\in{\mathcal{C}}({\mathbb{S}}^{1})$. Let
$p:{\mathbb{R}}\to{\mathbb{S}}^{1},\,p(t)=e^{t\mathbf{i}}$ be the universal
covering projection of ${\mathbb{R}}$ over of the circle. Since $[0,\infty)$
is contractible, $\beta$ can be lifted i.e., there exists a unique continuous
map $\tilde{\beta}:[0,\infty)\to{\mathbb{R}}$ such that
$\tilde{\beta}(0)\in[0,2\pi)$ and $\beta(t)=e^{{\mathbf{i}}\tilde{\beta}(t)}$.
Then, $\tilde{\beta}$ measures the angular change of $\beta$ i.e., the total
change of angle of $\beta$ in the interval $[0,T]$ is
$\tilde{\beta}(T)-\tilde{\beta}(0)=\tilde{\beta}(T)$. This is also the winding
number around the origin of the closed path in the unit circle, in the complex
plane, determined by $\beta$, when $\beta(0)=\beta(T)$.

Let $\beta\in{\mathcal{C}}({\mathbb{S}}^{1})$. Let $T\geq 0$, and let
$\beta_{T}$ be the restriction of $\beta$ to the interval $[0,T]$. Let
$C_{\beta}(T)\in{\text{H}_{1}({\mathbb{S}}^{1},{\mathbb{R}})}$ be the element
in $\text{H}_{1}({\mathbb{S}}^{1},{\mathbb{R}})\simeq{\mathbb{R}}$,
represented by the cycle on the circle obtained by closing $\beta_{T}$ by an
arc on the circle, of length $\leq\pi$, joining the end-points $\beta_{T}(0)$
with $\beta_{T}(T)$. Then, as shown by Spitzer [25] (with a shorter proof by
Durrett [5], see also [14] 7.15 page 270):

$\lim_{t\to\infty}\,P_{\bullet}\left(\frac{2{\tilde{\beta}(t)}}{\log(t)}\leq{y}\right)=\frac{1}{\pi}\int^{y}_{-\infty}\frac{dx}{1+x^{2}},$

which implies (because the short paths joining the endpoints of $\beta_{t}$
have uniformly bounded length) that:

$\lim_{t\to\infty}\,P_{\bullet}\left(\frac{{2C_{\beta}(t)}}{\log(t)}\leq{y}\right)=\frac{1}{\pi}\int^{y}_{-\infty}\frac{dx}{1+x^{2}},$

As a consequence, one has the following proposition:

######  Proposition 1 (Spitzer’s law,[25],[5] and [14] 7.15 page 270).

Let, for $t\geq 0$,

$A_{t}:{\mathcal{C}}({\mathbb{S}}^{1})\to\text{H}_{1}({\mathbb{S}}^{1},{\mathbb{R}})\simeq{\mathbb{R}}$

be the measurable function given by the formula
$A_{t}(\beta)=\frac{2{C_{\beta}(t)}}{\log(t)}$. Then, as $t\to\infty$, $A_{t}$
converges in distribution to a Cauchy distribution with parameter 1, with
respect to the Wiener measure $P_{\bullet}$ on the circle. More precisely:

$\lim_{t\to\infty}P_{\bullet}\\{A_{t}\leq{y}\\}=\frac{1}{\pi}\int_{-\infty}^{y}\frac{dx}{1+x^{2}}$

###### End of the proof of Theorem 1.

Since a linear combination of Cauchy-distributed random variables has also a
Cauchy distribution it is enough to prove the theorem for
$\alpha\in\text{H}^{1}(M,{\mathbb{Z}})$ represented by a closed 1-form
$\omega$ with integer periods because we can choose a basis of
$\text{H}^{1}(M,{\mathbb{R}})$ with integer periods. Such a basis generates
$\text{H}^{1}(M,{\mathbb{Z}})$. The group $\text{H}^{1}(M;{\mathbb{Z}})$
isomorphic to the abelian group $[M,{\mathbb{S}}^{1}]$ of homotopy classes
$[f]$ of continuous maps $f:M\to{\mathbb{S}}^{1}$ to the circle, with its
natural multiplication $\star$ given by the rule $[f]\star{[g]}=[f\cdot{q}]$ ,
where $f\cdot{g}$ is pointwise multiplication. We recall that
$\text{H}_{1}(M;{\mathbb{Z}})\subset\text{H}_{1}(M;{\mathbb{R}})$ is a lattice
of maximal rank. By Hodge theorem, given the class
$[\alpha]\in\text{H}^{1}(M;{\mathbb{R}})$, there exists a harmonic 1-form
$\omega$ representing $\alpha$ in the de Rham cohomology (see [9], sect.
(7.1)). If $\omega$ has integer periods (i.e.,
$\omega\in{\text{H}^{1}(M;{\mathbb{Z}})}$), then there exist a harmonic map
$\phi:M\to{\mathbb{S}}^{1}$ such that $\omega=\phi^{*}(d\theta)$ and such that
its homotopy class represents $\alpha$, here $d\theta$ is the canonical
nowhere-vanishing 1-form on the circle such that
$\int_{{\mathbb{S}}^{1}}d\theta=2\pi$ (see [27]). Then, $\phi$ is unique up to
post-composition with a rotation of the circle (Example 7, 3.3.8 in [2]). In
fact, $\phi$ is a harmonic morphism. Recall that a harmonic morphism is a map
$\phi$ that is harmonic and horizontally conformal, that is, at any point
$p\in M$ at which $d\phi_{p}\neq 0$, the restriction of $d\phi_{p}$ to the
orthogonal complement of $Ker(d\phi_{p})$ in $T_{p}(M)$ is conformal and
surjective. In our case, the image is of dimension one therefore, it is
immediately conformal. Since $\phi$ is a harmonic morphism, the image of
Brownian motion on $M$ is Brownian motion on ${\mathbb{S}}^{1}$ up to a random
bounded time change. This is carefully explained in the Lemma and Theorem C of
section 8 of the article by R. W. R. Darling [3], (see also[15], [10] Pag.
395, sect. 2.44).

The family of cycles $\\{\Gamma(B,T)\\}_{\\{B\in{\mathcal{C}}(M),T\geq 0\\}}$
given in Definition 3 are mapped under the harmonic map $\phi$ onto the family
of cycles $\\{\phi(\Gamma(B,T)\\})$ in ${\mathbb{S}}^{1}$. These cycles are
parametrized by $\phi\circ{F_{(B,T)}}$ as in definition 3. We can assume that
the collection of small closing curves in $M$ goes to the collection of small
closing curves in ${\mathbb{S}}^{1}$ so that $\phi|_{\Gamma(B,T)}$ is the
cycle in $S^{1}$ representing
$\phi_{*}([\Gamma(B,T)])\in\text{H}_{1}({\mathbb{S}}^{1},{\mathbb{Z}})$
corresponding to Brownian motion with a time change. That is, the image of the
cycle associated with a Brownian path in $M$ is a cycle corresponding to a
Brownian path in $S^{1}$. Then, by Spitzer’s theorem one has:

$\lim_{t\to\infty}\,P_{\bullet}\left(\frac{\phi_{*}([\Gamma(B,t)])}{\log{t}}\leq{y}\right)=\frac{1}{\pi}\int^{y}_{-\infty}\frac{dx}{1+x^{2}}.$

This finishes the proof since

$\phi_{*}([\Gamma(B,t)])=A_{(B,T)}([w])=\underset{\Gamma(B,T)}{\int}w\in{\mathbb{R}},$

and the change of time is a reparametrization of the curves $\Gamma(B,T)$
which does not affect the integration of closed 1-forms over these cycles. ∎

##  4\. Lifting to Holonomy

In this section, we review the construction of parallel transport along a
Brownian path for a principal bundle $P\rightarrow M$ with structure group a
compact Lie group $G$. We will follow closely Itô’s construction [12] of
parallel transport for the Riemannian connection on $M$.

Assume we have a principal bundle $P\rightarrow M$ with structure group $G$ a
compact connected Lie group. We suppose we have a connection, represented by a
Lie algebra valued 1-form $A\in\mathcal{A}^{1}(M,\mathfrak{g})$, the 1-forms
on $P$ with values in the Lie algebra $\mathfrak{g}$ of $G$. Given any closed,
$C^{1}$ curve $\gamma:[0,T]\mapsto M$ on $M$ starting and ending at a point
$p\in M$, we have the holonomy or parallel transport along $\gamma$ with
respect to the connection $A$ given by the time ordered (or path-ordered)
exponential integral as:

$P\,exp\int_{\gamma}A=$ $I+\int\limits_{0\leq s\leq
T}A(\gamma(s)){\gamma}^{\prime}(s)ds+\iint\limits_{0\leq s_{1}\leq s_{2}\leq
T}A(\gamma(s_{1}))A(\gamma(s_{2})){\gamma}^{\prime}(s_{2}){\gamma}^{\prime}(s_{1})ds_{2}ds_{1}+\cdots$

We observe that the ordering is important because in general the evaluation of
the connection form at different places of the curve give matrices that do not
commute in the multiple integrands, so we must take an ordering of the
multiplications of the connection forms evaluated under the multiple
integrals, and that causes the factorial constants to disappear (see Nelson
[20] pp. 15-16). In the case of an abelian Lie group, it is not necessary to
specify an ordering and the path-ordered exponential is an honest exponential.

We follow Itô’s construction closely. Consider a Brownian trajectory $B(t)$.
For each time $T$, we use the collection of closing curves to obtain a closed
cycle as above. Again, with a partition ${\mathcal{P}}$ of the interval
$[0,T]$, we approximate the closed cycle by a piecewise smooth closed curve,
obtaining the cycle $\Gamma^{{\mathcal{P}}}(B,T)$ as in the previous section.
Given any connection on the principal bundle, it is possible to construct the
holonomy along this cycle.

######  Theorem 2.

With the data as in the above paragraph, one can construct the holonomy along
the asymptotic cycle defined by Brownian motion. We denote it by $Hol_{A}(B)$.
It is a measurable function $Hol_{A}:{\mathcal{C}}_{p}\longrightarrow G.$

###### Proof.

The proof is the same as in Itô [12], except that we now use any connection on
a principal bundle, instead of the Riemannian connection. Another way to prove
it is as follows. We need to take the time-ordered exponential integral of the
connection form $A$ along the cycle, so we use the approximations of the
asymptotic cycle $\Gamma(B)$ by the piecewise smooth cycle
$\Gamma^{{\mathcal{P}}}(B,T)$ as described above. Now consider the path-
ordered exponential integral of the connection form $A$ over the piecewise
smooth cycle,

$Pexp\int\limits_{\Gamma^{{\mathcal{P}}}(B,T)}A.$

Each term of the path-ordered exponential integral can be uniformly bounded by
a bound of the value of the connection form on the compact manifold $M$.
Dividing the cycle by $T$, we get:

$Pexp\int_{\frac{1}{T}\Gamma^{{\mathcal{P}}}(B,T)}A$

where this is interpreted as dividing the $k^{th}$ term of the series of
integrals written above by $1/{T^{k}}$ (counting $k$ from $0$).

Take the limit as the mesh of ${\mathcal{P}}$ goes to zero, and the limit as
$T\rightarrow\infty$

$\lim_{T\to\infty}\left(\lim_{mesh{({\mathcal{P}})}\to
0}Pexp\int_{\frac{1}{T}\Gamma^{{\mathcal{P}}}(B,T)}A\right).$

The limits exist for each summand in the expression of the path-ordered
exponential integral with probability one in ${\mathcal{C}}$, by the lemmas of
the previous section. Note that we need the subadditivity of the cycles (in
the argument $T$), implying that each summand of the path-ordered exponential
defines a subadditive function in $T$ also. Since the convergence of the
series is uniform, this proves the theorem. It is obvious that this limit
coincides with the limit of the parallel transport on the closed cycles as
considered in [12]. ∎

We recall that the Wilson loop associated to a closed curve $\gamma$ in $M$,
is the trace of the holonomy determined by the curve,
$Tr_{\rho}(Hol_{A}(\gamma))$, where the trace $Tr_{\rho}$ depends on a
representation $\rho$ of the structure group $G$ of the principal bundle. It
is independent of a choice of a base point. Then we have:

######  Corollary 2.

For a Brownian trajectory on a compact Riemannian manifold, the asymptotic
Wilson loop for Brownian paths $Tr_{\rho}(Hol_{A}(B))$, is well defined.

######  Remark 4.

We note that J. Eells has given a construction of Brownian motion on the
fundamental group of a Riemannian manifold and as extension, a construction of
Brownian motion in the holonomy group corresponding to the Riemannian
connection [6].

Acknowledgements. The first author benefited from a grant IN103324 of PAPIIT,
DGAPA, Universidad Nacional Autónoma de México, Proyecto Laminaciones y
foliaciones reales y complejas y su interacción con la geometría, topología,
teoría de números y sistemas dinámicos. The second author by CONAHCYT research
grant A1-S-45886.

## References

  * [1] V. I. Arnold, B. A. Khesin. Topological Methods in Hydrodynamics. 2nd. Ed 2021\. Springer Verlag
  * [2] P. Baird, J. C.Wood, Harmonic morphisms between Riemannian manifolds. London Mathematical Society Monographs. New Series, 29. The Clarendon Press, Oxford University Press, Oxford, 2003
  * [3] R. W. R. Darling, Martingales in manifolds. Definition, examples and behaviour under maps Séminaire de probabilités (Strasbourg), tome S16 (1982), p. 217–236.
  * [4] G. de Rham, Georges.Differentiable manifolds. Forms, currents, harmonic forms. Grundlehren der mathematischen Wissenschaften. Vol. 266. Translated by Smith, F. R. With an introduction by S. S. Chern. (Translation of 1955 French original ed.). Berlin: Springer-Verlag
  * [5] R. Durret. A new proof of Spitzer’s result on the winding of two-dimensional Brownian motion Ann. of Prob. 10 1 244–246.
  * [6] J. Eells. Random Walk on the Fundamental Group. Proceed. of Symp. in Pure Math. AMS, 27, 1975, 211–217
  * [7] J. Eells, K. D. Elworthy. Wiener Integration on Certain Manifolds. Problems in Non-linear Analysis (CIME), Edizioni Cremonese, 1971, 67–94.
  * [8] J. Eells, K. D. Elworthy. Stochastic Dynamical Systems. Control Theory and Topics in Functional Analysis, International Seminar ICTP; Trieste III 1974, 179–185.
  * [9] J. Eells, L. Lemaire. A report on harmonic maps. Bull. London Math Soc. 10, 1978, 1–68.
  * [10] J. Eells, L. Lemaire. Another report on harmonic maps. Bull. London Math. Soc. 20, 1988, 285–524.
  * [11] E. Hsu. Stochastic Analysis on Manifolds. 2002\. GSM 138. AMS.
  * [12] K. Itô. The Brownian motion and tensor fields on a Riemannian manifold. Proc. Intern. Congr. Math. Stokholm, 1963, 536–539.
  * [13] K. Itô. Stochastic differential equations in a differentiable manifold. Nagoya Math. J., 1 (1950), 35–47.
  * [14] K. Itô, H. P. McKean Jr. Diffusion processes and their sample paths. Second printing, corrected. Die Grundlehren der mathematischen Wissenschaften, Band 125. Springer-Verlag, Berlin-New York, 1974. xv+321 pp.
  * [15] W. S. Kendall. Stochastic Differential Geometry: An Introduction. Acta Appl. Mathematicæ, 9, 1987, 29–60.
  * [16] F. Ledrappier Central Limit Theorem in Negative Curvature Ann. Probab. 23(3): 1219–1233 (July, 1995).
  * [17] P. Lévy. Processus Stochasitques et Mouvement Brownien, 1965, Gauthier-Villars
  * [18] H. P. McKean. Stochastic Integrals . Academic Press, 1969
  * [19] J. W. Milnor. Infinite cyclic coverings. Conference on the Topology of Manifolds (Michigan State Univ., E. Lansing, Mich., 1967), pp. 115–133.
  * [20] E. Nelson. Topics in Dynamics I. Flows. Princeton Univ. Press, 1969
  * [21] M. A. Pinsky. Isotropic Transport Process on a Riemannian Manifold. Trans. of the AMS, 218, 1976, 353–360.
  * [22] J. J. Prat. Étude asymptotique du mouvement brownien sur une variété riemannienne à courbure négative. C. R. Acad. Sci. Paris Sér. A-B 272 (1971), A1586–A1589.
  * [23] J. Pitman, M. Yor. Asymptotic Laws of Planar Brownian motion. Ann. of Prob. 11 (3), 1986, 733–779.
  * [24] S. Schwartzman. Asymptotic cycles. Ann. Math. 66, 2, 1957, 270–284
  * [25] F. Spitzer. Some theorems concerning two-dimensional Brownian motion. Trans. AMS, 87, 1958, 187–197.
  * [26] D. W. Stroock. Gaussian Measures in traditional and not so traditional settings. Bull. AMS, 33, 2, 1996, 135–255.
  * [27] D. Tischler, On fibering certain foliated manifolds over ${\mathbb{S}}^{1}$. Topology 9 (1970), 153–154.
  * [28] P. Walters. An Introduction to Ergodic Theory. Springer Verlag, 1982


# Truncation technique for variational quantum eigensolver for Molecular
Hamiltonians

Qidong Xu  Kanav Setia1<EMAIL_ADDRESS>1qBraid Co., 111 S Wacker Dr.,
Chicago, IL 60606, USA

###### Abstract

The variational quantum eigensolver (VQE) is one of the most promising quantum
algorithms for the near-term noisy intermediate-scale quantum (NISQ) devices.
The VQE typically involves finding minimum energy of a quantum Hamiltonian
through classical optimization of a parametrized quantum ansatz. One of the
bottlenecks in VQEs is the number of quantum circuits to be measured. In this
work, we propose a physically intuitive truncation technique that starts the
optimization procedure with a truncated Hamiltonian and then gradually
transitions to the optimization for the original Hamiltonian via an operator
classification method. This strategy allows us to reduce the required number
of evaluations for the expectation value of Hamiltonian on a quantum computer.
The reduction in required quantum resources for our strategy is substantial
and likely scales with the system size. With numerical simulations, we
demonstrate our method for various molecular systems.

##  I Introduction

Recent years have seen growing interest in variational quantum algorithms
(VQA) Peruzzo  _et al._ (2014) for applications in diverse fields Cerezo  _et
al._ (2021a); Amaro  _et al._ (2022); Bravo-Prieto  _et al._ (2023); Blunt
_et al._ (2022); Landman  _et al._ (2022), as they are particularly suitable
to study the potential advantage of quantum computation in the current noisy
intermediate-scale quantum (NISQ) computers stage. VQAs are typically used to
calculate the ground state energy of a Hamiltonian dictated by the problem of
interest. For such problems, it is often efficient to calculate the
expectation value of a given observable of the Hamiltonian for a parameterized
quantum state on a quantum computer. The expectation value of the Hamiltonian
is then used to calculate the updated parameters with the help of an optimizer
on a classical computer. The evaluation of the expectation value of the
Hamiltonian, followed by the calculation of updated parameters for the ansatz
on a classical computer, forms a loop that is run until this procedure
converges to a value.

Variational Quantum Eigensolver, a subclass of VQAs has been a subject of
extensive exploration for quantum chemistry McClean  _et al._ (2016); Kandala
_et al._ (2017); Cao  _et al._ (2019); Bauer  _et al._ (2020); Lee  _et al._
(2021), condensed matter Bravo-Prieto  _et al._ (2020); Uvarov  _et al._
(2020), and many other areas Cao  _et al._ (2018); Banuls  _et al._ (2020);
Liu  _et al._ (2022). This is primarily because there are many known
transformations from fermionic and bosonic operators to the qubit operators
Jordan and Wigner (1993); Somma (2005), and many physics-inspired ansatzes
could be constructed easily on a quantum computer O’Malley  _et al._ (2016);
Wecker  _et al._ (2015). Further, many physics-inspired ansatzes could be
efficiently prepared on a quantum computer. For recent reviews on VQE and best
practices, see Cerezo  _et al._ (2021b); Tilly  _et al._ (2022) and references
therein.

Since the proposal of VQE Peruzzo  _et al._ (2014), there have been many
extensive studies investigating the best ansatz Grimsley  _et al._ (2019);
Mehendale  _et al._ (2023), the classical optimizers used Arrasmith  _et al._
(2020); Lavrijsen  _et al._ (2020), and the convergence Cerezo  _et al._
(2021c); Arrasmith  _et al._ (2021) for VQE. The choice of the ansatz to build
parameterized quantum circuits is of particular importance as it typically can
only span a subspace of the full Hilbert space and, therefore, essentially
sets the upper limit for the distance between the final output minimal value
found by the optimizer to the true minimal cost function values corresponding
to the ground energy of the Hamiltonian. It also impacts the convergence speed
for optimization. Deeper (more complicated) quantum circuits allow for more
precise final results as they can explore larger subspaces, but the fact that
we are still in the NISQ age with noisy qubits also means that deeper circuits
are more vulnerable to noise-induced errors. There have been many recent
exploration in this area, trying to achieve the optimal balance between these
two issues Lee  _et al._ (2018); Ryabinkin  _et al._ (2018); Grimsley  _et
al._ (2019); Lang  _et al._ (2020); Ryabinkin  _et al._ (2020); Tang  _et al._
(2021); Huang  _et al._ (2022). For the classical optimization component,
studies have been focusing on how to improve the performance of the classical
optimizer, including achieving faster convergence and robustness against the
sampling and gate noise in quantum calculation Spall (1992); Kandala  _et al._
(2017); Nakanishi  _et al._ (2020); Sung  _et al._ (2020); Koczor and Benjamin
(2022); Kübler  _et al._ (2020); Sim  _et al._ (2021). Since there already
exists many optimizers in the long study of optimization field, several works
have also been devoted to comparing the performance of these classical
optimizersNannicini (2019); Lavrijsen  _et al._ (2020); Pellow-Jarman  _et
al._ (2021); Bonet-Monroig  _et al._ (2023).

One area where there has been some work but more effort is needed is the
number of terms to be measured for the Hamiltonian. Previous works in this
area have focused on strategies where the number of circuits to be run is
constrained. In this work, we propose a new strategy for reducing the number
of measurements for the circuits to be run in the optimization routine in VQE.
Our strategy does not depend on the optimization routine and is compatible
with many optimizers. It recognizes that many quantum chemistry Fermionic
Hamiltonians, when transformed to qubit operators, result in a large number of
qubit operator terms with small coefficients. This allows us to engineer an
adaptive strategy, wherein we start with a truncated Hamiltonian, which
contains a much smaller number of terms, and gradually increase the number of
terms in the qubit operator Hamiltonian during the optimization procedure. In
this work, we propose two methods for obtaining the truncated Hamiltonian and
discuss strategies for transitioning from truncated Hamiltonian to full
Hamiltonian to reduce the number of terms to be measured during the
optimization procedure. The first one depends on simple cutoffs which serves
as a motivational illustration purpose, and in the second method, we propose a
systematic truncation strategy via an operator classification with clear
physics intuition.

This paper is organized as follows. In Sec.II, we lay out the necessary
backgrounds for the VQE. In Sec.III, we describe two algorithms based on the
idea of Hamiltonian truncation for VQE. We start with a naive algorithm for
illustrative and motivational purposes, and we discuss its effectiveness with
numerical results validating a simple case of the algorithm. We then analyze
the shortcomings of our naive algorithm and present the main result in this
work: a truncation algorithm based on the operator classification method,
where we also provide numerical simulation results to demonstrate its
performance in comparison with the standard VQE algorithm. We conclude in
section IV and discuss some future directions.

##  II Background

We consider a Fermionic Hamiltonian in the second quantization formalism for a
given system:

$\displaystyle
H_{f}=\sum_{ij}^{M}h_{ij}a_{i}^{\dagger}a_{j}+\frac{1}{2}\sum_{ijkl}^{M}h_{ijkl}a_{i}^{\dagger}a_{j}^{\dagger}a_{k}a_{l}$
(1)

where $M$ is the number of fermionic modes in the system, $h_{ij}$ and
$h_{ijkl}$ are one-electron and two-electron integrals of the operators
projected into a given basis set, and $a_{i}^{\dagger}$, and $a_{i}$ are
fermionic creation and annihilation operators satisfying the following anti-
commutation relations:

$\displaystyle\left\\{a_{i},a_{j}^{\dagger}\right\\}=\delta_{ij}\mathbf{1}.$
(2)

One of the main problems of relevance in quantum chemistry is the calculation
of the energy of the ground-state and excited states of the above Hamiltonian
in Eq. (1). Finding these energies on classical computers is generally
intractable because the Fock space associated with the Hamiltonian above
scales exponentially, requiring exponential computational resources to store
the general state of fermions and calculate the energies of various quantum
states. Whereas on a quantum computer, the number of qubits required to
represent the fermionic state scale linearly.

The Fermionic Hamiltonian, Eq. (1), could be transformed into a qubit operator
Hamiltonian using various fermion-to-qubit transformations in the literature,
such as Jordan-Wigner Jordan and Wigner (1993), Bravyi-Kitaev Bravyi and
Kitaev (2002), Generalized Superfast Encoding Setia  _et al._ (2019), etc. The
transformed Hamiltonian is a weighted sum of qubit operators, which are
strings of Paulis, given as follows:

$\displaystyle H_{q}=\sum_{k=1}^{K}c_{k}\hat{P}_{k},$ (3)

where $c_{k}$ is the coefficient of the Pauli string $P_{k}$, and the Pauli
string is the tensor product of $N$ Pauli matrices that act on $N$ number of
qubits separately:

$\displaystyle\hat{P}_{k}=\hat{p}_{1k}\otimes\hat{p}_{2k}...\otimes\hat{p}_{Nk},$

where, $\hat{p}_{ik}\in\\{I,\sigma_{x},\sigma_{y},\sigma_{z}\\}$ with $I$ the
identity matrix, and $\sigma_{x},\sigma_{y}$, and $\sigma_{z}$ Pauli matrices.
The number of Pauli strings $K$ and the number of qubits $N$ depend on the
dimensionality of the original $H_{f}$ and the fermion-to-qubit mapping used.
With the fermionic Hamiltonian $H_{f}$ transformed to qubit operator
Hamiltonian $H_{q}$, the VQE problem reduces to finding the ground state
energy of the qubits Hamiltonian $H_{q}$. In a variational quantum
eigensolver, one prepares a parameterized ansatz quantum state that can be
constructed using various quantum gates. Starting from the initial state of
all qubits being in the $|0\rangle$ ground state, the parameterized quantum
circuit can be constructed by applying a unitary $U(\bm{\theta})$ to the
initial state. The expectation value of the Hamiltonian in Eq. (3) with
respect to the prepared parameterized quantum state is obtained by measuring
the expectation value of each Pauli string and taking their weighted sum:

$\displaystyle\langle
0|U^{\dagger}(\bm{\theta})H_{q}U(\bm{\theta})|0\rangle=\sum_{k=1}^{K}c_{k}\langle
0|U^{\dagger}(\bm{\theta})\hat{P}_{k}U(\bm{\theta})|0\rangle.$ (4)

Upon the measurement of the expectation value of the Hamiltonian, the
parameters are updated by an optimization routine which is run until it
converges to a minimum value.

$\displaystyle\bm{\theta}^{*}=\min_{\bm{\theta}}\sum_{k=1}^{K}c_{k}\langle
0|U^{\dagger}(\bm{\theta})\hat{P}_{k}U(\bm{\theta})|0\rangle.$ (5)

The optimization algorithm is run on a classical computer. We note that since
the parametrized state usually can only span a subspace of the Hilbert space
for the qubits, the minimized expectation value can only be regarded as an
upper bound for the ground energy; finding an efficient ansatz for the
parametrized state that can generate final result close to the true eigenvalue
with fast convergence speed is of key importance when designing the VQE
experiment.

##  III Truncation technique for qubits Hamiltonian

As the size of the system grows, one of the thing that is expected to slow the
optimization procedure is the number of evaluations (measurements) on a
quantum computer. Note that, from Eq. (1), it is apparent that the number of
terms grows as $O(M^{4})$. This implies that the cost function, defined in Eq.
(5), also contains $O(M^{4})$ Pauli strings. In recent years, many studies
have aimed to reduce the depth of the ansatz and reduce the number of
measurements required for each iteration of the VQE Arrasmith  _et al._
(2020). Further, there have also been studies that have looked at the
performance of the different classical optimizers for the cost function
defined in Eq. (5) Nakanishi  _et al._ (2020); Sung  _et al._ (2020); Koczor
and Benjamin (2022); Kübler  _et al._ (2020); Sim  _et al._ (2021). In this
work, we present techniques that complement the previous techniques.

Our strategy is based on the observation that even though the cost function
usually contains a large number of Pauli strings, many of the coefficients are
quite small. Fig. 1 shows the plot of coefficient ranges with the number of
terms in the Hamiltonian. Further, if we sum the norms of the Pauli terms
within the coefficients range, it can be seen in Fig. 2 that even though the
largest terms were very few, they still are dominant, and in fact, the plot
for the sum of norms of the coefficients is skewed opposite to the
distribution of the number of Pauli-terms. The plot presents the data for
small molecules, but this trend gets stronger for larger molecules. This
observation suggests that for the Hamiltonians plotted in Fig. 1, the terms
with coefficients greater than $0.1$ likely contribute the most towards the
ground state energy of the system. So, in our algorithm, we start with a qubit
Hamiltonian that contains Pauli terms with coefficients above a threshold and
start the optimization routine. As the routine progresses, we gradually add
more terms with smaller coefficients. Algorithm III describes our algorithm
formally.

Figure 1: Coefficient range vs. the number of terms for molecular Hamiltonians
in STO-3G basis at equilibrium geometry. Figure 2: Coefficient range vs. the
sum of coefficients within the coefficients range for molecular Hamiltonians
in STO-3G basis at equilibrium geometry.

[h] Basic Truncation algorithm

1:function trunc_ham($\mathcal{H}_{q}$, cut-off)

2: $\mathcal{H}_{list}\leftarrow\\{\mathcal{H}_{q}\\}$

3: for co in cut-off do

4: $\mathcal{H}_{i}$={$c_{k}P_{k}$} s.t.
($\left|c_{k}\right|P_{k}\in\mathcal{H}$ and $c_{k}\geq co$)

5: Append $\mathcal{H}_{i}$ to $\mathcal{H}_{list}$

6: Reverse $\mathcal{H}_{list}$

7: return $\mathcal{H}_{list}$

8:function Find-Min($\mathcal{H}$, $\ket{\Psi(\hat{\theta})}$)

9: cut-off $\leftarrow$ {$c_{off1},c_{off2},..c_{off(k)}$}

10: $\mathcal{H}_{list}\leftarrow$
$\\{\mathcal{H}_{k},..\mathcal{H}_{1},\mathcal{H}_{q}\\}$ =
trunc_ham($\mathcal{H}$, cut-off)

11: $\hat{\theta}_{0}\leftarrow$ random($len(\hat{\theta})$)

12: for $\mathcal{H}_{i}$ in $\mathcal{H}_{list}$ do

13: E, $\hat{\theta_{i}}$ =
VQE($\mathcal{H}_{i}$,$\ket{\psi(\hat{\theta}_{i-1})}$)

14: if $\mathcal{H}_{i}==\mathcal{H}_{q}$ then

15: $\hat{\theta}_{opti}=\hat{\theta}_{i}$

16: return $\text{E},\hat{\theta}_{i}$

To formalize the above ideas, assume that we have a qubit Hamiltonian, $H_{q}$
given in Eq. (3) for a VQE problem. We can construct a sequence of
Hamiltonians of length $k+1$ derived from $H_{q}$ such that each Hamiltonian
$H_{i}$ is contained within $H_{i-1}$ and $H_{1}$ is contained within $H_{q}$:

$\displaystyle H_{k}\subset H_{k-1}...\subset H_{1}\subset H_{q}$ (6)

During the optimization stage, we can optimize the sequence of the
Hamiltonians one by one, and when transitioning to the next Hamiltonian, we
shall use the previously found best parameters as initial parameters for the
next round. For such a strategy to work, the initial Hamiltonian to be
optimized $H_{k}$ needs to contain the Pauli strings in $H_{q}$ whose absolute
value of the coefficients are relatively large compared to other Pauli strings
in $H_{q}$ so that the optimization starts out in the correct direction
instead of randomly exploring the parameter space. As more terms that are
small enough to be considered perturbations are added, eigenvalues and
eigenvectors are stable Tao and Vu (2010). This is why the ansatz parameters
for the low energy state of the newly perturbed Hamiltonian could be found in
the vicinity of the parameters obtained when the optimization routine
converges for the truncated Hamiltonian.

###  III.1 Naive two-step hard cutoff strategy

In this subsection, we shall consider a naive hard-cutoff truncation
technique. The idea is to choose a hard cutoff $c_{\Lambda}$ and only pick
Pauli strings in the original Hamiltonian with coefficients whose absolute
value is larger than $c_{\Lambda}$ to obtain a truncated Hamiltonian $H_{1}$.
So, we start with $H_{q}$ given in Eq. (3) and obtain the truncated
Hamiltonian as

$\displaystyle H_{1}=\sum_{k=1}^{N}\Theta(|c_{k}|-c_{\Lambda})c_{k}P_{k},$ (7)

where $\Theta$ is the Heaviside step function. The VQE algorithm is started
with this truncated Hamiltonian and is allowed to run for a certain number of
iterations. Then, without terminating the VQE routine, the truncated
Hamiltonian is swapped with the full Hamiltonian, and the VQE routine is
continued until convergence or until the max iterations. This two-step
truncation technique is a subset of Algorithm 1, and it allows us to test the
algorithm and understand performance gains. Our numerical experiments show
that this naive strategy delivers considerable improvement in measurement
count over the standard VQE algorithm while producing similar final results.
It is quite general and does not depend on the choice of the optimizer,
ansatz, or a chemical system. Further, this strategy can easily be updated to
a more complex strategy with multi-step optimization, given in the previous
section, where the truncation thresholds, number of iterations in each step,
and convergence criteria all act as hyperparameters.

In our numerical simulations, we tested various molecular systems that
required up to 16 qubits. For each system, we used the variational quantum
eigensolver (VQE) algorithm, modified according to our strategy, to calculate
the ground energy. We used the Simultaneous Perturbation Stochastic
Approximation (SPSA) Spall (1997) as our optimizer for the VQE algorithm, as
many recent studies have shown that it provides stable and fast performance
Sung  _et al._ (2020); Lavrijsen  _et al._ (2020). It should be noted that our
strategy is not optimizer-dependent and should work with any of the
optimizers. We explored two different ansatzes for our calculations: Two Local
ansatzes and the Unitary Coupled Cluster Singles and Doubles (UCCSD) ansatz
Qiskit contributors (2023); Kandala  _et al._ (2017). Two-local ansatz is a
heuristic ansatz that can take into account the constraints imposed by the
quantum hardware but ignores the physics of the problem, thereby making it
hardware efficient, whereas the UCCSD ansatz is an ansatz that captures the
physics of the system. There have been many recent studies that strike a
balance between Two-Local and UCCSD in the sense that the proposed ansatz
tries to capture the physics and the constraints of the hardware
simultaneously Ryabinkin  _et al._ (2018); Grimsley  _et al._ (2019). Since
our strategy works for the two extremes, we expect it to be compatible with
other ansatz proposed in the literature as well.

Table 1 provides the summary of our numerical simulations and demonstrates
that even with our naive strategy, an improvement of over 50% is possible, and
the improvement $S$ is calculated with the formula:

$\displaystyle S=\sum_{n}k_{n}*i_{n}/(K*I),$ (8)

where $k_{n}$ and $i_{n}$ are the number of terms in the truncated Hamiltonian
being used and the number of iterations, respectively, and $K$ is the total
terms in the original Hamiltonian $H_{q}$ with $I=\sum_{n}i_{n}$ the maximum
iterations used for the standard VQE. For each molecule, we ran 100
simulations to compare the results between our strategy and standard VQE. The
performance of VQE depends on various aspects such as variational ansatz,
choice of optimizer, etc. So, for our study, we focused on how the results
from our strategy compare with the standard VQE. We found that our strategy
performs in line with results from standard VQE while delivering significant
improvement for the measurement counts. The convergence plots and the
distribution of the final minimum eigenvalues for the numerical simulation we
performed are given in the appendix A.

Molecule |  $H_{1}$/$H_{q}$ terms | Improvement
---|---|---
$H_{2}$ | 11/15 | 14%
$H_{4}$ | 31/361 | 48%
$H_{6}$ | 55/1819 | 51%
$BeH_{2}$ | 53/666 | 49%
$H_{2}O$ | 130/1086 | 47%
$LiH$ | 18/631 | 51%
$NH_{3}$ | 149/2941 | 50%
Table 1:  For each of the molecular systems, we provide the improvement in the
operator measurement count using the naive truncation strategy. The
improvement is calculated using Eq. 8. The maximum iteration for each case is
800, and for the naive hard cutoff strategy, the iterations for $H_{1}$ and
$H_{q}$ are both 400. The second column gives the number of terms for each
optimization stage.

###  III.2 Operator classification truncation strategy

In the previous subsection, we see that the truncation strategy can save a
significant amount of computation resources and lead to faster optimization
while achieving the same level of performance when compared to the standard
VQE. When adopting the naive hard cutoff strategy, one has to choose a hard
cutoff for the coefficients of the Pauli strings (0.1 in the case of the
molecular systems discussed in the previous subsection) to generate the
truncated Hamiltonian for the first stage of the optimization. However, such a
cutoff value is not available a priori; if we choose the cutoff to be too
small, we lose the potential benefit of the truncation strategy as we evaluate
more Pauli strings in the first stage, while if we choose the cutoff to be too
large, the parameters after the first stage might not be good initial
parameters for the second stage (this actually can be seen from several
subplots in Fig. 4 in the appendix A where a sudden jump appears when
transitioning to the second stage optimization). As we discussed, one way to
choose the cutoff value is to plot the distribution of the coefficients in the
original Hamiltonian $H_{0}$, from which one can make a reasonable guess. On
the other hand, it would be more desirable to have a strategic approach to
pick the suitable cutoff. In this section, we propose such an approach to
generate the truncated Hamiltonian via operator classification with clear
physics motivation.

(a)

(b)

Figure 3: For each of the molecular Hamiltonians in STO-3G basis we plot (a)
Number of terms vs each fermionic operator kind (b) Sum of norm of
coefficients vs fermionic operator kind.

In Figure 1, we analyzed the distribution of the coefficients of Pauli terms.
Instead, if we analyze the coefficients of fermionic terms in Eq. (1), we
observe that the biggest contribution comes from the number-operator terms
($a_{i}^{\dagger}a_{i}$). This is because the rest of the terms contain an
overlap integral, which, depending on the locality of the orbitals in the
basis set used, could be relatively small. This motivates segmenting the
Hamiltonian into five fermionic operator types: number-operators
($a_{i}^{\dagger}a_{i}$), excitation-operators ($a_{i}^{\dagger}a_{j}$),
number-excitation operators ($a_{i}^{\dagger}a_{k}^{\dagger}a_{l}a_{i}$),
coulomb operators ($a_{i}^{\dagger}a_{j}^{\dagger}a_{j}a_{i}$), and double-
excitation operators ($a_{i}^{\dagger}a_{j}^{\dagger}a_{k}a_{l}$) Setia and
Whitfield (2018). It is observed that the greatest contribution to the
Hamiltonian comes in the following order (Figure 3):

$\displaystyle H_{num}>H_{cou}>H_{exc}>H_{nex}>H_{dex},$

where $H_{num}=\sum^{M}_{i}a_{i}^{\dagger}a_{i}$,
$H_{cou}=\sum^{M}_{ij}a_{i}^{\dagger}a_{j}^{\dagger}a_{j}a_{i}$,
$H_{exc}=\sum^{M}_{ij}a_{i}^{\dagger}a_{j}$,
$H_{nex}=\sum^{M}_{ikl}a_{i}^{\dagger}a_{k}^{\dagger}a_{l}a_{i}$, and
$H_{dex}=\sum^{M}_{ijkl}a_{i}^{\dagger}a_{j}^{\dagger}a_{k}a_{l}$

So, instead of defining a truncation threshold and running the optimization,
we can make the truncation criteria to be the exclusion of certain fermionic
operators from the Hamiltonian. We start the optimization routine with just
the number-operator terms, followed by Coulomb operators, excitation
operators, number excitation operators, and double-excitation operators.
Furthermore, we note that the two leading Hamiltonins in the contribution
sequence, the number operator and Coulomb operator terms when transformed to
qubit operators using Jordan-Wigner transformation, lead to just $Z-$string
qubit operators, indicating all these terms commute and can be measured in a
single measurement, we can then start the optimization routine with both of
these terms included. The details of the algorithms are given in the Algorithm
III.2.

This operator-classification approach not only provides a good way of deciding
the truncation but also provides some physical meaning to how the Hamiltonian
is ‘evolved’ throughout the optimization routine. In the initial phases, the
optimization is performed for a Hamiltonian that only includes the potential
energy and electronic repulsion terms. Then, in the following steps, the
quantum mechanical terms are added. In our testing, we found that this
operator classification technique works better than the standard VQE and the
truncation technique proposed in the previous section. This is because in
first phase of optimization involves only a single measurement.

[h] Physics inspired truncation algorithm

1:function trunc_ham($\mathcal{H}$, cut-off)

2: $\mathcal{H}_{num}$={$h_{pq}a_{p}^{\dagger}a_{q}$} s.t.
($h_{pq}a_{p}^{\dagger}a_{q}$ $\in$ $\mathcal{H}$, p=q)

3: $\mathcal{H}_{cou}$={$h_{pqrs}a_{p}^{\dagger}a_{q}^{\dagger}a_{r}a_{s}$}
s.t. ($h_{pqrs}a_{p}^{\dagger}a_{q}^{\dagger}a_{r}a_{s}$ $\in$ $\mathcal{H}$,
len(set{p,q,r,s})=2)

4: $\mathcal{H}_{exc}$={$h_{pq}a_{p}^{\dagger}a_{q}$} s.t.
($h_{pq}a_{p}^{\dagger}a_{q}$ $\in$ $\mathcal{H}$, p$\neq$q)

5: $\mathcal{H}_{nex}$={$h_{pqrs}a_{p}^{\dagger}a_{q}^{\dagger}a_{r}a_{s}$}
s.t. ($h_{pqrs}a_{p}^{\dagger}a_{q}^{\dagger}a_{r}a_{s}$ $\in$ $\mathcal{H}$,
len(set{p,q,r,s})=3)

6: $\mathcal{H}_{dex}$={$h_{pqrs}a_{p}^{\dagger}a_{q}^{\dagger}a_{r}a_{s}$}
s.t. ($h_{pqrs}a_{p}^{\dagger}a_{q}^{\dagger}a_{r}a_{s}$ $\in$ $\mathcal{H}$,
len(set{p,q,r,s})=4)

7:
$\mathcal{H}_{list}=\\{\mathcal{H}_{num},\mathcal{H}_{cou},\mathcal{H}_{exc},\mathcal{H}_{nex},\mathcal{H}_{dex}\\}$

8: return $\mathcal{H}_{list}$

9:function Find-Min($\mathcal{H}_{q}$, $\ket{\Psi(\hat{\theta})}$)

10: $\mathcal{H}_{segments}\leftarrow$
$\\{\mathcal{H}_{num},\mathcal{H}_{cou},\mathcal{H}_{exc},\mathcal{H}_{nex},\mathcal{H}_{dex}\\}$
= segments($\mathcal{H}_{f}$)

11: $\mathcal{H}_{3}=\mathcal{H}_{num}+\mathcal{H}_{cou}$,

12: $\mathcal{H}_{2}=\mathcal{H}_{3}+\mathcal{H}_{exc}$,

13: $\mathcal{H}_{1}=\mathcal{H}_{2}+\mathcal{H}_{nex}$,

14: $\mathcal{H}_{q}=\mathcal{H}_{1}+\mathcal{H}_{dex}$

15:
$\mathcal{H}_{list}=\\{\mathcal{H}_{3},\mathcal{H}_{2},\mathcal{H}_{1},\mathcal{H}_{q}\\}$

16: $\hat{\theta}_{0}\leftarrow$ random($len(\hat{\theta})$)

17: for $\mathcal{H}_{i}\in\mathcal{H}_{list}$ do

18: E, $\hat{\theta_{i}}$ = VQE(jordan-
wigner$(\mathcal{H}_{i})$,$\ket{\Psi(\hat{\theta}_{i-1})}$)

19: if $\mathcal{H}_{i}==\mathcal{H}_{f}$ then

20: $\hat{\theta}_{opti}=\hat{\theta}_{i}$

21: return $\text{E},\hat{\theta}_{i}$

Molecule |  $H_{3}/H_{2}/H_{1}/H_{q}$ terms | Improvement
---|---|---
$H_{2}$ | 11/ 11/ 11/ 15 | 57%
$H_{4}$ | 37/ 61/ 205/ 361 | 68%
$H_{6}$ | 79/ 139/ 739/ 1819 | 72%
$BeH_{2}$ | 106/122/314/666, | 70%
$H_{2}O$ | 106 /134 /470 /1086 | 71%
$LiH$ | 79/ 103/ 343/ 631 | 69%
$NH_{3}$ | 137/ 201/ 1085/ 2941 | 73%
Table 2:  For each of the molecular systems, we provide the improvement in the
operator measurement count using the operator classification truncation
strategy. The improvement is calculated using Eq. 8. The maximum iteration for
each case is 1000, and the iterations for $H_{3}$, $H_{2}$, $H_{1}$, and
$H_{q}$ are 500, 100, 200 and 200 respectively. The second column gives the
number of terms for each optimization stage.

For our testing, we tested both the hardware-efficient ansatz and UCCSD. Our
proposed strategy worked for both, but we conjecture that our strategy will
work particularly best with UCCSD as the system size grows. Table 2 summarizes
our numerical simulations. As is evident from the data, an improvement of over
70% is achieved for many molecules along with substantial improvement for the
rest. From Eq. (8), notice that the improvement over standard VQE is a
function of the number of terms in each step of our new algorithm and the
total iterations used for that particular step. In particular, we note that
the number of terms in the initial steps is much smaller than the total number
of terms that are used in the final step. Therefore, the improvement is also
constrained by the number of iterations being used in the final step of the
algorithm. As the molecule size increases, we expect the number of iterations
required for convergence also to increase, implying that further improvement
might be possible.

##  IV Conclusion and outlook

One of the main concerns around the variational quantum eigensolver is the
number of quantum circuit measurements required for its implementation. In
this work, we proposed two different strategies for reducing the measurement
counts and we see both strategies deliver substantial performance gains,
indicating a solid step toward addressing the measurement problem in VQE. For
future directions, we reckon there is a way to combine both of our strategies
for further improvement. The improvement for our physics-inspired truncation
strategy is bound by the number of measurements in the last step of the
algorithm. In principle, it should be possible to use Algorithm 1 as a
subroutine for the last step of Algorithm 2. In this work, we considered the
operator classification method for the class of molecular Hamiltonians; it
would also be interesting to explore the possibility of applying such kind of
strategy to other types of Hamiltonians containing perturbative terms.

###### Acknowledgements.

This work is supported by Wellcome Leap as part of the Quantum for Bio
Program.

## References

  * Peruzzo  _et al._ (2014) A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’brien, Nature communications 5, 4213 (2014).
  * Cerezo  _et al._ (2021a) M. Cerezo, A. Arrasmith, R. Babbush, S. C. Benjamin, S. Endo, K. Fujii, J. R. McClean, K. Mitarai, X. Yuan, L. Cincio, and P. J. Coles, Nature Reviews Physics 3, 625 (2021a).
  * Amaro  _et al._ (2022) D. Amaro, C. Modica, M. Rosenkranz, M. Fiorentini, M. Benedetti, and M. Lubasch, Quantum Science and Technology 7, 015021 (2022).
  * Bravo-Prieto  _et al._ (2023) C. Bravo-Prieto, R. LaRose, M. Cerezo, Y. Subasi, L. Cincio, and P. J. Coles, Quantum 7, 1188 (2023).
  * Blunt  _et al._ (2022) N. S. Blunt, J. Camps, O. Crawford, R. Izsák, S. Leontica, A. Mirani, A. E. Moylett, S. A. Scivier, C. Sünderhauf, P. Schopf, J. M. Taylor, and N. Holzmann, Journal of Chemical Theory and Computation 18, 7001–7023 (2022).
  * Landman  _et al._ (2022) J. Landman, N. Mathur, Y. Y. Li, M. Strahm, S. Kazdaghli, A. Prakash, and I. Kerenidis, Quantum 6, 881 (2022).
  * McClean  _et al._ (2016) J. R. McClean, J. Romero, R. Babbush, and A. Aspuru-Guzik, New Journal of Physics 18, 023023 (2016).
  * Kandala  _et al._ (2017) A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta, nature 549, 242 (2017).
  * Cao  _et al._ (2019) Y. Cao, J. Romero, J. P. Olson, M. Degroote, P. D. Johnson, M. Kieferová, I. D. Kivlichan, T. Menke, B. Peropadre, N. P. Sawaya, _et al._ , Chemical reviews 119, 10856 (2019).
  * Bauer  _et al._ (2020) B. Bauer, S. Bravyi, M. Motta, and G. K.-L. Chan, Chemical Reviews 120, 12685 (2020).
  * Lee  _et al._ (2021) J. Lee, D. W. Berry, C. Gidney, W. J. Huggins, J. R. McClean, N. Wiebe, and R. Babbush, PRX Quantum 2, 030305 (2021).
  * Bravo-Prieto  _et al._ (2020) C. Bravo-Prieto, J. Lumbreras-Zarapico, L. Tagliacozzo, and J. I. Latorre, Quantum 4, 272 (2020).
  * Uvarov  _et al._ (2020) A. Uvarov, J. D. Biamonte, and D. Yudin, Physical Review B 102, 075104 (2020).
  * Cao  _et al._ (2018) Y. Cao, J. Romero, and A. Aspuru-Guzik, IBM Journal of Research and Development 62, 6 (2018).
  * Banuls  _et al._ (2020) M. C. Banuls, R. Blatt, J. Catani, A. Celi, J. I. Cirac, M. Dalmonte, L. Fallani, K. Jansen, M. Lewenstein, S. Montangero, _et al._ , The European physical journal D 74, 1 (2020).
  * Liu  _et al._ (2022) X. Liu, A. Angone, R. Shaydulin, I. Safro, Y. Alexeev, and L. Cincio, IEEE Transactions on Quantum Engineering 3, 1 (2022).
  * Jordan and Wigner (1993) P. Jordan and E. P. Wigner,  _Über das paulische äquivalenzverbot_ (Springer, 1993).
  * Somma (2005) R. D. Somma, “Quantum computation, complexity, and many-body physics,” (2005), arXiv:quant-ph/0512209 [quant-ph] .
  * O’Malley  _et al._ (2016) P. J. J. O’Malley, R. Babbush, I. D. Kivlichan, J. Romero, J. R. McClean, R. Barends, J. Kelly, P. Roushan, A. Tranter, N. Ding, B. Campbell, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, A. G. Fowler, E. Jeffrey, E. Lucero, A. Megrant, J. Y. Mutus, M. Neeley, C. Neill, C. Quintana, D. Sank, A. Vainsencher, J. Wenner, T. C. White, P. V. Coveney, P. J. Love, H. Neven, A. Aspuru-Guzik, and J. M. Martinis, Phys. Rev. X 6, 031007 (2016).
  * Wecker  _et al._ (2015) D. Wecker, M. B. Hastings, and M. Troyer, Phys. Rev. A 92, 042303 (2015).
  * Cerezo  _et al._ (2021b) M. Cerezo, A. Arrasmith, R. Babbush, S. C. Benjamin, S. Endo, K. Fujii, J. R. McClean, K. Mitarai, X. Yuan, L. Cincio, _et al._ , Nature Reviews Physics 3, 625 (2021b).
  * Tilly  _et al._ (2022) J. Tilly, H. Chen, S. Cao, D. Picozzi, K. Setia, Y. Li, E. Grant, L. Wossnig, I. Rungger, G. H. Booth, _et al._ , Physics Reports 986, 1 (2022).
  * Grimsley  _et al._ (2019) H. R. Grimsley, S. E. Economou, E. Barnes, and N. J. Mayhall, Nature communications 10, 3007 (2019).
  * Mehendale  _et al._ (2023) S. G. Mehendale, B. Peng, N. Govind, and Y. Alexeev, The Journal of Physical Chemistry A 127, 4526–4537 (2023).
  * Arrasmith  _et al._ (2020) A. Arrasmith, L. Cincio, R. D. Somma, and P. J. Coles, “Operator Sampling for Shot-frugal Optimization in Variational Algorithms,”  (2020), arXiv:2004.06252 [quant-ph].
  * Lavrijsen  _et al._ (2020) W. Lavrijsen, A. Tudor, J. Müller, C. Iancu, and W. De Jong, in  _2020 IEEE international conference on quantum computing and engineering (QCE)_ (IEEE, 2020) pp. 267–277.
  * Cerezo  _et al._ (2021c) M. Cerezo, A. Sone, T. Volkoff, L. Cincio, and P. J. Coles, Nature Communications 12 (2021c), 10.1038/s41467-021-21728-w.
  * Arrasmith  _et al._ (2021) A. Arrasmith, M. Cerezo, P. Czarnik, L. Cincio, and P. J. Coles, Quantum 5, 558 (2021).
  * Lee  _et al._ (2018) J. Lee, W. J. Huggins, M. Head-Gordon, and K. B. Whaley, Journal of chemical theory and computation 15, 311 (2018).
  * Ryabinkin  _et al._ (2018) I. G. Ryabinkin, T.-C. Yen, S. N. Genin, and A. F. Izmaylov, Journal of chemical theory and computation 14, 6317 (2018).
  * Lang  _et al._ (2020) R. A. Lang, I. G. Ryabinkin, and A. F. Izmaylov, Journal of Chemical Theory and Computation 17, 66 (2020).
  * Ryabinkin  _et al._ (2020) I. G. Ryabinkin, R. A. Lang, S. N. Genin, and A. F. Izmaylov, Journal of chemical theory and computation 16, 1055 (2020).
  * Tang  _et al._ (2021) H. L. Tang, V. Shkolnikov, G. S. Barron, H. R. Grimsley, N. J. Mayhall, E. Barnes, and S. E. Economou, PRX Quantum 2, 020310 (2021).
  * Huang  _et al._ (2022) Y. Huang, Q. Li, X. Hou, R. Wu, M.-H. Yung, A. Bayat, and X. Wang, Physical Review A 105, 052414 (2022).
  * Spall (1992) J. C. Spall, IEEE transactions on automatic control 37, 332 (1992).
  * Nakanishi  _et al._ (2020) K. M. Nakanishi, K. Fujii, and S. Todo, Physical Review Research 2, 043158 (2020).
  * Sung  _et al._ (2020) K. J. Sung, J. Yao, M. P. Harrigan, N. C. Rubin, Z. Jiang, L. Lin, R. Babbush, and J. R. McClean, Quantum Science and Technology 5, 044008 (2020).
  * Koczor and Benjamin (2022) B. Koczor and S. C. Benjamin, Physical Review Research 4, 023017 (2022).
  * Kübler  _et al._ (2020) J. M. Kübler, A. Arrasmith, L. Cincio, and P. J. Coles, Quantum 4, 263 (2020).
  * Sim  _et al._ (2021) S. Sim, J. Romero, J. F. Gonthier, and A. A. Kunitsa, Quantum Science and Technology 6, 025019 (2021).
  * Nannicini (2019) G. Nannicini, Physical Review E 99, 013304 (2019).
  * Pellow-Jarman  _et al._ (2021) A. Pellow-Jarman, I. Sinayskiy, A. Pillay, and F. Petruccione, Quantum Information Processing 20, 202 (2021).
  * Bonet-Monroig  _et al._ (2023) X. Bonet-Monroig, H. Wang, D. Vermetten, B. Senjean, C. Moussa, T. Bäck, V. Dunjko, and T. E. O’Brien, Physical Review A 107, 032407 (2023).
  * Bravyi and Kitaev (2002) S. B. Bravyi and A. Y. Kitaev, Annals of Physics 298, 210–226 (2002).
  * Setia  _et al._ (2019) K. Setia, S. Bravyi, A. Mezzacapo, and J. D. Whitfield, Phys. Rev. Res. 1, 033033 (2019).
  * Tao and Vu (2010) T. Tao and V. Vu, “Random matrices: Universality of local eigenvalue statistics,” (2010), arXiv:0906.0510 [math.PR] .
  * Spall (1997) J. Spall, in  _Proceedings of the 36th IEEE Conference on Decision and Control_, Vol. 2 (1997) pp. 1417–1424 vol.2.
  * Qiskit contributors (2023) Qiskit contributors, “Qiskit: An open-source framework for quantum computing,”  (2023).
  * Setia and Whitfield (2018) K. Setia and J. D. Whitfield, The Journal of Chemical Physics 148 (2018), 10.1063/1.5019371.
  * Sun  _et al._ (2020) Q. Sun, X. Zhang, S. Banerjee, P. Bao, M. Barbry, N. S. Blunt, N. A. Bogdanov, G. H. Booth, J. Chen, Z.-H. Cui, J. J. Eriksen, Y. Gao, S. Guo, J. Hermann, M. R. Hermes, K. Koh, P. Koval, S. Lehtola, Z. Li, J. Liu, N. Mardirossian, J. D. McClain, M. Motta, B. Mussard, H. Q. Pham, A. Pulkin, W. Purwanto, P. J. Robinson, E. Ronca, E. R. Sayfutyarova, M. Scheurer, H. F. Schurkus, J. E. T. Smith, C. Sun, S.-N. Sun, S. Upadhyay, L. K. Wagner, X. Wang, A. White, J. D. Whitfield, M. J. Williamson, S. Wouters, J. Yang, J. M. Yu, T. Zhu, T. C. Berkelbach, S. Sharma, A. Y. Sokolov, and G. K.-L. Chan, The Journal of Chemical Physics 153 (2020), 10.1063/5.0006074.

##  Appendix A Convergence plots for naive and the operator classification
truncation strategy

In this section, we provide the details of our numerical simulation setup. We
explored seven molecular systems, $H_{2}$, $H_{4}$, $H_{6}$, $H_{2}O$, $LiH$,
$H_{2}O$, $BeH_{2}$. We used qiskit Qiskit contributors (2023) along with the
pyscf Sun  _et al._ (2020) driver to run the Hartree-Fock method in STO-3G
basis and at equilibrium geometries to obtain the one-body and two-body terms.
The geometries are given in the following table:

Molecular system | Geometry ($\AA$)
---|---
$H_{2}$ | H 0 0 0; H 0 0 0.735
$H_{4}$ | H 0 0 0; H 0 0 0.735; H 0 0 1.535; H 0 0 2.135
$H_{6}$ | ”H 0 0 0; H 0 0 0.735; H 0 0 1.535; H 0 0 2.135; H 0 0 2.835; H 0 0 3.57
$LiH$ | Li .0 .0 .0; H .0 .0 1.5949
$NH_{3}$ | N1 0.0000 0.0000 0.0000;H2 0.0000 -0.9377 -0.3816;H3 0.8121 0.4689 -0.3816;H4 -0.8121 0.4689 -0.3816
$BeH_{2}$ | Be1 0.0000 0.0000 0.0000; H2 0.0000 0.0000 1.3264; H3 0.0000 0.0000 -1.3264
$H_{2}O$ | O1 0.0000 0.0000 0.1173; H2 0.0000 0.7572 -0.4692; H3 0.0000 -0.7572 -0.4692
Table 3: Geometries for various molecular systems. The units for the
coordinates are in angstrom.

For all our numerical simulations, we used the Jordan-Wigner transformation to
map the fermionic Hamiltonian to the qubit operator Hamiltonian. For
simulating the quantum circuits, we use the AerEstimator with statevector
simulator available in qiskit.

In Fig. 4 we compare the averaged convergence plots for the naive two-step
truncation strategy and the standard VQE for all seven molecules, and we show
the comparison of the final minimum eigenvalues found by both strategies in
the inset box plot. For both strategies, we ran 100 simulations and took the
average for the plots. For simulations in subplots (a) to (g) in Fig. 4, we
used the TwoLocal() ansatz with $R_{y}$ rotation blocks and $CZ$ entangling
gates, and the number of repetitions is three; for simulations in subplots (h)
and (i), we used the UCCSD ansatz. The classical optimizer used is SPSA, with
the maximum iterations set to 800 for the standard VQE and 400 for both
Hamiltonians in the naive two-step truncation strategy. We note that the SPSA
optimizer has a default initial calibration phase with 50 iterations, we see
that there are in total 850 evaluations in the plots. We also note that our
convention is slightly different from the hyperparameters ‘iteration’ defined
in the SPSA optimizer; one iteration in our convention corresponds to one
measurement while one iteration in SPSA corresponds to two evaluations
(measurements) of the cost function. From the subplots, we see that the naive
two-step method achieves similar convergence speed with a similar distribution
of the final minimum eigenvalues when compared with the standard VQE. some
kinks can be observed in the subplots because the energy to the subsequent
truncated Hamiltonian is higher (or lower), indicating the importance of a
proper choice of cutoff value for the original Hamiltonian.

In Fig. 5 we compare the averaged convergence plots for the operator
classification truncation strategy and the standard VQE for all seven
molecules with subplots (a) to (g) using the TwoLocal() ansatz with the same
settings as in Fig. 4, and subplot (h) and (i) using the UCCSD ansatz. All
subplots are obtained by averaging 100 simulations, except for subplot (h)
which is obtained by averaging 10 simulations due to its much longer clock run
time. The classical optimizer used is SPSA, with the maximum iterations set to
1000 for the standard VQE, and 500, 100, 200, and 200 iterations for $H_{3}$,
$H_{2}$, $H_{1}$, and $H_{q}$ respectively. It is of interest to note that the
kinks are mostly absent in the operator classification truncation method,
indicating that this is a generic algorithm suitable for the study of quantum
chemistry VQE problems.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

Figure 4: Convergence comparison plots between the naive two-step strategy and
the standard VQE for different molecules. The inset box plots display the
distribution for the final minimum eigenvalues. All subplots are obtained by
averaging 100 runs of simulations.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

Figure 5: Convergence comparison plots between the operator classification
truncation strategy and the standard VQE for different molecules. The inset
box plots display the distribution for the final minimum eigenvalues. All
subplots are obtained by averaging 100 runs of simulations except for subplot
(h), which is obtained by averaging 10 runs of simulations.