# Variational theory combining number-projected BCS and coupled-cluster
doubles

V.V. Baran1,2<EMAIL_ADDRESS>J. Dukelsky3<EMAIL_ADDRESS>1
Faculty of Physics, University of Bucharest, 405 Atomiştilor, POB MG-11,
Bucharest-Măgurele, RO-077125, Romania
2”Horia Hulubei” National Institute of Physics and Nuclear Engineering,
30 Reactorului, RO-077125, Bucharest-Măgurele, Romania
3 Instituto de Estructura de la Materia, CSIC, Serrano 123, E-28006 Madrid,
Spain

###### Abstract

The ground state pairing correlations in finite fermionic systems are
described with a high degree of accuracy within a variational approach based
on a combined coupled-cluster and particle-number-projected BCS ansatz. The
flexibility of this symmetry-preserving wavefunction enables a unified picture
valid from weak to strong coupling, both in small and large systems. The
present variational approach consistently yields an energy upper bound while
operating at the same level of precision of the non-variational particle-
number projected Bogoliubov-coupled-cluster theory [Phys. Rev. C 99, 044301
(2019)].

##  I Introduction

Pairing Hamiltonians are ubiquitous in quantum many-body physics. Starting
from their variational treatment in the microscopic theory of
superconductivity given by Bardeen, Cooper and Schrieffer (BCS) [1], they were
soon exported to nuclear physics for the descripton of the large gaps observed
in even-even nuclei [2]. However, the violation of particle number in the BCS
theory, which is negligible for macroscopic systems, represents a major
drawback when applied to finite systems. Therefore, techniques to implement
number projection on top of the BCS wavefunction (PBCS) were developed in
nuclear structure [3, 4] and more recently in quantum chemistry [5, 6] where
the PBCS wavefunction is known as the antisymmetrized geminal power (AGP).
PBCS improves over the BCS theory in finite systems, specially in the strong
coupling limit where superconductivity is well established, but it still fails
in the weak coupling limit dominated by pairing fluctuations, and all along
the transitional region [7, 8]. This fact was made evident in ultrasmall
superconducting grains, where PBCS predicted an abrupt metal-superconductor
transition as a function of the grain size [9] while the exact solution showed
a smooth crossover dominated by large fluctuations [10]. It was precisely in
the field of ultrasmall superconducting grains that the exact solution of the
constant pairing Hamiltonian given by Richardson in the sixties [11] was
recovered [12] and intensively used as a natural benchmark model for
superconducting theories beyond BCS [13, 14, 15, 16, 17, 18, 19, 20, 21].

In the extreme weak coupling limit pair coupled-cluster doubles (pCCD)
describes correctly the pairing fluctuating regime but it quickly overbinds
due to the non-variational character of the theory based on the left
projection to a Hartree-Fock Slater determinant [22, 23, 8]. It seems,
therefore, that a combination of pCCD and quasiparticle BCS would be able to
approach both limits correctly. Indeed, the extension of pCCD to BCS
quasiparticles, the BCS-CCD method [15], gave the correct behavior in the weak
coupling limit, but still suffers from large deviations ($\thicksim$ 10 % for
sizes of $\thicksim$ 100 particles) across the transition region. There were
attempts to interpolate between pCCD and PBCS [16, 8] or to diagonalize the
pairing Hamiltonian in a subspace defined by the reference PBCS state and
different two and four number projected quasiparticle states [24]. Other
possible ways to add correlations to the PBCS state were explored in [25, 26].

Perhaps, the most successful theory beyond BCS-CCD is simply its number
projected version coined as particle-number projected BCS coupled-cluster
doubles (PBCS-CCD) theory [27]. The theory is not Ritz variational and
therefore, it cannot assure an upper bound for the ground state energy.
However it has an affordable computational cost that scales polynomially with
the system size. Moreover, it gives excellent numerical results both in the
weak and strong coupling limits, as well as in the transitional region.

The aim of our paper is to design a Ritz variational method with a trial
wavefunction that combines pCCD and PBCS and produces numerical results with
the same level of precision of PBCS-CCD and with a similar computational cost.
We will benchmark our variational theory with the exact solution of the
Richardson model [28] and with PBCS-CCD results [27] where available.

##  II Theoretical background

###  II.1 Pairing Hamiltonian

We consider the generic pairing Hamiltonian

$H=\sum_{i=1}^{L}\epsilon_{i}\,(c^{\dagger}_{i}c_{i}+c^{\dagger}_{\bar{i}}c_{\bar{i}})+\sum_{i,j=1}^{L}V_{i,j}c^{\dagger}_{i}c^{\dagger}_{\bar{i}}c_{\bar{j}}c_{j}~{},$
(1)

where $i$ and $\bar{i}$ indicate one of the $L$ pairs of conjugated degenerate
single particle levels with energy $\epsilon_{i}=\epsilon_{\bar{i}}$.

The Hamiltonian (1) preserves seniority, and for simplicity we restrict
ourselves to the seniority zero ($v=0$) subspace . Then, in the absence of the
interaction term the ground state is given by the Hartree-Fock product state

$|\text{HF}\rangle=\prod_{i=1}^{M}c^{\dagger}_{i}c^{\dagger}_{\bar{i}}\,|0\rangle~{},$
(2)

where $M$ is the number of pairs in the system.

In preparation for the discussion that follows, we pass to the particle-hole
(ph) representation and denote the particle levels with $p>M$ and the hole
levels with $h\leq M$. We define the particle and hole pair and number
operators

$\displaystyle P^{\dagger}_{p}$
$\displaystyle=c^{\dagger}_{p}c^{\dagger}_{\bar{p}}~{},~{}P_{p}=c_{\bar{p}}c_{p}~{},~{}N_{p}=c^{\dagger}_{p}c_{p}+c^{\dagger}_{\bar{p}}c_{\bar{p}}~{},$
(3) $\displaystyle P^{\dagger}_{h}$
$\displaystyle=c_{\bar{h}}c_{h}~{},~{}P_{h}=c^{\dagger}_{h}c^{\dagger}_{\bar{h}}~{},~{}N_{h}=2-c^{\dagger}_{h}c_{h}-c^{\dagger}_{\bar{h}}c_{\bar{h}}~{},$

such that the Hartree-Fock state (2) is the vacuum to the $P$ and $N$
operators, $P_{p}|\text{HF}\rangle=P_{h}|\text{HF}\rangle=0$,
$N_{p}|\text{HF}\rangle=N_{h}|\text{HF}\rangle=0$. The pairing Hamiltonian (1)
is then expressed as

$\displaystyle H=$
$\displaystyle~{}E_{\text{HF}}+\sum_{p=M+1}^{L}\epsilon_{p}\,N_{p}+\sum_{h=1}^{M}(-\epsilon_{h}-V_{hh})\,{N}_{h}$
(4)
$\displaystyle+\sum_{p,p^{\prime}=M+1}^{L}V_{pp^{\prime}}P^{\dagger}_{p}P_{p^{\prime}}+\sum_{h,h^{\prime}=1}^{M}V_{hh^{\prime}}{P}^{\dagger}_{h}{P}_{h^{\prime}}$
$\displaystyle+\sum_{p=M+1}^{L}\sum_{h=1}^{M}V_{ph}\left(P^{\dagger}_{p}{P}^{\dagger}_{h}+{P}_{h}P_{p}\right)~{},$

with
$E_{\text{HF}}=\langle\text{HF}|H|\text{HF}\rangle=\sum_{h=1}^{M}\left(2\epsilon_{h}+V_{hh}\right)$
being the energy of the Hartree-Fock state (2).

###  II.2 Mean-field theory and its symmetry restoration

The standard description of the pairing correlations induced by the
Hamiltonian (1) is given within the BCS approximation [1] in terms of the pair
condensate

$|\text{BCS}\rangle=\exp[\Gamma^{\dagger}(x)]|0\rangle,\quad\Gamma^{\dagger}(x)\equiv\sum_{i=1}^{L}x_{i}c^{\dagger}_{i}c^{\dagger}_{\bar{i}}~{},$
(5)

which explicitly breaks the $U(1)$ gauge symmetry associated with particle
number conservation. For macroscopic systems the symmetry broken picture is
exact and the particle number fluctuations are negligible. However, for finite
systems like atomic nuclei or small superconducting grains one speaks only of
obscured or emergent symmetry breaking [29, 30, 31, 32], in which case the
quantum fluctuations inevitably lift any degeneracy associated with the broken
symmetry.

Much effort is thus devoted to restore the symmetry of the mean-field ansatz
with the help of projection techniques [33, 34, 35, 36, 37, 38, 39]. In the
BCS case upon particle number restoration we obtain the so-called number-
projected BCS (PBCS) [3] or antisymmetrized geminal power (AGP) in the context
of quantum chemistry [40]

$\displaystyle|\text{PBCS}(x)\rangle$
$\displaystyle=\mathcal{P}_{M}|\text{BCS}\rangle$ (6)
$\displaystyle=\frac{1}{M!}[\Gamma^{\dagger}(x)]^{M}|0\rangle$
$\displaystyle=\frac{1}{2\pi}\int_{0}^{2\pi}\text{d}\theta~{}e^{-i\theta
M}\exp[\Gamma^{\dagger}(e^{i\theta}x)]~{},$

where $\mathcal{P}_{M}$ is the projector onto the state of $M$ pairs. While
PBCS describes well the properties of superfluid nuclei with a small number of
valence nucleons [36], it cannot account for the weak pairing correlations
that develop within larger spaces, e.g. those considered in the large-scale
energy density functional treatments of finite nuclei or in the study of small
superconducting grains [41]. For a working description of the weak pairing
regime one usually turns to RPA [22] or coupled-cluster [15, 23, 8]
approaches. Generalizations of the PBCS ansatz have also been considered based
on its structural similarity with a particular coupled-cluster ansatz [8].
Specifically, the PBCS representation in the ph-basis is obtained as

$\displaystyle|\mathrm{PBCS}\rangle\propto\sum_{\ell=0}^{M}\frac{1}{\ell!^{2}}[\Gamma_{P}^{\dagger}(x)\Gamma^{\dagger}_{H}(1/x)]^{\ell}\,|\text{HF}\rangle$
(7)
$\displaystyle=\frac{1}{2\pi}\int_{0}^{2\pi}\text{d}\theta~{}\exp[\Gamma^{\dagger}_{P}(e^{i\theta}x)]~{}\exp[\Gamma^{\dagger}_{H}(e^{-i\theta}/x)]\,|\text{HF}\rangle~{},$

in terms of the particle and hole components of the collective pairs

$\Gamma_{P}^{\dagger}(x)=\sum_{p=M+1}^{L}x_{p}P_{p}^{\dagger},\quad\Gamma_{H}^{\dagger}(x)=\sum_{h=1}^{M}x_{h}P_{h}^{\dagger}~{}.$
(8)

The structure of the PBCS state is then defined by the inverse squared
factorials appearing as expansion coefficients in the collective ph-pair
basis.

###  II.3 Coupled-cluster theory

Analogously, by using a slightly modified expansion involving plain factorials
one obtains a separable pair coupled-cluster doubles variational ansatz
(vCCD${}_{\text{sep}}$)

