# Extension of the Trotterized Unitary Coupled Cluster to Triple Excitations

Mohammad Haidar  Sorbonne Université, Laboratoire de Chimie
Théorique(UMR-7616-CNRS), 4 place Jussieu-75005 Paris, France  Marko J. Rančić
TotalEnergies, Tour Coupole La Défense, 2 Pl. Jean Millier, 92078 Paris,
France  Yvon Maday  Sorbonne Université, CNRS, Université Paris Cité,
Laboratoire Jacques Louis Lions (LJLL), 4 place Jussieu-75005 Paris, France
Jean-Philip Piquemal  Sorbonne Université, Laboratoire de Chimie
Théorique(UMR-7616-CNRS), 4 place Jussieu-75005 Paris, France
Mohammadhaidar2016<EMAIL_ADDRESS>
###### Abstract

The Trotterized Unitary Coupled Cluster Single and Double (UCCSD) ansatz has
recently attracted interest due to its use in Variation Quantum Eigensolver
(VQE) molecular simulations on quantum computers. However, when the size of
molecules increases, UCCSD becomes less interesting as it cannot achieve
sufficient accuracy. Including higher-order excitations is therefore mandatory
to recover the UCC’s missing correlation effects. In this Letter, we extend
the Trotterized UCC approach via the addition of (true) Triple T excitations
introducing UCCSDT. We also include both spin and orbital symmetries. Indeed,
in practice, these later help to reduce unnecessarily circuit excitations and
thus accelerate the optimization process enabling to tackle larger molecules.
Our initial numerical tests (12-14 qubits) show that UCCSDT improves the
overall accuracy by at least two-orders of magnitudes with respect to standard
UCCSD. Overall, the UCCSDT ansatz is shown to reach chemical accuracy and to
be competitive with the CCSD(T) gold-standard classical method of quantum
chemistry.

###### keywords:

American Chemical Society, LaTeX

Sorbonne Université, CNRS, Université Paris Cité, Laboratoire Jacques Louis
Lions (LJLL), 4 place Jussieu-75005 Paris, France TotalEnergies, Tour Coupole
La Défense, 2 Pl. Jean Millier, 92078 Paris, France Institut Universitaire de
France, Paris, France IR,NMR,UV

##  1 Introduction

The use of quantum computers is a promising strategy to overcome challenges
occuring in the area of quantum chemistry.1, 2, 3, 4, 5. Such hardware
potentially facilitate the encoding of the full configuration interaction
(FCI) wavefunction of the many-electron molecular system thanks to entangled
quantum bits or qubits. In classical computers, the FCI includes a number of
determinants which scales exponentially with the number of electrons denoted
$n_{e}$, as roughly $O(n_{o}^{n_{e}})$, where $n_{o}$ is the number of spin-
orbitals, which makes the manipulation and storage of the wavefunction
inefficient6, 7, 8. However, using a quantum computer, we can instead store
the (FCI) wavefunction by using only $n_{o}$ orbitals which corresponds to
$n_{o}$ qubits1. This potential quantum advantage has recently excited both
hardware and software communities generating rapid progresses in the field.
Several quantum algorithms have been developed, and among them, the
Variational Quantum Eigensolver (VQE)9, 10, 11, 12 appears well suited for its
practical implementation on present Noisy Intermediate Scaled Quantum (NISQ)
devices 13, 14. VQE simulations have been performed numerically for molecules
using various virtual noiseless simulators15, 16, 17, 18 and have been tested
experimentally on actual NISQ devices9, 19, 20, 21, 22. An important component
in VQE is the parametrized ansatz, which represents the trial wavefunction and
is implemented as a quantum circuit, composed of unitary quantum gates, that
measures the energy expectation value of the trial wavefunction and then its
parameters are updated in a classical optimization loop. The Unitary Coupled-
Cluster (UCC) ansatz was used when VQE was initially proposed, and has granted
lots of attention since then. Indeed, UCC includes several attractive
features: (i) it is chemistry-inspired since it is a unitary version of the
classical coupled-cluster method, which is among the most accurate classical
quantum chemistry methods for many-body simulation (see exhaustive recent
reviews in 23, 24); (ii) A variant of the UCC theory including single (S) and
double (D) excitations called UCCSD is well suited for a direct use on the
current NISQ circuits (see comprehensive reviews of UCCSD in 25, 26).
Moreover, the Suzuki-Trotter approximation27 allows separating out the ansatz
into a product form of individual exponentials. This version of the ansatz is
then denoted Trotterized UCC. Since it is unitary, it can be easily
transformed into well-defined unitary gates and therefore its implementation
is possible on NISQ devices.

The Trotterized UCCSD ansatz is well documented 25, 26, 28, 29, 30 and has
been successfully tested for small molecules31, 32, 33. However, UCCSD
exhibits a limited accuracy in strongly correlated systems similarly as the
classical CCSD ansatz29, 34. Moreover, when the size of the molecules
increases, many unnecessary excitations can appear in UCCSD which do not
improve the chemical accuracy while adding extra-noise due to the extra
circuit-depth.35 This reveals the need for more accurate and compact
representations of the wavefunction. In that context, various alternative
versions of the UCC ansatz, also truncated at the double excitations level
have been proposed by the community and have brought interesting upgrades30,
36, 37, 38, 39. A possible direction also is to use the contracted quantum
eigensolvers (CQE) methods which provides an ansatz that can generalize CCSD
towards FCI40, 41. However, a possible direction towards improving accuracy is
to extend the Trotterized UCCSD through adding three-body clusters that
represent the (true) Triple T excitations. For example, the gold-standard
CCSD(T) method42 of (classical) quantum chemistry involves perturbative triple
excitations that play a key role in the inclusion of correlation effects 43.
The objective of this Letter is to introduce the Trotterized UCCSDT approach
and to analyze the behaviour of triple excitations on a set of molecules
compared to the initial UCCSD. Of course, an expected downside of the UCCSDT
approach is that it involves a large number of parameters to optimize and
therefore leads to generate circuits that are excessively long and therefore
not suited for their practical use on current NISQ devices. A possible way to
cope with this potential inefficiency of UCCSDT (i.e allowing to reduce the
required number of parameters) is to use the power of the total spin and point
group symmetries. Indeed, taking into account these symmetries enable to use
spin factorization and orbital symmetry techniques23 that have been shown to
enable simplifications in the classical coupled cluster equations through
eliminating unnecessarily single and double excitations. Very recently, the
point group symmetry constraint has been implemented with the UCCSD ansatz26,
18. The authors demonstrated a significant reduction of the number of
parameters without affecting the accuracy.
This article is organized as follows: the next section is devoted to a
description of the theory of Unitary Coupled Cluster method for electronic
structure calculations. It includes the detailed formalism of triple
excitations. It is given in its simplified form after application of both spin
and orbital symmetries. The subsequent section presents the computational
materials used to describe the UCCSDT-VQE approach and the last section
illustrates the application of UCCSDT-VQE by showing numerical results through
testing several molecules. Comparisons were made with UCCSD as well as the
classical methods such as CCSD, CCSD(T), CCSDT-full and FCI.

##  2 Theory

In this work we have based ourselves on the Born-Oppenheimer clamped-nuclei
Hamiltonian formalism:

$\hat{H}=\sum_{p,q}h_{pq}\hat{\textrm{a}}_{p}^{\dagger}\hat{\textrm{a}}_{q}+\frac{1}{2}\sum_{p,q,r,s}h_{pqrs}\hat{\textrm{a}}_{p}^{\dagger}\hat{\textrm{a}}_{q}^{\dagger}\hat{\textrm{a}}_{r}\hat{\textrm{a}}_{s}$
(1)

where $a_{p}^{\dagger}$ ($a_{q}$) are anti-commuting operators that create
(annihilate) electrons in molecular spin-orbital $p(q)$, respectively. The
symbols $h_{pq}$ and $h_{pqrs}$ denote the one- and two-body integrals of the
corresponding operators and spin-orbitals in Dirac notation, respectively.
These integrals can be easily computed on a classical computer.
The main goal is to find the ground state energy $E_{0}$ of Hamiltonian
$\hat{H}$ given in eq. (1). The Variational Quantum Eigensolver (VQE), a
hybrid quantum classical algorithm was designed to solve this problem (see
Figure 1 in 9). The principal of VQE is that the quantum computer prepares and
measures, for any given parameter $\vec{\theta}$, the parametrized quantum
state $|\Psi(\vec{\theta})\rangle$ while the classical computer is used for
optimization to propose and update $\vec{\theta}$, in order to minimize the
variational energy $\langle\Psi(\vec{\theta})|H|\Psi(\vec{\theta})\rangle\geq
E_{0}$.
We focus on the Unitary Coupled Cluster (UCC) wavefunction, which, as
mentioned in the introduction, is chemically inspired from the classical
Coupled Cluster (CC) method, that has yielded high-accuracy results in
classical quantum chemistry. The UCC theory is important in our study, since
its wavefunction is unitary and therefore can be efficiently implemented on
the quantum computer unlike the regular CC method. The UCC wavefunction can be
expressed as follows:

$|\Psi(\vec{\theta})\rangle=e^{\hat{T}(\vec{\theta})-\hat{T}^{\dagger}(\vec{\theta})}|\psi_{\mathrm{HF}}\rangle$
(2)

where $\hat{T}$ is the so-called cluster excitation operator and
$\hat{T}^{\dagger}$ is its Hermitian conjugate. $|\psi_{\mathrm{HF}}\rangle$
is the reference Hartree-Fock state. Since we are interested in investigating
the effect of Triple excitations on UCC’s accuracy, the coupled-cluster
excitations $\displaystyle\hat{T}(\vec{\theta})$ are truncated to single,
double and triple excitations to introduce UCCSDT,

$\hat{T}(\vec{\theta})=\hat{T_{1}}(\vec{\theta})+\hat{T_{2}}(\vec{\theta})+\hat{T_{3}}(\vec{\theta})=\sum_{a,i}\theta_{i}^{a}\,\hat{\textrm{a}}_{a}^{\dagger}\hat{\textrm{a}}_{i}+\sum_{a,b,i,j}\theta_{i,j}^{a,b}\,\hat{\textrm{a}}_{a}^{\dagger}\hat{\textrm{a}}_{b}^{\dagger}\hat{\textrm{a}}_{i}\hat{\textrm{a}}_{j}+\sum_{a,b,c,i,j,k}\theta_{i,j,k}^{a,b,c}\,\hat{\textrm{a}}_{a}^{\dagger}\hat{\textrm{a}}_{b}^{\dagger}\hat{\textrm{a}}_{c}^{\dagger}\hat{\textrm{a}}_{i}\hat{\textrm{a}}_{j}\hat{\textrm{a}}_{k}$
(3)

