#  $\pi$-Phase shift across stripes in a charge density wave system

T. Ying<EMAIL_ADDRESS>School of Physics, Harbin Institute of
Technology, Harbin 150001, China  R.T. Scalettar
<EMAIL_ADDRESS>Department of Physics, University of California,
Davis, CA 95616, USA  R. Mondaini<EMAIL_ADDRESS>Beijing Computational
Science Research Center, Beijing 100193, China

###### Abstract

Many strongly correlated materials are characterized by deeply intertwined
charge and spin order. Besides their high superconducting transition
temperatures, one of the central features of these complex patterns in
cuprates is a phase shift which occurs across lines of decreased hole density.
That is, when doped away from their AF phase, the additional charge is not
distributed uniformly, but rather in ‘stripes’. The sublattices preferentially
occupied by up and down spin are reversed across these stripes, a phenomonenon
referred to as a ‘$\pi$-phase shift’. Many of the spin-charge patterns,
including the $\pi$-phase shift, are reproduced by Density Matrix
Renormalization Group and Quantum Monte Carlo calculations of simplified tight
binding (repulsive Hubbard) models. In this paper we demonstrate that this
sublattice reversal is generic by considering the corresponding phenomenon in
the attractive Hubbard Hamiltonian, where a charge density wave phase forms at
half-filling. We introduce charge stripes via an appropriate local chemical
potential; measurements of charge correlation across the resulting lines of
lowered density reveal a clear $\pi$ phase.

##  I Introduction

Inhomogeneous phases are a central feature of strongly correlated materials,
and especially of oxide systems [1]. The manganites are one example [2, 3, 4,
5]. They have ferromagnetic and antiferromagnetic (AF) states in close
proximity in energy, and when a small quenched disorder is included, extended
glassy regions emerge in which these phases coexist. This regime is very
sensitive to external perturbation, e.g. the application of a small magnetic
field, leading to the phenomenon of colossal magnetoresistance. The cuprates
are another instance. Here a wide body of experiments, including transport
[6], nuclear magnetic resonance [7, 8], X-ray scattering [9, 10, 11], and
scanning tunneling microscopy [12], has indicated that complex patterns of
charge and spin develop upon doping [13], and that these inhomogeneous
structures are also present in the pairing gaps [14, 11], thereby suggesting a
possible connection to their high superconducting transition temperatures.

Quite remarkably, many of the intricate details of these structures are
reproduced in calculations on simple Hamiltonians, even when the models are
translation invariant [15]. Indeed, one of the earliest indications of stripe
physics, in which doped holes arrange themselves in linear patterns, came out
of mean-field calculations [16, 17, 18], pre-dating much of the experimental
work. These observations have been confirmed by a wide variety of methods like
functional renormalization group [19], the Density Matrix Renormalization
Group (DMRG) [20, 21, 22, 23], exact diagonalization (ED) [24], dynamic mean
field theory [25], Auxiliary Field Quantum Monte Carlo (AFQMC) [26, 27, 28,
29], infinite projected entangled pair states [30], and density matrix
embedding theory [31], all of which treat many-body effects more exactly.

One central feature of these striped phases in the repulsive Hubbard model is
their mixed charge and spin character, and, in particular a ‘$\pi$-phase
shift’ in the spin correlations which is found to exist across a linear
depletion of charge: The sublattice which holds the surplus of up spin
character on one side of a stripe instead holds a surplus of down spin as the
stripe is traversed [32].

In this paper we examine whether such $\pi$-phase shifts exist in a charge
density wave (CDW) phase across a (charge) stripe. That is, we address the
question of whether the mixed character of the stripe, and the type of order
being established across it, is essential to the existence of a phase shift.
We address this question by using Determinant Quantum Monte Carlo (DQMC)
simulations of the attractive Hubbard model in which stripes are imposed
externally via raising the local site energy along several rows of the
lattice. While this method does not establish the spontaneous formation of
inhomogeneities, it does allow an exploration of sublattice order reversal.
Significantly, because of the absence of a sign problem [33, 34], we are able
to examine low temperatures and especially the effect of stripes and the
$\pi$-phase shift on ($s$-wave) pairing correlations. Addressing such issues
is much more challenging in the repulsive Hubbard model [27] because of the
sign problem.

Our key results are that the introduction of a charge stripe of sufficient
depth does cause the CDW domains on opposite sides of the stripe to develop a
$\pi$-phase shift, so that their high and low density sublattices are
interchanged. In this way the structure of the charge order across a density
stripe mimics the well-known behavior of the spin order in the repulsive
Hubbard model. We also show that this $\pi$-phase shift appears to be
detrimental to pairing order.

This paper is organized as follows. In Sec. II we define the Hamiltonian and
correlation functions we will investigate, and give a brief description of the
DQMC and ED methodologies. Section III presents our results for charge and
pairing correlations, focusing especially on the issue of a $\pi$-phase shift.
Finally, Sec. IV summarizes our conclusions.

##  II Hamiltonian and methodology

We study a two dimensional square lattice attractive Hubbard Hamiltonian in
which stripes are introduced externally via a raised site energy $V_{0}$ on a
set of rows of period ${\cal P}$, ${\bf i}=(i_{x},i_{y})$ with ${\rm
mod}(i_{y},{\cal P})=0$.

$\displaystyle\hat{\mathcal{H}}=-t\sum_{\langle{\bf
i\,j}\rangle\,\sigma}(\hat{c}_{{\bf i}\sigma}^{\dagger}\hat{c}_{{\bf
j}\sigma}^{\phantom{\dagger}}+\hat{c}_{{\bf j}\sigma}^{\dagger}\hat{c}_{{\bf
i}\sigma}^{\phantom{\dagger}})+U\sum_{\bf i}\hat{n}_{{\bf
i}\uparrow}\hat{n}_{{\bf i}\downarrow}$ $\displaystyle-\mu\sum_{{\bf
i}}(\hat{n}_{{\bf i}\uparrow}+\hat{n}_{{\bf
i}\downarrow})+V_{0}\sum_{i_{y}\in{\cal P}}(\hat{n}_{{\bf
i}\uparrow}+\hat{n}_{{\bf i}\downarrow}).$ (1)

We choose ${\cal P}=4$, so that each stripe (row with $V_{0}$ term active,
blue spheres in Fig. 1) is separated by three rows where the $V_{0}$ term is
not present (light-colored spheres in Fig. 1).

Figure 1: (Color online) Illustration of the Hamiltonian in an $8\times 8$
lattice with period ${\cal P}=4$; sites with $V_{0}$ active are depicted in
blue, whereas the interstripe sites are white. The other relevant parameters,
hopping and attractive interactions are annotated.

Our primary methodology is the DQMC approach [35, 36]. We first express the
partition function associated with the Hamiltonian of Eq. 1, ${\cal Z}={\rm
Tr}\,e^{-\beta\hat{\cal H}}$, as a path integral by discretizing the imaginary
time $\beta=L_{\tau}\Delta\tau$. This allows us to use the Trotter
approximation $e^{-\Delta\tau\hat{\cal H}}\approx e^{-\Delta\tau\hat{\cal
K}}e^{-\Delta\tau\hat{\cal V}}$ where $\hat{\cal K}$ and $\hat{\cal V}$ are
the single- and two-particle terms of $\hat{\cal H}$ respectively. In this
work we use $\Delta\tau$ = 0.1; an analysis of Trotter errors in observables
is given in Appendix B. With $\hat{\cal V}$ isolated, we can introduce a
(discrete) Hubbard-Stratonovich field (HSF) [37] to decouple the interactions.
The trace which gives ${\cal Z}$ now contains only exponentials of quadratic
forms of fermion operators and can be done analytically. The resulting
expression for ${\cal Z}$ is a sum over the configurations of the HSF, which
is done stochastically. We will present results mainly for $8\times 8$
($12\times 12$ and $16\times 16$-lattice results are shown in Appendix C),
i.e. linear extent $L=8$, and supplement the DQMC calculations with ED to
confirm and further understand the physics. Our ED calculations, presented in
the Appendix E, focus on a $2\,\times\,4$ lattice, which has a manageable
Hilbert space size yet is large enough to hold a stripe and examine the spin
and charge patterns at finite temperature $T$.

For $V_{0}\neq 0$, the density is inhomogeneous. We denote by $\rho$ the
overall density, averaged over the entire lattice. $\rho_{\rm str}$ is the
density on the stripes, the sites with $V_{0}\neq 0$, and $\rho_{\rm dom}$ is
the density on the domains between the stripes. For ${\cal P}=4$ these are
related by $\rho=\frac{1}{4}\rho_{\rm str}+\frac{3}{4}\rho_{\rm dom}$.

We focus attention on the equal-time, real space density-density correlation
function

$\displaystyle c({\bf r})=\big{\langle}\,(\hat{n}_{\bf i}-1)(\hat{n}_{\bf
i+r}-1)\,\big{\rangle}$ (2)

between two sites ${\bf i}$ and ${\bf j}$. The subtraction is a convention,
useful for this paper where much of the lattice has $\langle\hat{n}_{\bf
i}\rangle=1$, since it then enables a focus on fluctuations away from the
average density. In a homogeneous CDW phase, this correlation function is
long-ranged, taking positive values for pairs of sites on the same sublattice
and negative values for sites on opposite sublattices. We also examine the
real space, equal-time $s$-wave pairing correlators $p({\bf r})$, and their
associated structure factor $P_{s}$ at ${\bf q}=0$:

$\displaystyle p({\bf r})$
$\displaystyle=\big{\langle}\,\hat{\Delta}^{\phantom{\dagger}}_{\bf
i+r}\hat{\Delta}^{\dagger}_{\bf i}\,\big{\rangle}\hskip
25.29494pt\hat{\Delta}^{\dagger}_{\bf i}\equiv\hat{c}^{\dagger}_{{\bf
i}\uparrow}\hat{c}^{\dagger}_{{\bf i}\downarrow}$ $\displaystyle P_{s}$
$\displaystyle=\sum_{\bf r}\,p({\bf r})$ (3)

which signal the formation of a superconducting phase. In the presence of off-
diagonal long range order $p({\bf r})$ approaches a non-zero value as ${\bf
r}\rightarrow\infty$, and $P_{s}$ grows linearly with the lattice size
$L^{2}$.

In the absence of $V_{0}$, and at half-filling $\mu=U/2$, the attractive
Hubbard model has simultaneous CDW and $s$-wave order. This enlarged symmetry
of possible ordered phases implies that long range order (LRO) is possible
only at $T=0$ [38, 39]. On finite lattices, LRO can be observed as long as the
temperature is lowered to a value for which the correlation length $\xi$
exceeds the linear lattice size $L$. The introduction of stripes ($V_{0}\neq
0$) breaks this degeneracy between CDW and pairing. As we will discuss below,
the resulting anisotropy favors pairing order, if small in strength, an
observation which will be confirmed by our DQMC simulations at small values of
$V_{0}$. The qualitative argument is similar to that for uniform doping of the
attractive Hubbard model [38, 39].

##  III Results

We first focus separately on the density and pairing correlations, and then
comment on how they are related, interpreting the enhancement of pairing order
for weak stripes in terms of a particle-hole transformation to the repulsive
Hubbard model.

###  III.1 Density-density correlations: Demonstration of Stripe Formation

We begin with DQMC results for the density-density correlations. The central
feature of the $\pi$-phase shift is the interchange of the sublattice ordering
pattern across the stripe. A spin (or charge) correlation connecting sites on
the same sublattice, which would be positive in the absence of a stripe,
becomes negative. Likewise, a spin (or charge) correlation connecting sites on
the different sublattices, which would be negative in the absence of a stripe,
becomes positive. The pairs of sites connected by blue and red arrows in the
inset of Fig. 2(a) show these two cases.

Figure 2: (a) Density-density correlator $c\big{(}{\bf r}=(0,2)\big{)}$ (blue
circles and blue arrow in the inset) and $c\big{(}{\bf r}=(1,2)\big{)}$ (red
squares, and red arrow in the inset), as functions of stripe strength $V_{0}$.
The emergence of $\pi$-phase shift at $V_{0}\sim t$ is signalled by the
reversal in sign of the correlation functions. (b) The total electron density
$\rho$ and the density on the stripe $\rho_{\rm str}$, as functions of stripe
strength $V_{0}$. The electron density off stripe $\rho_{\rm dom}$ is fixed to
be $1$ by adjusting the global chemical potential $\mu$. Here the lattice size
is $8\,\times\,8$, the on-site attraction $U=-4\,t$ is half the non-
interacting bandwidth, and the inverse temperature $\beta\,t=10$.

We perform a set of simulations in which $V_{0}$ is increased, gradually
forming well-defined stripes of reduced fermion occupation. We simultaneously
adjust the global chemical potential $\mu$ to keep the density on the sites of
the domains between the stripes $\rho_{\rm dom}=1$, a value which is optimal
for CDW order (other densities are investigated in Appendix D). Figure 2(a) is
one of the central results of our paper. It demonstrates that as stripes are
introduced, a $\pi$-phase shift emerges. For the parameters shown, $U=-4\,t$
and $T=t/10$, the crossover occurs at $V_{0}\sim t$.

The corresponding density evolution is shown in Fig. 2(b). The blue curve
shows the rapid depletion of fermion occupation on the stripes. At $V_{0}\sim
t$ where the $\pi$ phase shift emerges, $\rho_{\rm str}\sim 0.25$. The black
curve shows the overall fermion occupation $\rho$, which approaches 0.75 in
the large $V_{0}$ limit since 1/4 of the sites (the stripes) have been driven
to empty with the remaining 3/4 of the sites remaining at unit density. The
two curves are not independent, being related, as noted earlier, by
$\rho=\frac{3}{4}+\frac{1}{4}\rho_{\rm str}$; both are shown for clarity,
however.

If one instead investigates the repulsive Hubbard model ($U>0$), see Appendix
A, a qualitatively similar picture follows: Sufficiently large stripe energies
lead to the emergence of a magnetization reversal, i.e., a magnetic
$\pi$-phase shift appears, whose onset, however, occurs at much larger
$V_{0}$’s ($V_{0}\simeq 4t$ for $U=+4t$).

Turning back to the attractive case, to characterize the nature of the charge
pattern more precisely, Fig. 3(a) shows the density-density correlation
function $c\big{(}\,{\bf r}=(x,0)\,\big{)}$ between sites running immediately
parallel to the charge stripe. These correlations exhibit an interesting non-
monotonicity with $V_{0}$. As might be expected, the oscillating CDW pattern
is most robust at $V_{0}=0$ when one has a pristine half-filled lattice. For
weak stripe potentials $V_{0}=0.5,1.0$, there is only short range order, and
$c\big{(}\,{\bf r}=(4,0)\,\big{)}\sim 0$. However, as $V_{0}$ is further
increased, CDW order is recovered. For $V_{0}=3$, $c\big{(}\,{\bf
r}=(4,0)\,\big{)}$ is nearly as large as its $V_{0}=0$ value. Fig. 3(b)
emphasizes this behavior by showing $c\big{(}\,{\bf r}=(4,0)\,\big{)}$ as a
function of $V_{0}$. As a comparison, we also show $c\big{(}\,{\bf
r}=(4,0)\,\big{)}$ along the middle line in the domain as a function of
$V_{0}$.

Figure 3:  (a) Density-density correlations $c\big{(}{\bf r}=(x,0)\big{)}$ on
a row parallel to, and immediately neighboring, the stripe, as indicated in
the lattice of the inset. Results for different values of the stripe strength
$V_{0}$ are shown. (b) Density-density correlations $c\big{(}\,{\bf
r}=(4,0)\,\big{)}$ on a row parallel to the stripe, for both immediately
neighboring and in the middle of two stripes, as indicated in the lattice of
the inset. Results for different values of the stripe strength $V_{0}$ are
shown. In both panels, the density in the inter-stripe domains is fixed at
$\rho_{\rm dom}=1$ by adjusting $\mu$. Here the lattice size is
$8\,\times\,8$, the on-site attraction $U=-4\,t$, and the inverse temperature
$\beta\,t=10$.

The preceding results are for a coupling $U=-4\,t$ which is one-half the non-
interacting bandwidth. The low temperature spin and charge correlations of the
Hubbard Hamiltonian are usually strongest at somewhat larger interaction
strengths. In Fig. 4(a) we examine $c({\bf r})$ across a stripe for $U=-6\,t$.
The $\pi$-phase shift forms at slightly smaller values of $V_{0}$. A
compilation of the critical energies $V_{0}$ that trigger a $\pi$-phase shift
as a function of $U$ is given in Fig. 4(b), showing both the critical
$V_{0c}$, and the corresponding fermion occupation on the stripe $\rho_{\rm
str}$. The trend above described is confirmed, that a smaller stripe energy
$V_{0}$ with systematically larger interactions (in magnitude) is sufficient
for the onset of a $\pi$-phase shift. Even more remarkably, the electronic
density at the stripes where the transition occurs is largely $U$-independent
in the range investigated, pinned at filling $\rho_{\rm str}=0.25$. This
corresponds to the total density $\rho=0.8125$ [$\rho_{\rm dom}=1$].

In the Appendix C we show that these results for the cross-stripe spin
correlations are converged in both inverse temperature $\beta$ and spatial
lattice size $L$.

Figure 4: (a) Similar to Fig. 2 (a), but with the interaction strength
increased to $U=-6\,t$. The formation of $\pi$-phase shift occurs at a
slightly smaller stripe potential. (b) Compilation of the critical values
$V_{0c}$ where the $\pi$-phase shift appears with different interaction
strengths, and the corresponding electron density on the stripe $\rho_{\rm
str}$. The inverse temperature is set at $\beta t=10$.

###  III.2 $s$-wave pairing correlations

We next look at how the $s$-wave pairing is affected by the stripe strength.
We begin with the pair structure factor $P_{s}$, focusing, as before, on
fermionic ‘domain’ occupation $\rho_{\rm dom}=1$. Figure 5 shows the
dependence of $P_{s}$ on $V_{0}$. Initially (small $V_{0}$), the presence of
an extrinsic density modulation results in enhanced $s$-wave pairing. This
increase exists roughly up to the value of $V_{0}=V_{0c}(U)$ where the charge
$\pi$-phase shift is formed. Subsequently, a further increase in the stripe
strength is detrimental to the pairing. This behavior is consistent with
trends at different values of the interaction $U$ [See Figs. 5(a) and 5(b)].

Figure 5: The $s$-wave pairing structure factor $P_{s}$ (left y-axes) as a
function of $V_{0}$ and the corresponding pair susceptibility $\chi_{s}$
(right y-axes) for $\rho_{\rm dom}=1$, in an $8\times 8$ lattice with $U/t=-4$
(a) and $-8$ (b). Apart from the very different scales, the two quantities
qualitatively follow each other. Vertical dashed lines give the critical
stripe strength for $\pi$-phase shift formation $V_{0c}(U)$. The inverse
temperature $\beta t=10$.