$\displaystyle|\text{vCCD}_{\text{sep}}\rangle$
$\displaystyle=\sum_{\ell=0}^{M}\frac{1}{\ell!}[\Gamma_{P}^{\dagger}(y)\Gamma^{\dagger}_{H}(y)]^{\ell}\,|\text{HF}\rangle$
(9)
$\displaystyle=\exp[\Gamma_{P}^{\dagger}(y)\Gamma^{\dagger}_{H}(y)]\,|\text{HF}\rangle$
$\displaystyle=\exp\left[\sum_{p,h}y_{p}y_{h}P^{\dagger}_{p}P^{\dagger}_{h}\right]\,|\text{HF}\rangle~{}.$

This is a particular case of a pair coupled-cluster doubles variational
wavefunction (vCCD)

$|\text{vCCD}\rangle=\exp\left[\sum_{p,h}z_{ph}P^{\dagger}_{p}P^{\dagger}_{h}\right]\,|\text{HF}\rangle~{},$
(10)

involving the most general double excitations that do not break pairs through
a fully non-separable structure matrix $z_{ph}$ [22]. The separable case
vCCD${}_{\text{sep}}$ is thus recovered for $z_{ph}=y_{p}\,y_{h}$.

On the one hand, the full freedom in the structure matrix of the CCD ansatz
allows for an excellent description of the weak pairing regime (see also Fig.
1 below), due to the additional four-body correlations accounted for relative
to the separable case. Note, for example, that the Richardson solution
involving complex-conjugated pairs may be expressed exactly as a product of
four-body quartet structures [42]. On the other hand, the computational
complexity of evaluating operator matrix elements exactly in the non-separable
case grows exponentially with the size of the system. The usual approximation
of coupled cluster involves a left projection onto the subspaces of zero and
two ph pairs, thus breaking the Ritz variational principle. Furthermore, even
with a fully nonseparable structure matrix, the validity of the vCCD ansatz
breaks down around the critical value of the pairing strenght (see also Fig. 1
below). The choice of a BCS mean-field reference state does improve on this
aspect at the cost of effectively breaking the particle number symmetry [15].

Given the success of vCCD in the weak pairing regime and that of PBCS in the
strong pairing regime, it is then natural to combine them for a precise
unified description of all regimes.

###  II.4 Combining coupled-cluster and symmetry-restored mean-field theories

The symmetry restoration of broken-symmetry coupled-cluster theories has been
only recently considered [43, 44, 27, 45]. For the schematic pairing
Hamiltonian with the breaking and restoration of the particle number symmetry,
the so-called particle-number projected Bogoliubov-coupled-cluster theory
yields highly accurate results [18]. In this context, a set of differential
equations is set up for obtaining the gauge-angle-dependent excitation
operator. The practical need for truncating this set of ODEs implies however
an approximate action of the projection operator and a violation of the Ritz
variational principle.

It is the purpose of this work to explore a physically transparent alternative
in the form of a variational CCD-PBCS combined approach. Ideally the ground
state would involve the CCD excitations built directly on top of the particle-
number-projected BCS as

$|\text{vCCD-
PBCS}\rangle=\exp\left[\sum_{p,h}z_{ph}P^{\dagger}_{p}P^{\dagger}_{h}\right]|\text{PBCS}(x)\rangle~{},$
(11)

where the $L$ mixing amplitudes $x_{i}$ of the PBCS state and the $M(L-M)$
structure matrix elements $z_{ph}$ of the CCD excitations are to be treated as
free variational parameters. Indeed, the energy minimization procedure

$E_{\text{gs}}=\text{min}_{\psi}\frac{\langle\psi|H|\psi\rangle}{\langle\psi|\psi\rangle}$
(12)

yields extremely precise results for the few small systems where this ansatz
may actually be applied.

Computational access to larger systems is enabled upon various
simplifications. Let us consider first the separability approximation
$z_{ph}=y_{p}y_{h}$ for the structure matrix of the CCD excitations in Eq.
(11), yielding the $\text{vCCD}_{\text{sep}}\text{PBCS}$ ansatz

$\displaystyle|\text{vCCD}_{\text{sep}}\text{PBCS}\rangle$
$\displaystyle=\exp[\Gamma_{P}^{\dagger}(y)\Gamma^{\dagger}_{H}(y)]\,|\text{PBCS}(x)\rangle~{},$
(13)

which involves $2L$ independent variational parameters $x_{i},y_{i}$.
Remarkably, while the PBCS ansatz (6) and the $\text{vCCD}_{\text{sep}}$
ansatz (9) each individually fails in the weak pairing regime due to their
separable structure matrices, the combined
$\text{vCCD}_{\text{sep}}\text{PBCS}$ wavefunction will turn out to be quite
accurate. This may be easily understood by noting that for weak pairing the
above ansatz effectively involves a non-separable structure matrix in the form

$|\text{vCCD}_{\text{sep}}\text{PBCS}\rangle\approx[1+\sum_{p,h}(x_{p}x_{h}+y_{p}y_{h})P^{\dagger}_{p}P^{\dagger}_{h}]|\text{HF}\rangle,$
(14)

obtained after taking into account the PBCS (6) and $\text{vCCD}_{\text{sep}}$
(9) expansions and redefining $x_{h}\rightarrow 1/x_{h}$ in (7). The accuracy
of this doubly-separable $2L$-dimensional parametrization of the fully non-
separable $M(L-M)$-dimensional structure matrix naturally degrades with
increasing system size, the actual rate being numerically determined in the
next section.

An alternative approximation scheme would then involve a fully non-separable
CCD excitations limited to a relatively small finite window around the Fermi
level. We thus consider the combination

$\displaystyle|\text{vCCD}^{(w)}\text{PBCS}(x,z)\rangle=\text{vCCD}^{(w)}(z)|\text{PBCS}(x)\rangle$
(15)

with

$\displaystyle\text{vCCD}^{(w)}(z)=\exp\left[\sum_{p=M+1}^{M+w/2}\sum_{h=M+1-w/2}^{M}z_{ph}P^{\dagger}_{p}P^{\dagger}_{h}\right]~{},$
(16)

where $w$ denotes the size of the truncation window. As will be detailed in
the next section, the quality of the results obtained within this approach
will turn out to be inferior to that of the above
$\text{vCCD}_{\text{sep}}\text{PBCS}$ ansatz of Eq. (13) involving two sets of
global parameters.

Overall the optimal compromise between the computational complexity and the
accuracy of the results is found by combining the above two approximation
schemes (9) and (15) into the ansatz

$\displaystyle|\text{vCCD}^{(w)}_{\text{sep}}\text{PBCS}(x,y,z)\rangle=\text{vCCD}^{(w)}_{\text{sep}}(y,z)|\text{PBCS}(x)\rangle~{},$
(17)
$\displaystyle\text{vCCD}^{(w)}_{\text{sep}}=\text{vCCD}^{(w)}(z)\,\text{vCCD}_{\text{sep}}(y)~{},$

involving a total of $2L+(w/2)^{2}$ free variational parameters and leading to
highly accurate results comparable to those of Ref. [18], to be discussed
later on. Next, we shortly review the actual computational strategy for the
above mentioned wavefunctions.

###  II.5 Computational aspects

In this section we propose a novel efficient algorithm for the evaluation of
expectation values on the combined $\text{vCCD}^{(w)}_{\text{sep}}\text{PBCS}$
wavefunction (17). It is based on a representation of the PBCS (7),
vCCD${}_{\text{sep}}$ (19) and vCCD(w) (16) terms as disentangled particle and
hole gauge-angle-rotated BCS states.

We start from the discrete exact representation for the particle-number
projection operation and decompose the PBCS ansatz using the definitions of
Eqs. (7) and (8)

$\displaystyle|\text{PBCS}\rangle=\sum_{n=0}^{L}\exp[\Gamma^{\dagger}_{H}(e^{-i\theta_{n}}x)]\exp[\Gamma^{\dagger}_{P}(e^{i\theta_{n}}x)]\,|\text{HF}\rangle~{},$
(18)

where $\theta_{n}=2\pi n/(L+1)$. Note that we neglect the irrelevant constant
normalization factors throughout this section.

The collective pairs (for the particle and hole subspaces) appearing in the
vCCD${}_{\text{sep}}$ operator

$\text{vCCD}_{\text{sep}}=\sum_{\ell=0}^{M}\frac{1}{\ell!}[\Gamma_{P}^{\dagger}(y)\Gamma^{\dagger}_{H}(y)]^{\ell}$
(19)

are also expanded as superpositions of gauge-angle-rotated BCS operators

$\displaystyle\frac{[\Gamma_{P}^{\dagger}(y)]^{\ell}}{\ell!}$
$\displaystyle=\sum_{k=0}^{L-M}e^{-i\phi_{k}\ell}\exp[\Gamma_{P}^{\dagger}(e^{i\phi_{k}}y)]~{},$
(20) $\displaystyle\frac{[\Gamma_{H}^{\dagger}(y)]^{\ell}}{\ell!}$
$\displaystyle=\sum_{k=0}^{M}e^{-i\varphi_{k}\ell}\exp[\Gamma_{H}^{\dagger}(e^{i\varphi_{k}}y)]~{},$

with $\phi_{k}=2\pi k/(L-M+1)$ and $\varphi_{k}=2\pi k/(M+1)$. Finally we
obtain the representation

$\displaystyle|\text{vCCD}$
$\displaystyle{}_{\text{sep}}\text{PBCS}\rangle=\sum_{k_{p}=0}^{L-M}\sum_{k_{h}=0}^{M}\sum_{n=0}^{L}g_{k_{p}k_{h}}\,\times$
(21)
$\displaystyle\text{BCS}_{P}(k_{p},n)\,\text{BCS}_{H}(k_{h},n)\,|\text{HF}\rangle$

in terms of the gauge-angle-rotated BCS operators

$\displaystyle\text{BCS}_{P}(k_{p},n)$
$\displaystyle=\exp[\Gamma_{P}^{\dagger}(e^{i\theta_{n}}x+e^{i\phi_{k_{p}}}y)]~{},$
(22) $\displaystyle\text{BCS}_{H}(k_{h},n)$
$\displaystyle=\exp[\Gamma_{H}^{\dagger}(e^{-i\theta_{n}}x+e^{i\phi_{k_{p}}}y)]~{},$

with expansion coefficients

$g_{k_{p}k_{h}}=\sum_{\ell=0}^{M}\ell!\,\exp[-i(\phi_{k_{p}}+\varphi_{k_{h}})\ell]~{}.$
(23)

We are interested in the expectation values of generic particle-number-
conserving ph-factorized operators
$\mathcal{O}=\mathcal{O}^{(P)}\mathcal{O}^{(H)}$ on the
$\text{vCCD}_{\text{sep}}\text{PBCS}$ state (21). As the particle-number
projection operation needs only to be performed once on the mean-field BCS
wavefunction, we obtain

$\displaystyle\langle\mathcal{O}\rangle=\,\langle\text{vCCD}_{\text{sep}}\text{BCS}|\mathcal{O}^{(P)}\mathcal{O}^{(H)}|\text{vCCD}_{\text{sep}}\text{PBCS}\rangle$
(24)
$\displaystyle=\sum_{k_{p},k_{p}^{\prime}=0}^{L-M}\sum_{k_{h},k_{h}^{\prime}=0}^{M}g^{*}_{k_{p}k_{h}}\,g_{k_{p}^{\prime}k_{h}^{\prime}}\sum_{n=0}^{L}{O}^{(P)}_{k_{p},k_{p}^{\prime},n}\,\mathcal{O}^{(H)}_{k_{h},k_{h}^{\prime},n}~{},$