where $a,b,c\in virt,i,j,k\in occ$, $\theta_{i}^{a}$, $\theta_{i,j}^{a,b}$ and
$\theta_{i,j,k}^{a,b,c}$ are the variational parameters. The triple excitation
term $\hat{T_{3}}(\vec{\theta})$ (given in equation (3)) increases drastically
the number of parameters as $n_{o}$ increases. This makes the optimization
process very difficult and slow to reach convergence. However, as discussed in
the introduction, the many excitations appearing in the UCCSDT wavefunction
(single, double and triple) are not all important and will not affect all the
correlation energy. Therefore, determining, a priori, those unimportant
amplitudes, and thus not including them in equation (3) above before
compilation and execution of a variational optimization of the parameters, is
highly significant to save resources. In this work we focus only on neutral
closed-shell molecular species allow to make the electronic calculations on
such molecules very efficient by exploiting simultaneously both spin and
orbital symmetries in the UCCSDT ansatz. Indeed, since in the neutral closed-
shell molecules exhibit an even number of electrons, one can maintain only
excitations that can satisfy spin symmetry through keeping the balance between
spin-up($\uparrow$) and spin-down($\downarrow$) electrons. Therefore it is
straightforward to include the spin symmetry constraint in equation (3). The
inclusion of orbital symmetry constraints requires more precise explanation
for its implementation for each molecule subject to a non-trivial point group
symmetry 44, 45. For example the H2 and N2 molecules belongs to the $D_{\infty
h}$ point group whereas LiH and H2O belong to $C_{\infty v}$. Each such group
is characterized by its irreducible representations (irreps) which are listed
in 46. Let us now describe how to impose the orbital symmetry constraint from
each irreps group into equation (3) following the rules that are given in
chapter 2 of reference 23, and that were recently reviewed in 26. First, we
use the orbital symmetry operator $\hat{s}_{e}$, noted in 26, which maps

  * •

each orbital to the irreps of that orbital

  * •

each pair of orbitals to the tensor product of the two irreps of the two
orbitals

  * •

each triple of orbitals to the tensor product of the three irreps of the three
orbitals

Assuming e.g. that the $D_{2h}$ group is used, and taking for example the
triple excitation case, if the irreps of $i,j,k$ are
$i=A_{u},j=B_{1g},k=B_{2u}$ , then we have : $\hat{s}_{e}(i)=A_{u}$,
$\hat{s}_{e}(j)=B_{1g}$, $\hat{s}_{e}(k)=B_{2u}$ and
$\hat{s}_{e}(i,j,k)=B_{3g}$ according to the multiplication table for
$D_{2h}$. Then, the complexity reduction comes from the rule that single
excitations must only occur between orbitals that belong to the same irreps.
Double and Triple excitations must only occur between occupied and virtual
orbitals whose tensor product irreps are similar. Furthermore, the products of
irreps that are partaken in each excitation must contain the fully symmetric
irreps.

Merging these constraints with the spin symmetry, we obtain the operator for
Singles (S)

$\hat{T}_{1}(\vec{\theta})-\hat{T}^{\dagger}_{1}(\vec{\theta})=\sum_{\begin{subarray}{c}i,a\in\\{\uparrow\\}\\\
\hat{s}_{e}(i)=\hat{s}_{e}(a)\end{subarray}}\theta_{i}^{a}(\hat{\textrm{a}}_{a}^{\dagger}\hat{\textrm{a}}_{i}-\hat{\textrm{a}}_{i}^{\dagger}\hat{\textrm{a}}_{a})+\sum_{\begin{subarray}{c}\bar{i},\bar{a}\in\\{\downarrow\\}\\\
\hat{s}_{e}(\bar{i})=\hat{s}_{e}(\bar{a})\end{subarray}}\theta_{\bar{i}}^{\bar{a}}(\hat{\textrm{a}}_{\bar{a}}^{\dagger}\hat{\textrm{a}}_{\bar{i}}-\hat{\textrm{a}}_{\bar{i}}^{\dagger}\hat{\textrm{a}}_{\bar{a}})$
(4)

For Doubles (D),

$\begin{split}\hat{T}_{2}(\vec{\theta})-\hat{T}^{\dagger}_{2}(\vec{\theta})&=\sum_{\begin{subarray}{c}i<j,a<b\in\\{\uparrow\\}\\\
\hat{s}_{e}(i,j)=\hat{s}_{e}(a,b)\end{subarray}}\theta_{i,j}^{a,b}(\hat{\textrm{a}}_{a}^{\dagger}\hat{\textrm{a}}_{b}^{\dagger}\hat{\textrm{a}}_{j}\hat{\textrm{a}}_{i}-\hat{\textrm{a}}_{i}^{\dagger}\hat{\textrm{a}}_{j}^{\dagger}\hat{\textrm{a}}_{b}\hat{\textrm{a}}_{a})+\sum_{\begin{subarray}{c}\bar{i}<\bar{j},\bar{a}<\bar{b}\in\\{\downarrow\\}\\\
\hat{s}_{e}(\bar{i},\bar{j})=\hat{s}_{e}(\bar{a},\bar{b})\end{subarray}}\theta_{\bar{i},\bar{j}}^{\bar{a},\bar{b}}(\hat{\textrm{a}}_{\bar{a}}^{\dagger}\hat{\textrm{a}}_{\bar{b}}^{\dagger}\hat{\textrm{a}}_{\bar{j}}\hat{\textrm{a}}_{\bar{i}}-\hat{\textrm{a}}_{\bar{i}}^{\dagger}\hat{\textrm{a}}_{\bar{j}}^{\dagger}\hat{\textrm{a}}_{\bar{b}}\hat{\textrm{a}}_{\bar{a}})\\\
&+\sum_{\begin{subarray}{c}i,\bar{j},a,\bar{b}\\\
i,a\in\\{\uparrow\\};\bar{j},\bar{b}\in\\{\downarrow\\}\\\
\hat{s}_{e}(i,\bar{j})=\hat{s}_{e}(a,\bar{b})\end{subarray}}\theta_{i,\bar{j}}^{a,\bar{b}}(\hat{\textrm{a}}_{a}^{\dagger}\hat{\textrm{a}}_{\bar{b}}^{\dagger}\hat{\textrm{a}}_{\bar{j}}\hat{\textrm{a}}_{i}-\hat{\textrm{a}}_{i}^{\dagger}\hat{\textrm{a}}_{\bar{j}}^{\dagger}\hat{\textrm{a}}_{\bar{b}}\hat{\textrm{a}}_{a})\end{split}$
(5)

For Triples (T),

$\begin{split}\hat{T}_{3}(\vec{\theta})-\hat{T}^{\dagger}_{3}(\vec{\theta})&=\sum_{\begin{subarray}{c}i<j<k\in\\{\uparrow\\}\\\
a<b<c\in\\{\uparrow\\}\\\
\hat{s}_{e}(i,j,k)=\hat{s}_{e}(a,b,c)\end{subarray}}\theta_{i,j,k}^{a,b,c}(\hat{\textrm{a}}_{a}^{\dagger}\hat{\textrm{a}}_{b}^{\dagger}\hat{\textrm{a}}_{c}^{\dagger}\hat{\textrm{a}}_{k}\hat{\textrm{a}}_{j}\hat{\textrm{a}}_{i}-\hat{\textrm{a}}_{i}^{\dagger}\hat{\textrm{a}}_{j}^{\dagger}\hat{\textrm{a}}_{k}^{\dagger}\hat{\textrm{a}}_{c}\hat{\textrm{a}}_{b}\hat{\textrm{a}}_{a})+\sum_{\begin{subarray}{c}\bar{i}<\bar{j}<\bar{k}\in\\{\downarrow\\}\\\
\bar{a}<\bar{b}<\bar{c}\in\\{\downarrow\\}\\\
\hat{s}_{e}(\bar{i},\bar{j},\bar{k})=\hat{s}_{e}(\bar{a},\bar{b},\bar{c})\end{subarray}}\theta_{\bar{i},\bar{j},\bar{k}}^{\bar{a},\bar{b},\bar{c}}(\hat{\textrm{a}}_{\bar{a}}^{\dagger}\hat{\textrm{a}}_{\bar{b}}^{\dagger}\hat{\textrm{a}}_{\bar{c}}^{\dagger}\hat{\textrm{a}}_{\bar{k}}\hat{\textrm{a}}_{\bar{j}}\hat{\textrm{a}}_{\bar{i}}-\hat{\textrm{a}}_{\bar{i}}^{\dagger}\hat{\textrm{a}}_{\bar{j}}^{\dagger}\hat{\textrm{a}}_{\bar{k}}^{\dagger}\hat{\textrm{a}}_{\bar{c}}\hat{\textrm{a}}_{\bar{b}}\hat{\textrm{a}}_{\bar{a}})\\\
&+\sum_{\begin{subarray}{c}i<j,a<b\in\\{\uparrow\\}\\\
\bar{k},\bar{c}\in\\{\downarrow\\}\\\
\hat{s}_{e}(i,j,\bar{k})=\hat{s}_{e}(a,b,\bar{c})\end{subarray}}\theta_{i,j,\bar{k}}^{a,b,\bar{c}}(\hat{\textrm{a}}_{a}^{\dagger}\hat{\textrm{a}}_{b}^{\dagger}\hat{\textrm{a}}_{\bar{c}}^{\dagger}\hat{\textrm{a}}_{\bar{k}}\hat{\textrm{a}}_{j}\hat{\textrm{a}}_{i}-\hat{\textrm{a}}_{i}^{\dagger}\hat{\textrm{a}}_{j}^{\dagger}\hat{\textrm{a}}_{\bar{k}}^{\dagger}\hat{\textrm{a}}_{\bar{c}}\hat{\textrm{a}}_{b}\hat{\textrm{a}}_{a})+\sum_{\begin{subarray}{c}\bar{i}<\bar{j},\bar{a}<\bar{b}\in\\{\downarrow\\}\\\
k,c\in\\{\uparrow\\}\\\
\hat{s}_{e}(i,j,\bar{k})=\hat{s}_{e}(a,b,\bar{c})\end{subarray}}\theta_{\bar{i},\bar{j},k}^{\bar{a},\bar{b},c}(\hat{\textrm{a}}_{\bar{a}}^{\dagger}\hat{\textrm{a}}_{\bar{b}}^{\dagger}\hat{\textrm{a}}_{c}^{\dagger}\hat{\textrm{a}}_{k}\hat{\textrm{a}}_{\bar{j}}\hat{\textrm{a}}_{\bar{i}}-\hat{\textrm{a}}_{\bar{i}}^{\dagger}\hat{\textrm{a}}_{\bar{j}}^{\dagger}\hat{\textrm{a}}_{k}^{\dagger}\hat{\textrm{a}}_{c}\hat{\textrm{a}}_{\bar{b}}\hat{\textrm{a}}_{\bar{a}})\end{split}$
(6)

With these excitations we construct the symmetric UCCSDT, which is named
throughout the text as sym-UCCSDT.

##  3 Computational Procedure

The implementation of an exponential operator
$e^{\hat{T}(\vec{\theta})-\hat{T}^{\dagger}(\vec{\theta})}$ given in eq. (2)
cannot be done directly on quantum hardware. Therefore applying Trotterization
using the Suzuki–Trotter expansion27 is mandatory to break up the exponential
of a sum as a product of individual exponentials. Thus we used this approach
to write the Trotterized
$e^{\hat{T}(\vec{\theta})-\hat{T}^{\dagger}(\vec{\theta})}$ with
$T(\theta)-T(\theta)^{\dagger}$ (eq. (3)) in the form:

$e^{\hat{T}(\vec{\theta})-\hat{T}^{\dagger}(\vec{\theta})}=\bigg{(}\prod_{\rho}e^{\frac{\theta_{\rho}}{t}(T_{\rho}-T_{\rho}^{\dagger})}\bigg{)}^{t}+\mathcal{O}\,(\frac{1}{t})$
(7)

where $t$ is the number of Trotter steps and $\rho$ corresponds to the
elements of excitations introduced in equation (3). The essential tools needed
to perform the VQE process (such as fermionic second quantization of $\hat{H}$
equation (1), fermion-qubit transforms using Jordan-Wigner representation47,
Trotterization, etc) are implemented in the myQLM-fermion package48, 49. We
use the UCC family module from our OpenVQE package50, 49 to implement UCCSDT
excitations (equations (4) – (6)). We also used PySCF package51 to evaluate
the molecular orbital integrals $h_{pq}$ and $h_{pqrs}$ which are given in
equation (1). Using PySCF package, we determine automatically the irreps of
each molecular orbital as well as the direct tensor product of the irrep of
the molecular spin-orbitals by using the product table46. Finally, we used a
gradient-based optimization method (BFGS) from scipy.optimize library52 for
optimizing the variational parameters $\theta_{i}^{a}$, $\theta_{i,j}^{a,b}$
and $\theta_{i,j,k}^{a,b,c}$. The gradient norm used is equal to 10-10.
For the calculation of classical chemistry methods CCSD and CCSD(T), we use
PYSCF package. However for the full CCSDT, we use CCT353 plugin that is
provided in Psi4 package54. Given that CC(t;3) method is supposed to yield
practically as accurate as full CCSDT for total electronic energies. The
energy convergence threshold used is 10-10 Ha.

##  4 Numerical Results and Discussions

We tested our UCCSDT-VQE (eq. (3)) and sym-UCCSDT-VQE (eqs. (4) – (6)) methods
on three sets of molecules: LiH, BeH2 and H2O. For sym-UCCSD and sym-UCCSDT-
VQE calculations, we chose the highest possible Abelian point groups that
belong to each of the systems i.e by refereeing to 55 we have $C_{2v}$ point
group for LiH and H2O but D2h point group for BeH2. The minimal STO-3G basis56
is used in order to minimize the computational costs. Within the STO-3G basis
set, LiH has a Hilbert space spanned by 12 Hartree-Fock (HF) orbitals (4
occupied, 8 virtual). BeH2 has 14 HF orbitals (6 occupied, 8 virtual), H2O has
14 HF orbitals (10 occupied, 4 virtual). In the multi-qubit register, we
assign each qubit to a spin-orbital, thus LiH (BeH2 and H2O) is (are)
represented by 12 (14) qubits. We used the Quantum Learning Machine49 (QLM)
simulator that is able to simulate up to 40 qubits (with large memory), and we
used the noiseless simulation mode (i.e. such ideal simulator mimics the
behavior of performing the experiment with an infinite number of shots) in
order to minimize computational run time and memory demands. First, we tested
the two methods on the LiH molecule at two bond lengths: $R=1.0$ Å and $R=3.0$
Å. As shown in Table 1, we observe that UCCSD and sym-UCCSD brings similar
numerical results. A similar observation can be made for the UCCSDT and sym-
UCCSDT results. This demonstrates that the presence of spin and orbital
symmetries in sym-UCC-VQE reaches the same result as UCC-VQE but in a more
efficient way and without influencing the overall accuracy. For example, Table
2 shows that the number of parameters is reduced from 188 to 58 for the case
of LiH, which in practice accelerates the optimization due to the fact that a
lower dimensional optimization problem is tackled. We confirm the convergence
of sym-UCCSD and sym-UCCSDT energies at gradient norm equal to 10-10 inside
BFGS optimizer (see the convergence as a function of optimization steps in
Figure 1 of the Supplementary Materials57).

Method | R=1.0 Å | R=3.0 Å
---|---|---
Ener. (Ha) | Err. (Ha) | Ener. (Ha) | Err. (Ha)
UCCSD | $\textbf{-7.7844}5508682$ | $5.19\times 10^{-6}$ | $\textbf{-7.798}7523587$ | $9.08\times 10^{-5}$
| Sym.UCCSD | $\textbf{-7.7844}5509065$ | $5.18\times 10^{-6}$ | $\textbf{-7.798}75240726$ | $9.07\times 10^{-5}$
UCCSDT(this work) | $\textbf{-7.7844602}5863$ | $2.14\times 10^{-8}$ | $\textbf{-7.798843}08319$ | $7.63\times 10^{-8}$
| Sym.UCCSDT | $\textbf{-7.7844602}5863$ | $2.13\times 10^{-8}$ | $\textbf{-7.798843}08642$ | $7.30\times 10^{-8}$
FCI (STO-3G) | $-7.78446028003$ | $-7.79884315950$
Table 1: UCC-VQE simulations for LiH at various bond lengths: $R=1.0$ Å and
$R=3.0$ Å. The obtained UCC-energies are in Hartree (Ha). The error (Err) is
the difference between the estimated UCC-energy and the FCI (the reference
result is given in the last row).

As mentioned in the literature26, 58, 25, 30, 49, testing the Trotterization
effect in the UCC ansatz is important since, as seen in equation (7), using
$t>1$ might in principle, decreases the Trotterization error. Besides
Trotterization study, a different method called disentangled UCC ansatz59 has
been introduced, and in 60 the ordering of operators in this ansatz has been
discussed. In our case we are interested in focusing on Trotterized ansatz
truncated at Triple excitations. To test this in our symmetric operators, we
use the sym-UCCSD and sym-UCCSDT methods with different Trotter steps $t$. It
should be noted that when we increased it beyond 1, the same ansatz (and
associated parameters) was chosen first and then we extended it. We observe in
Table 2 of the Supplementary Materials57 that increasing $t$, from $t=1$ up to
$t=15$ steps does not impact the results. A similar conclusion is observed in
the sym-UCCSDT case when we increased $t$ from $t=1$ up to $t=3$ steps. We
have thus carefully checked that most of the operators that are involved in
the symmetric version actually do not commute (see Table 3, Supplementary
Materials57 ) so that – a priori – there should be, at least, a slight
difference in the optimized discrete energies with different values of $t$,
but also there should be a difference on the coefficients in front of the
operators (results not shown : the coefficients for $t$ are those obtained for
$t=1$ divided by $t$). The reason we found is that – a posteriori, i.e. after
computation – we verified that, in both cases the operators associated to non-
negligible coefficients do commute (i.e. the coefficients in front of non
commuting operators are negligible). We do not see reasons for this but it
needs to be further analyzed.

In the remaining part of this Letter, we therefore use a UCC ansatz truncated
at the Single, Double, and Triple excitations using a single Trotter step.
When applying the point group symmetry method, the number of optimization
parameters was reduced between 14.2% (case of BeH2 having a D2h point group
symmetry) and 30.8% (case of LiH having a C2v point group symmetry), as shown
in Table 2. This consequently will significantly reduce the gate cost of
implementing triple excitations on quantum computers. Despite this reduction
in parameters in the sym-UCCSDT circuit-ansatz, this circuit depth can be
still considered as high compared to current NISQ devices capabilities. This
is linked to the fact that we are using the CNOT staircase4, 61 method that
maps the exponential of each fermionic excitation into one- and two-qubit
gates. Therefore, the number of CNOT (2-qubit gate) increases with the level
of excitations: 12 for single, 48 for double and 320 for Triples. Then the
circuit depths are 21, 84, 464, respectively. This can be illustrated for
example with the BeH2 case, as shown in Table 2. It contains 92 parameters,
associated to the following excitations: 8 singles, 26 doubles, 52 triples,
leading to a total CNOT gates in Sym-UCCSDT circuit equal to 12$\times$ 8 + 48
$\times$ 26 + 320 $\times$ 52 = 17984 . Attempts for reducing the CNOT counts
in triple excitation operators and higher-orders have been recently proposed
62. However, even with this reduction in CNOT counts, UCCSDT ansatz are not
yet amendable to current NISQ devices.

| Qbits | Sym. | Para Bef.(Af.) | $\%$ |  $\Delta E_{SymUCCSD}$ 18 | $\Delta E_{SymUCCSD}^{*}$ | $\Delta E_{SymUCCSDT}^{*}$
---|---|---|---|---|---|---|---
LiH | 12 | C2v | 188 (58) | 30.8% | $1.09\times 10^{-5}$ | $1.06\times 10^{-5}$ | $2.16\times 10^{-8}$
H2O | 14 | C2v | 340 (104) | 30.5% | $1.09\times 10^{-4}$ | $1.0\times 10^{-4}$ | $2.11\times 10^{-6}$
BeH2 | 14 | D2h | 644 (92) | 14.2% | $3.82\times 10^{-4}$ | $3.81\times 10^{-4}$ | $6.63\times 10^{-6}$
Table 2: UCC-VQE simulations for LiH, H2O and BeH2 calculated at equilibrium.
The geometry data are taken from CCCBDB-NIST Database55. Columns 2 and 3
corresponds to the number of qubits and to the point group used in each of the
molecules, respectively. The original parameters (and the parameters after
reduction) are given in the 4th column and the percent of the operators used
(i.e single, double and triple excitations) after parameters reduction are in
column 5. The UCC-VQE energy differences compared with the FCI energy (Ha) are
shown from the 6th to 8th columns. Data of $\Delta E_{SymUCCSD}$ are taken
from ref18. Symbol(*) represents this work.

Now let us discuss about sym-UCCSD and sym-UCCSDT energy errors, one of the
main focus in the article. It is important to note that we will not consider
the ordering of operators in the ansatz since it is believed to bring only
little accuracy improvements.30 However, it will, of course, remain of
potential interest to further study its impact. We believe that the inclusion
of triple excitations should generate a sufficient improvement of the ansatz
quality. The energy error we study is the difference between the predicted UCC
energy and FCI energy at the STO-3G basis set level. When the molecules are
chosen at their equilibrium geometry (see Table 2), the obtained sym-UCCSD
error is similar to the one obtained in reference 18, with only a very small
difference of $\approx$ 10-7 Ha. The calculated sym-UCCSDT improves the
accuracy with respect to sym-UCCSD by about 3 orders of magnitude in LiH and
by about 2-order of magnitudes in H2O and BeH2 molecules. This illustrates the
capabilities of Triple excitations to recover the correlation energy missed by
the sym-UCCSD ansatz at equilibrium. Moreover we calculate the sym-UCCSD and
sym-UCCSDT energies for LiH and BeH2 at several bond lengths $R$ ranging
between 0.5 Å to 3.5 Å, see Figure 1. As shown in the Figure, in the LiH case,
the sym-UCCSD and sym-UCCSDT errors remain below the chemical accuracy for all
bond lengths and the difference between sym-UCCSDT and sym-UCCSD error is
again about 3 orders of magnitude smaller. In the BeH2 case, when R is 2.0 Å,
sym-UCCSD error starts to get above the chemical accuracy (1.6 $\times$ 10-3
Hartree) and it gets further away as R increases. Overall, in term of error,
Sym-UCCSDT is better than sym-UCCSD ones by about 1-2 order of magnitudes at
all bond lengths. It remains below the chemical accuracy, except when R is 2.5
Å and beyond, where the error starts to approach the blue line delimiting the
chemical accuracy. When R is equal to 2.75 Å, a discontinuity is observed. In
fact this has been also seen in a previous work, see reference 63, around this
bond length in BeH2 calculated at STO-3G basis set level. There is one
hypotheses which illustrate this anomaly. Indeed, the importance of dynamic
correlation starts to increase when the bond length is around 2.7Å. At this
point, according to figure 5 in 63 the occupation number of the natural
orbitals starts to differ from zero or two. Hence many natural orbitals would
have significant occupation numbers (i.e different from zero), so the
difference between occupied and virtual orbitals tend to disappear. This
probably leads to the inability of our sym-UCCSDT approach to capture some
important excitations that might contribute strongly to the dynamical
correlation, and thus the sym-UCCSDT errors could increase. In general, this
limitation from triple excitations at large bond lengths in BeH2 can be
probably improved by adding higher-order excitations but is consistent with
classical couple cluster behaviours and the limits of single determinant
approaches. The data which were used to plot Figure 1 are provided in Tables 4
and 5 from the Supplementary Informations57.