The pair-field susceptibility,

$\displaystyle\chi_{s}=\sum_{\bf r}\int d\tau\,\,\langle e^{\tau\hat{\cal
H}}\,\hat{\Delta}_{\bf i+r}\,e^{-\tau\hat{\cal
H}}\,\hat{\Delta}^{\dagger}_{\bf i}\rangle,$ (4)

provides a more sensitive measurement of superconductivity (SC), since it also
samples over the correlations in imaginary time. $\chi_{s}$ as a function of
$V_{0}$ for the density $\rho_{\rm dom}=1$ is shown (red line) in Fig. 5. The
pair structure factor $P_{s}$ (black line of Fig. 5) and $\chi_{s}$ maintain
the same ratio independent of $V_{0}$: $\chi_{s}/P_{s}\sim\beta$. Indeed,
since $\chi_{s}$ includes $\int_{0}^{\beta}d\tau$ (see Eq. 4), a factor of
$\beta$ is the ‘expected’ ratio, assuming there is very little decay in
imaginary time. The data of Fig. 5 therefore are indicative of long range pair
correlations in $\tau$.

###  III.3 Intertwined Order

The preceding results suggest that in this model pairing is initially enhanced
when charge stripes are introduced, but that the subsequent development of a
$\pi$-phase shift between the CDW domains then rapidly suppresses the SC. This
is a non-trivial result. One might intuitively expect that turning on $V_{0}$
would promote SC by doping the stripes away from half-filling, but the
subsequent diminution by the emergence of CDW sublattice reversal does not
have an obvious origin.

We can get a more nuanced view of the interplay of the charge patterns and
pairing by recalling that, on a bipartite lattice, the Hamiltonian of Eq. 1
with $U<0$ has a well known mapping which reverses the sign to $U>0$. This is
accomplished via a partial particle-hole transformation (PHT) on just one of
the spin species: $c_{{\bf i}\uparrow}\rightarrow(-1)^{\bf i}c_{{\bf
i}\uparrow}$, leaving $c_{{\bf i}\downarrow}$ unchanged. Here $(-1)^{\bf
i}=+1(-1)$ on sublattice A(B). The SC correlations in the attractive case map
onto AF correlations in the $xy$ plane of the repulsive model, while CDW
correlations map onto AF correlations in the $z$ direction of the repulsive
model.

Under this PHT, the additional on-site term $V_{0}$, and global chemical
potential $\mu$, couple to the $z$ component of spin $S^{z}_{\bf i}=n_{{\bf
i}\uparrow}-n_{{\bf i}\downarrow}$. As has previously been noted [38, 39], the
full symmetry of the spin correlators of a Heisenberg antiferromagnet is
broken in a Zeeman field. $B_{z}$ will preferentially favor order in the $xy$
plane. Returning to the $-U$ model, we conclude that the $V_{0}$ and $\mu$
terms are both likely to favor SC over CDW formation, as we indeed observe for
$V_{0}\lesssim t$. This analysis, while very useful in the $V_{0}=0$ case,
does not appear to lend additional insight into why the enhancement of SC is
limited to stripes which are not accompanied by a $\pi$-phase shift.

Finally, note that this familiar mapping of attractive to repulsive Hubbard
model does not allow one to infer stripe formation in the attractive case from
existing results for $U>0$. (A full PHT preserves the sign of $U$, and instead
relates the model at fixed $U$ and $V_{0}$ to the model at the same $U$ but
reversed $-V_{0}$.)

##  IV Conclusions

In this paper we have shown that the well-known behavior of the AF spin order
across a charge stripe in the repulsive Hubbard model, in which the
magnetization of the sublattices on opposite sides of a stripe is
interchanged, also occurs in CDW order in the attractive Hubbard model. Thus
the fact that the magnetic order in the inter-stripe domain is distinct from
the charge pattern of the stripe appears to be irrelevant to the occurrence of
a $\pi$-phase shift. That phenomenon occurs even if the inter-stripe domain
order and stripe physics are both associated with charge degrees of freedom.
We note that in this work we have exclusively consider stripes which are
imposed externally via the potential $V_{0}$. We have not considered
spontaneous stripe formation, which is a considerably more challenging
calculation. In the much-studied repulsive Hubbard mode ($U>0)$ case,
spontaneous stripes were found early on in inhomogeneous Hartree-Fock
calculations [16, 17], but are very challenging to see in DQMC [27, 28].

Although we have found here an analogous $\pi$-phase shift between CDW across
a stripe in the attractive Hubbard model, similar to what is known for SDW
domains in the repulsive case, there is also a potential difference between
the two situations. DQMC studies of the $+U$ Hamiltonian [40] found an
increase of the stripe strength and the $\pi$-phase shift keeps enhancing
pairing monotonically. In contrast, we find here that the evolution of
superconducting correlations with stripe stength is non-monotonic, and that
when the stripes are sufficiently robust to support a $\pi$-phase shift,
pairing is suppressed. A DCA treatment [41] of the repulsive case indicates an
optimal stripe strength for pairing, but did not examine the $\pi$ phase
shift.

##  V Acknowledgments

T.Y. was supported by the joint guiding project of Natural Science Foundation
of Heilongjiang Province (Grant No. LH2019A011). R.T.S. was supported by the
grant DE‐SC0014671 funded by the U.S. Department of Energy, Office of Science.
R.M. acknowledges support from the National Natural Science Foundation of
China (NSFC) Grants No. U1930402, 12050410263, 12111530010 and No. 11974039.

##  VI Appendix

These Appendices are divided as follows. Appendix A contrasts some of the
results to the case with repulsive interactions. Appendices B and C consider
the systematic Trotter, finite temperature, and finite lattice size errors in
our DQMC calculations, showing that they do not affect our conclusions.
Appendix D analyzes the charge patterns when $\mu$ is not tuned to keep the
density in the regions between the stripes at the commensurate filling which
is optimal for CDW, $\rho_{\rm dom}=1$. Finally, Appendix E compares our DQMC
results with ED, showing that the qualitative physics is unchanged.

##  Appendix A The repulsive case: magnetic $\pi$-phase shift

A direct parallel to the results in the main text is given by the repulsive
Hubbard model ($U>0$). These can display a magnetization reversal across a
hole-rich region, leading to a magnetic $\pi$-phase shift [40], much like what
is observed in certain classes of cuprates with static stripe formation [13,
42]. Figure 6 gives the equivalent of Fig. 2, but for $U/t=+4$ instead, and
similarly tuning $\rho_{\rm dom}=1$. Here, we display the spin correlations,

$\displaystyle c_{s}({\bf r})=\big{\langle}\hat{S}_{\bf i}^{z}\hat{S}_{\bf
i+r}^{z}\big{\rangle},$ (5)

connecting the same [different] sublattice $c_{s}\big{(}{\bf r}=(0,2)\big{)}$
[$c\big{(}{\bf r}=(1,2)\big{)}$] across a stripe line. Due to the presence of
the sign problem, we focus instead at higher temperatures, $\beta t=5$.
Although qualitatively similar, the value of the stripe energy that leads to
the AF $\pi$-phase shift is much larger ($V_{0c}\simeq 4t$) than the one that
leads to the density $\pi$-phase shift ($V_{0c}\simeq t$), both of which with
$|U|/t=4$; yet, the stripe electronic filling at which the magnetization
reversal takes place is similar, $\rho_{\rm str}\simeq 0.23$.

The fact that larger $V_{0}$ is required to lower the density on a stripe in
the repulsive case is due to the presence of a Mott gap. Single occupancy of
sites is energetically preferred, and $V_{0}$ must overcome that tendency.
Indeed, when thermal and quantum fluctuations are turned off, $T=t=0$,
$V_{0c}=U$ is precisely the critical value for stripe formation. This can be
seen by considering a two site system with $V_{0}$ on one site: The
configuration $|\uparrow\,\,\,\,\downarrow\,\,\rangle$ with two singly
occupied sites has energy $-\frac{U}{4}+\big{(}-\frac{U}{4}+V_{0}\big{)}$. The
configuration $|\uparrow\downarrow\,\,\,\,*\,\,\rangle$ with a doubly occupied
and an empty site has energy $+\frac{U}{4}+\frac{U}{4}$. These become
degenerate at $V_{0c}=U$, in agreement with Fig. 6.

Meanwhile, for the attractive case, double occupation is already favored (on
alternating sites). A potential $V_{0}$ on a linear set of sites needs only to
overcome the charge alternation, an energy scale $\sim 4t^{2}/U$. This
suggests $V_{0c}\sim t$, again, in rough agreement with Fig. 2. In any case,
the difference between the values of $V_{0c}$ stresses that the two models are
not connected by PHT.

Figure 6: (a) Spin-spin correlator $c_{s}\big{(}{\bf r}=(0,2)\big{)}$ and
$c_{s}\big{(}{\bf r}=(1,2)\big{)}$, as functions of stripe strength $V_{0}$
for the repulsive Hubbard model. Magnetic $\pi$-phase shifts are formed at
$V_{0}\sim 4t$, where the reversal in sign of the spin correlation functions
takes place. (b) The total electron density $\rho$ and the density on the
stripe $\rho_{\rm str}$, as functions of stripe strength $V_{0}$. As in Fig.
2, the off-stripe electron density, $\rho_{\rm dom}$, is fixed to be $1$ by
adjusting the global chemical potential $\mu$. Here the lattice size is
$8\,\times\,8$, the on-site repulsive interactions is $U=+4\,t$, and the
inverse temperature $\beta\,t=5$.