involving the matrix elements between particle and hole BCS states (22)

$\displaystyle{O}^{(P)}_{k_{p},k_{p}^{\prime},n}$
$\displaystyle=\langle\text{BCS}_{P}(k_{p},0)|{O}^{(P)}|\text{BCS}_{P}(k_{p}^{\prime},n)\rangle~{},$
(25) $\displaystyle\mathcal{O}^{(H)}_{k_{h},k_{h}^{\prime},n}$
$\displaystyle=\langle\text{BCS}_{H}(k_{h},0)|{O}^{(H)}|\text{BCS}_{H}(k_{h}^{\prime},n)\rangle~{}.$

By considering all terms in Eq. (4), the expression (24) may be directly
employed to compute the energy expectation value $E=\langle H\rangle/\langle
I\rangle$ on the vCCD${}_{{\text{sep}}}$PBCS wavefunction.

To include the effects of the non-separable vCCD(w) excitations of Eq. (17) we
use a Hubbard-Stratonovich transformation [46] and pass to a disentangled BCS
representation of the vCCD operator

$\displaystyle\text{vCCD}(z)=\exp\left(\sum_{p,h}z_{ph}P^{\dagger}_{p}P^{\dagger}_{h}\right)$
(26)
$\displaystyle=\int\text{d}^{L}\xi\,\exp\left(-\frac{1}{2}\,{\xi}^{T}\,Z^{-1}\,{\xi}\right)\,\exp[\Gamma^{\dagger}_{P}(\xi)+\Gamma^{\dagger}_{H}(\xi)],$

with

$Z=\begin{pmatrix}0&z\\\ z^{T}&0\end{pmatrix}~{}.$ (27)

This BCS representation allows us to generalize the above Eq. (24) to the
expectation values of ph-factorized operators on the
vCCD${}^{(w)}_{\text{sep}}$PBCS wavefunction (17). In practice the matrix
elements on the resulting BCS wavefunctions which generalize Eq. (25) are to
be treated as polynomials in the integration variables $\xi$. Only specific
terms corresponding to the various nonzero Wick contractions are to be
selected according to

$\displaystyle\int\text{d}^{L}\xi$
$\displaystyle\,\exp\left(-\frac{1}{2}\,\xi^{T}\,Z^{-1}\,\xi\right)\,\xi_{i}\xi_{j}\dots\xi_{k}\xi_{l}=$
(28) $\displaystyle\sum_{\text{Wick }}Z_{ab}\dots Z_{cd}~{},$

the set of indices $\\{a,b,...,c,d\\}$ representing a permutation of
$\\{i,j,...,k,l\\}$ [47].

With the non-separable CCD excitations restricted to a small window of $w$
levels around the Fermi level as in Eq. (15), the matrix elements of the
relevant operators on generic vCCD(w)BCS states may be computed analytically
(the notebook in Cadabra2 [48] is available upon request from the authors).
Their subsequent coupling to the separable subspace (i.e. the complement of
the $w$-level window) leads to the final vCCD${}^{(w)}_{\text{sep}}$PBCS
results presented in the next section.

##  III Numerical results

In this section we benchmark the variational calculations for the various
wavefunctions presented above against the exact solution [28] for a constant
pairing Hamiltonian ($V_{i,j}=-G$ in Eq. (1)) of equally spaced single
particle levels with energies $\epsilon_{k}=k\epsilon,\,k=1,\dots,L$. In order
to avoid the singularities in the equations that solve exactly the picket
fence model we use the same algorithm proposed by Richardson in [28].

The numerical code used to compute the expectation value of the Hamiltonian
(4) on the various derivatives of the vCCD-PBCS wavefunction is freely
available upon request from the authors. The minimization procedure for the
energy function (12) is performed using the e04ucf routine of the NAG library
[49].

Figure 1: Error in the correlation energy (29) relative to its exact value (in
percentages) versus the pairing strength $G$ (in units of the level spacing
$\epsilon$) for the fully variational vCCD (10), PBCS (6), and
vCCD${}_{\text{sep}}$PBCS (13) wavefunctions in the case of an $L=12$ level
system at half filling. The “vCCD${}_{\text{sep}}$+PBCS” procedure indicated
by a red dashed line involves an energy minimization with respect to the
vCCD${}_{\text{sep}}$ amplitudes in the presence of the fixed optimal PBCS
reference. Within the “PBCS+vCCD${}_{\text{sep}}$” procedure indicated by a
green dashed line the PBCS amplitudes are varied in the presence of the fixed
optimal vCCD${}_{\text{sep}}$ structure.

We present in Figure (1) the errors for the correlation energy

$E_{c}=\frac{\langle\psi|H|\psi\rangle}{{\langle\psi|\psi\rangle}}-E_{\text{HF}}~{},$
(29)

relative to its exact value for $L=12$ at half filling. On the one hand,
notice how PBCS (6) and CCD${}_{\text{sep}}$ (9) are limited in the weak
pairing regime by their common separable structure, both reducing to
$(1+\Gamma^{\dagger}_{P}\Gamma^{\dagger}_{H})|\text{HF}\rangle$ at
$G/\epsilon\ll 1$. On the other hand, the full generality of the structure
matrix of the vCCD wavefunction (10) allows it to capture precisely all the
correlations in this regime. However, beyond the critical value
$G_{\text{cr}}$ of the HF to BCS transition the vCCD ansatz quickly loses its
ability to describe the stronger correlations, becoming indistinguishable from
its separable version vCCD${}_{\text{sep}}$ for $G/\epsilon>1$ (see also the
discussion around Fig. 4 of Ref. [8]).

As indicated by the thick green line next to the horizontal axis of Fig. (1),
the combined vCCD${}_{\text{sep}}$PBCS wavefunction (13) leads to very
accurate energetics across all regimes with the relative errors in the
correlation energy not exceeding $3\cdot 10^{-4}$. Note however that it is
essential to enable all parameters to vary freely in the minimization process
as to retain the full flexibility of the wavefunction and thus recover all
available dynamical correlations. Indeed, by considering a variational
CCD${}_{\text{sep}}$ on top of the frozen optimal PBCS reference (or
viceversa) we find significant improvements only in the weak pairing regime,
as indicated by the dashed lines in Fig. (1). In particular, by varying in
this way just the PBCS amplitudes within vCCD${}_{\text{sep}}$PBCS while
keeping fixed the optimal CCD${}_{\text{sep}}$ structure (computed beforehand)
we find no improvement over the full vCCD results for $G>G_{\text{cr}}$. At
$G/\epsilon\ll 1$, the results differ due to the limitations discussed around
Eq. (14).

Figure 2: Error in the correlation energy (29) relative to its exact value (in
percentages) versus the pairing strength $G$ (in units of the level spacing
$\epsilon$) for the fully variational PBCS (6), vCCD${}_{\text{sep}}$PBCS
(13), vCCD(w=6)PBCS (15) and vCCD${}^{(w=4,6)}_{\text{sep}}$PBCS (17)
wavefunctions in the case of an $L=20$ level system at half filling. Figure 3:
Error in the correlation energy (29) relative to its exact value (in
percentages) versus the pairing strength $G$ (in units of the level spacing
$\epsilon$) for the fully variational vCCD${}_{\text{sep}}$PBCS (13) and
vCCD${}^{(w)}_{\text{sep}}$PBCS (17) wavefunctions, in the case of an $L=60$
level system (left panel) and for an $L=100$ level system (right panel) at
half filling. The red squares indicate the results of the particle-number
projected Bogoliubov-coupled-cluster theory of Ref. [18]. Note also that for
$L=100$, the errors for the pure PBCS (6) reach a maximum of about 22% at
$G/\epsilon=0.24$ and decrease down to 6% at very small $G$ (see also Fig. 1
of Ref. [18]).

The relative errors for the correlation energy (29) relative to its exact
value are shown in logarithmic scale in Fig. (2) for a slightly larger system
of $L=20$ levels at half-filling. Note that the vCCD${}_{\text{sep}}$PBCS
errors have slighly increased with respect to the $L=12$ case, but are still
respectably accurate when compared to PBCS. The limitations of the
vCCD${}_{\text{sep}}$PBCS wavefunction in the weak pairing regime originate in
the approximate form of its structure matrix, as remarked around Eq. (14).

This weak pairing behaviour may be improved at a reasonable computational cost
by including additional non-separable CCD excitations on top of
vCCD${}_{\text{sep}}$PBCS. The vCCD${}^{(w)}_{\text{sep}}$PBCS ansatz (17)
already provides more than an order of magnitude lower errors with respect to
vCCD${}_{\text{sep}}$PBCS at weak pairing even for the modest $w=4$. A larger
non-separable CCD excitation window $w=6$ naturally accounts for an additional
amount of correlations, further reducing the errors at weak pairing. Beyond
the critical value $G_{\text{cr}}$ however there is no significant benefit of
the supplementary non-separable excitations. This is also the situation for
the vCCD(w)PBCS ansatz (15) which displays the same strong pairing behaviour
as the plain PBCS. With the Fermi sea being washed out at strong pairing, only
a global deformation of PBCS such as vCCD${}_{\text{sep}}$PBCS is able to
bring substantial improvement to the energetics, as opposed to any local
deformation such as vCCD(w)PBCS (15).

The numerical values for the optimal non-separable structure matrix elements
$z_{ph}$ are typically found to be very small, of the order $10^{-2}-10^{-3}$
upon the full vCCD${}^{(w)}_{\text{sep}}$PBCS energy minimization. In
practice, within the vCCD${}^{(w)}_{\text{sep}}$PBCS approach we chose to
limit the non-separable CCD excitations to linear order

$\displaystyle\text{vCCD}^{(w)}(z)=1+\sum_{p=M+1}^{M+w/2}\sum_{h=M+1-w/2}^{M}z_{ph}P^{\dagger}_{p}P^{\dagger}_{h}~{},$
(30)

which allows for a decrease in the computational complexity without any
noticeable loss in precision with respect with the full form of Eq. (15).

An additional computational speed-up is enabled by limiting the action of the
particle number projection operations in Eqs. (18) and (20) to a reduced
subspace of particle and hole collective pairs. For the moderate values of the
pairing strength $G$ considered here (2 to 3 times the value of $G_{cr}$) this
still allows for an exact particle-number conservation (within the numerical
accuracy). This key computational aspect is ensured by the negligible
contributions of the high order terms within the PBCS (6) and
vCCD${}_{\text{sep}}$ (19) expansions, due to the small (subunitary) numerical
values of their corresponding collective pair amplitudes. All results
presented below for $L=60$ and $L=100$ were obtained upon projecting the
particle-number only within a $\mathcal{C}_{1}=15$ ph-pair subspace for PBCS
in Eq. (18) and within a $\mathcal{C}_{2}=10$ ph-pair subspace for
CCD${}_{\text{sep}}$ in Eq. (20). The gauge angles were adjusted accordingly,
i.e. $\theta_{n}=2\pi n/(\mathcal{C}_{1}+1)$, $\phi_{k}=\varphi_{k}=2\pi
k/(\mathcal{C}_{2}+1)$. We also note that for all considered systems it was
sufficient to restrict the amount of CCD${}_{\text{sep}}$ excitations by
imposing an upper bound at $\ell_{\text{max}}=7$ in the $\ell$-sum of Eq.
(23), as to avoid the numerical errors originating from the combination of
large factorials and rapidly oscillating phases and, at the same time, but
preserving the numerical precision.