(a) LiH: 12 qubits and C2v point group

(b) BeH2: 14 qubits and D2h point group

Figure 1: Error of the UCCSD and UCCSDT simulations using both spin and point
group symmetries for geometries LiH and BeH2 at different bond lengths. The
horizontal blue line represents the chemical accuracy, i.e. 1.6 $\times
10^{-3}$ Ha.

In order to compare the classical methods with UCC-VQE simulations, we
calculated CCSD, CCSD(T) and CCSDT-full errors for LiH, H2O and BeH2 at
equilibrium. The results are shown in Table 3. CCSD and sym-UCCSD errors are
very similar. This is also observed in 34 where Orbital optimized (OO)-UCCSD
method has been tested. The CCSD(T) errors deviates clearly from sym-UCCSDT
errors by at least one order of magnitudes. sym-UCCSDT also surpasses CCSDT-
full(i.e Coupled cluster with a full treatment of singles,doubles and triples)
in the three molecules, which is the most remarkable outcome here (see also
same behavior in Figures 2 (a) and 2(b) of the supplementary materials57 at
some bond lengths in LiH and BeH2, respectively). Such behaviour has been
noticed in 37, where it is always CCSDT that has lower quality than UCCSDT,
which was calculated analytically.

| $\Delta E_{CCSD}$ | $\Delta E_{SymUCCSD}^{*}$ | $\Delta E_{CCSD(T)}$ | $\Delta E_{CCSDT-full}$ | $\Delta E_{SymUCCSDT}^{*}$ | FCI
---|---|---|---|---|---|---
LiH | $1.05\times 10^{-5}$ | $1.06\times 10^{-5}$ | $2.11\times 10^{-6}$ | $2.79\times 10^{-8}$ | $2.16\times 10^{-8}$ | $-7.882403410335502$
H2O | $1.17\times 10^{-4}$ | $1.0\times 10^{-4}$ | $4.91\times 10^{-5}$ | $2.17\times 10^{-5}$ | $2.11\times 10^{-6}$ | $-75.01257824109094$
BeH2 | $3.94\times 10^{-4}$ | $3.81\times 10^{-4}$ | $1.84\times 10^{-4}$ | $4.22\times 10^{-5}$ | $6.63\times 10^{-6}$ | $-15.595176868923053$
Table 3: Error comparison between Classical methods (CCSD, CCSD(T) and CCSDT-
full) and UCC-VQE (sym.UCCSD and sym.UCCSDT) simulations for LiH, H2O and BeH2
calculated at equilibrium. The geometry data are taken from CCCBDB-NIST
Database. The UCC-VQE results are referred from Table 2, and those of CCSD,
CCSD(T), CCSDT-full are from Table 5 of the Supplementary Material57. The last
column represents the FCI energies at STO-3G basis set. The error is the
difference between the obtained energy and FCI in (Ha).

##  5 Conclusion

In conclusion, in this article, we presented the Trotterized UCCSDT ansatz. We
simplified the unitary coupled cluster triple excitation equations by using
the spin and point group symmetries. By testing several molecules on QLM
simulator, we demonstrated that the efficiency of the symmetric UCCSDT method
can greatly improved. Indeed, reducing the non-important single, double and
triple excitations from the circuit ansatz, one can accelerate the convergence
process as well as reducing the circuit depth. Moreover, by performing several
numerical simulations, we showed that the symmetric UCCSDT is able to reach a
superior accuracy compared to UCCSD (by at least two orders of magnitudes). In
the future, as quantum simulators and Quantum Processing Units improve and
enable to handle more gates, the Trotterized UCCSDT should enable to perform
accurate calculations with large chemical basis sets (i.e beyond the minimal
STO-3G). Overall, the Trotterized UCCSDT can bring interesting results for the
general quantum chemistry community since it is shown to be competitive with
”gold-standard” CCSD(T) classical methods.

M.R. acknowledges funding from European Union’s Horizon 2020 research and
innovation program, more specifically the $\langle$NE$|$AS$|$QC$\rangle$
project under grant agreement No. 951821. This work has also received funding
from the European Research Council (ERC) under the European Union’s Horizon
2020 research and innovation program (grant agreement No 810367), project EMC2
(JPP, YM). Support from the PEPR EPiQ and HQI programs is acknowledged.

In the Supplementary Materials57, we present all data related to the computed
total energies for sym-UCCSD and sym-UCCSDT found in Tables 1, 2 and 3.
Additionally, a detailed sketch algorithm that describe the triple excitations
within spin and point group symmetries is presented and coded in our openVQE
package50.