##  Appendix B Trotter Errors

In Fig. 7, we re-compute the results for $c({\bf r})$ of Fig. 2 at a smaller
$\Delta\tau=0.0625$. The values are not significantly shifted, and conclusions
of the main text are unaltered.

Figure 7: Density-density correlator $c\big{(}{\bf r}=(0,2)\big{)}$ (blue
symbols and blue arrow in the inset) and $c\big{(}{\bf r}=(1,2)\big{)}$ (red
markers, and red arrow in the inset), as functions of stripe strength $V_{0}$,
for two different imaginary time discretizations $\Delta\tau=0.1$ (empty
markers), the value presented in the main text, and $\Delta\tau=0.0625$ (cross
markers). The results are the same to within the statistical error bars.

##  Appendix C Finite temperature and finite size effects

Here we explore the robustness of our results to lowering the temperature
further and to increasing the lattice size. Figure 8(a) compares the results
for the density correlations across a stripe of Fig. 2(a) (at $\beta\,t=10$)
with DQMC simulations at $\beta\,t=20$. The $\pi$-phase shift is still
observed, and occurs at the same $V_{0c}\sim t$.

To quantify finite size effects, we increase the lattice size to $12\times
12$, leaving the other parameters unchanged. The $\pi$-phase shift is still
observed, as shown in Fig. 8(b). There is some reduction in the density
correlations in going from $8\times 8$ to $12\times 12$ at large $V_{0}$.
However, the $16\times 16$ lattice results lie on top of those for $12\times
12$, indicating convergence to a non-zero $\pi$-phase shifted value in the
thermodynamic limit.

Figure 8: (a) Density-density correlations $c\big{(}{\bf r}=(0,2)\big{)}$
(blue markers, and blue arrow in the inset) and $c\big{(}{\bf r}=(1,2)\big{)}$
(red markers, and red arrow in the inset), as functions of stripe strength
$V_{0}$, for two different temperatures $\beta=10$ (empty markers), the value
presented in the main text, and $\beta=20$ (cross markers) for an $8\times 8$
lattice. The results are the same to within the statistical error bars. (b)
Same correlations as functions of stripe strength $V_{0}$, but now for three
different lattice sizes $8\times 8$ (empty markers), $12\times 12$ (thin cross
markers) and $16\times 16$ (full cross markers) at $\beta t=10$. Data
parameters in both panels are $U/t=-4$ and $\rho_{\rm dom}=1$.

##  Appendix D Results for general $\rho_{\rm dom}$

In Fig. 9(a), we show two density-density correlations traversing the stripe
$c\big{(}{\bf r}=(1,2)\big{)}$ and $c\big{(}{\bf r}=(0,2)\big{)}$ for
simulations in which $\mu$ is not tuned to keep $\rho_{\rm dom}=1$ but instead
the total density $\rho=1$. Since $\rho=\frac{1}{4}\rho_{\rm
str}+\frac{3}{4}\rho_{\rm dom}$ and $\rho_{\rm str}\rightarrow 0$ as $V_{0}$
becomes large, the density in the inter-stripe domains interpolates between
$\rho_{\rm dom}=1$ and $\rho_{\rm dom}=\frac{4}{3}$ as $V_{0}$ goes from
$V_{0}=0$ to $V_{0}=\infty$. This is seen in Fig. 9(b). Although $c\big{(}{\bf
r}=(0,2)\big{)}$ initially decreases from its positive value at $V_{0}=0$, it
recovers and always remains positive. There is no clear signature of
$\pi$-phase shift.

Figure 9: (a) Similar to Fig. 2(a), but with the total electron density $\rho$
fixed to be $1$. Although $c\big{(}{\bf r}=(1,2)\big{)}$ exhibits a sign
change, $c\big{(}{\bf r}=(0,2)\big{)}$ does not. When the interstripe domains
are not pinned at half-filling, complete sublattice reversal across a stripe
does not occur. (b) The electron density on and off the stripe, as functions
of stripe strength $V_{0}$, with the total electron density fixed to be $1$,
and for parameters listed in the figure.

We next performed QMC calculations fixing the total electron density to be
$0.75$, as shown in Fig. 10(a). We find that, with the increase of $V_{0}$,
$c\big{(}{\bf r}=(0,2)\big{)}$ goes from positive to negative values as should
occur for a $\pi$-phase shift. Meanwhile $c\big{(}{\bf r}=(1,2)\big{)}$ begins
already at $V_{0}=0$ with a positive value appropriate to a $\pi$ phase shift,
and remains so as $V_{0}$ grows. Thus at large $V_{0}$ both correlation
functions are consistent with sublattice reversal. The corresponding densities
are shown in Fig. 10(b).

Figure 10: Same as Fig. 9 but with the total electron density fixed to be
$0.75$, and for the parameters listed in the figure.

In Fig. 11, we show the $s$-wave pairing structure factor $P_{s}$ and the
$s$-wave pairing susceptibility $\chi_{s}$, separately, as functions of
$V_{0}$. The two panels are similar to Fig. 5, but for $\rho=1$ and
$\rho=0.75$. For $\rho=1$ there are peaks in $P_{s}$ and $\chi_{s}$ similar to
those occurring for fixed $\rho_{\rm dom}=1$. However $\rho=0.75$ shows no
enhancement of SC with the imposition of stripes.

Figure 11: The $s$-wave pairing structure factor $P_{s}$, left axis (pairing
susceptibility $\chi_{s}$, right axis), as a function of $V_{0}$, similar to
Fig. 5, but for $\rho=1$ (a) and $\rho=0.75$ (b).

##  Appendix E Exact diagonalization results

As a complement to the QMC data investigated, we also performed ED
calculations on a $2\times 4$ lattice, with the total electron density fixed
to be $1$, as shown in Fig. 12. The density correlators $c\big{(}{\bf
r}=(0,2)\big{)}$ and $c\big{(}{\bf r}=(1,2)\big{)}$ behave in a qualitatively
similar way to the QMC data in Fig. 9(a) of Appendix D.

Figure 12: ED results of the density-density correlator $c\big{(}{\bf
r}=(0,2)\big{)}$ (blue circles and blue arrow in the inset) and $c\big{(}{\bf
r}=(1,2)\big{)}$ (red squares, and red arrow in the inset), as functions of
stripe strength $V_{0}$. The total electron density is fixed to be $1$.