We present in Fig. (3) the relative errors in the correlation energy (29)
relative to its exact values for $L=60$ (left panel) and $L=100$ (right panel)
at half filling. The vCCD${}_{\text{sep}}$PBCS wavefunction remains
surprisingly accurate given that its doubly-separable structure matrix (14) at
weak coupling involves only $L$ independent parameters (down from a total of
$2L$ parameters due to the present ph-symmetry) as compared to the fully non-
separable CCD structure matrix that requires $L^{2}/8$ independent parameters
(still for the present ph-symmetric systems). More precisely, the
vCCD${}_{\text{sep}}$PBCS errors at weak coupling are 0.3% and 0.5% for $L=60$
and $L=100$ respectively, increasing slightly until $G\sim G_{\text{cr}}$ and
then rapidly decreasing at stronger couplings.

While the inclusion of non-separable local CCD excitations within the
vCCD${}^{(w)}_{\text{sep}}$PBCS ansatz (17) significantly improves the error
at weak coupling even for relatively very small values of $w$, perfectly
accurate energetics would still require a set of global ph-excitations on top
of vCCD${}_{\text{sep}}$PBCS. Nevertheless, the overall quality of the
vCCD${}^{(w=6)}_{\text{sep}}$PBCS results is at the level of the more involved
particle-number projected Bogoliubov-coupled-cluster theory of Ref. [18],
while consistently providing an energy upper bound across all regimes.

Figure 4: Error in the correlation energy (29) relative to its exact value (in
percentages) versus the filling fraction $M/L$ of an $L=100$ level system for
the fully variational PBCS (6) and vCCD${}_{\text{sep}}$PBCS (13)
wavefunctions. The pairing strength is $G=1.5G_{\text{cr}}$, with
$G_{\text{cr}}$ the critical strength at half filling. The blue diamonds
indicate the results of the particle-number projected Bogoliubov-coupled-
cluster theory of Ref. [18].

The situation is similar away from half-filling, as shown in Fig. (4). To
provide a comparison with the results of Ref. [18], we consider the $L=100$
case at an interaction strength of $G=1.5G_{\text{cr}}$. For this value the
particle-number symmetry is broken for all filling fractions. Both PBCS and
vCCD${}_{\text{sep}}$PBCS are most accurate for small (or large) filling
fractions, but their errors live on different scales (0.5% vs 20% at half
filling). As seen in Fig. (3), in the considered regime the non-separable
local CCD excitations bring no improvement over the vCCD${}_{\text{sep}}$PBCS
results, which are well matched against those of Ref. [18].

Finally, we consider the behaviour of the canonical gap

$\Delta=G\sum_{i=1}^{L}\sqrt{n_{i}(1-n_{i})}~{},$ (31)

where $n_{i}=\langle
c^{\dagger}_{i}c_{i}+c^{\dagger}_{\bar{i}}c_{\bar{i}}\rangle/2$ indicates the
occupation probability of each level $i$. This quantity exhibits a more
pronounced sensitivity to the structure of the wavefunction than the
correlation energy due to its dependence on the occupation probabilities.

We show in Fig. (5) the error for the canonical gap $\Delta$ relative to its
exact value for $L=60$ at half filling. Improving on both individual PBCS and
vCCD${}_{\text{sep}}$ wavefunctions, their combination
vCCD${}_{\text{sep}}$PBCS only shows visible $\sim 1\%$ errors in the
$G<G_{\text{cr}}$ region. These are further reduced within the
vCCD${}^{(w)}_{\text{sep}}$PBCS approach which exhibits highly accurate
occupations across all regimes.

##  IV SUMMARY AND CONCLUDING REMARKS

In this work, we considered a variational approach for the ground state of
finite paired systems based on a combined coupled-cluster and particle-number-
projected BCS wavefunction. We benchmarked our results against the exact
solution for a picket-fence model involving a pure pairing force acting on a
space of doubly degenerate, equally distanced levels. The analyzed systems
systems range from small ($L=M=12$) to relatively large ($L=M=100$).

We confirmed within the variational context that the combination of symmetry-
restored mean field theory and coupled-cluster theory leads to a wavefunction
that is significantly better than either of the two taken separately, which is
also the main conclusion of Ref. [18]. By incorporating pure four-body
correlations, our vCCD${}^{(w)}_{\text{sep}}$PBCS (17) is able to reproduce
well the physics at weak pairing while also offering orders of magnitude
smaller errors relative to the standard pair-only PBCS (6) across all other
regimes. Our results match the high level of precision of the particle-number
projected Bogoliubov-coupled-cluster theory of Ref. [18], while consistently
providing an energy upper bound.

Computational limitations include the restriction of the pure quartet
correlations to a relatively small window around the Fermi level (discussed in
detail in the main text) and also the need of a restricted space for particle-
number projection. While the latter does not spoil the exact particle-number
conservation (within the numerical accuracy) for moderate values of the
pairing strength, alternative approaches need to be considered for the very
strong pairing regime $G\gg G_{\text{cr}}$. One possibility would involve
limiting the separable CCD excitations to linear order within the
vCCD${}_{\text{sep}}$PBCS approach. Computations with this simpler
wavefunction could then be performed efficiently without resorting to any
approximations for an improvement over the already very good PBCS results for
this regime.

The vCCD${}^{(w)}_{\text{sep}}$PBCS (17) wavefunction was shown to be quite
flexible but it is still affected by structural limitations leading to a
visible maximum in the energy error around $G\sim 1.5G_{\text{cr}}$ in all
analyzed cases. Attempts at mitigating these effects have included treating as
independent variational parameters the factorials in Eq. (23) originating from
the vCCD${}_{\text{sep}}$ expansion (9), with only marginal benefits.

Possible avenues to explore could involve the variational treatment of beyond-
ph CCD excitations on top of the symmetry-restored mean-field state as in

$|\Psi\rangle=\exp\left(\sum_{i,j=1}^{L}z_{ij}\,c^{\dagger}_{i}c^{\dagger}_{\bar{i}}\,c_{\bar{j}}c_{j}\right)|\text{PBCS}\rangle~{},$
(32)

which would act as a more natural choice (albeit more computationally
challenging) for a regime lacking a well-defined Fermi sea. As a first step
one could envision building a multi-separable beyond-ph CCD excitation
operator which could be optimized with the very recently developed methods for
constructing linearly-independent PBCS states [19].

Finally, we leave to future studies extensions of the theory presented in this
work that could incorporate the effect of seniority breaking terms within a
generalized variational wavefunction.

Figure 5: Error in the canonical gap (31) relative to its exact value (in
percentages) versus the pairing strength $G$ (in units of the level spacing
$\epsilon$) for the fully variational PBCS (6), vCCD${}_{\text{sep}}$ (9),
vCCD${}_{\text{sep}}$PBCS (13) and vCCD${}^{(w=6)}_{\text{sep}}$PBCS (17)
wavefunctions in the case of an $L=60$ level system at half filling.

###### Acknowledgements.

This work was supported by a grant of the Romanian Ministry of Education and
Research, CNCS - UEFISCDI, project number PN-III-P1-1.1-PD-2019-0346, within
PNCDI III, and PN-19060101/2019-2022; and by the Spanish Ministerio de Ciencia
e Innovación, and the European regional development fund (FEDER) project Nº
PGC2018-094180-B-I00. We thank the authors of Ref. [18] for sharing with us
their numerical results.