## References

  * Aspuru-Guzik et al. 2005 Aspuru-Guzik, A.; Dutoi, A. D.; Love, P. J.; Head-Gordon, M. Simulated quantum computation of molecular energies. _Science_ 2005, _309_ , 1704–1707
  * Reiher et al. 2017 Reiher, M.; Wiebe, N.; Svore, K. M.; Wecker, D.; Troyer, M. Elucidating reaction mechanisms on quantum computers. _Proceedings of the national academy of sciences_ 2017, _114_ , 7555–7560
  * Bauer et al. 2020 Bauer, B.; Bravyi, S.; Motta, M.; Chan, G. K.-L. Quantum algorithms for quantum chemistry and quantum materials science. _Chemical Reviews_ 2020, _120_ , 12685–12717
  * McArdle et al. 2020 McArdle, S.; Endo, S.; Aspuru-Guzik, A.; Benjamin, S. C.; Yuan, X. Quantum computational chemistry. _Reviews of Modern Physics_ 2020, _92_ , 015003
  * Bassman et al. 2021 Bassman, L.; Urbanek, M.; Metcalf, M.; Carter, J.; Kemper, A. F.; de Jong, W. A. Simulating quantum materials with digital quantum computers. _Quantum Science and Technology_ 2021, _6_ , 043002
  * Gan et al. 2006 Gan, Z.; Grant, D. J.; Harrison, R. J.; Dixon, D. A. The lowest energy states of the group-IIIA–group-VA heteronuclear diatomics: BN, BP, AlN, and AlP from full configuration interaction calculations. _The Journal of chemical physics_ 2006, _125_ , 124311
  * Lehtola et al. 2017 Lehtola, S.; Tubman, N. M.; Whaley, K. B.; Head-Gordon, M. Cluster decomposition of full configuration interaction wave functions: A tool for chemical interpretation of systems with strong correlation. _The Journal of chemical physics_ 2017, _147_ , 154105
  * Vogiatzis et al. 2017 Vogiatzis, K. D.; Ma, D.; Olsen, J.; Gagliardi, L.; De Jong, W. A. Pushing configuration-interaction to the limit: Towards massively parallel MCSCF calculations. _The Journal of chemical physics_ 2017, _147_ , 184111
  * Peruzzo et al. 2014 Peruzzo, A.; McClean, J.; Shadbolt, P.; Yung, M.-H.; Zhou, X.-Q.; Love, P. J.; Aspuru-Guzik, A.; O’brien, J. L. A variational eigenvalue solver on a photonic quantum processor. _Nature communications_ 2014, _5_ , 1–7
  * McClean et al. 2016 McClean, J. R.; Romero, J.; Babbush, R.; Aspuru-Guzik, A. The theory of variational hybrid quantum-classical algorithms. _New Journal of Physics_ 2016, _18_ , 023023
  * Cerezo et al. 2021 Cerezo, M.; Arrasmith, A.; Babbush, R.; Benjamin, S. C.; Endo, S.; Fujii, K.; McClean, J. R.; Mitarai, K.; Yuan, X.; Cincio, L., et al. Variational quantum algorithms. _Nature Reviews Physics_ 2021, _3_ , 625–644
  * Li et al. 2022 Li, W.; Huang, Z.; Cao, C.; Huang, Y.; Shuai, Z.; Sun, X.; Sun, J.; Yuan, X.; Lv, D. Toward practical quantum embedding simulation of realistic chemical systems on near-term quantum computers. _Chemical science_ 2022, _13_ , 8953–8962
  * O’Malley et al. 2016 O’Malley, P. J.; Babbush, R.; Kivlichan, I. D.; Romero, J.; McClean, J. R.; Barends, R.; Kelly, J.; Roushan, P.; Tranter, A.; Ding, N., et al. Scalable quantum simulation of molecular energies. _Physical Review X_ 2016, _6_ , 031007
  * Bharti et al. 2022 Bharti, K.; Cervera-Lierta, A.; Kyaw, T. H.; Haug, T.; Alperin-Lea, S.; Anand, A.; Degroote, M.; Heimonen, H.; Kottmann, J. S.; Menke, T., et al. Noisy intermediate-scale quantum algorithms. _Reviews of Modern Physics_ 2022, _94_ , 015004
  * Kühn et al. 2019 Kühn, M.; Zanker, S.; Deglmann, P.; Marthaler, M.; Weiß, H. Accuracy and resource estimations for quantum chemistry on a near-term quantum computer. _Journal of chemical theory and computation_ 2019, _15_ , 4764–4780
  * Lolur et al. 2021 Lolur, P.; Rahm, M.; Skogh, M.; García-Álvarez, L.; Wendin, G. Benchmarking the variational quantum eigensolver through simulation of the ground state energy of prebiotic molecules on high-performance computers. AIP Conference Proceedings. 2021; p 030005
  * Yeter-Aydeniz et al. 2021 Yeter-Aydeniz, K.; Gard, B. T.; Jakowski, J.; Majumder, S.; Barron, G. S.; Siopsis, G.; Humble, T. S.; Pooser, R. C. Benchmarking quantum chemistry computations with variational, imaginary time evolution, and Krylov space solver algorithms. _Advanced Quantum Technologies_ 2021, _4_ , 2100012
  * Cao et al. 2022 Cao, C.; Hu, J.; Zhang, W.; Xu, X.; Chen, D.; Yu, F.; Li, J.; Hu, H.-S.; Lv, D.; Yung, M.-H. Progress toward larger molecular simulation on a quantum computer: Simulating a system with up to 28 qubits accelerated by point-group symmetry. _Phys. Rev. A_ 2022, _105_ , 062452
  * Quantum et al. 2020 Quantum, G. A.; Collaborators*†,; Arute, F.; Arya, K.; Babbush, R.; Bacon, D.; Bardin, J. C.; Barends, R.; Boixo, S.; Broughton, M.; Buckley, B. B., et al. Hartree-Fock on a superconducting qubit quantum computer. _Science_ 2020, _369_ , 1084–1089
  * Kandala et al. 2017 Kandala, A.; Mezzacapo, A.; Temme, K.; Takita, M.; Brink, M.; Chow, J. M.; Gambetta, J. M. Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets. _Nature_ 2017, _549_ , 242–246
  * Hempel et al. 2018 Hempel, C.; Maier, C.; Romero, J.; McClean, J.; Monz, T.; Shen, H.; Jurcevic, P.; Lanyon, B. P.; Love, P.; Babbush, R., et al. Quantum chemistry calculations on a trapped-ion quantum simulator. _Physical Review X_ 2018, _8_ , 031022
  * Tazhigulov et al. 2022 Tazhigulov, R. N.; Sun, S.-N.; Haghshenas, R.; Zhai, H.; Tan, A. T.; Rubin, N. C.; Babbush, R.; Minnich, A. J.; Chan, G. K.-L. Simulating Models of Challenging Correlated Molecules and Materials on the Sycamore Quantum Processor. _PRX Quantum_ 2022, _3_ , 040318
  * Crawford and Schaefer III 2007 Crawford, T. D.; Schaefer III, H. F. An introduction to coupled cluster theory for computational chemists. _Reviews in computational chemistry_ 2007, _14_ , 33–136
  * Bartlett and Musiał 2007 Bartlett, R. J.; Musiał, M. Coupled-cluster theory in quantum chemistry. _Reviews of Modern Physics_ 2007, _79_ , 291
  * Romero et al. 2018 Romero, J.; Babbush, R.; McClean, J. R.; Hempel, C.; Love, P. J.; Aspuru-Guzik, A. Strategies for quantum computing molecular energies using the unitary coupled cluster ansatz. _Quantum Science and Technology_ 2018, _4_ , 014008
  * Anand et al. 2022 Anand, A.; Schleich, P.; Alperin-Lea, S.; Jensen, P. W.; Sim, S.; Díaz-Tinoco, M.; Kottmann, J. S.; Degroote, M.; Izmaylov, A. F.; Aspuru-Guzik, A. A quantum computing view on unitary coupled cluster theory. _Chemical Society Reviews_ 2022,
  * Hatano and Suzuki 2005 Hatano, N.; Suzuki, M. _Quantum annealing and other optimization methods_ ; Springer, 2005; pp 37–68
  * Sokolov et al. 2020 Sokolov, I. O.; Barkoutsos, P. K.; Ollitrault, P. J.; Greenberg, D.; Rice, J.; Pistoia, M.; Tavernelli, I. Quantum orbital-optimized unitary coupled cluster methods in the strongly correlated regime: Can quantum algorithms outperform their classical equivalents? _The Journal of chemical physics_ 2020, _152_ , 124107
  * Xia and Kais 2020 Xia, R.; Kais, S. Qubit coupled cluster singles and doubles variational quantum eigensolver ansatz for electronic structure calculations. _Quantum Science and Technology_ 2020, _6_ , 015001
  * Grimsley et al. 2019 Grimsley, H. R.; Claudino, D.; Economou, S. E.; Barnes, E.; Mayhall, N. J. Is the trotterized uccsd ansatz chemically well-defined? _Journal of chemical theory and computation_ 2019, _16_ , 1–6
  * Yung et al. 2014 Yung, M.-H.; Casanova, J.; Mezzacapo, A.; Mcclean, J.; Lamata, L.; Aspuru-Guzik, A.; Solano, E. From transistor to trapped-ion computers for quantum chemistry. _Scientific reports_ 2014, _4_ , 1–7
  * Shen et al. 2017 Shen, Y.; Zhang, X.; Zhang, S.; Zhang, J.-N.; Yung, M.-H.; Kim, K. Quantum implementation of the unitary coupled cluster for simulating molecular electronic structure. _Physical Review A_ 2017, _95_ , 020501
  * O’Brien et al. 2022 O’Brien, T. E.; Anselmetti, G.; Gkritsis, F.; Elfving, V.; Polla, S.; Huggins, W. J.; Oumarou, O.; Kechedzhi, K.; Abanin, D.; Acharya, R., et al. Purification-based quantum error mitigation of pair-correlated electron simulations. _arXiv preprint arXiv:2210.10799_ 2022,
  * Mizukami et al. 2020 Mizukami, W.; Mitarai, K.; Nakagawa, Y. O.; Yamamoto, T.; Yan, T.; Ohnishi, Y.-y. Orbital optimized unitary coupled cluster theory for quantum computer. _Physical Review Research_ 2020, _2_ , 033421
  * Sennane et al. 2023 Sennane, W.; Piquemal, J.-P.; Rančić, M. J. Calculating the ground-state energy of benzene under spatial deformations with noisy quantum computing. _Phys. Rev. A_ 2023, _107_ , 012416
  * Lee et al. 2018 Lee, J.; Huggins, W. J.; Head-Gordon, M.; Whaley, K. B. Generalized unitary coupled cluster wave functions for quantum computation. _Journal of chemical theory and computation_ 2018, _15_ , 311–324
  * Köhn and Olsen 2022 Köhn, A.; Olsen, J. Capabilities and limits of the unitary coupled-cluster approach with generalized two-body cluster operators. _The Journal of Chemical Physics_ 2022, _157_ , 124110
  * Grimsley et al. 2019 Grimsley, H. R.; Economou, S. E.; Barnes, E.; Mayhall, N. J. An adaptive variational algorithm for exact molecular simulations on a quantum computer. _Nature communications_ 2019, _10_ , 1–9
  * Fedorov et al. 2022 Fedorov, D. A.; Alexeev, Y.; Gray, S. K.; Otten, M. Unitary Selective Coupled-Cluster Method. _Quantum_ 2022, _6_ , 703
  * Smart et al. 2022 Smart, S. E.; Boyn, J.-N.; Mazziotti, D. A. Resolving correlated states of benzyne with an error-mitigated contracted quantum eigensolver. _Physical Review A_ 2022, _105_ , 022405
  * Smart and Mazziotti 2021 Smart, S. E.; Mazziotti, D. A. Quantum solver of contracted eigenvalue equations for scalable molecular simulations on quantum computing devices. _Physical Review Letters_ 2021, _126_ , 070504
  * Raghavachari 2013 Raghavachari, K. Historical perspective on: A fifth-order perturbation comparison of electron correlation theories [Volume 157, Issue 6, 26 May 1989, Pages 479–483]. _Chem. Phys. Lett_ 2013, _589_ , 35–36
  * Ramabhadran and Raghavachari 2013 Ramabhadran, R. O.; Raghavachari, K. Extrapolation to the gold-standard in quantum chemistry: computationally efficient and accurate CCSD (T) energies for large molecules using an automated thermochemical hierarchy. _Journal of chemical theory and computation_ 2013, _9_ , 3986–3994
  * Dresselhaus et al. 2007 Dresselhaus, M. S.; Dresselhaus, G.; Jorio, A. _Group theory: application to the physics of condensed matter_ ; Springer Science & Business Media, 2007
  * Tinkham 1964 Tinkham, M. Group Theory and Quantum Mechanics McGraw-Hill: New York. 1964
  * Atkins et al. 1970 Atkins, P. W.; Child, M. S.; Phillips, C. S. G. _Tables for group theory_ ; Oxford University Press Oxford, 1970; Vol. 6
  * Fradkin 1989 Fradkin, E. Jordan-Wigner transformation for quantum-spin systems in two dimensions and fractional statistics. _Physical review letters_ 1989, _63_ , 322
  * 48 myqlm-fermion. https://github.com/myQLM/myqlm-fermionhttps://myqlm.github.io/qat-fermion.htmlRepository: https://github.com/myQLM/myqlm-fermion and documentation: https://myqlm.github.io/qat-fermion.html
  * Haidar et al. 2022 Haidar, M.; Rančić, M. J.; Ayral, T.; Maday, Y.; Piquemal, J.-P. Open Source Variational Quantum Eigensolver Extension of the Quantum Learning Machine (QLM) for Quantum Chemistry. _arXiv preprint arXiv:2206.08798_ 2022,
  * 50 Open-VQE package. https://github.com/OpenVQE/OpenVQE.githttps://openvqe.github.io/OpenVQE/Repository: https://github.com/OpenVQE/OpenVQE.git and documentation:https://openvqe.github.io/OpenVQE/
  * Sun et al. 2018 Sun, Q.; Berkelbach, T. C.; Blunt, N. S.; Booth, G. H.; Guo, S.; Li, Z.; Liu, J.; McClain, J. D.; Sayfutyarova, E. R.; Sharma, S., et al. PySCF: the Python-based simulations of chemistry framework. _Wiley Interdisciplinary Reviews: Computational Molecular Science_ 2018, _8_ , e1340
  * Virtanen et al. 2020 Virtanen, P. et al. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. _Nature Methods_ 2020, _17_ , 261–272
  * 53 CCT3: A PSI4 plugin which performs active-space coupled-cluster CCSDt calculations and which can determine noniterative corrections to CCSDt defining the CC(t;3) approach. https://github.com/piecuch-group/psi4_cct3
  * Smith et al. 2020 Smith, D. G.; Burns, L. A.; Simmonett, A. C.; Parrish, R. M.; Schieber, M. C.; Galvelis, R.; Kraus, P.; Kruse, H.; Di Remigio, R.; Alenaizan, A., et al. PSI4 1.4: Open-source software for high-throughput quantum chemistry. _The Journal of chemical physics_ 2020, _152_ , 184108
  * Johnson 2022 Johnson, N. R. Computational Chemistry Comparison and Benchmark DataBase. 2022,
  * Hehre et al. 1970 Hehre, W.; Ditchfield, R.; Stewart, R.; Pople, J. self-consistent molecular orbital methods. iv. use of Gaussian expansions of Slater-type orbitals. Extension to second-row molecules. _The Journal of Chemical Physics_ 1970, _52_ , 2769–2773
  * 57 Haidar, M.; Rančić, M. J.; Maday, Y.; Piquemal, J.-P. See Supplementary Materials for Triple excitation algorithmic code and the total UCC energies.
  * Barkoutsos et al. 2018 Barkoutsos, P. K.; Gonthier, J. F.; Sokolov, I.; Moll, N.; Salis, G.; Fuhrer, A.; Ganzhorn, M.; Egger, D. J.; Troyer, M.; Mezzacapo, A., et al. Quantum algorithms for electronic structure calculations: Particle-hole Hamiltonian and optimized wave-function expansions. _Physical Review A_ 2018, _98_ , 022322
  * Evangelista et al. 2019 Evangelista, F. A.; Chan, G. K.-L.; Scuseria, G. E. Exact parameterization of fermionic wave functions via unitary coupled cluster theory. _The Journal of Chemical Physics_ 2019, _151_ , 244112
  * Izmaylov et al. 2020 Izmaylov, A. F.; Díaz-Tinoco, M.; Lang, R. A. On the order problem in construction of unitary operators for the variational quantum eigensolver. _Phys. Chem. Chem. Phys._ 2020, _22_ , 12980–12986
  * Whitfield et al. 2011 Whitfield, J. D.; Biamonte, J.; Aspuru-Guzik, A. Simulation of electronic structure Hamiltonians using quantum computers. _Molecular Physics_ 2011, _109_ , 735–750
  * Magoulas and Evangelista 2022 Magoulas, I.; Evangelista, F. A. CNOT-Efficient Circuits for Arbitrary Rank Many-Body Fermionic and Qubit Excitations. _arXiv preprint arXiv:2210.05771_ 2022,
  * Van Raemdonck et al. 2015 Van Raemdonck, M.; Alcoba, D. R.; Poelmans, W.; De Baerdemacker, S.; Torre, A.; Lain, L.; Massaccesi, G. E.; Van Neck, D.; Bultinck, P. Polynomial scaling approximations and dynamic correlation corrections to doubly occupied configuration interaction wave functions. _The Journal of Chemical Physics_ 2015, _143_ , 104106