## References

  * Dagotto [2005] E. Dagotto, Complexity in strongly correlated electronic systems, Science 309, 257 (2005).
  * Dörr [2006] K. Dörr, Ferromagnetic manganites: spin-polarized conduction versus competing interactions, Journal of Physics D: Applied Physics 39, R125 (2006).
  * Eerenstein  _et al._ [2006] W. Eerenstein, N. Mathur, and J. F. Scott, Multiferroic and magnetoelectric materials, Nature 442, 759 (2006).
  * Israel  _et al._ [2007] C. Israel, M. J. Calderon, and N. D. Mathur, The current spin on manganites, Materials Today 10, 24 (2007).
  * Zhang  _et al._ [2019] H. Zhang, Y. Zhang, R. Li, J. Yu, W. Dong, C. Chen, K. Wang, X. Tang, and J. Chen, Room temperature hidden state in a manganite observed by time-resolved x-ray diffraction, npj Quantum Materials 4, 1 (2019).
  * Doiron-Leyraud  _et al._ [2007] N. Doiron-Leyraud, C. Proust, D. LeBoeuf, J. Levallois, J.-B. Bonnemaison, R. Liang, D. Bonn, W. Hardy, and L. Taillefer, Quantum oscillations and the Fermi surface in an underdoped high-${T}_{c}$ superconductor, Nature 447, 565 (2007).
  * Wu  _et al._ [2011] T. Wu, H. Mayaffre, S. Krämer, M. Horvatić, C. Berthier, W. Hardy, R. Liang, D. Bonn, and M.-H. Julien, Magnetic-field-induced charge-stripe order in the high-temperature superconductor YBa2Cu3O7+y, Nature 477, 191 (2011).
  * Wu  _et al._ [2015] T. Wu, H. Mayaffre, S. Krämer, M. Horvatić, C. Berthier, W. Hardy, R. Liang, D. Bonn, and M.-H. Julien, Incipient charge order observed by NMR in the normal state of YBa2Cu3O7+y, Nature Communications 6, 1 (2015).
  * Ghiringhelli  _et al._ [2012] G. Ghiringhelli, M. Le Tacon, M. Minola, S. Blanco-Canosa, C. Mazzoli, N. Brookes, G. De Luca, A. Frano, D. Hawthorn, F. He, _et al._ , Long-range incommensurate charge fluctuations in (Y, Nd) Ba2Cu3O6+x, Science 337, 821 (2012).
  * Chang  _et al._ [2012] J. Chang, E. Blackburn, A. Holmes, N. B. Christensen, J. Larsen, J. Mesot, R. Liang, D. Bonn, W. Hardy, A. Watenphul, _et al._ , Direct observation of competition between superconductivity and charge density wave order in YBa2Cu3O6.67, Nature Physics 8, 871 (2012).
  * Comin and Damascelli [2016] R. Comin and A. Damascelli, Resonant x-ray scattering studies of charge order in cuprates, Annual Review of Condensed Matter Physics 7, 369 (2016).
  * Hoffman  _et al._ [2002] J. Hoffman, E. W. Hudson, K. Lang, V. Madhavan, H. Eisaki, S. Uchida, and J. Davis, A four unit cell periodic pattern of quasi-particle states surrounding vortex cores in Bi2Sr2CaCu2O8+δ, Science 295, 466 (2002).
  * Tranquada  _et al._ [1995] J. Tranquada, B. Sternlieb, J. Axe, Y. Nakamura, and S. Uchida, Evidence for stripe correlations of spins and holes in copper oxide superconductors, Nature 375, 561 (1995).
  * Julien [2015] M.-H. Julien, Magnetic fields make waves in cuprates, Science 350, 914 (2015).
  * Zheng  _et al._ [2017] B.-X. Zheng, C.-M. Chung, P. Corboz, G. Ehlers, M.-P. Qin, R. M. Noack, H. Shi, S. R. White, S. Zhang, and G. K.-L. Chan, Stripe order in the underdoped region of the two-dimensional Hubbard model, Science 358, 1155 (2017).
  * Zaanen and Gunnarsson [1989] J. Zaanen and O. Gunnarsson, Charged magnetic domain lines and the magnetism of high-${T}_{c}$ oxides, Phys. Rev. B 40, 7391 (1989).
  * Poilblanc and Rice [1989] D. Poilblanc and T. M. Rice, Charged solitons in the Hartree-Fock approximation to the large-$U$ Hubbard model, Phys. Rev. B 39, 9749 (1989).
  * Kato  _et al._ [1990] M. Kato, K. Machida, H. Nakanishi, and M. Fujita, Soliton lattice modulation of incommensurate spin density wave in two dimensional Hubbard model -a mean field study-, J. Phys. Soc. Jpn. 59, 1047 (1990).
  * Yamase  _et al._ [2016] H. Yamase, A. Eberlein, and W. Metzner, Coexistence of incommensurate magnetism and superconductivity in the two-dimensional Hubbard model, Phys. Rev. Lett. 116, 096402 (2016).
  * White and Scalapino [1998] S. R. White and D. J. Scalapino, Density matrix renormalization group study of the striped phase in the 2D t-J model, Phys. Rev. Lett. 80, 1272 (1998).
  * White and Scalapino [2003] S. R. White and D. Scalapino, Stripes on a 6-leg Hubbard ladder, Phys. Rev. Lett. 91, 136403 (2003).
  * Hager  _et al._ [2005] G. Hager, G. Wellein, E. Jeckelmann, and H. Fehske, Stripe formation in doped Hubbard ladders, Phys. Rev. B 71, 075108 (2005).
  * [23] H.-C. Jiang and S. A. Kivelson, Stripe order enhanced superconductivity in the Hubbard model, arXiv:2105.07048 .
  * Fleck  _et al._ [2000] M. Fleck, A. I. Lichtenstein, E. Pavarini, and A. M. Oleś, One-dimensional metallic behavior of the stripe phase in La2-xSrxCuO4, Phys. Rev. Lett. 84, 4962 (2000).
  * Vanhala and Törmä [2018] T. I. Vanhala and P. Törmä, Dynamical mean-field theory study of stripe order and $d$-wave superconductivity in thetwo-dimensional Hubbard model, Phys. Rev. B 97, 075112 (2018).
  * Chang and Zhang [2010] C.-C. Chang and S. Zhang, Spin and charge order in the doped Hubbard model: Long-wavelength collective modes, Phys. Rev. Lett. 104, 116402 (2010).
  * Huang  _et al._ [2017] E. W. Huang, C. B. Mendl, S. Liu, S. Johnston, H.-C. Jiang, B. Moritz, and T. P. Devereaux, Numerical evidence of fluctuating stripes in the normal state of high-${T}_{c}$ cuprate superconductors, Science 358, 1161 (2017).
  * Huang  _et al._ [2018] E. W. Huang, C. B. Mendl, H.-C. Jiang, B. Moritz, and T. P. Devereaux, Stripe order from the perspective of the Hubbard model, npj Quantum Materials 3, 1 (2018).
  * Qin  _et al._ [2020] M. P. Qin, C. M. Chung, H. Shi, E. Vitali, C. Hubig, U. Schollwöck, S. R. White, and S. W. Zhang, Absence of superconductivity in the pure two-dimensional hubbard model, Phys. Rev. X 10, 031016 (2020).
  * Corboz  _et al._ [2014] P. Corboz, T. M. Rice, and M. Troyer, Competing states in the t-J model: Uniform $d$-wave state versus stripe state, Phys. Rev. Lett. 113, 046402 (2014).
  * Zheng and Chan [2016] B.-X. Zheng and G. K.-L. Chan, Ground-state phase diagram of the square lattice Hubbard model from density matrix embedding theory, Phys. Rev. B 93, 035126 (2016).
  * [32] M. Qin, T. Schäfer, S. Andergassen, P. Corboz, and E. Gull, The Hubbard model: A computational perspective, arXiv:2104.00064 .
  * Loh  _et al._ [1990] E. Loh, J. Gubernatis, R. Scalettar, S. White, D. Scalapino, and R. Sugar, Sign problem in the numerical simulation of many-electron systems, Phys. Rev. B 41, 9301 (1990).
  * Troyer and Wiese [2005] M. Troyer and U.-J. Wiese, Computational complexity and fundamental limitations to fermionic quantum Monte Carlo simulations, Phys. Rev. Lett. 94, 170201 (2005).
  * Blankenbecler  _et al._ [1981] R. Blankenbecler, D. J. Scalapino, and R. L. Sugar, Monte carlo calculations of coupled boson-fermion systems. i, Phys. Rev. D 24, 2278 (1981).
  * White  _et al._ [1989] S. R. White, D. J. Scalapino, R. L. Sugar, E. Y. Loh, J. E. Gubernatis, and R. T. Scalettar, Numerical study of the two-dimensional hubbard model, Phys. Rev. B 40, 506 (1989).
  * Hirsch [1983] J. E. Hirsch, Discrete Hubbard-Stratonovich transformation for fermion lattice models, Phys. Rev. B 28, 4059 (1983).
  * Scalettar  _et al._ [1989] R. Scalettar, D. Scalapino, and N. Bickers, Competition of pairing and Peierls–charge-density-wave correlations in a 2–d electron-phonon model, Phys. Rev. B 40, 197 (1989).
  * Moreo and Scalapino [1991] A. Moreo and D. J. Scalapino, Two-dimensional negative-U Hubbard model, Phys. Rev. Lett. 66, 946 (1991).
  * Mondaini  _et al._ [2012] R. Mondaini, T. Ying, T. Paiva, and R. T. Scalettar, Determinant quantum Monte Carlo study of the enhancement of d-wave pairing by charge inhomogeneity, Phys. Rev. B 86, 184506 (2012).
  * Maier  _et al._ [2010] T. A. Maier, G. Alvarez, M. Summers, and T. C. Schulthess, Dynamic cluster quantum Monte Carlo simulations of a two-dimensional Hubbard model with stripelike charge-density-wave modulations: Interplay between inhomogeneities and the superconducting state, Phys. Rev. Lett. 104, 247001 (2010).
  * Tranquada  _et al._ [2004] J. M. Tranquada, H. Woo, T. G. Perring, H. Goka, G. D. Gu, G. Xu, M. Fujita, and K. Yamada, Quantum magnetic excitations from stripes in copper oxide superconductors, Nature 429, 534 (2004).


# Chain-structured neural architecture search for financial time series
forecasting

Denis Levchenko
Efstratios Rappos*
Shabnam Ataee*
Biagio Nigro
Stephan Robert* School of Engineering and Management Vaud (HEIG-VD),
University of Applied Sciences and Arts Western Switzerland (HES-SO), Yverdon-
les-Bains, SwitzerlandPredictive Layer SA, Rolle, Switzerland

###### Abstract

We compare three popular neural architecture search strategies on chain-
structured search spaces: Bayesian optimization, the hyperband method, and
reinforcement learning in the context of financial time series forecasting.

##  1 Introduction

Deep neural networks have been very successful in a wide variety of tasks over
the last two decades. In large part their success is attributed to their
ability to perform very well without major manual feature engineering required
when compared to more classical techniques. The hierarchical structure of the
networks weighs and extracts important features automatically from the data
during the learning process [GBC16]. However, the exact architecture of the
neural network still has to be prescribed manually by the user. This led to
the development of so-called auto-ML techniques that aim to automate this
process. In the context of deep neural networks, auto-ML has a very large
overlap with neural architecture search (NAS), itself having a large overlap
with hyperparameter optimization.

A lot of research has been done in NAS in recent years, see [WSS+23] for an
overview and insights from over 1000 papers. However, most work focused on
computer vision or natural language applications, with less investigation into
architectures for analyzing time series data. In this work, we attempt to
bridge this gap by evaluating and comparing the performance of three popular
simple NAS strategies on 3 distinct yet similar time series datasets.

A seminal survey article on NAS [EMH19] breaks down the search into 3
components: search spaces, search strategies and performance estimation
strategies. Once a search space is set up, an algorithm following the search
strategy explores the search space looking for the best neural network
configuration, evaluated and eventually selected by the performance estimation
strategy. In Section 2, we discuss the exact practical problem from project
partner Predictive Layer SA that our neural networks were trying to solve. We
discuss the search spaces we set up for NAS in Section 3, the search
strategies in Section 5 and the performance estimation together with the
general challenges of our problem in Section 4. Our results are summarized and
discussed in Section 6 and Section 7 respectively. Finally, we outline
possible future directions for this research in Section 8.