## References

  * Bardeen  _et al._ [1957] J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Theory of superconductivity, Phys. Rev. 108, 1175 (1957).
  * Bohr  _et al._ [1958] A. Bohr, B. R. Mottelson, and D. Pines, Possible analogy between the excitation spectra of nuclei and those of the superconducting metallic state, Phys. Rev. 110, 936 (1958).
  * Dietrich  _et al._ [1964] K. Dietrich, H. J. Mang, and J. H. Pradal, Conservation of particle number in the nuclear pairing model, Phys. Rev. 135, B22 (1964).
  * Sheikh and Ring [2000] J. A. Sheikh and P. Ring, Symmetry projected hartree fock bogoliubov equations, Nuclear Physics A 665, 71 (2000).
  * Scuseria  _et al._ [2011] G. E. Scuseria, C. A. Jimenez-Hoyos, T. M. Henderson, K. Samanta, and J. K. Ellis, Projected quasiparticle theory for molecular electronic structure, The Journal of Chemical Physics 135, 124108 (2011).
  * Khamoshi  _et al._ [2019] A. Khamoshi, T. M. Henderson, and G. E. Scuseria, Efficient evaluation of agp reduced density matrices, The Journal of Chemical Physics 151, 184103 (2019), https://doi.org/10.1063/1.5127850 .
  * Sandulescu and Bertsch [2008] N. Sandulescu and G. F. Bertsch, Accuracy of bcs-based approximations for pairing in small fermi systems, Phys. Rev. C 78, 064318 (2008).
  * Dukelsky  _et al._ [2016] J. Dukelsky, S. Pittel, and C. Esebbag, Structure of the number-projected bcs wave function, Phys. Rev. C 93, 034313 (2016).
  * Braun and von Delft [1998] F. Braun and J. von Delft, Fixed-n superconductivity: The crossover from the bulk to the few-electron limit, Phys. Rev. Lett. 81, 4712 (1998).
  * Dukelsky and Sierra [1999a] J. Dukelsky and G. Sierra, Density matrix renormalization group study of ultrasmall superconducting grains, Phys. Rev. Lett. 83, 172 (1999a).
  * Richardson [1963] R. Richardson, A restricted class of exact eigenstates of the pairing-force hamiltonian, Physics Letters 3, 277 (1963).
  * Sierra  _et al._ [2000] G. Sierra, J. Dukelsky, G. G. Dussel, J. von Delft, and F. Braun, Exact study of the effect of level statistics in ultrasmall superconducting grains, Phys. Rev. B 61, R11890 (2000).
  * Hirsch  _et al._ [2002] J. G. Hirsch, A. Mariano, J. Dukelsky, and P. Schuck, Fully self-consistent rpa description of the many level pairing model, Annals of Physics 296, 187 (2002).
  * Dukelsky  _et al._ [2003a] J. Dukelsky, G. Dussel, J. Hirsch, and P. Schuck, Comparison between exact and approximate treatments of the pairing interaction for finite fermi systems, Nuclear Physics A 714, 63 (2003a).
  * Henderson  _et al._ [2014] T. M. Henderson, G. E. Scuseria, J. Dukelsky, A. Signoracci, and T. Duguet, Quasiparticle coupled cluster theory for pairing interactions, Phys. Rev. C 89, 054305 (2014).
  * Degroote  _et al._ [2016] M. Degroote, T. M. Henderson, J. Zhao, J. Dukelsky, and G. E. Scuseria, Polynomial similarity transformation theory: A smooth interpolation between coupled cluster doubles and projected bcs applied to the reduced bcs hamiltonian, Phys. Rev. B 93, 125124 (2016).
  * Rubio-García  _et al._ [2018] A. Rubio-García, D. R. Alcoba, P. Capuzzi, and J. Dukelsky, Benchmarking the variational reduced density matrix theory in the doubly occupied configuration interaction space with integrable pairing models, Journal of Chemical Theory and Computation 14, 4183 (2018), pMID: 29906104.
  * Qiu  _et al._ [2019] Y. Qiu, T. M. Henderson, T. Duguet, and G. E. Scuseria, Particle-number projected bogoliubov-coupled-cluster theory: Application to the pairing hamiltonian, Phys. Rev. C 99, 044301 (2019).
  * Dutta  _et al._ [2021] R. Dutta, G. P. Chen, T. M. Henderson, and G. E. Scuseria, Construction of linearly independent non-orthogonal agp states (2021), arXiv:2101.08911 [physics.chem-ph] .
  * Henderson and Scuseria [2019] T. M. Henderson and G. E. Scuseria, Geminal-based configuration interaction, The Journal of Chemical Physics 151, 051101 (2019), https://doi.org/10.1063/1.5116715 .
  * Henderson and Scuseria [2020a] T. M. Henderson and G. E. Scuseria, Correlating the antisymmetrized geminal power wave function, The Journal of Chemical Physics 153, 084111 (2020a), https://doi.org/10.1063/5.0021144 .
  * Dukelsky  _et al._ [2003b] J. Dukelsky, G. Dussel, J. Hirsch, and P. Schuck, Comparison between exact and approximate treatments of the pairing interaction for finite fermi systems, Nuclear Physics A 714, 63 (2003b).
  * Henderson  _et al._ [2015] T. M. Henderson, I. W. Bulik, and G. E. Scuseria, Pair extended coupled cluster doubles, The Journal of Chemical Physics 142, 214116 (2015), https://doi.org/10.1063/1.4921986 .
  * Ripoche  _et al._ [2017] J. Ripoche, D. Lacroix, D. Gambacurta, J.-P. Ebran, and T. Duguet, Combining symmetry breaking and restoration with configuration interaction: A highly accurate many-body scheme applied to the pairing hamiltonian, Phys. Rev. C 95, 014326 (2017).
  * Henderson and Scuseria [2020b] T. M. Henderson and G. E. Scuseria, Correlating the antisymmetrized geminal power wave function, The Journal of Chemical Physics 153, 084111 (2020b), https://doi.org/10.1063/5.0021144 .
  * Khamoshi  _et al._ [2021] A. Khamoshi, G. P. Chen, T. M. Henderson, and G. E. Scuseria, Exploring non-linear correlators on agp, (2021), arXiv:2012.02377 [cond-mat.str-el] .
  * Qiu  _et al._ [2017] Y. Qiu, T. M. Henderson, J. Zhao, and G. E. Scuseria, Projected coupled cluster theory, The Journal of Chemical Physics 147, 064111 (2017), https://doi.org/10.1063/1.4991020 .
  * Richardson [1966] R. W. Richardson, Numerical study of the 8-32-particle eigenstates of the pairing hamiltonian, Phys. Rev. 141, 949 (1966).
  * Ui and Takeda [1983] H. Ui and G. Takeda, A Class of Simple Hamiltonians with Degenerate Ground State. II A Model of Nuclear Rotation: Spontaneous Breakdown of Rotation Symmetry and Goldstone Theorem for Finite Dimensional System, Progress of Theoretical Physics 70, 176 (1983), https://academic.oup.com/ptp/article-pdf/70/1/176/5367313/70-1-176.pdf .
  * Koma and Tasaki [1994] T. Koma and H. Tasaki, Symmetry breaking and finite-size effects in quantum many-body systems, J. Stat. Phys 76, 745 (1994).
  * Yannouleas and Landman [2007] C. Yannouleas and U. Landman, Symmetry breaking and quantum correlations in finite systems: studies of quantum dots and ultracold bose gases and related nuclear and chemical methods, Reports on Progress in Physics 70, 2067 (2007).
  * Papenbrock and Weidenmüller [2014] T. Papenbrock and H. A. Weidenmüller, Effective field theory for finite systems with spontaneously broken symmetry, Phys. Rev. C 89, 014334 (2014).
  * Löwdin [1955] P.-O. Löwdin, Quantum theory of many-particle systems. iii. extension of the hartree-fock scheme to include degenerate systems and correlation effects, Phys. Rev. 97, 1509 (1955).
  * Mayer [1980] I. Mayer, The spin-projected extended hartree-fock method,  Advances in Quantum Chemistry, 12, 189 (1980).
  * Bender  _et al._ [2003] M. Bender, P.-H. Heenen, and P.-G. Reinhard, Self-consistent mean-field models for nuclear structure, Rev. Mod. Phys. 75, 121 (2003).
  * Schmid [2004] K. Schmid, On the use of general symmetry-projected hartree-fock-bogoliubov configurations in variational approaches to the nuclear many-body problem, Progress in Particle and Nuclear Physics 52, 565 (2004).
  * Jiménez-Hoyos  _et al._ [2012] C. A. Jiménez-Hoyos, T. M. Henderson, T. Tsuchimochi, and G. E. Scuseria, Projected hartree-fock theory, The Journal of Chemical Physics 136, 164109 (2012), https://doi.org/10.1063/1.4705280 .
  * Ring and Schuck [1980] P. Ring and P. Schuck,  _The Nuclear Many-Body Problem_ (Springer-Verlag, Berlin Heidelberg, 1980).
  * Blaizot and Ripka [1985] J.-P. Blaizot and G. Ripka,  _Quantum Theory of Finite Systems_ (The MIT Press, Cambridge, MA, 1985).
  * Coleman [1965] A. J. Coleman, Structure of fermion density matrices. ii. antisymmetrized geminal powers, Journal of Mathematical Physics 6, 1425 (1965), https://doi.org/10.1063/1.1704794 .
  * Dukelsky and Sierra [1999b] J. Dukelsky and G. Sierra, Density matrix renormalization group study of ultrasmall superconducting grains, Phys. Rev. Lett. 83, 172 (1999b).
  * Sambataro and Sandulescu [2013] M. Sambataro and N. Sandulescu, Treatment of the like-particle pairing of quartets, J. Phys. G: Nucl. Part. Phys. 40, 055107 (2013).
  * Duguet [2015] T. Duguet, Symmetry broken and restored coupled-cluster theory: I. rotational symmetry and angular momentum, J. Phys. G: Nucl. Part. Phys. 42, 025107 (2015).
  * Duguet and Signoracci [2017] T. Duguet and A. Signoracci, Symmetry broken and restored coupled-cluster theory: Ii. global gauge symmetry and particle number, J. Phys. G: Nucl. Part. Phys. 44, 015103 (2017).
  * Tsuchimochi and Ten-no [2017] T. Tsuchimochi and S. Ten-no, Bridging single- and multireference domains for electron correlation: Spin-extended coupled electron pair approximation, Journal of Chemical Theory and Computation 13, 1667 (2017), pMID: 28240896, https://doi.org/10.1021/acs.jctc.7b00073 .
  * Baran and Delion [2020] V. V. Baran and D. S. Delion, A quartet bcs-like theory, Physics Letters B 805, 135462 (2020).
  * Zee [2010] A. Zee,  _Quantum field theory in a nutshell_ (Princeton University Press, 2010).
  * Peeters [2018] K. Peeters, Cadabra2: computer algebra for field theory revisited, Journal of Open Source Software 3, 1118 (2018).
  * [49] The nag library, the numerical algorithms group (nag), oxford, united kingdom www.nag.com.,  .


# Neural Network Exemplar Parallelization with Go

Georg Wiesinger  Faculty of Computer Science
University of Vienna
Vienna, Austria
<EMAIL_ADDRESS>Erich Schikuta  Faculty of Computer Science
University of Vienna
Vienna, Austria
<EMAIL_ADDRESS>
###### Abstract

This paper presents a case for exemplar parallelism of neural networks using
Go as parallelization framework. Further it is shown that also limited multi-
core hardware systems are feasible for these parallelization tasks, as
notebooks and single board computer systems. The main question was how much
speedup can be generated when using concurrent Go goroutines specifically. A
simple concurrent feedforward network for MNIST[1] digit recognition with the
programming language Go[2, 3, 4, 5] was created to find the answer. The first
findings when using a notebook (Lenovo Yoga 2) showed a speedup of 252% when
utilizing 4 goroutines. Testing a single board computer (Banana Pi M3)
delivered more convincing results: 320% with 4 goroutines, and 432% with 8
goroutines.

###### Index Terms:

Backpropagation, Exemplar Parallelization, Go Programming Language, MNIST

##  I Introduction

Neural networks and artificial intelligence are becoming more and more
important not only in research, but also in daily used technology. Due to
higher amounts of data these artificial agents have to analyze there is the
need for a larger throughput and highly efficient neural networks. The
programming language Go looks promising for developing such a highly efficient
agent, as the language itself has been made not only for highly efficient
parallelization but also with fast development in mind. The main question is
if Go is suitable for a highly efficient parallelization of neural networks.
The main objective is the creation of an efficient parallelized neural
network. There is a possibility that Go could lead to a higher parallelization
efficiency/speedup than other programming languages. As Go is a young
programming language, the literature about this specific topic is very sparse
to almost nonexistent. There are tertiary sources like websites comparing the
general throughput of Go in comparison to web languages like NodeJS, PHP and
Java111https://www.toptal.com/back-end/server-side-io-performance-node-php-
java-go. Other literature is related to parallelization speedup. There are
also some neural networks realized in Go. No sources for a better comparison
of parallelization techniques has been found. The scope of this work is to
find out the speedup when using multiple goroutines with a neural network
while maintaining a high and sustainable classification accuracy. A working
MNIST digit recognition system has been created for testing the speedup with
up to sixteen cores. The network and parameters have been optimized, but due
to only negligible improvements with more than 100 hidden layer nodes this
amount has not been exceeded. The execution time for one epoch has been sped
up from 856.57-1271.73 (median 1005.80) seconds with 1 goroutine to only
171.50-221.38 (median 201.82) seconds with 16 goroutines with a Banana Pi M3.
The Lenovo Yoga 2 showed a less significant speedup with 137.29-146.01 (median
142.33) for 1 goroutine to 55.10-64.62 (median 56.49) with 4 goroutines.
Additional goroutines exceeding the maximum thread limit brought further
speedup due to pipelining with the Banana Pi, but a negligible speed loss for
the Lenovo Yoga.

##  II Related Work and Baseline Research

Artificial neural networks and their parallel simulation gained high attention
in the scientific community. Parallelization is a classic approach for
speeding up execution times and exploiting the full potential of modern
processors. Still, not every algorithm can profit from parallelization, as the
concurrent execution might add a non-negligible overhead. This can also be the
case for data parallel neural networks, where accuracy problems usually occur,
as the results have to be merged.

In the literature a huge number of papers on parallelizing neural networks can
be found. An excellent source of references is the survey by Tal Ben-Nun and
Torsten Hoefler [7]. However, only few research was done on using Golang in
this endeavour.