# Deep neural networks from the perspective of ergodic theory

Fan Zhang  Gravitational Wave and Cosmology Laboratory, Department of
Astronomy, Beijing Normal University, Beijing 100875, China  Advanced
Institute of Natural Sciences, Beijing Normal University at Zhuhai 519087,
China

###### Abstract

The design of deep neural networks remains somewhat of an art rather than
precise science. By tentatively adopting ergodic theory considerations on top
of viewing the network as the time evolution of a dynamical system, with each
layer corresponding to a temporal instance, we show that some rules of thumb,
which might otherwise appear mysterious, can be attributed heuristics.

##  I Introduction and motivation

Artificial neural networks have demonstrated great potential in their ability
to learn existing knowledge, and interpolate or even slightly extrapolate to
new situations. They however lack the ability to understand causations and
other logical relations that indicate general intelligence. No matter, much of
human activities are experience-based, and so the present incarnation of
artificial intelligence algorithms are sufficient to revolutionize society. In
particular, one can highlight medicine as one area where vast amounts
experiential enigmas persist, and a doctor’s effectiveness is largely
restricted by their capacity to learn rather than the ability to understand.

Furthermore, human bodies are complex systems whereby heterogeneity creates
self-organization, which an artificial neural network, as also a complex
dynamical system, may be able to (intentionally or unintentionally) emulate.
Thereby, a network’s learning power may well be tunable to arise out of
emergent phenomena that share traits with the processes underlying human body
functions. In this way, a network could serve as a crude simulation, and
subsequently stands a better chance at accurately and efficiently learning and
storing medical knowledge, providing better interpolations and extrapolations
than the usually naive and linear manner with which humans try to carry out
such tasks.

Despite such great potential and practical successes, a theory of deep
learning is still lacking however, so we are in want of general guidelines as
to what architecture might perform well. In such a pursuit, it is beneficial
to be able to examine deep neural networks from as many perspectives as
possible, thereby acquiring a more complete picture. In this brief note, we
advocate also adopting the ergodic theory approach by showing how it could
offer simple intuitive (although arguably hand-wavy at its present state of
deployment) explanations to some properties that we have observed in the
behavior of the deep neural networks. We begin by relating some understood
aspects of the networks to ergodicity concepts, which could then serve as the
conduits connecting the two disciplines.

###  I.1 As fitting functions to training data

We are just now beginning to understand some aspects of the behavior of deep
neural networks. For example, one’s experiences with lower dimensional non-
convex optimization problems might suggest that there would be an enormous
number of local minima, on the cost function surface in parameter space, that
trap our optimization procedure and keep it away from the optimal solution.
This is in fact not a problem Pascanu et al. (2014), because the local minima
are replaced with saddle points in very high dimensions. The Hessian at
critical points have many eigenvalues in high parameter dimensions, and it is
more likely that they take on both positive and negative values giving us
saddle points, so we are less likely to be trapped, and could instead slip out
along the downwardly curving directions. But more fundamental than not being
trapped, is that a good-enough solution is at all on that surface to begin
with. In other words,
$\mathcal{C}_{1}:$ _the network needs to be sufficiently flexible so the
surface extends a large range of cost values._

Another aspect of neural networks that people had expected to be a problem is
statistical sample complexity, i.e., having too many parameters, even more
than the training data sets, should lead to overfitting, so while fitting to
the (possibly noisy) training data may be perfect, interpolation, let alone
extrapolation, would not work as the overfitted function may oscillate wildly
in-between the training data points on which they are pinned. This behavior is
not observed in reality. One explanation is proposed by Bartlett et al.
(2020), which rigorously proved that having many parameters being unimportant
for the fitting quality is key to the effectiveness of overfitted linear
regression. An intuitive understanding being offered in literature (see e.g.,
Sejnowski (2020)) is that there are then many different solutions that work
rather well, differing only in the unimportant parameters, thus it becomes
much easier to find them. However, we note that for any fitting problem, one
could always introduce completely spurious parameters that don’t appear in the
actual fitting procedure (we just vacuously append them to the formal
parameter set to be fitted), thus are unimportant to the extreme, but this
should not change at all the quality of output of that fitting procedure. So
this explanation would make sense only if the original procedure is not
overfitting to begin with, or in other words, so many of the fitting
parameters in the linear regression are unimportant or spurious that the
surviving important ones are so few in numbers that there is no overfitting in
reality. This is thus a rather trivial explanation – the overfitting problem
is not problematic if for the problem at hand the fitting architecture is
nothing but superficially over-parameterized. This does not appear to be the
case for deep neural networks though, for changing the network structure,
especially the number of neurons, appears to alter the fitting quality.
Therefore, an alternative explanation for why overfitting does not present a
catastrophic snag may be needed.

To this end, we first note that condition $\mathcal{C}^{1}$ already helps, as
overfitting often occurs in regression problems when the fitting curve is just
sufficiently flexible to hit all the training data points if we strong-arm it,
but not enough to do so smoothly if noises are present or the choice of basis
functions are not appropriate. This rigidity means the fitting curves have to
take a large detour in order to make the correct directional changes in order
to hit the next training point needing to be fitted, much like how a fast
aircraft needs to make a very large divagation in order to make a turn. When
we do interpolation or extrapolation, we have to sample from within those
large detours and thus yield terrible results. When the parameter number goes
much larger than training point population however, the curve becomes
extremely flexible, so it could possibly effortlessly (smoothly) transit
between training data points, without having to take up some very contorted
shape that shoots off to large extremes in the intervening values. This is the
case with deep neural networks that are complex and flexible (see e.g., Hornik
(1991); Elbrächter et al. (2019) for neutral networks as universal
approximators). This is not enough though. If the fitting quality to the
training data is the only criteria in the cost function, then there is in
principle no guarantee that the flexible fitting curve won’t still whip around
like crazy. That is, the smooth curves exist, but are not necessarily the ones
we get if we do not deliberately look for them. To cure this, one typically
introduces these so-called regularizers into the cost functions to penalize
undesired parameters that could cause instabilities of the output, as one
rather contrived approach to achieve
$\mathcal{C}_{2}:$ _the network needs to be sufficiently wellposed so the
output doesn’t depend so extremely sensitively on initial data that
infinitesimally separated initial data points leads to wildly diverging
outcomes._

###  I.2 As dynamical systems

The conditions $\mathcal{C}_{1}$ and $\mathcal{C}_{2}$ are qualitative and
difficult to translate into some enforceable quantitative criteria on network
design. To make further progress, we enlist the dynamical systems view of
neural networks (see e.g., Furuhata et al. (2021)), whereby each layer of
neurons corresponds to a time instance of a dynamical system, so updates
propagating through the network according to some given initial data behaves
like a time evolution of the system within a state space having the same
dimension as the number of neurons in each layer, where the number of layers
becomes the number of discrete111Previous literature would typically go to the
continuous limit and adopt differential equation results in order to glean
some insight into the behavior and properties of the neural networks (see in
particular Belozyorov and Volkova (2022), where some criteria for neural ODEs
to develop chaotic behavior is worked out). This strategy is not suitable for
us however, since the chaotic nature or lack thereof, differs quite
drastically for discrete and continuous systems, not least because the
requirement for no-self-intersection or uniqueness of trajectories is so very
much relaxed for discrete evolutions. The (spatial) differentiability also
effectively restricts $K_{\alpha}^{[j]\beta}$ of Eq. (3) below to tridiagonal
or other nearly diagonal (depending on the order of derivatives and the finite
difference scheme adopted) forms that limit the dependence of a neuron’s time
evolution to only its nearest neighbours. time steps. This dynamical systems
perspective had already yielded important insight, e.g, the Lemma 1 of Haber
and Ruthotto (2018), which is a well-established result applicable to the
numerical methods for solving differential equation, essentially imposes a
necessary condition for $\mathcal{C}_{2}$. On the other hand, the same paper
also recognizes the importance for the forward propagation to not be too lossy
(i.e., the Lyapunov exponents cannot all be too negative, or the state space
volume shrinks very quickly, losing the ability to distinguish information in
the initial data), essentially offering a caution on what might definitely
violate $\mathcal{C}_{1}$.

In this note, we further enlist the powerful mathematical tool of ergodic
theory222The usual ergodic theory studies measure-invariant or state space
volume preserving evolutions, or in other words, there is no dissipation in
the system, which may not be the case for deep neural networks. But
dissipation may not matter directly, because the sum of _all_ Lyapunov
exponents determines whether we have dissipation, but chaos ($\exists$ at
least one positive exponent), ergodicity, weakly-mixing (only constant
functions are the eigenfunctions corresponding to the exponents as
eigenvalues) and entropy (sum of all positive exponents) are really more about
the positive ones. Also, as discussed in the main text in the last paragraph,
a well-designed neural network should not be too dissipative. to help tighten
the discussion, inching a little more towards necessary and sufficient
conditions. In this language, $\mathcal{C}_{1}$ translates to a desirability
to have ergodicity333No invariant subsets that trap orbits thus prevent
effective migration in e.g., the classification problem. Beginning from any
set of initial data, the bundle of trajectories will eventually cover almost
all of the allowed state space (besides perhaps sets of zero measure; these
don’t matter for statistical considerations, but if we are aiming to train the
neural network to search for very special properties that nearly never occurs
in a big data set, these will be relevant and so neural networks cannot
produce master pieces of exceptional qualities, they are however efficient at
churning out mediocrity). We can see this from the celebrated Birkhoff’s
ergodicity theorem that underlies statistical mechanics, which states that
spatial average equals temporal average along the dynamical evolution, so any
macrostate with non-vanishing measure must eventually be reached., while
$\mathcal{C}_{2}$ means we should avoid mixing444Stronger than ergodicity,
requiring that memory of the initial data be completely lost as we take many
steps, in the sense that the trajectories emerging from any initial data set
will have to spread out over the entire state space completely randomly
according to the probability distribution of the overall state space, with
conditional probability conditioning on the initial data adding no further
information. For mixing systems, trajectories from any two initially disjoint
initial data sets get meshed together thoroughly and permanently after
sufficiently long evolution. Note mere ergodicity could also mesh two bundles
of trajectories out of two disjoint sets of initial data, but the meshing can
be transitory, meaning the two bundles of trajectories can separate again, and
then repeat in the meshing and separating cycle (cf., the claim that merely
ergodic systems are not even apparently random Frigg et al. (2020)). With
mixing though, the meshing is permanent with no subsequent re-separation. In a
way, one can visualize the two bundles as meeting (if they meet at all)
transversely over and over again with an ergodic-but-not-mixing system, but
they intersect and merge tangentially with a mixing system.. This preference
to wedge in-between ergodicity and mixing is easiest to understand in regard
to classification tasks, where we need ergodicity to be able to move any
initial point almost anywhere else in state space to achieve segregation (see
e.g., Fig. 3 in Furuhata et al. (2021)), but not mixing that causes
neighbourhoods to all get mangled together that there is no segregation
possible, and classification ends up having to be done point-wise while
interpolation becomes impossible.