##  2 Data and problem formulation

This project was a collaboration with Predictive Layer SA, who provided real-
world customer data and prediction requirement. The data under consideration
are financial multivariate daily time series in tabular form, each row
corresponding to one day and each column a separate feature. The task consists
of predicting whether the target feature will increase or decrease 5 days (or
10 days, depending on the dataset) in the future. Thus the problem is
essentially binary classification. The models output a number between $0$ and
$1$ as output, interpreted as probability that the target feature will grow on
a given day.

The challenge comes with financial markets being notoriously hard to predict
[HH19]. Our three source datasets that we report on in Section 6 are for
Japanese, German, and US bonds and come with an additional challenge of very
high number of features (around 1000 each), comparable with the number of
training examples (under 4000), which is very low in a deep learning setting
[GBC16], especially relatively to the number of features. To circumvent the
high dimensionality, the number of features was reduced approximately by a
factor of $3$ by removing features in the prepared datasets ‘time-derived’
from original features, such as by lagging original features by fixed time
steps or taking the mean over the same time steps. Such time-derived features
are standard in the industry [HH19], but are fundamentally used because the
algorithms applied to the data can only take in one row of data as input. All
the more sophisticated neural network architectures take in a sequence of rows
of data, making time-derived features redundant. See Section 3 for more
details. Additionally, principal component analysis (PCA) was applied to
further reduce the number of features. Based on the scree plot coming from PCA
analysis, two breakpoints were tested: the point of highest curvature (the
usual choice) and the point where the graph flattens, retaining most of the
information. Thus for each of the 3 source datasets, we obtained 4 variations:
no modification, time-derived features removed and 2 PCA transformations,
resulting in 12 datasets total.

##  3 Architecture types and their search spaces

In NAS, the search spaces and search strategies go hand in hand and are
designed together essentially at the same time. The simplest among the search
spaces are chain-structured. Search spaces with a chain structure feature a
straightforward architectural topology: a sequential chain of operational
layers. These configurations frequently utilize cutting-edge manually designed
architectures as their foundational framework, trying to find the best
configuration for the specific dataset at hand by essentially hyperparameter
optimization. Despite their simplicity, chain-structured search spaces often
yield very good results [WSS+23]. The biggest downsides of this approach when
compared to the more complicated cell-based and one-shot approaches are its
limited flexibility and scalability when applied to very large datasets. While
being limited to a pre-defined overall structure prevents discovering truly
novel architectures, in practice well-designed and well-understood
architectures usually perform best [WSS+23]. Scalability becomes important
when dealing with large image and natural language datasets that most research
on NAS focuses on, however in the case of time series and especially our small
datasets, our neural networks train in seconds. Therefore scalability was not
an issue for us, and thus we focused on chain-structured search spaces.

For the foundational frameworks the search spaces were built upon, we
concentrated on simple well-known architecture types, like feedforward (FFNN,
also known as multi-layer perceptron), convolutional (CNN), and recurrent
(RNN) networks. This was done both for their proven perforamnce in time-series
context and due to them having fewer parameters than most cutting-edge
architectures, which is more appropriate when dealing with very little data as
in our case. However, we also examined the state of the art Temporal Fusion
Transformer (TFT) architecture for comparison (see below). Each of these
architectures was adapted to the problem at hand and for each of them a search
space of hyperparameters was set up to determine the exact configuration. The
search strategies discussed in Section 5 then optimized each of these
separately, and the best-performing ones were compared. This approach is
similar to what [ZL17] did with a single type of CNN for image classification.

###  3.1 Feedforward networks

Feedforward networks have the advantage of being the simplest type of deep
neural networks. The main disadvantage for our problem comes from the fact
that they cannot take the time-dimension into account explicitly. Feedforward
networks take a row of numbers as input, one row of our datasets corresponding
to one day information for various time series. Implicitly, the time
information is taken into account in the form of delayed (or averaged)
features – a column in the dataset that is essentially the same as another
column with an added delay of a couple days (or averaged over a couple days,
respectively), a standard practice in the field [HH19].

As far as NAS is concerned, our hyperparameters for optimization are the
number of hidden layers and the number of units per hidden layer. Dropout
layers were used for regularization and the dropout rate was another
hyperparameter optimized. Although not part of the architecture itself, the
learning rate is a crucial hyperparameter for any learning algorithm that was
also optimized.

###  3.2 Convolutional neural networks

Although designed for and used mainly with images, convolutional neural
networks (CNNs) are suitable for any kind of grid-like data, especially when
there is a natural fixed notion of distance between neighboring data points,
such as the case for time series [GBC16].

Now, a ‘chunk’ of our dataset – a subtable of fixed length (in time) can be
passed as a single input to our model. Evaluating the history from the last,
say, $15$ days explicitly, the model is making a prediction for the following
day as before.

Both one- and two-dimensional CNNs can be used for our problem. One-
dimensional (1D) CNNs only perform the convolution operation over the time
dimension, treating separate features as parallel input channels, just like
the red, green and blue channels are considered when working with images. Two-
dimensional (2D) CNNs take each chunk as a single-channel ‘image’, performing
the convolution operation over both the time dimension and across the various
features. The latter seems counterintuitive as there is no fixed distance
between the features, but [HH19] found success with this technique. After
experimenting with 2D CNNs on our data, we found them a lot more
computationally expensive than others, while often performing worse. For this
reason 2D CNNs were later dropped from our consideration.

The length in time of the chunks passed as individual examples (‘chunk
length’) was a crucial hyperparameter to optimize for all of the CNNs. The
number of convolutional layers (including pooling and activation functions)
was fixed at $3$ for all CNN architectures to reduce the size of the search
spaces, but the kernel sizes were optimized. The number of convolution filters
per layer was optimized too, but it works differently for different
implementations, as discussed below. Dropout layers were used once again for
regularization and the dropout rate was optimized too, together with the
learning rate.

####  3.2.1 Depthwise 1D CNNs

Among 1D CNNs, there are two distinct ways to perform the convolution
operation. In _depthwise_ convolution, each channel (input feature) is passed
separately through its own convolutional and pooling layers. Each of these
then gives a fixed number (the number of separate convolution filters) of
output channels per input channel, interpreted as a summary of what happened
to this feature in the time span considered. The convolution and pooling
operations will shrink the length (in time) of each output channel, possibly
to $1$ (otherwise flattening is applied), and all of these are then passed to
a single dense layer that makes a prediction based on these summaries.

Although easy to interpret, this method had the problem of blowing up the
dimensionality of the problem, as each of the features gets its own set of
convolution parameters. This made computation infeasible and more work is
needed to be carried out for this architecture.

####  3.2.2 Ordinary 1D CNNs

In basic (we call them ‘ordinary’) 1D convolution, all input channels are
mixed already in the first (and subsequent) convolution layers. The output of
convolution is flattened and passed to a single dense layer for making the
final prediction.

###  3.3 Recurrent neural networks

Recurrent neural networks (RNNs) were designed for and are very successful at
treating sequential data [GBC16] and are thus the natural choice for our
problem. Both simple RNNs and more sophisticated long short-term memory (LSTM)
RNNs were implemented and optimized separately, but they have the same search
spaces of hyperparameters. In both cases, _stacked_ RNNs were used, where the
output of one RNN is taken as input to the next one. The number of such
stacked layers was a hyperparameter to optimize. Each RNN layer also has its
number of hidden units, which was also optimized. The chunk length, here
naturally interpreted as the length of the input sequence, was once again
optimized, together with dropout and learning rates as before.

LSTMs provide many advantages over basic RNNs [GBC16] without too much
additional computational cost, especially in the case of our small datasets.
For this reason after initial trials, we proceeded with just LSTMs for the
final tests.

###  3.4 Temporal fusion transformer

Forecasting across multiple time horizons often involves a intricate
amalgamation of inputs, encompassing static (i.e., time-invariant) covariates,
known future inputs, and additional exogenous time series observed solely in
the past. The challenge lies in handling this complexity without prior
knowledge of how these inputs interact with the target variable. While various
deep learning approaches have been proposed, they often manifest as ‘black-
box’ models, lacking transparency regarding their utilization of the diverse
input types encountered in practical scenarios.

In a recent work [LALP21], the Temporal Fusion Transformer (TFT) was
introduced as an innovative attention-based architecture addressing this
issue. TFT not only achieves high-performance multi-horizon forecasting but
also provides interpretable insights into temporal dynamics. In order to
capture temporal relationships at different scales, TFT incorporates recurrent
layers for local processing and interpretable self-attention layers for
modelling long-term dependencies. Specialized components within TFT are
employed to select relevant features, and a series of gating layers
effectively suppress unnecessary components, resulting in high performance
across a broad spectrum of scenarios.

[LALP21] demonstrated significant performance enhancements over existing
benchmarks across various real-world datasets. The study also highlighted
three practical use cases illustrating the interpretability of TFT, showcasing
its efficacy in shedding light on the decision-making process. The PyTorch
Forecasting package [PyT23] has on open-source implementation of the model. It
comes with its own optimizer, selecting the optimal number of attention heads,
the network’s hidden size, and learning and dropout rates.