In the following only specific references are listed, which influenced the
presented approach directly. The authors of [8] presented a parallel
backpropagation algorithm dealing with the accuracy problem only by using a
MapReduce and Cascading model. In the course of our work on parallel and
distributed systems [9, 10, 11] we developed several approaches for the
parallelization of neural networks. In [12], two novel parallel training
approaches were presented for face recognizing backpropagation neural
networks. The authors use the OpenMP environment for classic CPU
multithreading and CUDA for parallelization on GPU architectures. Aside from
that, they differentiated between topological data parallelism and structural
data parallelism [13], where the latter is focus of the presented approach
here. [14] gave a comparison of different parallelization approaches on a
cluster computer. The results differed depending on the network size, data set
sizes and number of processors. Using Go as parallelization tool was already
analyzed for the Single-Program-Multiple-Data approach [15, 16] and showed
promising results. However, in this paper we focus on exemplar
parallelization. Besides parallelizing the backpropagation algorithm for
training speed-up, alternative training algorithms like the Resilient
Backpropagation described in [17] might lead to faster convergence. One major
difference to standard backpropagation is that every weight and bias has a
different and variable learning rate. A detailed comparison of both network
training algorithms was given in [18] in the case of spam classification.

##  III Data and Methods

The following data and methods have been used to gain insight on the speedup
possibilities.

###  III-A Choosing the data and parallelization method

Different approaches of which data could be used have been evaluated. Weather
data, crime rates, etc. all seemed to be a good fit, but with the possibility
of very inconclusive outputs. Finally the "Hello World!" of neural networks
has been chosen: The MNIST dataset[1]. With this ready to use dataset the
development process sped up as the convolutional part was already done.

Exemplar parallelism[19] has been chosen as parallelization technique. Within
the workers the learning method is stochastic gradient descent[20], but due to
combining the data in the main connectome, it behaves like a mini-batch
update[21].

###  III-B Basic structure

First a functional code for a basic neural network has been prepared. With
this code it is also possible to define simple multi-layer feedforward
networks. From that stable basis more functionality has been added (i.E.
different activation functions) to ease up the future development. Then the
parallelization of the neural network has been implemented. There were
additional challenges with avoiding possible race conditions. Go was very
helpful with its built in race detector which can be utilized with the "-race"
flag. It was easy to spot any race conditions and therefore the development
sped up in the area deemed to take the most time. Afterwards the possibility
to input and compute large datasets has been implemented. A batch file
functionality for ease of testing as well as data output functionality have
been added too. Afterwards the code and neural network have been optimized for
a balance of speed, memory usage and training quality. Shuffling of the
training data has been implemented to prevent any unwanted behavior that comes
from repeated data. The elastic net regularization[22] has been chosen to get
better results and more stability for the neural network.

###  III-C The math

Different activation functions have been tested to get a high accuracy,
although this is not the purpose of this work. At the end it has been
concluded that the best way is to start without any data normalization for the
data to be put into the input layer. But before any activation function runs
over the layer, the data is normalized by dividing each value by the size of
the layer (including the bias) to minimize the risk of exploding
gradients[23].

The following activation functions have been used:

  * •

Input layer: Identity

  * •

Hidden layer: ELU[24]

  * •

Output layer: SoftMax

Variables are aas followed:

  * •

$\eta$ = learning rate

  * •

t = target

  * •

x = neuron value before activation

  * •

$\varphi$ = activation function

  * •

$\delta$ = error

  * •

w = weight

####  III-C1 Activation functions

Identity
The simplest activation function is the identity function. It simply states
that the value stays the same, as in equation (1).

$\varphi(x)=x$ (1)

So the derivation simply is 1 as shown in equation (2).

$\varphi^{\prime}(x)=1$ (2)

Exponential Linear Unit
In comparison to ReLU and leaky ReLU, ELU "[…] speeds up learning in deep
neural networks and leads to higher classification accuracies."[24] and
therefore has been chosen over the other options. Equation (3) shows the math,
where alpha is a positive value that can be freely chosen.

$\varphi(x)=\begin{cases}x&\text{if x $\geq$ 0}\\\
\alpha*(e\textsuperscript{x}-1)&\text{if x < 0}\end{cases}$ (3)

The derivation for training is shown in equation (4).

$\varphi^{\prime}(x)=\begin{cases}1&\text{if x $\geq$ 0}\\\
\varphi(x)+\alpha&\text{if x < 0}\end{cases}$ (4)

SoftMax
The SoftMax function gives us a classification of the likelihood that the
current input represents a certain number. The math is straightforward and
shown in equation (5).

$\varphi\textsubscript{i}(\overrightarrow{x})=\frac{e\textsuperscript{x\textsubscript{i}}}{\sum\limits_{j=1}^{J}e\textsuperscript{x\textsubscript{j}}}$
(5)

But with this equation, exploding or vanishing gradients[25] can become a
problem due to the high likelihood of numbers exceeding.For the SoftMax
activation there is a little "trick". It is possible to add a scalar, as shown
in (6), without changing the value of the softmax function[26].

$\varphi\textsubscript{i}(\overrightarrow{x})=\frac{e\textsuperscript{x\textsubscript{i}
+ S}}{\sum\limits_{j=1}^{J}e\textsuperscript{x\textsubscript{j} + S}}$ (6)

So, instead of using softmax(x), softmax(z) - with a scalar value of the
negative x maximum - has been used, as in equation (7).

$z\textsubscript{i}=(x\textsubscript{i}-max\textsubscript{i}(x\textsubscript{i}))$
(7)

If we use the maximum, we push the calculation into the negative number
spectrum. So, instead of having values ranging over ]-$\infty$ , $\infty$[,
they’ve been shifted to ]0, 1], as in (8).

$\varphi\textsubscript{i}(\overrightarrow{x})=\frac{e\textsuperscript{z\textsubscript{i}}}{\sum\limits_{j=1}^{J}e\textsuperscript{z\textsubscript{i}}}$
(8)

Equation (9) derivation for training is a little bit more complicated.

$\varphi^{\prime}\textsubscript{i}(\overrightarrow{x})=\frac{\partial\varphi\textsubscript{i}(\overrightarrow{x})}{\partial
x\textsubscript{j}}=\begin{cases}\varphi\textsubscript{i}(\overrightarrow{x})*(1-\varphi\textsubscript{j}(\overrightarrow{x}))&\text{i
= j}\\\
\varphi\textsubscript{i}(\overrightarrow{x})*(0-\varphi\textsubscript{j}(\overrightarrow{x}))&\text{i
$\neq$ j}\end{cases}$ (9)

Mathematicans use (10) to shorten the equation to (11).

$\delta\textsubscript{ij}=\begin{cases}1&\text{i = j}\\\ 0&\text{i $\neq$
j}\end{cases}$ (10)

$\varphi^{\prime}\textsubscript{i}(\overrightarrow{x})=\frac{\partial\varphi\textsubscript{i}(\overrightarrow{x})}{\partial
x\textsubscript{j}}=\varphi\textsubscript{i}(\overrightarrow{x})*(\delta\textsubscript{ij}-\varphi\textsubscript{j}(\overrightarrow{x}))$
(11)

But in the end it comes down to the same function as the logistic derivation,
shown in equation (12). As the result of the derivation is a diagonal
matrix[27] there is no need to calculate the whole matrix.

$\varphi^{\prime}\textsubscript{i}(\overrightarrow{x})=(1-\varphi\textsubscript{i}(\overrightarrow{x}))*\varphi\textsubscript{i}(\overrightarrow{x})=(1-x)*x$
(12)

####  III-C2 Elastic Net Regularization

The elastic net regularization[22] has been used for weight updates.

It’s a combination of the lasso regression (13) and the ridge regression (14).

$L\textsuperscript{1}=\lambda|w|$ (13)

$L\textsuperscript{2}=\lambda w\textsuperscript{2}$ (14)

The elastic net (15) is simple.

$ElasticNet=\lambda|w|+\lambda w\textsuperscript{2}$ (15)

Computational optimization (16)

$ElasticNet=\lambda(|w|+w\textsuperscript{2})$ (16)

For the derivation (17) the signum function (18) is needed.

$|w|=w*sgn(w)$ (17)

$sgn(w)=\begin{cases}1&\text{w > 0}\\\ 0&\text{w = 0}\\\ -1&\text{w <
0}\end{cases}$ (18)

Which leads to (19).

$ElasticNet^{\prime}=\lambda(sgn(w)+2w)$ (19)

####  III-C3 Loss function

Quadratic loss has been chosen as loss function. Although the classification
of handwritten digits - as the name says - is a classification problem and
therefore cross entropy loss should show better results. Different loss
functions will be implemented in the future of this work.

$L=-\frac{1}{2}\sum\limits_{i}^{nclass}\Bigl{(}t\textsubscript{i}-\varphi(x\textsubscript{i})\Bigr{)}\textsuperscript{2}+\lambda\frac{1}{2}\sum\limits_{i}^{k}\Bigl{(}|w\textsubscript{i}|+w\textsubscript{i}\textsuperscript{2}\Bigr{)}$
(20)

The derivative is the logistic equation, therefore the loss function is
equation (21).

$L^{\prime}=-\sum\limits_{i}^{nclass}\Bigl{(}t\textsubscript{i}-\varphi(x\textsubscript{i})\Bigr{)}+\lambda\sum\limits_{i}^{k}\Bigl{(}\frac{1}{2}sgn(w\textsubscript{i})+w\textsubscript{i})\Bigr{)}$
(21)

####  III-C4 Forward and backward propagation

All the previous information is needed to understand the forward and backward
propagation methods.

Forward
After setting the inputs and targets, the first layer of the neural network
gets activated. Then for each layer, the next neuron gets excited with the
product of the activated value and the weight of the connection between the
neurons (22).

$x\textsubscript{j}=\varphi(x\textsubscript{i})*w\textsubscript{ij}$ (22)

Backward
When learning, equation (23) is used to calculate the error.

$\delta=t-\varphi(x)$ (23)

The formula for the weight update, with learning rate and the regularization
in (24).

$\Delta w=\eta*(\delta*\varphi(x)+\lambda(sgn(w)+w))$ (24)

The final equation in (25).

$w\textsubscript{ij}\textsuperscript{+}=w\textsubscript{ij}-\Delta w,$ (25)

###  III-D Choosing the parameters

ELU alpha = 0.5 (currently hardcoded)
The alpha value for ELU has been hardcoded as there was no incentive to do
otherwise in the current software iteration.

workerBatch: 100
The worker batch has been chosen to merge the single neural networks as often
as possible, but without losing too much performance due to context switching.

Minimum/Maximum weight (starting weights): [-0.1; 0.1]
As the weights are usually getting smaller, when learning occurs, the starting
values have to be chosen to be 0.1 instead of 1. This led to the best outcome.

LearningRate: 0.8
The learning rate has been set to 0.8, as this led to the best outcome.

Lambda: 0.0000001
The multiplicator for the elastic net has been set to this value, as it
provided the highest accuracy for the training and test set. As it is hard to
tell if either L1 or L2 regularization is the best, there is only one lambda
for setting both methods, to achieve a balance between the two methods.

##  IV Results/Evaluation