We will hereafter refer to this dedicate balanced wildness of the dynamical
system’s orbits as being on the edge of chaos, taking the view that strong
mixing is a hallmark for chaotic systems555Even though it does not necessarily
imply exponential divergence between trajectories, so mixing systems are
sometimes referred to as being weakly chaotic. Werndl (2009); Frigg et al.
(2020), while merely ergodic systems are usually not chaotic. We caution that
this relationship between chaos and ergodicity concepts is not rigorous, and
is but a pragmatism that is useful for practical applications. We will adopt
it in this note understanding this lack of rigor. We will similarly be loose
with terminologies in the interest of brevity, especially when migrating
concepts defined for infinite time system to finite ones, and those defined
for measure-preserving systems to more general cases.

##  II Network spectroscopy

As always with dynamical systems, turning to spectral considerations tends to
simplify computations. The relevant quantities for our considerations, in the
case of deep neural networks of finite depth, are the finite time Lyapunov
exponents, which measure the tendency of the dynamics to locally drive nearby
trajectories apart (note the qualifying “finite time” does not mean they are
gauges of the cumulative divergence across the entire network depth, they are
still “per layer” quantities as the division by $[j]$ in Eq. 5 below shows),
and can be seen as the average local Lyapunov exponents across the available
depth of the network. Once we have these numbers, the various dynamical system
qualities can be assessed, e.g., whether the system is dissipative (volume
preserving in the Liouville sense) depends on whether the sum of the exponents
is negative, and particularly relevant for us, the system is chaotic and
likely practically mixing if there exists at least one positive exponent
(becomes hyperchaotic if there are more than one).

To this end, we first recast the neural network into the dynamical systems
language (see e.g., Benning et al. (2019)). For definitiveness, we assume a
basic multi-layer perceptron style network. If we see any particular
configuration of a layer (labelled by $[j]\in 0,\cdots N-1$) in the deep
neural network as being a point in a state variable space, with each neuron
within occupying a dimension (labelled by $\alpha$), then the value being
carried by that neuron $y^{[j]\alpha}_{i}$, as corresponding to some input
state $y^{[0]\beta}_{i}$ ($i$ indexes the training set), gives the $\alpha$th
coordinate of that point in the state variable space. We can subsequently
regard the layers as time steps in an evolution in this state space, where the
discrete evolution is written as

$\displaystyle y^{[j+1]\alpha}_{i}=$
$\displaystyle\tilde{f}^{[j]\alpha}\left(y^{[j]\beta}_{i}\right)$ (1)
$\displaystyle=$ $\displaystyle
y^{[j]\alpha}_{i}+f^{\alpha}\left(y^{[j]\beta}_{i},{\bf u}^{[j]}\right)\Delta
t\,,$ (2)

where usually

$\displaystyle f^{\alpha}\left(y^{[j]\beta}_{i},{\bf
u}^{[j]}\right)=\sigma_{\beta}^{\alpha}\left(K_{\delta}^{[j]\beta}y_{i}^{[j]\delta}+\xi^{[j]\beta}\right)\,,$
(3)

with ${\bf u}^{[j]}$ representing the weights $K_{\alpha}^{[j]\beta}$ and bias
$\xi^{[j]\alpha}$ connecting the $[j]$th layer to the $[j+1]$th (note these
parameters become independent of the input once training is complete). The
$\sigma_{\beta}^{\alpha}$ in Eq. (3) is an activation function for the neurons
that’s usually of a diagonal form. The expression (1) is a more abstract
representation of the dynamical process, which absorbs the layer dependent
variations in the parameters ${\bf u}^{[j]}$ into the $[j]$ label of
$\tilde{f}^{[j]\alpha}$.

Assuming for simplicity that the width of the layers do not change, then one
obtains, for each reference trajectory (e.g., starting from the input of a
training set labelled by $i$), and an end time $j$ (usually $N-1$ but we keep
the definition general here; we always fix the initial time at $j=0$ however),
a square matrix

$\displaystyle\mathbb{M}^{[j]\alpha}_{i\,\beta}\equiv\frac{\partial
y_{i}^{[j]\alpha}}{\partial y_{i}^{[0]\beta}}\,.$ (4)

This matrix, when acting on a perturbation vector of the initial data $\delta
y^{[0]\beta}_{i}$, yields the leading order changes in the state
$y^{[j]\alpha}_{i}$ at the $[j]$th layer. Therefore the exponential rate of
growth of the perturbation would be related to the logarithm of its
eigenvalues, rescaled by the number of time steps $[j]$. However, there is no
guarantee that $\mathbb{M}^{[j]\alpha}_{i\,\beta}$ is diagonalizable or even
square in the more general cases, so one instead goes to the singular values
$\mu^{[j]\alpha}_{i}$ of $\mathbb{M}^{[j]\alpha}_{i\,\beta}$ in its singular
value decomposition. The finite time Lyapunov exponents are then defined
as666The singular values $\mu^{[j]\alpha}_{i}$ can always be chosen to be
positive, since multiplying a negative singular value and the corresponding
left singular vector both by minus one will preserve the validity of the
decomposition.

$\displaystyle\lambda^{[j]\alpha}_{i}\equiv\frac{\ln\mu^{[j]\alpha}_{i}}{[j]}\,.$
(5)

To compute such exponents, we begin with the most general expression (1),
which after iterations yield

$\displaystyle
y^{[j]\alpha}_{i}=\circ^{j-1}_{q=0}\tilde{f}^{[q]\alpha}\left(y^{[0]\beta}_{i}\right)\,.$
(6)

So when the $\tilde{f}^{[q]\alpha}$ are all differentiable (because of the
need to back-propagate errors in neural networks, its derivatives usually
exist, maybe aside from at isolated points like with the ReLU activation
function), we can write

$\displaystyle\frac{\partial y_{i}^{[j]\alpha}}{\partial
y_{i}^{[0]\beta}}=\delta^{\alpha}{}_{\gamma_{j}}\delta_{\beta}{}^{\gamma_{0}}\prod^{j-1}_{q=0}\tilde{J}^{[q]\gamma_{q+1}}_{i\gamma_{q}}\,,$
(7)

where the local Jacobians are

$\displaystyle\tilde{J}^{[q]\gamma_{q+1}}_{i\gamma_{q}}\equiv\frac{\partial\tilde{f}^{[q]\gamma_{q+1}}}{\partial
y^{[q]\gamma_{q}}_{i}}\,.$ (8)

If each layer is identical so all the Jacobians are the same, then the method
of powers for singular value decomposition suggests that
$\mathbb{M}^{[j]\alpha}_{i\,\beta}$ is essentially rank one, dominated by the
largest singular value when the neural network is deep. This low rank scenario
has minimal expressivity, so it is important that the Jacobians are varied. In
this case, there is no simple relationship between the singular values of the
individual Jacobians and the final $\mathbb{M}^{[j]\alpha}_{i\,\beta}$. At
best, it is possible to approximately almost-diagonalize
$\mathbb{M}^{[j]\alpha}_{i\,\beta}$ (taking it to a bidiagonal form) with a
Householder procedure that turns all the Jacobians upper triangular Golub et
al. (2000). Then, the bidiagonal form is further fully diagonalized according
to the “chasing the bulge” sweep of e.g. Demmel and Kahan (1990), using 2-D
rotations that suppress the superdiagonal entries. This entire procedure is
rather opaque and blends the entries of the Jacobians in nontrivial manners,
therefore, without specializing to specific neural networks, we cannot make
definitive statements regarding the finite time Lyapunov exponents.
Nevertheless, we could perhaps make some vaguely probabilistic statements on
how various aspects of the network architecture (e.g., the number of layers,
the width of each layer, the rank of each Jacobian etc) would _likely_ affect
the singular values of $\mathbb{M}^{[j]\alpha}_{i\,\beta}$.

###  II.1 Effect of network architectural traits

####  II.1.1 Depth of network vs. activation function

A major hurdle in applying ergodic theory to deep neural networks is that the
concepts of the former are defined in the infinite evolution time limit, or in
other words asymptotically, while the actual networks tend to be of finite
depth. For the same local spectral characteristics, the ergodicity (or
transitivity in the language of topological dynamics) and mixing properties
therefore manifest more fully in deeper networks. As a result, greater depth
is desirable for networks whose individual layers don’t tend to push nearby
trajectories apart, i.e., for dynamics that are on the regular side, meaning
barely ergodic and far from mixing. On the other hand, for networks whose
individual-layer-driven local dynamics have a strong tendency to cause
neighboring trajectories to diverge, or in other words, a system that is
highly chaotic and deeply in the mixing regime asymptotically, shallower
networks are beneficial, as transitivity may have had time to transpire while
mixing hasn’t (viewing these concepts through the lens of losing dependence on
initial data, then mixing is a more complete amnesia that tends to happen
further into the evolution than ergodicity, which is only partial loss of
memory).

With the prevailing network designs, because the activation functions (some of
which are really just distributions) are more or less binary (switching
between active or dormant states) by definition, the outputs of each neuron
for nearby trajectories that happen to lay on either side of the threshold
tend to be very different777This variability against initial data is a
signature of the nonlinearity introduced by the activation functions. A truly
linear activation function will lead to a constant $\mathbb{M}$ matrix, and
the entire network becomes a linear regression, which may not be able to fit
the training data (e.g., when two inputs in the training set are related by a
simple $1/2$ rescaling, but their corresponding outputs differ not by the same
rescaling). On an intuitive level, such variability across finite and not
infinitesimal shifts in initial data may also be conducive to generate the
ergodic-but-not-mixing behavior discussed in footnote 4, because different
bundles of trajectories experiencing different local Jacobians tends to drive
the bundles, but not necessarily the individual trajectories within each
bundle, to move apart.. So we usually have the latter case, thus shallower
networks may often be desirable (see e.g., Liu and Zaidi (2016)). However, the
extent to which this is true is dependent on the activation functions, which
is supposed to bring in nonlinearity, thus should naturally relate to how
chaotic the network evolution is. More specifically:

  * •