Although very promising and successful for other datasets, we found
performance of TFT rather poor on our data: the models reduce to a trivial
binary predictor, always predicting $1$ or $0$ no matter the input test data.
We believe this was largely due to the lack of data. As previously mentioned,
novel complex architectures like TFT have tens if not hundreds of thousands of
parameters which need a lot of data to train well. It would be interesting to
revisit TFT and other models derived from it that have recently appeared in
the literature on other, bigger datasets.

##  4 Challenges and performance estimation

There are many evaluation metrics applicable to the binary classification
problem. Although common for other classification tasks, accuracy ($ACC$) is
not a good metric as a trivial predictor that always gives the same output no
matter the input can have high accuracy in the case of an unbalanced dataset.
Balanced accuracy ($bACC$), the average between the true negative and true
positive rates, mitigates this problem. Another popular alternative is the
$F_{1}$ score, the harmonic mean of precision and recall. The area under the
receiver operating characteristic (ROC) curve, or $AUC$ for short is another
good choice. It is hard to say which one is better, and a good strategy could
mix a combination of these.

Predicting financial markets is extremely challenging and this is reflected in
our test metrics always being very close to $0.5$, i.e. the predictors being
very close to random. This would be regarded as bad performance in many areas,
but for financial datasets it is quite the norm: many companies make
successful investments using algorithms (whether traditional or ML-based)
being only marginally better than a random predictor [HH19]. Further, neural
networks in particular in this context suffer from comparatively high variance
depending on the random seed chosen during training due to their stochastic
nature [HH19]. Even identical neural network architectures, trained in exactly
the same way but with different random seeds can have vastly different
performance on the test dataset, e.g. $0.56$ on a ‘good’ random seed and
$0.48$ (i.e. worse than random) on a ‘bad’ one. This is a problem for NAS as
what seemed like a good (or bad) neural network configuration might have just
been a lucky (or unlucky respectively) seed. To mitigate this problem, we
trained and tested each architecture configuration chosen $15$ times with
different random seeds, averaging the test metrics and using the averages as
our guide. At the end of optimization, the best performing configuration was
trained and tested $50$ times (again, with different random seeds every time)
to further eliminate random seed variation.

##  5 Search strategies

We started with $7$ different broad architecture types: FFNN, depthwise 1D
CNN, ordinary 1D CNN, 2D CNN, simple RNN, LSTM, and TFT, but for final
optimization studies only FFNN, ordinary 1D CNN, and LSTM were used, as
discussed in Section 3. Each of them have their own search space of
hyperparameters determining the exact network architecture.

Classical approaches for hyperparameter tuning in machine learning, such as
grid search and random search, have been widely used due to their simplicity
[BB12]. However, these methods can be computationally expensive and
inefficient, especially when dealing with high-dimensional search spaces or
expensive objective functions [BBBK11]. The survey on NAS [WSS+23] favourably
featured three other popular search strategies: Bayesian optimization,
reinforcement learning, and the hyperband method which are the three methods
we explored extensively.

Fundamentally, all of these suggest a specific configuration of
hyperparameters of a neural network. The network is then trained and tested
(in our case repeatedly so, see Section 4), the results saved, and the
algorithm suggests the next configuration to try. This continues for a set
number of trials (in our studies, we set it at $300$) and in the end the best
performing configuration is selected.

###  5.1 Bayesian optimization

Bayesian optimization has emerged as a popular alternative to grid and random
search for hyperparameter tuning, as it efficiently explores the search space
and intelligently guides the optimization process [SLA12]. One widely used
technique in Bayesian optimization is based on Gaussian Processes (GP), which
fit the objective function using Gaussian processes as probabilistic models
and leverages acquisition functions to balance exploration and exploitation
[RW06]. GPs provide a flexible, non-parametric tool for modeling complex
functions, allowing for uncertainty quantification and adaptation to new data
[SSW+16].

However, GP-based methods have some limitations, such as the inability to
handle categorical features or dependent parameters without ad hoc
modifications [BCdF10]. This has led to the development of alternative
methods, such as the Tree-structured Parzen Estimator (TPE), which has gained
attention for its scalability and effectiveness [BBBK11]. TPE models the
search space using hierarchical Parzen estimators, which adaptively partition
the space to model complex, high-dimensional functions [BBBK11].

TPE is particularly suitable for addressing the shortcomings of GP-based
methods, as it can naturally handle categorical features and dependent
parameters without requiring extensive adaptations. Performance comparisons
between TPE, GP-based methods, and classical search techniques have shown that
TPE often yields improved results across various tasks and datasets [EFH+13].

A notable implementation of TPE is available in the Optuna package, a
versatile optimization library that has demonstrated its utility in various
machine learning tasks [ASY+19].

###  5.2 Reinforcement learning approach

The reinforcement learning method consists of treating the neural architecture
search problem as a reinforcement learning problem. We construct an iterative
approach where the parameters (and structure) of the neural architecture are
initially chosen randomly, just as in the case of Bayesian optimization. After
the neural network is constructed, it is trained for a number of periods using
a train and validation time-series dataset, sometimes even a number of times
in order to correctly assess its performance in predicting the time series,
and exclude issues relating to the choice of initial random seed. At the end
of the run, a reward function is calculated, which represents the incremental
increase in performance. In this way choices of parameters which produce
“good” neural networks are rewarded, whereas those which produce badly-
performing neural networks are penalized. We let the system evolve based for a
number of iterations (reinforcement learning trials) and at the end we obtain
the best overall neural network configuration.

The specific method we used introduced in [ZL17] uses an RNN as a controller,
trained with REINFORCE. An open-source implementation of the controller RNN to
use with image data is available at [Maj18], which we adpated to work with our
time series data.

The reinforcement learning approach has several advantages compared to the
Bayesian optimization: it is able to capture all the dynamics and
interdependencies between the parameters, something that is difficult with
classic bayesian optimization without testing every possible combination of
parameters. Moreover we can more easily use customized reward functions (for
example, in one case the reward function used was the average of the AUC and
F1 scores) which can better capture the performance of the neural network in
the presence of imbalanced datasets.

###  5.3 Hyperband

Hyperband is a method introduced in 2018 by Li et al., bringing in a clever
tactic known as ‘bandit’ to make the hyperparameter optimization process very
effective [LJD+18]. It explores various combinations, paying more attention to
the promising ones and discarding the under-performers.

What makes Hyperband notable are its efficiency, scalability, and the ability
to use multiple GPUs simultaneously for faster optimization. It is compatible
with popular machine learning libraries like Scikit-learn and Keras, embedded
in the KerasTuner framework [OBL+19].

Hyperband is designed for the efficient exploration and optimization of
hyperparameter configurations. It selectively targets promising models in the
search space, optimizing computational resources and time. The algorithm
employs successive halving, a subprocess that involves eliminating less
promising models and further training the survivors.

In successive halving, you provide a list of models ($models$), the total
epochs needed for full training ($max\\_epochs$), the initial training epochs
for all models ($start\\_epochs$), and a factor determining the speed of
reducing models and increasing training epochs for the survivors ($factor$).
This process iterates until all models are fully trained, resulting in the
identification of the best model and its validation loss.

To achieve a balance between exploration and optimization, Hyperband executes
successive halving multiple times, labeling each iteration as a bracket. The
number of brackets is determined by a parameter, $s\\_max$, approximately
equal to $log(max\\_epochs,factor)$. In each bracket, the number of models
decreases by a factor, and new models are generated. Hyperband ensures a
robust exploration by randomly generating models for each bracket, preventing
too few models in a bracket. The initial training epochs for each bracket
strictly increase by a factor.

Practically implemented in KerasTuner as the Hyperband tuner, it is available
for single or multiple GPUs for parallel tuning. Users can specify $factor$
and $max\\_epochs$, along with the number of brackets using
$hyperband\\_iterations$. This allows control over the search time while
efficiently exploring the hyperparameter space [SJH22].

##  6 Results

The data was split into (in historical order): first $70\%$ for training, the
following $20\%$ for validation used in the evaluation during optimization
runs, and the last $10\%$ for final testing on unseen data. When the
optimization studies were completed, the architectures performing best on the
validation dataset were selected and evaluated on the test set. This involved
retraining and testing each network $50$ times, and average performance was
compared. This was done for each search strategy, architecture type and
dataset combination.

Our best results were obtained on the Germany dataset, with time-derived
redundant features removed but no PCA applied. The LSTM with parameters
selected by the hyperband method applied to the unseen test data achieved an
$AUC$ score of $0.56$ on average over $50$ test runs starting from different
random seeds (see Section 4), with $0.05$ standard deviation. The same
architecture showed balanced accuracy score of $0.54\pm 0.04$. 1D CNNs were
also able to achieve good performance on this dataset, with the best
architecture provided by Bayesian optimization giving an $AUC$ score of
$0.54\pm 0.05$ and a high $F1$ score of $0.6\pm 0.06$.

For the Japan dataset, the redundant time-derived features removed and PCA
applied retaining most of the information in the dataset while further cutting
the number of features approximately in half gave the best results. The best
performing architecture was a 1D CNN coming from Bayesian optimization,
achieving an $AUC$ score of $0.54\pm 0.03$ over $50$ test runs, with a high
$F1$ score of $0.65\pm 0.02$.

The US dataset was most challenging for us. Although tuned 1d CNNs gave an AUC
score of $0.6\pm 0.02$ on validation data (not used for training) in our
repeated testing, and an individual lucky random seed gave an $AUC$ of $0.58$
on test data, on average no architecture could achieve $AUC$ over $0.5$ on the
test dataset. We believe something very significant happened in the US market
in the last year or so of our data that changed market dynamics.