For testing the neural network, two available systems have been chosen: The
Lenovo Yoga 2 laptop, as it is a dual core consumer product which utilizes
threading and a turbo mode for higher workloads, with 64 Bit Linux. The Banana
Pi M3, as it is a well known octa core home server, with no threading, and
without data distortion due to turbo mode kicking in, and 32 Bit Linux. Both
systems have a standard CPU frequency of 1.80 GHz, although the minimum and
maximum values differ.

There are stark differences in computation speed as well as speedup between
the Intel and the ARM architecture. As RISC and CISC lost their meaning to
describe newer architectures, it is not possible to draw the conclusion here ,
although the main effect could come from the smaller - and therefore faster
access rates - of the Intel L1 and L2 caches, or the lack of an L3 cache in
the ARM architecture. Further research would be needed.

###  IV-A Benchmark

When using pprof for checking the total cpu usage of the code parts with
BenchBatch (it utilizes 4 cores, uses a worker batch of 100 lines, and
processes 1 training with 60.000 MNIST lines as well as 1 test with 10.000
MNIST lines), it can be seen in Figure 1 that thinking and training takes up
about 96% of the total time. Thinking takes about 40% of the time, training
takes about 56%. Thinking is the forward propagation, training is the backward
propagation. Due to that high amount of cpu usage heavy optimizations were
made in these code parts, as these had the greatest effects.

Figure 1: pprof worker CPU profile

The utility functions only play a marginal role. Even though they wouldn’t
need any optimization, they’ve been optimized for general code quality
reasons. In example the garbage collector (mallocgc) is hardly used in the
utility functions, and almost never in the main code part. As strings are only
converted when needed, these parts of the code - even though they’re not
impacting the measurements - have been highly optimized. Maybe there’s still
room for further optimization, but for the general purpose this goal has been
exceeded.

###  IV-B Test systems

The final tests were made with a "Lenovo Yoga 2 Pro Multimode Ultrabook" as
well as a "Banana Pi M3".

Specifications of the Lenovo: Intel(R) Core(TM) i7-4500U CPU @ 1.80GHz, Dual
Core (4 threads)
Min CPU: 800 MHz, Max CPU: 3.0 GHz
32 KiB L1 cache, 256 KiB L2 cache, 4 MiB L3 cache
2x4096 DIMM @ Clockspeed 1600 MHz
64 Bit Linux Ubuntu 18.04.1 LTS

Specifications of the Banana Pi M3:
A83T ARM Cortex-A7 octa-core CPU @ 1.80 GHz, Octa Core (8 threads), 4800
BogoMIPS
ARMv7 Processor rev 5 (v7l)
Min CPU: 480 MHz, Max CPU: 1.8 GHz
512 KiB L1 cache, 1 MiB L2 cache
2GB LPDDR3
32 bit (armv7l) Linux Ubuntu 16.04.5 LTS, MATE Desktop Environment 1.12.1

####  IV-B1 Lenovo Yoga 2

The 252% speedup generated with 4 goroutines on the Lenovo Yoga 2 when
utilizing more than 1 processor is clearly visible in Figure 2. It is also
visible that using more goroutines than processors slows the execution time
down only by an almost negligible amount.

Figure 2: Lenovo Yoga 2 speedup

Parallelization speedup comes at a price. Although very small, there is a
slight decrease in recognition rates when utilizing more goroutines as shown
in Figure 3.

Figure 3: Lenovo Yoga 2 accuracy

####  IV-B2 Banana Pi M3

When looking at the results of the Banana Pi M3 in Figure 4, it is apparent
that utilizing multiple cores leads to an even greater benefit than with the
Lenovo. It was possible to generate a 320% speedup with 4 goroutines, and -
due to pipelining - it was even possible to generate over 498% speedup when
using more goroutines than there were threads available.

Figure 4: Banana Pi M3 speedup

The training and test set accuracies look promising too. A 99.26% training set
accuracy and 97.14% test set accuracy with only one core has been
accomplished. The accuracy does not get lower when utilizing more cores, even
though quality differences in the recognition rate can occur. In Figure 5 it
is clearly visible that recognition rate drops can occur at any time.

Figure 5: Banana Pi M3 accuracy

###  IV-C Accuracy growth depending on goroutines

When only one goroutine is usedFigure 6 with the Banana Pi, the neural network
starts with a very high recognition accuracy after the first epoch and has a
very good learning rate.

Figure 6: Accuracy growth with 1 Goroutine

With 16 goroutinesFigure 7 the recognition accuracy starts lower and the
network takes longer to learn.

Figure 7: Accuracy growth with 16 goroutines

Recognition rates in a nutshell

1 Goroutine, accuracy > 90%/95%/99%:
93.24% accuracy after 1 epoch, 1040 seconds
95.22% accuracy after 2 epochs, 2006 seconds
99.03% accuracy after 15 epochs, 15608 seconds

16 goroutines, accuracy > 90%/95%/99%:
90.18% accuracy after 2 epochs, 392 seconds
95.12% accuracy after 8 epochs, 1583 seconds
99.02% accuracy after 49 epochs, 9628 seconds

To reach a higher accuracy with more goroutines more epochs and training
samples are needed. But the speedup allows to train it in shorter time - or,
to look at it from another perspective - to compute more inputs in a much
shorter timespan.

##  V Lessons Learned

The final part of this work is to look at what has been learned about the
"do’s and don’ts of implementing neural networks", Go as a language, the drawn
conclusion, and possible future work.

###  V-A "To do, or not to do?" of implementing neural networks

There are certain roads to victory and many paths to development hell. The
latter leads to a steeper learning curve and should therefore be preferred
when trying to understand the implications of certain design decisions - but
under normal circumstances the beaten path is the quicker route. These
recommendations for other coders shall make implementing neural networks a
little bit easier and shine a light on which thought processes are good and
which are impractical to do.

####  V-A1 Arrays instead of structs

Do not use one struct instance per neuron, connection, etc. as it has a large
overhead. The compiler is able to optimize the usage of arrays. The first
iteration of the neural network took hours for just one epoch on the Lenovo,
while the array version takes less than a minute.

####  V-A2 Only save when necessary

Only save and load data when needed. In the context of the neural network:
Save either batchwise or after every epoch. Try to hold the data in the memory
as long as possible.

####  V-A3 Machine readable is better than human readable

The conversion of data to XML, JSON, or any other human readable format takes
a higher amount of computation time, memory, and disk space, than machine
readable formats. If a human readable format is needed, it should only be
created, if a human wants to read it and there is a need for them to do so.
Sifting through millions of weights and updates is not something a human
should do. But, depending on the use case, the human readable format can be
created, when

  * •

The process is finished and the results shall be shown.

  * •

An error occurs and the data is necessary to fix it.

If human entities want to access data while the process is running (in real
time, or stepwise for debugging) there are different approaches:

  * •

Create only one file every few epochs which can be accessed by multiple human
entities. Do NOT create it for every entity that accesses the file.

  * •

Duplicate the machine readable results and parse them on a different system.
For snapshots a simple ID can be given to every file.

####  V-A4 Parallelization and context switches

It takes time to store states of threads. Data has to be shoved around between
CPU caches. If applicable give a worker as much data as possible, with one
drawback in mind: More data merges mean higher fluctuations and slower
computation - less data merges can lead to a more stable convergence and
faster computation[28], as well as a higher level of generalization[23]. All
while being able to perform online learning due to the singular workers
performing stochastic gradient descent.

####  V-A5 Struct packing

In Go it is possible to pack structs. That means organizing the data types in
a way so that they waste the least amount of memory. The principle for this
work was "Memory is cheap, but even though students lack the money to buy
some, there is no need to overdo it". But one late evening (at 7 o’clock in
the morning) these principles had been thrown over board. So structs have been
packed.

####  V-A6 Slices

Do not loop through lines of data to append it to a batch line by line. Use
the slice functionality of Go - which passes them by reference - if
applicable. The following data has been taken from the old model with using
the updated weights for the error calculation Figure 8.

Figure 8: Code changes from wrong to correct code

In example the code in Figure 9 takes 129.68 seconds for 1 training and 1
testing with the MNIST dataset, 4 cores, and a worker batch setting of 100, as
shown in Figure 10.

Figure 9: Old code, appending to a new slice Figure 10: CPU profile before
code changes

In comparison when utilizing slices instead of making a slice and appending
the data lines within core.run() to send the worker batches to the workers as
shown in Figure 11 saved about 17 seconds on the Lenovo, as is visible in
Figure 12.

Figure 11: New code, using the slice functionality Figure 12: CPU profile
after code changes

But when comparing the old speedup in Figure 13 with the new speedup in Figure
14, the computation speed is on average a little slower. The assumption is
that - due to code optimizations and therefore less workload - the CPU chose
to use a lower clock speed, so computation took longer. This assumption is
untested, but context switching can be ruled out with a high certainty as the
possible speedup of the new code has been tested with 4 cores: If context
switching would have any negative impact, it should have shown in this test.

Figure 13: Benchmark before code changes Figure 14: Benchmark after code
changes

####  V-A7 Never trust a single source

"Why is this important, when the objective is just to write code?" one may
ask. The answer is pretty straightforward and was also a pitfall in the early
stages of this work.

When implementing software to solve a certain problem, there are often
different sources available. It is paramount to not trust a single source. If
possible, look at the original studies, try to get the original data - but
watch out for personal interests and possible skewed or even fabricated data.
Therefore also look for meta analyses and systematic reviews, which look at a
broader spectrum, the methods used, and other vital data points to find
eventual outliers or problems with the data. Also don’t trust sources that are
given to you by a single source. Check it first. Everything else costs time
and nerves: Only assuming, but not knowing, that the source is correct can
lead to triple-checking code, math, and data, even though the error lies
somewhere else. Two examples follow.

Many different sources in the web, including Wikipedia, managed to quote the
correct source for linear unit activation functions, but used the old
equations from the first two versions of the paper. There is a small effect on
the derivative functions when x = 0. In example with ELU: The old and new
equation are only equal when $\alpha$ = 1. Otherwise, when x = 0, f’(x) should
be $\alpha$, not 1.

One trusted source - due to one University module handed out an excerpt
without citing the source which was the basis for an assignment to calculate a
forward, backward, and forward propagation by hand, and creating a simple
neural network, which both got graded - used the updated weights as basis for
the backpropagation.

There was a larger accuracy dropFigure 15 when using multiple goroutines.

Figure 15: Old Plot showing the accuracy decrease

Although it was possible to see when to tweak the parameters to gain a higher
accuracy with a single coreFigure 16, there has been found no practical use of
these values for the correct implementation - the hyperparameters vary widely,
so they have to be tuned differently.

Figure 16: Accuracy with wrong parameters

Also mathematical sources are often not the best source for calculations in
information systems. Math has been spared the problem of both errors due to
overflow and underflow (except when using calculators or information systems).
Also there is no need to optimize for memory or computation speed.

The problem with trusting the wrong data has been solved with further research
from different sources and consulting a mathematician to check if the partial
derivatives and all formulas have been implemented correctly. The error has
then been found very quick when checking against the standard reference[26].

###  V-B Comment on Go

Go is a wonderful language to write code. Implementation and testing of the
neural network seemed to be easier than with other programming languages. But
go also has some drawbacks (as does any language).