One category consists of the binary step, sigmoid, or tanh activation
functions, which all have small derivatives far from the transition region,
but a large derivative within. As a result, their Jacobians contain a delta-
function style spike that, at occasions, contributes to a very large Frobenius
norm to $\tilde{J}$ (we henceforth suppress unimportant indices for brevity)
and subsequently $\mathbb{M}$ (sans chance cancellations). Because the
Frobenius norm of $\mathbb{M}$ is the L2 norm of its singular values, this
then implies large singular values, and thus large (more positive) finite time
Lyapunov exponents, and consequently deeper protrusion into the mixing regime.
In summary, if one would like to adopt these $S$-shaped activation functions,
shallower networks are needed if one sets the transition region in these
functions to be very narrow.

  * •

The recently more popular ReLU, ELU or swish activation functions behave very
differently, as there is not spike, but only a (actual or almost)
discontinuity in the Jacobians, so there is no large Frobenius norm at
individual $\tilde{J}$ level, just jumps in their entries when the state of
the relevant layer changes. Once the local Jacobians multiply into the overall
$\mathbb{M}$, the entries in that matrix end up jumping frequently
(multiplication of many Heaviside functions jumping at different values).
Therefore, when we scan across all possibilities (varying $i$), the
$\mathbb{M}_{i}$ matrix changes often and $\lambda_{i}$ tends to explore large
ranges, thus have a high probability to hitting large positive values. This
effect is likely less pronounced than that of the previous item where large
$\lambda_{i}$ values are more definitively hit (one can alternatively think of
these call-option-payoff shaped functions as being less binary so nearby
trajectories don’t elicit very different activation function outcomes,
implying less chaotic propagations), so we predict that ReLU family functions
would be more suitable for networks of greater depth.

####  II.1.2 Width of layers vs. connectivity

The width $D$ of the layers of a neural network (number of neurons in each
layer, assuming to be a constant for the present discussion for brevity)
corresponds to the dimension of the state space of the dynamical system, and
$\mathbb{M}$ is a $D\times D$ matrix. How the finite time Lyapunov exponents
vary with $D$ is heavily dependent on how the connection matrix
$K^{\alpha}_{\beta}$ in Eq. (3), between neurons in adjacent layers, changes
when we scale up the state space dimension:

  * •

First consider increasing $D$ without changing the connectivity between the
layers (i.e., the percentage of neurons in the next layer that a particular
neuron is connected to, or in other words the percentage of non-vanishing
entries in each row of $K^{\alpha}_{\beta}$) or the coupling strength (the
typical size of the entries in the weighting matrix $K$), then the sparsity of
$\mathbb{M}$ and the typical amplitude of entries in it won’t change,
resulting in the Frobenius norm scaling as $D^{1}$ since the number of these
entries grow as $D^{2}$. The number of singular values on the other hand only
grows as $D^{1}$, so the average singular value would have to scale as
$\propto\sqrt{D}$. In other words, without changing other features like
connectivity and weight ranges etc, a wider neural network tends to possess
more positive finite time Lyapunov exponents, thus more likely to stray into
the undesirable deep mixing regime. Therefore, even without heeding
limitations on computational resources, wider networks should be made
shallower.

  * •

When we renormalize the coupling strengths so the weighting matrix becomes
more like a probability with weights summing up to a constant (e.g.,
$\sum_{\alpha}K^{\alpha}_{\beta}=1$), the situation changes. The Frobenius
norm now scales as $D^{0}$, and the average singular value now must scale as
$1/\sqrt{D}$, so the dynamics become less chaotic as the width of the layers
increases. With this approach, one could thus simultaneous increase both the
width and the depth of the network, without degrading performance. Although
there doesn’t seem to be a strong incentive to doing so, given the higher
drain on computational resources.

  * •

Another way to curtail chaos while increasing $D$ is to make the neurons in
the wider network more sparsely connected (setting a higher proportion of
entries in each row of $K$ to zero). A particularly interesting physical
interpretation that is relevant for a subset of such a strategy is related to
path dependence, which by definition has a tendency to help retain reliance on
initial data, thereby lessen mixing. Given a dynamical system, path dependence
can be folded into a particular type of higher dimensionality, with many new
auxiliary state variables that only depend on one other variable, since they
are supposed to just passively record the past states of that variable and
carry them forward without modifications. For example, if $y^{[j+1]\alpha}$
depend not only on $y^{[j]\beta}$ but also on $y^{[j-1]\gamma}$, then we can
define an additional set of variables $x^{[j]\beta}$ and simply let them be
updated by $x^{[j]\alpha}=y^{[j-1]\alpha}$ (note each $x$ node only depends on
one $y$ node in the previous step). This way, the $j+1$th step of the
$2D$-dimensional $y\oplus x$ combined system depends only on its step $j$
state, and path dependence is formally removed so a dimensionally expanded
version of Eq. (1) remains valid. Mapping such a dynamical system into a
neural network, we see that the width of the network can be increased in order
to simulate the effect of path dependence, so long as the newly added
dimensions don’t link up with too many of the existing ones. This could
explain why pruning networks sometimes helps when the network is otherwise
mixing (suffers from e.g., overfitting problems).

##  III Conclusion

In this brief note, we advocated for the enlistment of ergodic theory to help
intuit behaviors of deep neural networks. In particular, we argued that a
highly effective deep neural network would likely operate on the edge of
chaos. This is the easiest to intuit in the case of a classification problem.
The corresponding dynamical system evolution of the neural network rearranges
the data into orderly and well-segregated tiles (the specifics of the tiling
depends on the hypothesis function choice), each representing a class, and so
given any input data to be processed, the output will land in one of them. For
this to work, we need the flow to be sufficiently flexible to contort any
_contiguous_ (the problem needs to be interpolatable for training to be
useful) but possibly exotic-looking initial shape representing a class in the
initial data space, into a regular simply hypercubic tile (assuming the
simplest threshold-based hypothesis function). This means we would need
(quasi-)ergodicity, so any initial data point within that shape can eventually
arrive at the desired tile site (tiles are larger than infinitesimal
neighbourhoods of points, thus strictly speaking we don’t need full
ergodicity, and subsequently the prefix “quasi”). On the other hand, we must
also avoid the more strongly chaotic mixing behavior, otherwise the class
shape being convected by the flow will be shredded and thoroughly mixed up
with sibling shapes, making clean well-segregated tiles in output space
impossible. In other words, we need moderation in the wildness of the
dynamical system trajectories, and just-on-the-cusp of chaotic regime seems
ideal.

We then discussed some network architectural traits that may be advantageous
in terms of finding such ergodic-but-not-mixing niche. For our generic
deliberation, we have stayed with qualitative properties and intuitive
guidelines. However, after the implementation and training of an actual neural
network, it should not be prohibitively difficult to compute the actual
numerical values of the finite time Lyapunov exponents, as there exists mature
and efficient numerical routines for computing singular values.

The largest values of these could then serve as a quality control indicator
that informs on whether the result of the training is suitable for
interpolation and extrapolation. This is potentially an important application,
as the tunable parameters of a neural network are so very numerous, giving it
the ability to always fit to any training data we present to it, but real life
applications require the trained network to respond to new situations in a
moderate and controlled manner, rather than jerks around wildly (i.e., we want
to avoid the overfitting problem). Previous approaches to ensuring that this
is the case is largely empirical, by testing the trained network against
additional datasets not included in the training stage. This wastes precious
labelled data, and one can never be sure these tests properly cover all
plausible new situations. The finite time Lyapunov exponents and their ergodic
theory significances could possible provide a valuable alternative.

One could of course also compute the finite time Lyapunov exponents for the
purpose of debugging. For example, if it turns out that the trained network
lacks expressivity, then one can imagine that the problem might be that the
sum $\sum_{\alpha}\lambda^{[N-1]\alpha}_{i}$ is too negative, so the dynamics
ends up being highly dissipative and collapses onto a fixed point or an
attractor, occupying only a corner of the state space, thereby preventing the
network, as an approximator function, from taking up a chunk of the codomain.

###### Acknowledgements.

This work is supported by the National Natural Science Foundation of China
grants 12073005, 12021003.

## References

  * Pascanu et al. (2014) R. Pascanu, Y. N. Dauphin, S. Ganguli, and Y. Bengio, arXiv e-prints arXiv:1405.4604 (2014), eprint 1405.4604.
  * Bartlett et al. (2020) P. L. Bartlett, P. M. Long, G. Lugosi, and A. Tsigler, Proceedings of the National Academy of Science 117, 30063 (2020), eprint 1906.11300.
  * Sejnowski (2020) T. J. Sejnowski, Proceedings of the National Academy of Science 117, 30033 (2020), eprint 2002.04806.
  * Hornik (1991) K. Hornik, Neural Networks 4, 251 (1991), ISSN 0893-6080, URL https://www.sciencedirect.com/science/article/pii/089360809190009T.
  * Elbrächter et al. (2019) D. Elbrächter, D. Perekrestenko, P. Grohs, and H. Bölcskei, arXiv e-prints arXiv:1901.02220 (2019), eprint 1901.02220.
  * Furuhata et al. (2021) G. Furuhata, T. Niiyama, and S. Sunada, Physical Review Applied 15, 034092 (2021), eprint 2012.08761.
  * Belozyorov and Volkova (2022) V. Y. Belozyorov and S. A. Volkova, International Journal of Bifurcation and Chaos 32, 2250218-1583 (2022).
  * Haber and Ruthotto (2018) E. Haber and L. Ruthotto, Inverse Problems 34, 014004 (2018), eprint 1705.03341.
  * Frigg et al. (2020) R. Frigg, J. Berkovitz, and F. Kronz, in _The Stanford Encyclopedia of Philosophy_ , edited by E. N. Zalta (Metaphysics Research Lab, Stanford University, 2020), Fall 2020 ed.
  * Werndl (2009) C. Werndl, The British Journal for the Philosophy of Science 60, 195 (2009), eprint https://doi.org/10.1093/bjps/axn053, URL https://doi.org/10.1093/bjps/axn053.
  * Benning et al. (2019) M. Benning, E. Celledoni, M. J. Ehrhardt, B. Owren, and C.-B. Schönlieb, arXiv e-prints arXiv:1904.05657 (2019), eprint 1904.05657.
  * Golub et al. (2000) G. Golub, K. Solna, and P. Van Dooren, SIAM Journal on Matrix Analysis and Applications 22, 1 (2000), eprint https://doi.org/10.1137/S0895479897325578, URL https://doi.org/10.1137/S0895479897325578.
  * Demmel and Kahan (1990) J. Demmel and W. Kahan, SIAM Journal on Scientific and Statistical Computing 11, 873 (1990), eprint https://doi.org/10.1137/0911052, URL https://doi.org/10.1137/0911052.
  * Liu and Zaidi (2016) N. Liu and N. A. Zaidi, in _AI 2016: Advances in Artificial Intelligence_ , edited by B. H. Kang and Q. Bai (Springer International Publishing, Cham, 2016), pp. 535–541, ISBN 978-3-319-50127-7.
  * Melby et al. (2000) P. Melby, J. Kaidel, N. Weber, and A. Hübler, Phys. Rev. Lett. 84, 5991 (2000), eprint nlin/0007006.