###  6.1 Architectures and search strategies compared

Our results show LSTMs and 1D CNNs outperforming simple feedforward networks
for nearly all search strategies and datasets, as expected, but were neck in
neck between each other.

Figure 1: Best performing architectures selected by each search strategy on
each dataset. Every point represents average $AUC$ score and standard
deviation on test data after retraining the selected architecture $50$ times.
Type of the neural network chosen by the search strategy is displayed above
each point, search strategies are color-coded.

For search strategies, overall the hyperband method and bayesian optimization
showed better performance than the reinforcement learning based approach, but
all three were very close to each other, well within standard deviation as
shown in Figure 1. Optimization time was very similar between the three
methods, mostly determined by the architecture type, total number of trials,
repetitions per trial (to minimize variation due to different random seeds),
and the number of epochs to train each neural network. These were set manually
by us to 300, 15, and 80 respectively, for a total optimization time of about
12 hours per architecture type using a single Nvidia Quadro RTX 4000 GPU, with
feedforward networks training faster than CNNs and LSTMs due to their lower
complexity.

Practically, the reinforcement learning approach required the most work to
adapt from the open-source implementation of the method in [ZL17], and has
many parameters of its own to be optimized for best performance. In contrast,
both the hyperband implementation in KerasTuner and the Bayesian optimization
from Optuna are easy to use and well-supported. The Optuna implementation is
flexible in that nearly the same code can be used to optimize hyperparameters
of a different kind of solver other than neural networks, such as decision
trees-based algorithms popular for time series prediction. Meanwhile,
KerasTuner implementation of hyperband has advantage of easy parallelization
over several GPUs (although we didn’t have the hardware to make use of it).

##  7 Discussion

Overall, we didn’t observe significant convergence over time in the
hyperparameter optimization. An example history of optimization plot is shown
in Figure 2. This was a general problem present for all neural architecture
types, optimization strategies and datasets.

Figure 2: Bayesian optimization history for the ordinary 1D CNN architecture
on the US dataset. Each point represents average $AUC$ score on the validation
dataset over $15$ runs for the same network configuration.

We believe all search strategies were working as they should, but the main
problem of financial time series prediction on our datasets was just too
difficult to decisively improve upon through neural architecture search alone.
Many steps were taken to improve overall performance. In the beginning the
networks were heavily overfitting, where we struggled to get any metrics above
$0.5$ even on the validation data. The feature reduction measures such as
removing time-derived features and PCA discussed in Section 2 helped mitigate
this. Unfortunately, this still did not always translate to good performance
on the test dataset, even though we heavily used dropout layers for
regularization. We also tried converting the main binary clasification problem
to regression, but this did not result in substantial improvement.

As a consequence, simply taking the final suggested architecture by either
search strategy didn’t give the best results. The best performing ones were
also not always the best choice due to the random seed variance discussed in
Section 4. Even though each network configuration was trained and tested $15$
times during the optimization process, training and testing the best
architectures more (e.g., 50) times sometimes showed worse performance,
indicating that the issue of random seed variation was not fully resolved in
our approach.

Nevertheless, all search strategies were still successful in highlighting best
choices for network parameters, at least for the validation set, which is most
evident from slice plots as in Figure 3. For bayesian optimization, the best
parameters for each architecture were read off slice plots like this. However
we note that following this approach somewhat negates the main advantages in
going for more advanced search strategies as compared to simple random search.

Figure 3: Slice plot for bayesian optimization of LSTMs on the US dataset.
Each point represents average $AUC$ score on the validation dataset over $15$
runs for the same network configuration. It is clear that setting chunk_length
(i.e. the length in time of a sequence passed to the LSTM for individual
prediction) to $10$ gives the best results, while other parameters are less
impactful.

To combat the big variance between networks starting from different random
seeds (see Section 2 and Figure 1), we tried filtering the outputs,
disregarding those with probabilities close to $0.5$ (interpreted as those
where the network ‘was not sure’), a standard approach [GBC16]. Although this
improved overall performance a little on average, it instead greatly increased
the variance. For a fixed architecture, selecting networks trained from random
seeds giving good performance on validation set did not translate to better
performance on the test set. Our best suggestion to combat the variance is to
use an ensemble model built from many incarnations of a chosen network
architecture.

##  8 Future work

Performing the same studies on more and publicly available time-series
datasets, especially non-financial data would shed more light whether our
findings were specific to our data or a more general feature.

In addition, cell-based and hierarchical search spaces [WSS+23] could be
explored in the future in the context of time-series forecasting. These are
based on searching over structural blocks within a network piece by piece
instead of navigating through a grid-like space as in the standard approach we
took. Finally, one-shot NAS techniques such as supernet-based methods
(differentiable and not) and hypernetworks [WSS+23] could be explored. With
one-shot techniques, one trains a single (massive) super/hyper-network, which
can be used to subsample and evaluate many smaller networks to find the
optimal one, without re-training. This is opposed to the standard approach of
training and testing many architecture configurations that we took.

## References

  * [ASY+19] Takuya Akiba, Shotaro Sano, Toshihiko Yanase, Takeru Ohta, and Masanori Koyama.  Optuna: A next-generation hyperparameter optimization framework.  Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, 2019.
  * [BB12] James Bergstra and Yoshua Bengio.  Random search for hyper-parameter optimization.  J. Mach. Learn. Res., 13:281–305, 2012.
  * [BBBK11] James Bergstra, Rémi Bardenet, Yoshua Bengio, and Balázs Kégl.  Algorithms for hyper-parameter optimization.  In J. Shawe-Taylor, R. Zemel, P. Bartlett, F. Pereira, and K.Q. Weinberger, editors, Advances in Neural Information Processing Systems, volume 24. Curran Associates, Inc., 2011.
  * [BCdF10] Eric Brochu, Vlad M. Cora, and Nando de Freitas.  A tutorial on bayesian optimization of expensive cost functions, with application to active user modeling and hierarchical reinforcement learning.  ArXiv, abs/1012.2599, 2010.
  * [EFH+13] Katharina Eggensperger, Matthias Feurer, Frank Hutter, James Bergstra, Jasper Snoek, Holger H. Hoos, and Kevin Leyton-Brown.  Towards an empirical foundation for assessing bayesian optimization of hyperparameters.  2013\.
  * [EMH19] Thomas Elsken, Jan Hendrik Metzen, and Frank Hutter.  Neural architecture search: A survey.  Journal of Machine Learning Research, 20(55):1–21, 2019.
  * [GBC16] Ian J. Goodfellow, Yoshua Bengio, and Aaron Courville.  Deep Learning.  MIT Press, Cambridge, MA, USA, 2016.
  * [HH19] Ehsan Hoseinzade and Saman Haratizadeh.  Cnnpred: Cnn-based stock market prediction using a diverse set of variables.  Expert Systems with Applications, 129:273–285, 2019.  https://doi.org/10.1016/j.eswa.2019.03.029.
  * [LALP21] Bryan Lim, Sercan Ö. Arık, Nicolas Loeff, and Tomas Pfister.  Temporal fusion transformers for interpretable multi-horizon time series forecasting.  International Journal of Forecasting, 37(4):1748–1764, 2021.
  * [LJD+18] Lisha Li, Kevin Jamieson, Giulia DeSalvo, Afshin Rostamizadeh, and Ameet Talwalkar.  Hyperband: A novel bandit-based approach to hyperparameter optimization.  Journal of Machine Learning Research, 18(185):1–52, 2018.
  * [Maj18] Somshubra Majumdar.  neural-architecture-search.  https://github.com/titu1994/neural-architecture-search, 2018.  Accessed: 2023-02-21.
  * [OBL+19] Tom O’Malley, Elie Bursztein, James Long, François Chollet, Haifeng Jin, Luca Invernizzi, et al.  Kerastuner.  https://github.com/keras-team/keras-tuner, 2019.
  * [PyT23] PyTorch Forecasting Contributors.  Pytorch forecasting.  https://github.com/jdb78/pytorch-forecasting, 2023.
  * [RW06] Carl Edward Rasmussen and Christopher K. I. Williams.  Gaussian processes for machine learning. Adaptive computation and machine learning. MIT Press, 2006.
  * [SJH22] Q. Song, H. Jin, and X. Hu.  Automated Machine Learning in Action.  Manning, 2022.
  * [SLA12] Jasper Snoek, Hugo Larochelle, and Ryan P Adams.  Practical bayesian optimization of machine learning algorithms.  In F. Pereira, C.J. Burges, L. Bottou, and K.Q. Weinberger, editors, Advances in Neural Information Processing Systems, volume 25. Curran Associates, Inc., 2012.
  * [SSW+16] Bobak Shahriari, Kevin Swersky, Ziyu Wang, Ryan P. Adams, and Nando de Freitas.  Taking the human out of the loop: A review of bayesian optimization.  Proceedings of the IEEE, 104(1):148–175, 2016.
  * [WSS+23] Colin White, Mahmoud Safari, Rhea Sukthanker, Binxin Ru, Thomas Elsken, Arber Zela, Debadeepta Dey, and Frank Hutter.  Neural architecture search: Insights from 1000 papers, Jan 2023.
  * [ZL17] Barret Zoph and Quoc V. Le.  Neural architecture search with reinforcement learning.  In ICLR International Conference on Learning Representations, 2017\.