The main annoyance were the "unused imports" bugs. Sometimes only certain
outputs are needed for testing which will get dropped by the developer
immediately afterwards. It’s good that the Go compiler sees these oversights
as errors, even though they are a huge annoyance. A probably better way would
be if unused imports won’t be tested in a debug environment, only in
production. But this would have additional drawbacks when Go is used in
environments where code quality is not highly valued.

Another annoyance is the "sudden ’bad file descriptor’ of doom". Sometimes
it’s just a "data reader error: file already closed". It was not possible to
pin down what exactly causes the error, only that it affects the file as a
whole. Not even deleting and creating a new file with the same name helps to
overcome that error. Further testing is needed.

An additional observation that can ruin ones day is, that the Go compiler for
some reason accesses trashed files, at least under Linux. There is no problem
when files are overwritten by new files. But if a file gets deleted, and a new
one inserted instead, Go sometimes seems to try to compile the deleted files,
which can lead to hard to trace errors. If there’s, in example, an error where
the Go compiler expects an integer value, the code provides an integer value,
but recently a file with the same function expecting a double value had been
trashed, simply empty the trash bin.

Another "hard to debug except when you know it" part is: "panic: runtime
error: invalid memory address or nil pointer dereference". This error occurs
when the object has not been created with new(…). If it’s further up in the
code, i.E. some struct attribute, this error is not easy to find. When
starting with Go that panic tells almost nothing about its nature.

Circular dependencies are not allowed. They can happen while refactoring code
or when making some design mistake. It’s good that Go does not permit them as
they are a sign of bad software design.

The short variable declaration := is very handy. Go recognizes the type and
assigns the value to the left hand variable. The best part: It won’t break
type safety, which prevents weird behavior.

With the test coverage profiler it is easy to see the current code coverage.
There is also the possibility to create test heat maps and to show the test
coverage in the browser with highlighting good, poor, or not covered code
parts222Go test coverage and html heatmap: https://blog.golang.org/cover.

There are memory and cpu profilers, and even a profiling tool333Go profiling:
https://blog.golang.org/profiling-go-programs. It is easy to list the top cpu
or memory consumers or show a profiling web. Therefore memory issues can be
found easily, as well as slow code parts.

Go uses function inlining which is a great method for speeding up code.

Goroutines are very lightweight444Currently 2, 4, or 8 kB per Goroutine,
depending on the version, i.E. https://github.com/golang/go/issues/7514. As
they’re very efficient and only start to run if they get data from a channel,
there’s the probability of an application for parallelized neurons instead of
only parallelized networks.

It’s easy to find and fix race conditions with Go as it comes with its own
race detector.

##  VI Conclusion and Future work

It has been learned how to use the programming language Go and about its
parallel speedup possibilities. The main accomplishment of this work is to
have managed to create a stable and fast neural network. The hardest part was
to understand the mathematical concepts and ramifications behind neural
networks and how to implement them software wise.

The main focus of this work was to see how the parallel speedup of a neural
network behaves with the language Go. Due to time and resource restrictions
only little derivations from the main focus were made. There are still ways
left to make this neural network even more efficient, with higher accuracy,
and so on. The current version could have some possible memory leaks. They
will be fixed in a future version. As there will be further changes due to
development and additional insights, the code will probably be refined and
refactored in the future.

Some parts of the code are still untested - mainly file reading and writing.
As they work as intended no additional effort has been made to get 100% test
coverage in these areas. Here is room for improvement.

Optimization of the neural network would be the largest part of the future
work. Currently it is only a simple network with Bias. It would be possible to
implement momentum[29] and other artifacts to achieve higher accuracies. NADAM
and other stochastic gradient descent optimization algorithms[28] could be
implemented too.

Smaller changes will also include several options, in example if the user
wants bias nodes, which error severity to log, and to choose different lambdas
for the L1 and L2 Regularization in the elastic net. Adaptive learning
rates[29] would be of interest too. Different loss functions, especially Cross
Entropy Loss[30] will be implemented in the future.

There is an interest to look into Self-Normalizing Neural Networks[31].

## References

  * [1] C. J. B. Yann LeCun, Corinna Cortes, “THE MNIST DATABASE of handwritten digits,” http://yann.lecun.com/exdb/mnist/index.html, accessed: 2018-08-12.
  * [2] Google LLC, “Go website,” https://golang.org/, accessed: 2018-08-12.
  * [3] ——, “Go blog,” https://blog.golang.org/, accessed: 2018-08-12.
  * [4] ——, “Go tour,” https://tour.golang.org/, accessed: 2018-08-12.
  * [5] Multiple contributors, visible on https://github.com/golang/go/wiki/Home/_history, “Github go wiki,” https://github.com/golang/go/wiki/, accessed: 2018-08-12.
  * [6] Google LLC, “Go version 1 release,” https://blog.golang.org/go-version-1-is-released, accessed: 2018-08-12.
  * [7] T. Ben-Nun and T. Hoefler, “Demystifying parallel and distributed deep learning: An in-depth concurrency analysis,” _ACM Computing Surveys (CSUR)_ , vol. 52, no. 4, pp. 1–43, 2019.
  * [8] Y. Liu, W. Jing, and L. Xu, “Parallelizing backpropagation neural network using mapreduce and cascading model,” _Computational intelligence and neuroscience_ , vol. 2016, 2016.
  * [9] E. Schikuta, T. Fuerle, and H. Wanek, “Vipios: The vienna parallel input/output system,” in _Euro-Par’98 Parallel Processing_ , D. Pritchard and J. Reeve, Eds. Berlin, Heidelberg: Springer Berlin Heidelberg, 1998, pp. 953–958.
  * [10] P. Brezany, T. A. Mueck, and E. Schikuta, “A software architecture for massively parallel input-output,” in _Applied Parallel Computing Industrial Computation and Optimization_ , J. Waśniewski, J. Dongarra, K. Madsen, and D. Olesen, Eds. Berlin, Heidelberg: Springer Berlin Heidelberg, 1996, pp. 85–96.
  * [11] E. Schikuta and T. Weishaupl, “N2grid: neural networks in the grid,” in _2004 IEEE International Joint Conference on Neural Networks (IEEE Cat. No.04CH37541)_ , vol. 2, 2004, pp. 1409–1414 vol.2.
  * [12] A. A. Huqqani, E. Schikuta, S. Ye, and P. Chen, “Multicore and gpu parallelization of neural networks for face recognition,” _Procedia Computer Science_ , vol. 18, no. Supplement C, pp. 349 – 358, 2013, 2013 International Conference on Computational Science.
  * [13] E. Schikuta, “Structural data parallel neural network simulation,” in _Proceedings of 11th Annual International Symposium on High Performance Computing Systems (HPCS’97), Winnipeg, Canada_ , 1997.
  * [14] M. Pethick, M. Liddle, P. Werstein, and Z. Huang, “Parallelization of a backpropagation neural network on a cluster computer,” in _International conference on parallel and distributed computing and systems (PDCS 2003)_ , 2013.
  * [15] D. Kalwarowskyj and E. Schikuta, “Spmd-based neural network simulation with golang,” in _International Conference on Computational Science_. Springer, 2023, pp. 563–570.
  * [16] D. Turner and E. Schikuta, “Go for parallel neural networks,” in _Advances in Computational Intelligence: 15th International Work-Conference on Artificial Neural Networks, IWANN 2019, Gran Canaria, Spain, June 12-14, 2019, Proceedings, Part II 15_. Springer, 2019, pp. 833–844.
  * [17] M. Riedmiller and H. Braun, “A direct adaptive method for faster backpropagation learning: The rprop algorithm,” in _Neural Networks, 1993., IEEE International Conference on_. IEEE, 1993, pp. 586–591.
  * [18] N. Prasad, R. Singh, and S. P. Lal, “Comparison of back propagation and resilient propagation algorithm for spam classification,” in _2013 Fifth International Conference on Computational Intelligence, Modelling and Simulation_ , Sept 2013, pp. 29–34.
  * [19] R. Rogers and D. Skillicorn, “Strategies for parallelizing supervised and unsupervised learning in artificial neural networks using the bsp cost model,” _Queens University, Kingston, Ontario, Tech. Rep_ , 1997.
  * [20] L. Bottou, “Large-scale machine learning with stochastic gradient descent,” in _Proceedings of COMPSTAT’2010_. Springer, 2010, pp. 177–186.
  * [21] M. Zinkevich, M. Weimer, L. Li, and A. J. Smola, “Parallelized stochastic gradient descent,” in _Advances in Neural Information Processing Systems 23_ , J. D. Lafferty, C. K. I. Williams, J. Shawe-Taylor, R. S. Zemel, and A. Culotta, Eds. Curran Associates, Inc., 2010, pp. 2595–2603. [Online]. Available: http://papers.nips.cc/paper/4006-parallelized-stochastic-gradient-descent.pdf
  * [22] H. Zou and T. Hastie, “Regularization and variable selection via the elastic net,” _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , vol. 67, no. 2, pp. 301–320, 2005.
  * [23] S. Ioffe and C. Szegedy, “Batch normalization: Accelerating deep network training by reducing internal covariate shift,” _CoRR_ , vol. abs/1502.03167, 2015. [Online]. Available: http://arxiv.org/abs/1502.03167
  * [24] D. Clevert, T. Unterthiner, and S. Hochreiter, “Fast and accurate deep network learning by exponential linear units (elus),” _CoRR_ , vol. abs/1511.07289, 2015. [Online]. Available: http://arxiv.org/abs/1511.07289
  * [25] R. Pascanu, T. Mikolov, and Y. Bengio, “Understanding the exploding gradient problem,” _CoRR_ , vol. abs/1211.5063, 2012. [Online]. Available: http://arxiv.org/abs/1211.5063
  * [26] I. Goodfellow, Y. Bengio, and A. Courville, _Deep Learning_. MIT Press, 2016, http://www.deeplearningbook.org.
  * [27] I. M. Elfadel and J. L. Wyatt, Jr., “The "softmax" nonlinearity: Derivation using statistical mechanics and useful properties as a multiterminal analog circuit element,” in _Advances in Neural Information Processing Systems 6_ , J. D. Cowan, G. Tesauro, and J. Alspector, Eds. Morgan-Kaufmann, 1994, pp. 882–887.
  * [28] S. Ruder, “An overview of gradient descent optimization algorithms,” _CoRR_ , vol. abs/1609.04747, 2016. [Online]. Available: http://arxiv.org/abs/1609.04747

  * [29] Y. A. LeCun, L. Bottou, G. B. Orr, and K.-R. Müller, _Efficient BackProp_. Berlin, Heidelberg: Springer Berlin Heidelberg, 2012, pp. 9–48. [Online]. Available: https://doi.org/10.1007/978-3-642-35289-8_3
  * [30] P. Sadowski, “Notes on backpropagation,” _homepage: https://www. ics. uci. edu/~ pjsadows/notes. pdf (online)_ , 2016.
  * [31] G. Klambauer, T. Unterthiner, A. Mayr, and S. Hochreiter, “Self-normalizing neural networks,” _CoRR_ , vol. abs/1706.02515, 2017. [Online]. Available: http://arxiv.org/abs/1706.02515