# From laminar to chaotic flow via stochastic resonance in viscoelastic
channel flow

Yuke Li<EMAIL_ADDRESS>Department of Physics of Complex Systems,
Weizmann Institute of Science, Rehovot 7610001, Israel  Victor Steinberg
<EMAIL_ADDRESS>Department of Physics of Complex Systems,
Weizmann Institute of Science, Rehovot 7610001, Israel

###### Abstract

Recent research indicates that low-inertia viscoelastic channel flow
experiences supercritical non-normal mode elastic instability from laminar to
sustained chaotic flow due to finite-size perturbations. The challenge of this
study is to elucidate a realization of such a pathway when the intensity of
the elastic wave is too low to amplify velocity fluctuations above the
instability onset. The study identifies two subregions in the transition flow
regime at Weissenberg number $Wi>Wi_{c}$, the instability onset. In the lower
subregion at $Wi_{c}\leq Wi\leq 300$, we discover periodic spikes in the
streamwise velocity time series $u(t)$ that appear in the chaotic power
spectrum as low-frequency, high-intensity peaks resembling stochastic
resonance (SR). In contrast, the spanwise velocity power spectrum, $E_{w}$,
remains flat with low-intensity, noisy, and broad elastic wave peaks. The
spikes significantly distort the probability density function of $u$,
initiating and amplifying random streaks and wall-normal vorticity
fluctuations. The SR appearance is similar to dynamical systems where chaotic
attractor and limit cycle interact with external white noise. This similarity
is confirmed by presenting a phase portrait in two subregions of the
transition regime. In the upper subregion at $Wi>400$ the periodic spikes
disappear and $E_{w}$ becomes chaotic with a large intensity elastic wave
sufficient to self-organize and synchronize the streaks into cycles and to
amplify the wall normal vorticity according to a recently proposed mechanism.

††preprint: APS/123-QED

##  I Introduction

The transition from laminar to turbulent flow in Newtonian or to chaotic flow
in viscoelastic fluids follows two distinct pathways. One pathway, common to
both fluids, is attributed to a linear normal mode instability followed by a
sequence of linear instabilities leading to turbulent or chaotic flow [1].
Linear stability analysis has been successfully applied and validated for
various flow geometries with curved streamlines [1, 2, 3, 4]. It is important
to note that this approach is only valid for Hermitian equations (or
operators) characterized by normal eigenmodes, which are prone to linear
instability. The onset of instability is determined by the critical eigenvalue
of the most unstable normal eigenmode, which grows exponentially. Furthermore,
nonlinear effects stabilize the critical eigenmode at a sufficiently large
amplitude, thereby replacing a laminar flow that is stable at $Wi<Wi_{c}$.

In the second scenario, which occurs in both Newtonian and viscoelastic
parallel shear flows, turbulent or chaotic flows are induced by finite-sized
external perturbations just above the instability threshold and persist under
steady conditions. This type of instability, known as non-normal mode
instability, occurs in a flow described by non-Hermitian equations
characterized by both normal and non-normal modes, as first discussed for
Newtonian parallel shear flows [5, 6, 7]. The study of the linear instability
of Newtonian parallel shear flows uses a linearized version of the Navier-
Stokes equation, called the Orr-Sommerfeld equation, a classic example of non-
Hermitian equations. Historically, this equation has been used extensively in
its Hermitian approximation at Reynolds number $Re\gg 1$, leading to its
linear stability at all $Re$ [1]. More recently, its non-Hermitian form has
been used to study the instability of Newtonian flows. This approach shows
that while stable normal modes decay exponentially, unstable non-normal modes
grow algebraically and can be amplified by several orders of magnitude over
time [5, 6, 7].

A linear normal mode instability is driven by a deterministic instability
mechanism based on the competition between a destabilizing force and a
stabilizing factor, such as dissipation or external fields, especially when
the former predominates [8, 1]. Specifically, in polymer solution flows with
curved streamlines, linear instability arises due to “hoop stress”, an elastic
stress generated by polymers stretched along curved streamlines by the
velocity gradient. The hoop stress generates a destabilizing bulk force
directed toward the center of curvature, which competes with relaxation
dissipation of the elastic stress. The onset of linear elastic instability,
denoted by the Weissenberg critical number $Wi_{c}$, is determined by its
mismatch [2, 3, 4]. Here, $Wi\equiv U\lambda/d$ is the main control parameter
in inertia-less, viscoelastic fluid flow and denotes the ratio of elastic
stress, generated by polymer stretching due to its entropic elasticity to
stress dissipation due to its relaxation, where $U$ is the mean flow velocity,
$d$ is the characteristic vessel size, and $\lambda$ is the longest polymer
relaxation time. In the present experiment, we consider $Wi\gg 1$ and the
Reynolds number $Re\equiv Ud\rho/\eta\ll 1$, which corresponds to an extremely
large elasticity number $El=Wi/Re=\lambda\eta/\rho d^{2}>>1$, where $\rho$ and
$\eta$ are the solution density and dynamic viscosity, respectively (see
Materials and Methods for details).

However, the instability mechanism becomes ineffective at a zero-curvature
limit [2, 9]. Therefore, viscoelastic parallel shear flows have been shown to
be linearly stable over all $Wi$ and for $Re\ll 1$ [10, 11], which is
analogous to Newtonian parallel shear flows at all $Re$. For the latter,
however, the proven linear stability fails to explain Reynolds’ seminal
experiments, where a transition from laminar to turbulent flow was observed in
Newtonian pipe flow at finite $Re$ [1, 6, 12]. Thus, linear stability does not
necessarily imply global stability in either Newtonian [1, 6] or viscoelastic
[13, 14, 15, 16, 17, 12] parallel shear flows. Similar to Newtonian flows,
inertia-less viscoelastic parallel shear flows exhibit non-normal mode
instability at $Wi\gg 1$ and $Re\ll 1$ due to the non-Hermitian nature of the
linearized elastic stress equation [13, 14, 15, 16, 17, 12]. These papers
focus on theoretical investigations and numerical analyses of the non-normal
unstable modes in viscoelastic channel and plane Couette flows only during the
linear transient growth at $Wi\gg 1$ and $Re\ll 1$. During the transient
growth, the emergence of coherent structures (CSs), such as streamwise streaks
and rolls, is observed. However, the linear growth model fails to predict the
presence of nonlinearly stabilized sustained CSs.

The first experiments on viscoelastic parallel shear flows are carried out in
pipes [18, 19, 20], then extended to straight square microchannels with strong
pre-arranged perturbations at the inlet [21, 22, 23]. These experiments report
the finding of significant velocity fluctuations above $Wi_{c}$, which
contradicts the theoretically proven linear stability of such viscoelastic
flows [10, 11]. Subsequent experiments performed in our laboratory on
viscoelastic straight channel flows at $Wi\gg 1$ and $Re\ll 1$ with different
external perturbation intensities [24, 25, 26, 27, 28, 29], find a
supercritical instability at $Wi\geq Wi_{c}$ leading directly to a sustained
chaotic flow. This transition is characterized by the dependence of the
friction factor, the normalized root mean square (rms) streamwise velocity and
pressure fluctuations on $Wi-Wi_{c}$, with scaling exponents that deviate from
the 0.5 typical for normal mode supercritical bifurcation [1, 8]. Moreover,
$Wi_{c}$ depends on the amplitude and structure of the external perturbations,
in contrast to the normal mode bifurcation [1]. These two key features - the
direct transition from laminar to chaotic flow and the dependence of $Wi_{c}$
on the intensity and spectral characteristics of the external perturbations -
validate the non-modal elastic instability [14, 12, 30]. At $Wi>Wi_{c}$,
regardless of the strength and structure of the external perturbations, three
chaotic flow regimes emerge: transition, elastic turbulence (ET), and drag
reduction (DR), each accompanied by elastic waves just above $Wi_{c}$. The
experimental results show a universal scaling of the flow properties, a
dependence of the elastic wave propagation velocity on $Wi-Wi_{c}$, a power-
law decay of the velocity power spectra, and a predominant presence of streaks
[24, 25, 26, 27, 28, 29].

The existence of elastic waves in three flow regimes is the characteristic
feature of viscoelastic chaotic channel flow. The waves are predicted to
appear in the turbulent flow of viscoelastic fluid at $Re\gg 1$ and $Wi\approx
1$, and are expected to be suppressed by large viscous dissipation at $Re\ll
1$ and $Wi\gg 1$ [31]. Contrary to these predictions, the elastic waves were
first discovered in the inertia-less viscoelastic flow at $Re\ll 1$ and $Wi\gg
1$ in the ET flow regime. This discovery is characterized by distinct peaks in
the spanwise velocity power spectra, with frequency and intensity varying with
$Wi$, and in particular by their propagation velocity in the direction of the
streamwise velocity, which is characterized by a scaling law as a function of
$Wi-Wi_{c}$, as reported in Ref. [32]. Subsequent research has confirmed this
phenomenon and established its universality in channel flows with different
external noise intensities [4, 27, 29, 28, 24, 25, 26]. A simple physical
explanation for the appearance of elastic waves is based on the analogy of the
response of an elastic stress field to transverse perturbations, similar to
the response of an elastic string when plucked [4].

Here, the experimental study focuses on the immediate vicinity above the onset
of the supercritical elastic instability of inertia-less viscoelastic planar
channel flow subjected to finite-size perturbations induced by the unsmoothed
inlet and two small cavities in the upper wall at both channel ends. For
$Wi\geq Wi_{c}$, the instability leads to a transition from laminar to
sustained chaotic flow. Recent experimental evidence reveals the key role of
elastic waves in amplifying wall normal vorticity fluctuations through
resonant pumping of elastic wave energy withdrawn from the mean shear flow
into the fluctuating wall normal vorticity [28]. The higher(lower) the elastic
wave intensity, the higher(lower) the flow resistance and the rotational
vorticity fluctuations [28]. This mechanism also helps to explain the
transition from ET to DR first described in Ref. [33] and explained in Ref.
[34]. Given the key role of elastic waves in promoting sustained chaotic flow,
the main goal of the current experiment is to uncover the pathway to chaotic
flow, in particular when the energy of the elastic wave is insufficient to
initiate and amplify streaks and wall normal vorticity fluctuations above
$Wi_{c}$.

##  II Results

We performed the experiment in a long channel of dimensions 500 ($L$, $x$:
streamwise) $\times$ 0.5 ($h$, $y$: wall-normal) $\times$ 3.5 ($W$, $z$:
spanwise) mm3,shown in Fig. 1. The channel has an unsmoothed inlet and two
small holes of 0.5 mm diameter in the top plate near the inlet and outlet,
used for pressure measurements, which introduce finite-size perturbations,
much weaker than strong pre-arranged perturbations used in Refs. [24, 25] and
slightly weaker than in the channel discussed in Ref. [29]. The latter is
determined by the rms of streamwise velocity fluctuations, $u_{rms}$, at the
inlet. A dilute aqueous solution with 64% sucrose concentration and 80 ppm
high molecular weight polyacrylamide (PAAm) is the same as used in Refs. [29,
28]. Its longest relaxation time ($\lambda$) is 13 s. In the Materials and
Methods section, we discuss details of the experimental setup, preparation and
characterization of the polymer solution, pressure ($p$) fluctuations, and
Particle Image Velocimetry (PIV) measurement techniques.

Figure 1: Schematic of the experimental setup. Flow in the long straight
channel with the non-smoothed inlet and two small cavities in the top plate
close to the inlet and outlet, used for absolute pressure measurement, is
driven by compressed Nitrogen gas (from left to right in the schematic).
Pressure, velocity, and flow discharge are measured by absolute pressure
sensors, a high-speed and high spatial resolution camera with laser sheet
illumination, and a balance interfaced with a computer, correspondingly.

###  II.1 Characterization of the elastic transition at $Wi\geq Wi_{c}$

The onset of the elastic instability at $Wi_{c}=150\pm 10$ and the
characterization of the channel flow at $Wi\geq Wi_{c}$ in three subsequent
chaotic flow regimes, namely transition, ET, and DR, are described in Ref.
[28]. These three flow regimes are well characterized by the dependence of the
friction factor, normalized pressure, and streamwise velocity fluctuations on
$Wi/Wi_{c}-1$ at the downstream location $l/h=380$, as reported in the same
reference [28]. The corresponding exponents are found to be universal for the
studied viscoelastic channel flows with different external perturbations [24,
26, 29]. Examining $u_{rms}$ and rms pressure fluctuations, $p_{rms}$, versus
$Wi$ below and above $Wi_{c}$, one finds their continuous growth without any
clear indication of the instability onset as shown in Figs. 2(a) and 2(c).
Algebraic fits of the normalized streamwise velocity fluctuations,
$u_{rms}/u_{rms,lam}$, and pressure fluctuations, $p_{rms}/p_{rms,lam}$,
versus $Wi/Wi_{c}-1$ up to $Wi=280$ yield exponent values of 0.85 and 1.10,
respectively. These values are consistent with those obtained in the previous
studies, where the exponents were obtained for data taken in the entire
transition regime as shown in Figs. 2(b) and 2(d) [24, 26, 29, 27]. Here
$u_{rms,lam}$ and $p_{rms,lam}$ are rms velocity and pressure fluctuations in
laminar flow, respectively.

Figure 2: Characterization of elastic instability and transition regime as a
function of $Wi$ in a limited range above $Wi_{c}$. (a) rms stream-wise
velocity fluctuations, $u_{rms}$, (b) normalized rms stream-wise velocity
fluctuations,$u_{rms}/u_{rms,lam}$, (c) rms pressure fluctuations, $p_{rms}$,
(d) normalized rms pressure fluctuations, $p_{rms}/p_{rms,lam}$, versus $Wi$
in linear coordinates. The black dashed line is the power-law fit of
$u_{rms}/u_{rms,lam}$ and $p_{rms}/p_{rms,lam}$ versus
$(Wi/Wi_{c}-1)^{\epsilon}$+1 with $\epsilon=0.85$ and 1.10, respectively,
above and in close vicinity to $Wi_{c}$.

###  II.2 Temporal evolution of stream-wise velocity fluctuations and low-
frequency periodic spikes

By monitoring the temporal evolution of the velocity at the channel center at
$l/h=380$ by PIV, followed by a slight pressure increase at the inlet at $t=0$
s, the streamwise velocity $u(t)$ in laminar flow at $Wi=139$ starts to grow
algebraically up to $t\sim 100$ s, as shown in Fig. 3(a). This growth pattern
is in stark contrast to the exponential growth in a case of normal mode
instability [1]. At $t>100$ s, the growth rate of $u(t)$ decreases sharply,
with $u(t)$ reaching a saturation point after $t\approx 150$ s at
$Wi=169>Wi_{c}$. During the transient growth and subsequent saturation, the
streamwise velocity fluctuations, $u^{\prime}=u-u_{mean}$, show a slight
increase, in agreement with the observations in Figs. 4(a) and 4(b). In
particular, at $t\approx 100$ s, large and almost periodic negative spikes
appear and persist until the end of time series (Fig. 3(a)). Here $u$ and
$u_{mean}$ are the streamwise velocity and the mean streamwise velocity
averaged over a narrow spanwise extent around the channel centerline,
respectively.

To enhance the visibility of these spikes, $u(t)$ is partially averaged over
each 1 s interval (100 measurement samples), and the resulting partially
smoothed data is plotted as $u^{p}_{rms}(t)$ in Fig. 3(b). It can be seen that
$u^{p}_{rms}(t)$ shows only a marginal increase at $t>0$ s, with much more
pronounced periodic spikes at $t\approx 100$ s. In Figs. 3(c) and 3(d) of the
streamwise velocity power spectra, $E_{u}$, and pressure power spectra,
$E_{p}$, presented in lin-log coordinates, show narrow spike peaks in the low
frequency range for three values of $Wi$ up to 280. However, at $Wi>300$ the
energy of these spikes is barely detectable and eventually disappears. Thus,
the spikes exist only in the lower subrange of the transition flow regime from
$Wi_{c}$ up to $Wi\approx 300$. Finally, the dependence of the normalized
spike frequency, $\lambda f_{p}$, obtained from $E_{u}$, as a function of $Wi$
is shown in Fig. 3(e) along with a fitted relation $\lambda
f_{p}\sim(Wi/Wi_{c}-1)^{0.45}$.

Figure 3: Temporal evolution of the streamwise velocity $u(t)$, above
$Wi_{c}=150\pm 10$ and low-frequency periodic spikes. (a) Streamwise velocity
$u(t)$, and (b) its corresponding partially averaged velocity fluctuations
$u^{p}_{rms}(t)$ over each 1s at 100 fps at $l/h=380$ due to a pressure
increase at $t=0$ s. The flow varies from $Wi=139$ to $Wi=169$. (c) Streamwise
velocity, (d) pressure power spectra at low normalized frequencies, $\lambda
f$, in lin-log coordinates at three $Wi$ values. Note the sharp peaks of
velocity and pressure fluctuations. (e) $Wi$ dependence of the normalized
periodic spikes frequency, $\lambda f_{p}$, above $Wi_{c}$ to about the edge
of their existence. Figure 4: Time evolution of velocity fluctuations
witnessing a sudden change of inlet pressure. (a) $u(t)$ after pressure
increase at the inlet at $t=0$, resulting in $Wi$ growth from 140 to 204. The
grey line presents the raw data taken with 100 fps of stream-wise velocity as
a function of $t$, $u(t)$, the black presents $u(t)$ averaged over 1 s, and
the red dashed line is the fit $A\times t^{\epsilon}+u_{0}$, with fitting
parameters $A=0.0384$, $\epsilon$=0.0829, and $u_{0}$=6.52 mm/s. (b) Seven
cases of $u(t)$ include two scenarios: from laminar to laminar ($Wi$:
131$\rightarrow$116 and 132$\rightarrow$138) and from laminar to transition
flow regime at $Wi>Wi_{c}$ ($Wi$:139$\rightarrow$169, 148$\rightarrow$180,
143$\rightarrow$198, 140$\rightarrow$204, and 140$\rightarrow$224). (c)
streamwise velocity at a fixed location as a function of time, $u(t)$, during
50 s at different $Wi=169,180,198,204,224$.

###  II.3 Statistical properties of velocity fluctuations

Figure 5(a) shows the probability density functions (PDFs) of normalized
streamwise velocity fluctuations, $P((u-u_{mean})/u_{rms})$, over a range of
$Wi$ from 52 to 280 below and above $Wi_{c}$. In particular, the first
deviations from Gaussian PDFs appear for negative values of
$(u-u_{mean})/u_{rms}$ at $Wi=160$, just above $Wi_{c}$, with these deviations
becoming more pronounced up to $Wi=280$. This is highlighted in Fig. 5(c) by
showing the PDFs for $Wi\geq Wi_{c}$, with black arrows indicating the primary
changes from the Gaussian PDF as $Wi$ increases. As can be seen, the PDFs of
the positive fluctuations decrease, while the PDFs of the negative
fluctuations increase significantly. This is due to the spikes at lower $u(t)$
values, which are smaller than $u_{mean}$ at $Wi>Wi_{c}$ (Figs. 3(a) and 4).
The inset in Fig. 5(d) compares PDFs at $Wi=52<Wi_{c}$ and $Wi=253>Wi_{c}$,
where the latter is fitted by two Gaussian PDFs to emphasize the distortion in
the PDF. The main plot in Fig. 5(d) illustrates the abrupt changes in skewness
($S$) and flatness (kurtosis) ($F$) corresponding to the third and fourth
moments of $P((u-u_{mean})/u_{rms})$ at $Wi>Wi_{c}$. Remarkably, their
deviations from zero precisely determine the instability onset at
$Wi_{c}=150$. In contrast, the PDFs of the normalized spanwise velocity
fluctuations, $P(w/w_{rms})$, show only small deviations from the Gaussian
distribution within the error bars in Fig. 5(b) for the same $Wi$ range (Figs.
6(c) and 6(d)). Moreover, their $S$ and $F$ remain zero both below and above
$Wi_{c}$ (Fig. 5(d)).

Figure 5: PDFs of normalized stream-wise $P((u-u_{mean})/u_{rms})$ and span-
wise $P(w/w_{rms})$ velocity fluctuations at the channel center and various
$Wi$ values below and above $Wi_{c}$. (a) $P((u-u_{mean})/u_{rms})$ and (b)
$P(w/w_{rms})$ in the whole lower sub-range of $Wi$, (c)
$P((u-u_{mean})/u_{rms})$ only above $Wi_{c}$, and (d) skewness, $S$, and
flatness (kurtosis), $F$, of $P((u-u_{mean})/u_{rms})$ and $P(w/w_{rms})$,
respectively. Inset: $P((u-u_{mean})/u_{rms})$ for $Wi=52<Wi_{c}$ and
$Wi=253>Wi_{c}$. The dashed black line is a Gaussian fit to the data at
$Wi<Wi_{c}$, and dashed green lines are Gaussian partial fits to the data at
$Wi>Wi_{c}$ with extensions to negative and positive values of
$(u-u_{mean})/u_{rms}$, respectively.  Figure 6: PDFs of spanwise and
streamwise velocities at channel center. PDFs of $(u-u_{mean}/u_{rms})$ at (a)
$Wi<Wi_{c}$, and (b) $Wi>Wi_{c}$. PDFs of $w/w_{rms}$ at (c) $Wi<Wi_{c}$, and
(d) $Wi>Wi_{c}$.

###  II.4 Streamwise and spanwise velocity power spectra and elastic waves

At $Wi<Wi_{c}$, the streamwise velocity power spectrum, $E_{u}$, shows a
slight increase from nearly flat to slight growth toward lower normalized
frequencies, $\lambda f$, reminiscent of a white noise spectrum (Fig. 7(a)).
However, above $Wi_{c}$ up to $Wi=280$, $E_{u}$ shows a significant increase,
up to three orders of magnitude at $\lambda f\leq 10$ (Figs. 7(a) and 7(c)).
This increase is accompanied by a sharp power-law decay at $\lambda f\geq 10$,
with the exponent $\alpha$ ranging from -0.8 to -2.5 (see Fig. 7(d)). Here
$\alpha$ is denoted as the exponent of the power-law fit of $E_{u}$ to
$\lambda\cdot{f}$ in the decay regime.

In the lower subrange of $Wi>Wi_{c}$ up to $Wi=280$, the spanwise velocity
power spectra, $E_{w}$, plotted on a linear-logarithmic scale in Fig. 7(b) and
on a log-log scale in Fig. 8, remain flat and similar to those observed at
$Wi<Wi_{c}$, which characterizes white noise. However, at $Wi>Wi_{c}$,
pronounced broad and noisy peaks associated with the elastic waves appear at
the top of the flat power spectra (Fig. 7(b)) and are visible in Fig. 8.

Since the identification of the peak maximum location and the estimation of
its height become unambiguous near $Wi_{c}$ due to a drastic increase in peak
widths and fluctuations as $Wi$ approaches $Wi_{c}$, the resolution of both
$\lambda f_{el}$ and the normalized intensity of the elastic waves,
$I_{el}/w_{rms}^{2}$, is limited to about $\sim 5$ and $\approx 2\times
10^{-6}$, respectively (Figs. 7(e) and 7(e)(f)). Isolated single peaks in
$I_{el}/w_{rms}^{2}$ versus $\lambda f$ coordinates for $Wi>Wi_{c}$ provide
better resolution (Fig. 9). In addition, Fig. 7(f) includes a power-law fit of
the normalized elastic wave energy dependence on $Wi$ as
$I_{el}/w_{rms}^{2}\sim(Wi/Wi_{c}-1)^{\beta}$ with $\beta=0.20$. Remarkably,
in the upper subrange of the transition regime at $400\leq Wi<\approx 1000$,
and further in the ET and DR regimes, the elastic wave peaks in $E_{w}$ grow
by up to four orders of magnitude, as shown in Fig. 8. The location of the
peak maximum shifts more than tenfold to higher values of $\lambda f$, and a
power-law decay occurs with the scaling exponent increasing with $Wi$. This
observation reveals a clear distinction between two subranges in the
transition flow regime.

Figure 7: Energy spectra and properties of elastic waves. (a) $E_{u}$ and (b)
$E_{w}$ vs. $\lambda f$ at different $Wi$ from 52 to 280 in log-log and lin-
log scales, respectively. (c) $E_{u}$ vs. $\lambda f$ at $Wi>Wi_{c}$ with
power-law decays fitted by dashed lines with corresponding values of the decay
exponents $\alpha$ and colors. (d) Dependence of $\alpha$ on $Wi$. (e)
Dependence of the normalized elastic wave frequency, $\lambda f_{el}$, on
$Wi>Wi_{c}$. (f) Dependence of the normalized elastic wave energy,
$I_{el}/w_{rms}^{2}$, on $Wi$ at $Wi>Wi_{c}$. The dashed black line is the
fit, which gives a scaling relation:
$I_{el}/w_{rms}^{2}\sim(Wi/Wi_{c}-1)^{0.2}$. In both (e) and (f) plots the
last two points near the onset reach the experimental resolution limit.
Figure 8: Energy spectra of spanwise velocity, $E_{w}$, versus normalized
frequency, $\lambda\cdot f$ in log-log scales at $Wi<Wi_{c}$ in laminar flow
and up to $Wi=3051$ in transition, ET and DR flow regimes. It is clearly
notable that in the lower sub-range of the transition flow regime elastic wave
peaks are not noticeable and an order of magnitude smaller than in the upper
sub-range at $Wi=415$ in log-log coordinates.  Figure 9: Individual looks of
elastic waves. $I_{el}/w_{rms}^{2}$ at $Wi$=146, 160, 172, 190, 195, 253, and
280. Each dashed circle with the corresponding color indicates the appearance
of one significant peak.

###  II.5 Observation of random streaks

One of the striking findings in inertia-less viscoelastic channel flows,
especially when subjected to external perturbations of different amplitudes at
$Wi>Wi_{c}$, is the phenomenon of cyclic self-organized streaks synchronized
by elastic waves [29]. To verify the streak synchronization by the elastic
waves we explore the approach based on the velocity difference across the
counter-propagating streak interface as a function of normalized time,
$t^{*}=tf_{el}$, developed in our laboratory and used in several experiments
[24, 25, 26, 29, 27]. In all these experiments, at $Wi>Wi_{c}$, we discovered
cycling with $t^{*}\simeq 1$ in all three flow regimes.

On the contrary, in the lower subrange at $Wi>Wi_{c}$ up to $Wi=280$, we
reveal the emergence of random streaks shown in Fig. 10 in two rows of 5
consecutive images each, illustrating the streaks as a function of $t^{*}$ at
$Wi=193$ and 207. At $Wi\leq 280$ only random streaks appear, although the
number of streaks increases significantly with $Wi$ (Fig. 11). These streaks
lack the organized, cyclic nature observed at higher $Wi$ values. At $Wi\geq
400$, the cycling of self-organized streaks synchronized by the elastic waves,
whose energy increases by orders of magnitude, is in agreement with our
previous experiments [24, 25, 26, 29, 27, 28].

Figure 10: Instantaneous images of streamwise velocity fluctuations obtained
from the full streamwise velocity $u$ by first subtracting the mean streamwise
velocity profile, $u_{mean}(z)$, and then divided by it. Thus, the normalized
streamwise velocity fluctuations, $u^{\prime}(x,z)/u_{mean}(z)$, at $Wi=193$
(upper 5 images) and $Wi$=207 (lower 5 images) for 5 normalized times
$t*=tf_{el}$ are presented for each $Wi$. PIV measurements present images in
$\Delta z/W=0.67$ and $\Delta l/h=4$ spatial window.  Figure 11: Number of
streaks observed within 70 seconds as a function of $Wi$. The appearance of
streaks is identified by a contour plot of streamwise velocity fluctuations
with high resolution via PIV. The color bar of streaks is presented for
$Wi>Wi_{c}$ in the lower sub-range of transition flow regime at $Wi=193$ and
207.

##  III Discussion

Figure 12: Phase portraits of the instantaneous dynamics with increasing $Wi$.
(a) $(u^{p}_{rms})^{2}$ versus $(w^{p}_{rms})^{2}$ at six $Wi$ values
representing three states: (i) laminar basin ($Wi$=93); (ii) chaotic state of
$(u^{p}_{rms})^{2}$ along with periodic orbits and white noise of
$(w^{p}_{rms})^{2}$ at $Wi$=159, 190, and 280; and (iii) the pure chaotic
state of both $(u^{p}_{rms})^{2}$ and $(w^{p}_{rms})^{2}$ at $Wi$=415 and 750.
$(u^{p}_{rms})^{2}$ and $(w^{p}_{rms})^{2}$ are computed from squared velocity
fluctuations averaged over 1 second (100 points). Individual views of state
(ii) are shown separately at (b) $Wi$=280, (c) $Wi$=190, and (d) $Wi$=159.

This experimental study attempts to elucidate the mechanism behind the
supercritical, non-normal mode elastic instability in viscoelastic channel
flow with negligible inertia initiated by external finite-size perturbations
due to an unsmoothed inlet and two small holes, in the narrow range from
$Wi\sim 60$ to $Wi\sim 300$ close to $Wi_{c}=150$.

The first finding of this study is the existence of two subranges in the
transition flow regime above $Wi_{c}$ with different flow properties, as
detailed in the Results section. However, the most surprising discovery is the
observation of low-frequency periodic spikes in $u(t)$ and $p(t)$ detected
exclusively in the lower subrange at $Wi_{c}\leq Wi\leq 300$ (Figs. 3 and 4).
The periodic spikes appear as sharp low-frequency peaks in the streamwise
velocity, $E_{u}$, and pressure, $E_{p}$, power spectra presented in lin-log
coordinates in Figs. 3(c) and 3(d), and only $E_{u}$ in log-log coordinates in
Figs. 7(a) and 7(c). In the latter, the low-frequency $\lambda f<1$ sharp
peaks appear on top of the streamwise velocity fluctuation power spectrum
characterized by a power-law decay at $\lambda f>1$, indicating chaotic flow
in the streamwise direction. In contrast, the power spectrum of the spanwise
velocity fluctuations, $E_{w}$, slowly increases at $Wi>Wi_{c}$ in the lower
subrange up to $Wi\approx 300$, where $E_{w}$ still remains flat, similar to
that observed at $Wi<Wi_{c}$, indicating a white noise spectrum in the
spanwise direction (Figs. 7(b) and 8). Thus, in the lower subrange, the
chaotic streamwise velocity fluctuations coexist and interact with the
spanwise white noise velocity fluctuations in the presence of low-intensity
elastic waves. Furthermore, the periodic spikes presented and characterized in
Results are reminiscent of stochastic resonance (SR), a phenomenon considered
in autonomous dynamical chaotic systems interacting with external white noise
in the presence of weak periodic modulation. In dynamical chaotic systems, the
strange attractors found in practice are structurally unstable. This means
that such a quasi-attractor exists as several regular chaotic attractors
coexisting in the phase space of the system, which merge into one chaotic set
at the critical value of the control parameter, the so-called attractor crisis
[35]. Above the crisis, the merged attractor shows intermittent switching of
the “chaos-chaos” type, where the ”deterministic stochastic attractor” is
observed [36, 35]. The dynamical intermittent switching is a result of
additive noise depending on its intensity, which can be synchronized by weak
periodic modulation resulting in SR [36, 35]. However, the deterministic SR is
fundamentally different from the classical SR realized in a bistable system
driven simultaneously by noise and a weak periodic signal [37]. To further
substantiate the connection of the observations with deterministic SR, we
characterize the flow dynamics by plotting the phase portrait in coordinates
of the streamwise, $(u^{p}_{rms})^{2}$, versus spanwise, $(w^{p}_{rms})^{2}$,
velocity energies in the two subranges of the transition flow regime (Fig.
12(a)). In the lower subrange above $Wi_{c}$ we find a one-dimensional
streamwise velocity chaotic attractor interacting with spanwise velocity white
noise and weak elastic waves, resulting in a SR periodic orbit at $Wi=159$,
190, and 280 (Figs. 12(b)-12(d)). We note that in the lower subrange, the
energy of the external white noise perturbations $(w^{p}_{rms})^{2}$ at
$Wi\leq Wi_{c}$ increases slightly up to about 50%, while the chaotic
streamwise velocity energy, $(u^{p}_{rms})^{2}$, increases about tenfold. In
contrast, in the upper subrange we detect a two-dimensional chaotic attractor
of both $(u^{p}_{rms})^{2}$ and $(w^{p}_{rms})^{2}$ that grows by more than
two orders of magnitude compared to external perturbations in a laminar flow
(Fig. 12(a)).

Moreover, in the lower subrange above $Wi_{c}$, the presence of periodic
spikes leads to significant deviations from Gaussian distributions in the PDFs
of $(u-u_{mean})/u_{rms}$ at the channel center (Figs. 5(a), 5(c), 5(d), and
6). Other notable features in the lower sub-range include the abrupt change
from zero to high absolute values of the PDFs’ skewness, $S$, and flatness
(kurtosis), $F$, from $Wi_{c}=150$ to $Wi<300$; the appearance of random
streaks (Fig. 11); and the observation of wall-normal vortex fluctuations
already in the lower subrange (Fig. 13). However, the synchronization at
$f_{el}$ of the self-organized cycling streaks and the enhancement of the
wall-normal vortex fluctuations occur only in the upper subrange, at $Wi\geq
400$, where the elastic wave energy, $I_{el}$, increases by orders of
magnitude compared to the lower subrange (Figs. 7(b) and 8).

Figure 13: An instantaneous snapshot of the 2D image of streamlines of
velocity fluctuations field $(u^{\prime},w^{\prime})$ in $x$–$z$ plane at the
channel center overlapped with the field of the absolute value of the vortical
part of the wall-normal vortex fluctuations, $\lvert\Omega_{R}\rvert$, at
$Wi$=210. The fluid is seeded with red-fluorescent 0.2 $\mu$m traces of
$\sim$0.15% w/w concentration (Thermo Scientific). The resolution is
512$\times$512 pxl2 at 4000 fps. The PIV experimental setup is the same as
that in Fig. 2A in Ref. [28]. The color bar of the vortical part of vorticity
fluctuations is shown on the right.

It is also unexpected that the wall-normal vortex fluctuations are likely
generated and supported by strong periodic spikes in the lower subrange, while
the elastic waves take over their amplification only in the upper subrange at
$Wi\geq 400$ [28]. Remarkably, at $Wi>300$ the low-frequency periodic spikes
disappear, the strongly distorted negative tail of the PDF of
$(u-u_{mean})/u_{rms}$ becomes exponential, and the flow becomes chaotic in
both the $u$ and $w$ velocity fields. The flow dynamics at $Wi\geq 400$ are
characterized by a two-dimensional chaotic attractor, and the scaling of the
flow properties is found to be the same as in the transition flow regime of
viscoelastic channel flows with different perturbation intensities, where the
measurements are made with lower resolution and further away from $Wi_{c}$
[29, 28].

In conclusion, the presented experimental results reveal novel and unexpected
features of the non-modal instability evolution above and near $Wi_{c}$ and
provide strong evidence for the stochastic nature of the instability
mechanism. Above $Wi_{c}$, the low-frequency periodic spikes observed in the
streamwise velocity time series, characterized by the sharp peaks in $E_{u}$
in the lower subrange at $Wi_{c}\leq Wi\leq 300$ of the transitional flow
regime, are reminiscent of SR. Based on the experimental phase portrait of the
instantaneous dynamics, it is proposed that SR arises due to the interaction
of streamwise velocity chaotic flow with spanwise velocity white noise in the
presence of low-intensity elastic waves. This scenario is similar to the
deterministic SR found in autonomous chaotic nonlinear dynamical systems,
where deterministic SR is generated due to the interaction of a chaotic
attractor with external white noise in the presence of a weak periodic signal
at the fixed value of the control parameter [36, 35]. SR greatly increases the
probability of slow streamwise velocity fluctuations and promotes the
generation of wall-normal vorticity fluctuations. Moreover, since the elastic
waves have low intensity at $Wi_{c}\leq Wi\leq 300$, which is insufficient to
initiate and synchronize the streaks and amplify the wall-normal vorticity
fluctuations, SR takes over the role of the elastic waves in initiating the
random streaks and amplifying the wall-normal vorticity fluctuations.

The discovery of SR in the limited range of $Wi_{c}\leq Wi\leq 300$ is crucial
for understanding the pathway above supercritical non-modal elastic
instability to sustained chaotic flow. This phenomenon may have broader
implications for various flows, including those in stress fields, such as
magneto-hydrodynamic [38], viscoelastic electrolyte [39, 40], viscous solution
of rods [41] and suspension of flexible fibers [42], active nematics [43] and
bacterial and active turbulence [44] in parallel shear flows.

Finally, we would like to point out similarities and differences with the
widely studied theoretically, numerically, and experimentally Newtonian pipe
flow at high $Re$, especially in the last two decades [45, 46, 47, 7, 48, 49,
50, 51, 52, 53, 54]. While Newtonian pipe flows have long been shown to be
linearly stable, their global stability is not guaranteed. In fact, despite
their proven linear stability, Newtonian parallel shear flows become unstable
to finite-size external perturbations at finite $Re$ [1]. Only about four
decades ago it was realized that the discrepancy between the proven linear
stability and the experimental observations is explained by the non-Hermitian
Orr-Sommerfeld equation, the linearized Navier-Stokes equation for Newtonian
parallel shear flows [1]. At finite $Re$, this equation is prone to transient
algebraic time growth of non-normal modes, whereas linear eigenmodes decay
exponentially, since the flow is linearly stable [1, 5, 6, 7]. The stochastic
nature of the non-normal mode instability, recognized since Reynolds’ seminal
experiments, has only recently been experimentally demonstrated in
viscoelastic parallel shear flows [55, 48, 53]. The dependence of $Re_{c}$ on
the external perturbation intensity contradicts the normal mode linear
instability onset independent of perturbation intensity [1]. However, the
pathway from laminar flow to sustained turbulence via subcritical non-modal
instability in Newtonian pipe flow at $Re\gg 1$ by convective and absolute
instabilities is significantly different and more complicated [56, 46, 48, 50,
51, 52, 53, 57], than the pathway from laminar to sustained chaotic flow via
supercritical non-normal-mode elastic absolute instability in viscoelastic
channel flows at $Wi\gg 1$, $Re\ll 1$ and $El\gg 1$ [24, 26, 29, 28, 27]. As
suggested in Ref. [29], the surprising difference between the flows is due to
the presence of only one nonlinear advection term in the Navier-Stokes
equation, as opposed to three nonlinear terms in the elastic stress equation,
of which only one is an advection term.

###### Acknowledgements.

We thank Guy Han and Rostyslav Baron for their help with the experimental
setup. This work was supported in part by the Israel Science Foundation (ISF,
grant #784/19).

## Materials and Methods

###  .1 Experimental setup and flow discharge measurements

The experiments are conducted in a straight channel of 500(L) $\times$
3.5($W$) $\times$ 0.5(h) mm3 dimensions, shown in Fig. 1. The fluid is driven
by N2 gas at pressure up to 10psi. The fluid discharge is weighed
instantaneously, $m(t)$, by a PC-interfaced balance (BPS-1000-C2-V2, MRC) to
measure the time-averaged fluid discharge rate $Q=\langle\Delta m/\Delta
t\rangle$ to get the mean velocity $U=Q/\rho Wh$. Then $Wi=\lambda U/h$ and
$Re=\rho Uh/\eta$ vary in the ranges (30, 300) and (0.005, 0.045),
respectively. High resolution ($0.1\%$ of full scale) absolute pressure
sensors located near the inlet and outlet at $\Delta l/h\approx 10$ (HSC
series, Honeywell) with range up to 5 psi are used to measure the pressure
fluctuations. For Fig. 3 and Fig. 4, the container of the solution is
immediately switched to the other $N_{2}$ supply line with a slightly higher
or lower pressure at $t=0$ and just after switching, $Wi$ is calculated the
flow discharge when the flow system reaches equilibrium.

###  .2 Polymer solution preparation and characterization

As the working fluid, a dilute polymer solution of high molecular weight
polyacrylamide (Polysciences, $M_{w}=18MDa$) at a concentration of $c=80$ ppm
($c/c*\approx 0.4$ with the overlapping polymer concentration of $c*\approx
200$ ppm [58]) is prepared using a water-sucrose solvent with a sugar weight
fraction of $64\%$. The properties of the solution are as follows: the
solution density ($\rho$) is 1320 kg/$m^{3}$, the solvent viscosity
($\eta_{s}$) is 0.13 Pa$\cdot$s, and the total solution viscosity ($\eta$) is
0.17 Pa$\cdot$s. The ratio of solvent viscosity to total solution viscosity is
$\eta_{s}/(\eta_{s}+\eta_{p})=0.765$, where $\eta_{p}$ is the polymer
contribution to the solution viscosity. The longest polymer relaxation time
($\lambda$) is 13 seconds, obtained by the stress relaxation method [58]. The
result is $El=Wi/Re=\lambda\rho/\eta d^{2}=6.8\times 10^{3}$.

###  .3 Imaging system and PIV measurements

We perform velocity field measurements at various distances $l/h$ downstream
of the inlet using the Particle Image Velocimetry (PIV) method. The PIV setup
consists of 3.2$\mu$m latex fluorescent particle tracers of $\sim$1% w/w
concentration (Thermo Scientific) illuminated by a laser sheet of $\approx
50\mu$m thickness over the central channel plane, i.e. the $x-z$ plane. A high
speed camera (Mini WX100 FASTCAM, Photron) has a high spatial resolution and
images of tracer pairs are acquired at 200 to 2000 fps. The OpenPIV software
[59] is used to analyze $u(x,z,t)$ and $w(x,z,t)$ in the 2D $x$-$z$ plane to
record data for $\sim\mathcal{O}(15)$ minutes, or
$\sim\mathcal{O}(50\lambda)$, for each $Wi$ to obtain sufficient statistics.
For the velocity fluctuations in Figs. 3, 7, and 4, we use spatial averaging
over narrow windows with a resolution of 256 ($x$) $\times$96 ($z$) pxl2 with
a 4x objective, which serves as a single point velocity measurement at the
channel center. The window size for PIV is $32\times 32$ pxl2 with 50% overlap
and 200% search window size. For the velocity profile visualization in Fig.
10, we keep the same scale but increase the spatial resolution to
1280$\times$1280 pxl2.

## References

  * Drazin and Reid [2004] P. G. Drazin and W. H. Reid,  _Hydrodynamic stability_ (Cambridge university press, 2004).
  * Larson  _et al._ [1990] R. G. Larson, E. S. Shaqfeh, and S. J. Muller, A purely elastic instability in taylor–couette flow, J. Fluid Mech. 218, 573 (1990).
  * Shaqfeh [1996] E. S. Shaqfeh, Purely elastic instabilities in viscometric flows, Annu. Rev. Fluid Mech. 28, 129 (1996).
  * Steinberg [2021] V. Steinberg, Elastic turbulence: an experimental view on inertialess random flow, Annu. Rev. Fluid Mech. 53, 27 (2021).
  * Schmid [2007] P. J. Schmid, Nonmodal stability theory, Annu. Rev. Fluid Mech. 39, 129 (2007).
  * Trefethen  _et al._ [1993] L. N. Trefethen, A. E. Trefethen, S. C. Reddy, and T. A. Driscoll, Hydrodynamic stability without eigenvalues, Science 261, 578 (1993).
  * Kerswell [2018] R. Kerswell, Nonlinear nonmodal stability theory, Annu. Rev. Fluid Mech. 50, 319 (2018).
  * Cross [1993] P. C. Cross, M. C. & Hohenberg, Pattern formation outside of equilibrium, Rev. Mod. Phys. 65, 851 (1993).
  * Pakdel and McKinley [1996] P. Pakdel and G. H. McKinley, Elastic instability and curved streamlines, Phys. Rev. Lett. 77, 2459 (1996).
  * Gorodtsov and Leonov [1967] V. Gorodtsov and A. Leonov, On a linear instability of a plane parallel couette flow of viscoelastic fluid, J. Appl. Math. Mech. 31, 310 (1967).
  * Renardy and Renardy [1986] M. Renardy and Y. Renardy, Linear stability of plane couette flow of an upper convected maxwell fluid, J. Nonnewton. Fluid Mech. 22, 23 (1986).
  * Sánches  _et al._ [2022] H. A. C. Sánches, M. R. Jovanović, S. Kumar, A. Morozov, V. Shankar, G. Subramanian, and H. J. Wilson, Understanding viscoelastic flow instabilities: Oldroyd-b and beyond, J. Nonnewton. Fluid Mech. 302, 104742 (2022).
  * Jovanović and Kumar [2010] M. R. Jovanović and S. Kumar, Transient growth without inertia, Phys. Fluid 22, 023101 (2010).
  * Jovanović and Kumar [2011] M. R. Jovanović and S. Kumar, Nonmodal amplification of stochastic disturbances in strongly elastic channel flows, J. Nonnewton. Fluid Mech. 166, 755 (2011).
  * Page and Zaki [2014] J. Page and T. A. Zaki, Streak evolution in viscoelastic couette flow, J. Fluid Mech. 742, 520 (2014).
  * Hariharan  _et al._ [2021] G. Hariharan, M. R. Jovanović, and S. Kumar, Localized stress amplification in inertialess channel flows of viscoelastic fluids, J. Nonnewton. Fluid Mech. 291, 104514 (2021).
  * Lieu  _et al._ [2013] B. K. Lieu, M. R. Jovanović, and S. Kumar, Worst-case amplification of disturbances in inertialess couette flow of viscoelastic fluids, J. Fluid Mech. 723, 232 (2013).
  * Yesilata [2002] B. Yesilata, Nonlinear dynamics of a highly viscous and elastic fluid in pipe flow, Fluid Dyn. Res. 31, 41 (2002).
  * Yesilata [2009] B. Yesilata, Temporal nature of polymeric flows near circular pipe-exit, Polym-Plast. Technol. 48, 723 (2009).
  * Bonn  _et al._ [2011] D. Bonn, F. Ingremeau, Y. Amarouchene, and H. Kellay, Large velocity fluctuations in small-reynolds-number pipe flow of polymer solutions, Phys. Rev. E 84, 045301 (2011).
  * Pan  _et al._ [2013] L. Pan, A. Morozov, C. Wagner, and P. Arratia, Nonlinear elastic instability in channel flows at low reynolds numbers, Phys. Rev. Lett. 110, 174502 (2013).
  * Qin and Arratia [2017] B. Qin and P. E. Arratia, Characterizing elastic turbulence in channel flows at low reynolds number, Phys. Rev. Fluids 2, 083302 (2017).
  * Qin  _et al._ [2019] B. Qin, P. F. Salipante, S. D. Hudson, and P. E. Arratia, Flow resistance and structures in viscoelastic channel flows at low re, Phys. Rev. Lett. 123, 194501 (2019).
  * Jha and Steinberg [2020] N. K. Jha and V. Steinberg, Universal coherent structures of elastic turbulence in straight channel with viscoelastic fluid flow (2020), arXiv preprint https://arxiv.org/abs/2009.12258.
  * Jha and Steinberg [2021] N. K. Jha and V. Steinberg, Elastically driven kelvin–helmholtz-like instability in straight channel flow, Proceedings of the National Academy of Sciences 118, e2105211118 (2021).
  * Shnapp and Steinberg [2022] R. Shnapp and V. Steinberg, Nonmodal elastic instability and elastic waves in weakly perturbed channel flow, Phys. Rev. Fluids 7, 063901 (2022).
  * Steinberg [2022] V. Steinberg, New direction and perspectives in elastic instability and turbulence in various viscoelastic flow geometries without inertia, Low Temp. Phys. 48, 492 (2022).
  * Li and Steinberg [2023a] Y. Li and V. Steinberg, Mechanism of vorticity amplification by elastic waves in a viscoelastic channel flow, Proc. Nat. Acad. Sci. U. S. A. 120, e2305595120 (2023a).
  * Li and Steinberg [2023b] Y. Li and V. Steinberg, Elastic instability in a straight channel of viscoelastic flow without prearranged perturbations, Sci. Rep. 13, 1064 (2023b).
  * Datta  _et al._ [2022] S. S. Datta, A. M. Ardekani, P. E. Arratia, A. N. Beris, I. Bischofberger, G. H. McKinley, J. G. Eggers, J. E. López-Aguilar, S. M. Fielding, A. Frishman, M. D. Graham, J. S. Guasto, S. J. Haward, A. Q. Shen, S. Hormozi, A. Morozov, R. J. Poole, V. Shankar, E. S. G. Shaqfeh, H. Stark, V. Steinberg, G. Subramanian, and H. A. Stone, Perspectives on viscoelatic flow instabilities and elastic turbulence, Phys. Rev. Fluids 7, 080701 (2022).
  * Balkovsky  _et al._ [2001] E. Balkovsky, A. Fouxon, and V. Lebedev, Turbulence of polymer solutions, Phys. Rev. E 64, 056301 (2001).
  * Varshney and Steinberg [2019] A. Varshney and V. Steinberg, Elastic alfven waves in elastic turbulence, Nat. Commun. 10, 652 (2019).
  * Varshney and Steinberg [2018] A. Varshney and V. Steinberg, Drag enhancement and drag reduction in viscoelastic flow, Phys. Rev. Fluids 3, 103302 (2018).
  * Kumar  _et al._ [2022] M. V. Kumar, A. Varshney, D. Li, and V. Steinberg, Relaminarization of elastic turbulence, Phys. Rev. Fluids 7, L081301 (2022).
  * Anishchenko  _et al._ [1999] V. S. Anishchenko, A. B. Neiman, F. Moss, and L. Shimansky-Geier, Stochastic resonance: noise-enhanced order, Physics-Uspekhi 42, 7 (1999).
  * Anishchenko  _et al._ [1993] V. S. Anishchenko, A. B. Neiman, and M. A. Safanova, Stochastic resonance in chaotic systems, J. Stat. Phys. 70, 183 (1993).
  * Benzi  _et al._ [1981] R. Benzi, A. Sutera, and A. Vulpiani, The mechanism of stochastic resonance, J. Phys. A: Math. Gen. 14A, L451 (1981).
  * Schekochihin [2022] A. A. Schekochihin, Mhd turbulence: a biased review, J. Plasma Phys. 88, 155880501 (2022).
  * Li  _et al._ [2019] G. Li, L. A. Archer, and D. L. Koch, Electroconvection in a viscoelastic electrolyte, Phys. Rev. Lett. 122, 124501 (2019).
  * Das  _et al._ [2022] S. K. Das, G. Natale, and A. M. Benneker, Viscoelastic behavior of dilute polyelectrolyte solutions in complex geometries, J. Nonnewton. Fluid Mech. 309, 104920 (2022).
  * Emmanuel  _et al._ [2017] L. Emmanuel, S. Musacchio, and D. Vincenzi, Emergence of chaos in a viscous solution of rods, Phys. Rev. E 96, 053108 (2017).
  * Brouzet  _et al._ [2014] C. Brouzet, G. Verhille, and P. Le Gal, Flexible fiber in a turbulent flow: A macroscopic polymer, Phys. Rev. Lett. 112, 074501 (2014).
  * Doostmohammadi  _et al._ [2018] A. Doostmohammadi, J. Ignés-Mullol, J. M. Yeomans, and F. Sagués, Active nematics, Nat. Commun. 9, 3246 (2018).
  * Alert  _et al._ [2022] R. Alert, J. Casademunt, and J.-F. Joanny, Active turbulence, Annual Review of Condensed Matter Physics 13, 143 (2022).
  * Faisst and Eckhardt [2004] H. Faisst and B. Eckhardt, Sensitive dependence on initial conditions in transition to turbulence in pipe flow, J. Fluid Mech. 504, 343 (2004).
  * Wedin and Kerswell [2004] H. Wedin and R. R. Kerswell, Exact coherent structures in pipe flow: travelling wave solutions, J. Fluid Mech. 508, 333 (2004).
  * Kerswell [2005] R. Kerswell, Recent progress in understanding the transition to turbulence in a pipe, Nonlinearity 18, R17 (2005).
  * Peixinho and Mullin [2006] J. Peixinho and T. Mullin, Decay of turbulence in pipe flow, Phys. Rev. Lett. 96, 094501 (2006).
  * Mullin [2011] T. Mullin, Experimental studies of transition to turbulence in a pipe, Annu. Rev. Fluid Mech. 43, 1 (2011).
  * Hof  _et al._ [2006] B. Hof, J. Westerweel, T. M. Schneider, and B. Eckhardt, Finite lifetime of turbulence in shear flows, Nature (London) 443, 59 (2006).
  * Hof  _et al._ [2008] B. Hof, A. D. Lozar, D. J. Kuik, and J. Westerweel, Repeller or attractor? selecting the dynamical model for the onset of turbulence in pipe flow, Phys. Rev. Lett. 101, 214501 (2008).
  * Avila  _et al._ [2011] K. Avila, D. Moxey, A. De Lozar, M. Avila, D. Barkley, and B. Hof, The onset of turbulence in pipe flow, Science 333, 192 (2011).
  * Barkley [2016] D. Barkley, Theoretical perspective on the route to turbulence in a pipe, J. Fluid Mech. 803, P1 (2016).
  * Avila  _et al._ [2023] M. Avila, D. Barkley, and B. Hof, Transition to turbulence in pipe flow, Annu. Rev. Fluid Mech. 55, 575 (2023).
  * Hof  _et al._ [2003] B. Hof, A. Juel, and T. Mullin, Scaling of the threshold of pipe flow turbulence, Phys. Rev. Lett. 91, 244502 (2003).
  * Faisst and Eckhardt [2003] H. Faisst and B. Eckhardt, Traveling waves in pipe flow, Phys. Rev. Lett. 91, 224502 (2003).
  * Paranjape  _et al._ [2023] C. S. Paranjape, G. Yaln ız, Y. Duguet, N. B. Budanur, and B. Hof, Direct path from turbulence to time-periodic solutions, Phys. Rev. Lett. 131, 034002 (2023).
  * Liu  _et al._ [2009] Y. Liu, Y. Jun, and V. Steinberg, Concentration dependence of the longest relaxation times of dilute and semi-dilute polymer solutions, J. Rheol. 53, 1069 (2009).
  * Liberzon  _et al._ [2021] A. Liberzon, T. Käufer, A. Bauer, P. Vennemann, and E. Zimmer, Openpiv/openpiv-python: Openpiv-python v0.23.6 (2021).


# VISA: Variational Inference with Sequential Sample-Average Approximations

Heiko Zimmermann  Christian A. Naesseth  Jan-Willem van de Meent

###### Abstract

We present variational inference with sequential sample-average approximations
(VISA), a method for approximate inference in computationally intensive
models, such as those based on numerical simulations. VISA extends importance-
weighted forward-KL variational inference by employing a sequence of sample-
average approximations, which are considered valid inside a trust region. This
makes it possible to reuse model evaluations across multiple gradient steps,
thereby reducing computational cost. We perform experiments on high-
dimensional Gaussians, Lotka-Volterra dynamics, and a Pickover attractor,
which demonstrate that VISA can achieve comparable approximation accuracy to
standard importance-weighted forward-KL variational inference with
computational savings of a factor two or more for conservatively chosen
learning rates.

Variational Inference, Importance Sampling

##  1 Introduction

The development of general-purpose methods for Bayesian inference, such as
those provided by probabilistic programming systems (van de Meent et al.,
2018), has made it possible to apply Bayesian analysis to simulation-based
models in the natural and physical sciences. The reason that these methods are
so broadly applicable is that they generate samples by way of repeated
evaluation of the model. However, this generality often comes at significant
computational cost. Inference may require thousands of model evaluations,
where each evaluation may itself involve computationally intensive operations
such as numerical integration of a differential equation, or a call to a
numerical solver.

Gradient-based methods have become the workhorse for inference in simulation-
based models. When a model defines a fully differentiable density, methods
based on Hamiltonian Monte Carlo (Hoffman & Gelman, 2014) and reparameterized
variational inference (Kingma & Welling, 2013; Rezende et al., 2014) are often
considered the gold standard for generating high quality samples from the
posterior distribution. However, it is not always practical or possible to use
a differentiable model. The implementation of the simulator may not support
differentiation, or the model itself may not be differentiable, for example
because it employs discrete random variables or stochastic control flow. In
such cases, inference falls back on methods based on score-function estimators
(Glynn, 1990; Wingate & Weber, 2013; Ranganath et al., 2014) or importance-
weighted forward-KL variational inference (IWFVI), which derives from
reweighted wake-sleep methods (Bornschein & Bengio, 2014; Le et al., 2018).
These methods are less computationally efficient, but remain the most viable
option in a substantial number of use cases.

In this paper, we present VISA, a method that has the potential to
substantially improve the computational efficiency of variational inference
for models that are non-differentiable and computationally intensive. The key
idea is that evaluation of the variational approximation will typically be
cheap relative to that of the model. This means that we can save computation
by reusing model evaluations across multiple updates of the variational
posterior. To this end, we adapt IWFVI to employ a series of Sample-Average
Approximations (SAA), which use a fixed set of samples that defines a
deterministic surrogate to the objective, rather than generating a fresh set
of samples at each gradient step.

SAA methods were recently studied in the context of reparameterized black-box
variational inference (Giordano et al., 2024; Burroni et al., 2023) which
optimizes the reverse KL-divergence. These methods fix samples from a
parameter-free distribution, which are transformed to samples from the
approximate posterior using a differentiable map, whose parameters are
optimized to maximize the variational bound. VISA differs from these methods
in that it optimizes a forward-KL divergence and does not require a
differentiable model. Concretely, VISA fixes samples from a parameterized
variational distribution, rather than samples from a parameter-free
distribution. Since the variational distribution will change during
optimization, we construct a new SAA whenever the optimization leaves a trust
region, which we define in terms of the effective sample size. VISA is a drop-
in replacement for IWFVI, in which samples are re-used as much as possible,
thereby saving computation.

We evaluate VISA in the context of three experiments. We first consider high-
dimensional Gaussians, where the approximation error can be computed exactly.
We then consider inference in a Lotka-Volterra system and a Pickover
Attractor, where numerical integration is performed as part of the forward
simulation. Our results show that VISA with a conservative (i.e. smaller than
needed) step size can converge in a smaller number of model evaluations than
IWFVI with a more carefully tuned step size. These results come with the
caveat that VISA is more susceptible to bias than IWFVI, especially when used
with a low effective sample size threshold. In our experiments, savings of a
factor two or more are realizable with conservatively chosen learning rates,
while VISA performs on par with IWFVI for more carefully tuned step sizes.

##  2 Background

We first briefly review VI with SGD and SAAs, before we introduce VISA in
Section 3. Readers familiar with these topics can safely skip ahead.

###  2.1 Variational Inference

Variational Inference (VI) methods approximate an intractable target density
with a tractable variational distribution by solving an optimization problem.
The objective is typically to minimize a divergence measure $\mathcal{D}$
between the variational approximation $q_{\phi}$ with parameters $\phi\in\Phi$
and the target density $\pi$,

$\displaystyle\min_{\phi\in\Phi}\left\\{L(\phi):=\mathcal{D}(q_{\phi},\pi)+const.\right\\}.$
(1)

In this exposition, we assume that the target density is the posterior of a
probabilistic model $\pi(z)=p(z\mid y)$ for which we are able to point-wise
evaluate the joint density $p(y,z)$.

The two most common approaches to VI are to minimize the reverse or forward KL
divergence, for which objectives can be defined in terms of a lower and upper
bound on the log marginal $\log p(y)$,

$\displaystyle L_{\textsc{r}}(\phi)$
$\displaystyle=-\mathop{\mathbb{E}}_{q_{\phi}}\left[\log\frac{p(y,z)}{q_{\phi}(z)}\right]\\!=\text{KL}(q_{\phi}\,||\,p(\cdot|y))-\log
p(y),$ $\displaystyle L_{\textsc{f}}(\phi)$
$\displaystyle=\\!\\!\\!\\!\mathop{\mathbb{E}}_{p(\cdot\mid
y)}\\!\\!\left[\log\frac{p(y,z)}{q_{\phi}(z)}\right]\\!=\text{KL}(p(\cdot|y)\,||\,q_{\phi})+\log
p(y).$

We will briefly discuss standard reparameterized VI, which maximizes the lower
bound $-L_{\textsc{r}}$, and IWFVI, which minimizes the upper bound
$L_{\textsc{f}}$.

#### Reparameterized VI.

When maximizing a lower bound with stochastic gradient descent, we can either
employ score-function estimators (Ranganath et al., 2014; Paisley et al.,
2012), which tend to exhibit a high degree of variance, or make use of the
reparameterization trick, which is commonly employed in variational
autoencoders (Kingma & Welling, 2013; Rezende et al., 2014).

Reparameterization is a method for generating $z\sim q_{\phi}$ by way of a
parametric pushforward of random variables $\xi$ that are distributed
according to a parameter-free density $q_{\xi}$,

$\displaystyle z$ $\displaystyle=T_{\phi}(\xi),$ $\displaystyle\xi\sim
q_{\xi}.$

In VI, reparameterized samples can be used to compute an unbiased estimate of
the gradient,

$\displaystyle-\frac{d}{d\phi}L_{\textsc{r}}(\phi)$
$\displaystyle=\mathop{\mathbb{E}}_{q_{\xi}}\left[\frac{d}{d\phi}\log\frac{p(y,T_{\phi}(\xi))}{q_{\phi}(T_{\phi}(\xi))}\right]$
$\displaystyle\approx\frac{1}{N}\sum_{i=1}^{N}\frac{d}{d\phi}\log\frac{p(y,T_{\phi}(\xi^{(i)}))}{q_{\phi}(T_{\phi}(\xi^{(i)}))},$
$\displaystyle\xi^{(i)}\sim q_{\xi}.$

Reparameterized VI requires a model $p(y,z)$ that is differentiable with
respect to $z$ in order to compute the pathwise derivative (Rezende et al.,
2014),

$\displaystyle\frac{d}{d\phi}\log
p(y,T_{\phi}(\xi))=\left.\frac{\partial}{\partial z}\log
p(y,z)\right|_{z=T_{\phi}(\xi)}\frac{\partial T_{\phi}(\xi)}{\partial\phi}.$

This means that the model must support (automatic) differentiation and must
not make use of discrete random variables or stochastic control flow, such as
if statements that branch on random values (van de Meent et al., 2018).

#### Importance-Weighted Forward-KL VI.

When minimizing the forward KL divergence, we approximate the gradient

$\displaystyle-\frac{d}{d\phi}L_{\textsc{f}}(\phi)=\mathop{\mathbb{E}}_{p(\cdot\mid
y)}\\!\\!\left[\frac{d}{d\phi}\log q_{\phi}(z)\right].$

This does not require differentiability of the model $p(y,z)$ with respect to
$z$, but does require approximate inference to generate samples from the
posterior $p(z\mid y)$.

In IWFVI, the variational distribution $q_{\phi}$ is used as a proposal for a
self-normalized importance sampler. This approach makes use of the fact that
we can express the gradient as an expectation with respect to $q_{\phi}$,

$\displaystyle-\frac{d}{d\phi}L_{\textsc{f}}(\phi)=\mathop{\mathbb{E}}_{q_{\phi}}\\!\\!\left[w_{\phi}(z)\frac{d}{d\phi}\log
q_{\phi}(z)\right],~{}~{}w_{\phi}(z)=\frac{p(z\,|\,y)}{q_{\phi}(z)}.$

Here the ratio $w_{\phi}(z)$, known as the importance weight, needs to be
approximated, since we cannot compute $p(z\mid y)$. To do so, we define the
self-normalized weights

$\displaystyle\hat{w}_{\phi}^{(i)}$
$\displaystyle=\frac{\bar{w}_{\phi}^{(i)}}{\sum_{j=1}^{N}\bar{w}_{\phi}^{(j)}},$
$\displaystyle\bar{w}_{\phi}^{(i)}=\frac{p(y,z^{(i)})}{q_{\phi}(z^{(i)})},$

leading to the gradient estimate

$\displaystyle\frac{d}{d\phi}L_{\textsc{f}}(\phi)\simeq\sum_{i=1}^{N}\hat{w}^{(i)}\frac{d}{d\phi}\log
q_{\phi}(z^{(i)}),$ $\displaystyle z^{(i)}\sim q_{\phi}.$ (2)

The resulting estimate is biased but consistent, meaning that it converges to
the true gradient as $N\to\infty$ almost surely.

###  2.2 VI with Sample-Average Approximations

In the stochastic optimization literature, sample-average approximations are
used to approximate an expected loss with a surrogate loss in the form of a
Monte Carlo estimate (see Kim et al. (2015) for a review). Importantly, the
samples that the SAA is based on remain _fixed_ throughout the optimization
process. This means that the surrogate objective can be treated like any other
deterministic function, which can be optimized using second-order methods and
other standard optimization tools.

Concretely, a sample-average approximation applies to an optimization problem
of the form

$\displaystyle\min_{\phi\in\Phi}\left\\{L(\phi):=\mathop{\mathbb{E}}_{\rho}\left[\ell(z,\phi)\right]\right\\},$
(3)

in which the density $\rho(z)$ does not depend on the parameters $\phi$. This
means that we can compute a surrogate loss $\hat{L}(\phi)$ that is an unbiased
estimate of the original loss $L(\phi)$ by averaging over samples from $\rho$,

$\displaystyle\hat{L}(\phi)$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}l(z^{(i)},\phi),$ $\displaystyle
z^{(i)}$ $\displaystyle\sim\rho.$

Under mild conditions on $\ell$ and $\rho$, as the number of samples
$N\to\infty$, the minimizer $\hat{\phi}=\arg\min_{\phi}\hat{L}(\phi)$ and the
minimal value $\hat{L}(\hat{\phi})$ converge almost surely to the minima
$\phi^{*}=\arg\min{L}_{\phi}(\phi)$ and minimal value $L(\phi^{*})$ of the
original problem.

In the context of reparameterized VI, a sample-average approximation can be
constructed by fixing a set of samples $\\{\xi^{(i)}\sim q_{\xi}\\}_{i=1}^{N}$
from a parameter-free distribution,

$\displaystyle\hat{L}_{\textsc{r}}(\phi)$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}\log\frac{p\big{(}y,T_{\phi}(\xi^{(i)})\big{)}}{q_{\phi}\big{(}T_{\phi}(\xi^{(i)})\big{)}}.$

In an SAA-based approach to reparameterized VI (Giordano et al., 2024; Burroni
et al., 2023), optimization of the parameters $\phi$ will move the transformed
samples $z^{(i)}=T_{\phi}(\xi^{(i)})$ to match the posterior density, whilst
keeping the noise realizations $\xi^{(i)}$ fixed. Empirical evaluations show
that combining the SAA approximation with an off-the-shelf second-order
optimizer can result in substantial computational gains as well as more
reliable convergence to the optimum.

##  3 SAA for Forward-KL Variational Inference

The primary motivation behind existing SAA-based methods for reparameterized
VI (Giordano et al., 2024; Burroni et al., 2023) is that fixing the noise
realizations defines a completely deterministic surrogate objective
$\hat{L}_{\textsc{r}}(\phi)$, which can then be used with any number of
existing optimizers. The main requirement from an implementation point of view
is that the model density $p(y,z)$ is differentiable with respect to $z$. In
this setting it is also necessary to evaluate the model for every update,
since any change to $\phi$ also changes the values of the transformed samples
$z^{(i)}=T_{\phi}(\xi^{(i)})$.

In developing VISA, both our motivation and implementation requirements are
somewhat different. Our primary interest is in minimizing the total number of
model evaluations at convergence. We also wish to develop a method that is
applicable when the model density $p(y,z)$ is not differentiable, either
because the implementation simply does not support (automatic) derivatives, or
because the model incorporates discrete variables or stochastic control flow,
which introduce discontinuities in the density $p(y,z)$.

To this end, we propose a method that optimizes a forward KL with an
importance weighted objective that incorporates ideas from SAA-based
approaches. In a setting where we already have access to samples from the
posterior, we could trivially define a SAA for the upper bound
$L_{\textsc{f}}(\phi)$,

$\displaystyle\hat{L}_{\textsc{f}}(\phi)$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}\log\frac{p(y,z^{(i)})}{q_{\phi}(z^{(i)})},$
$\displaystyle z^{(i)}$ $\displaystyle\sim p(\cdot\mid y).$

In practice, this naive approach is unlikely to be useful in a setting where
evaluation of $p(y,z)$ is computationally expensive, since we would still need
to carry out approximate inference to generate a set of samples from the
posterior.

We therefore adopt the approach used in IWFVI, which uses the variational
distribution as a proposal in a self-normalized importance-sampler. To define
an SAA for the objective in this setting, we express the objective at
parameters $\phi$ in terms of an expectation with respect to a distribution
from the same family with fixed parameters $\tilde{\phi}$,

$\displaystyle L_{\textsc{f}}(\phi)$
$\displaystyle=\mathop{\mathbb{E}}_{p(\cdot|y)}\left[\log\frac{p(y,z)}{q_{\phi}(z)}\right]=\mathop{\mathbb{E}}_{q_{\tilde{\phi}}}\left[w_{\tilde{\phi}}(z)\log\frac{p(y,z)}{q_{\phi}(z)}\right].$

As in standard IWFVI, the importance weights $w_{\tilde{\phi}}(z)$ are
intractable, but we can define an SAA for the objective in terms of self-
normalized weights,

$\displaystyle\hat{L}_{\textsc{f}}(\phi\,;\tilde{\phi})$
$\displaystyle=\sum_{i=1}^{N}\hat{w}^{(i)}_{\tilde{\phi}}\log\frac{p(y,z^{(i)})}{q_{\phi}(z^{(i)})},\qquad
z^{(i)}\sim q_{\tilde{\phi}}.$

In this surrogate objective, which we will optimize with respect to $\phi$,
the quality of the approximation depends on how closely the proposal with
parameters $\tilde{\phi}$ matches the posterior. Since our approximation of
the posterior will improve during optimization, we will update $\tilde{\phi}$
to the current parameter values $\phi$ at some interval, resulting in an
approach that we will refer to as a _sequential_ sample-average approximation.

To determine when we need to generate a fresh SAA, we will define the notion
of a trust region. This defines an optimization process in which the SAA is
refreshed whenever the optimization trajectory leaves the current trust
region. This optimization process is illustrated in Figure 1 and described
schematically in Algorithm 1. We begin by setting the proposal parameters
$\tilde{\phi}=\phi_{0}$ to the initial variational parameters, generating a
set of samples $\mathcal{Z}=\\{z^{(i)}\\}_{i=1}^{N}$ and defining an SAA of
the objective $\hat{L}_{\textsc{f}}(\phi;\tilde{\phi})$ and a trust region
$S_{\mathcal{Z},\alpha}(\tilde{\phi})$ based on these samples. We then
repeatedly update $\phi_{t}$ using an optimizer until the value $\phi_{t}$ no
longer lies in the trust region. At this point, we update the proposal
parameters $\tilde{\phi}=\phi_{t}$ and generate a fresh sample set, which we
then use to update the SAA and the trust region.

Algorithm 1 VISA

0: Initial param. $\phi_{0}$, trust region threshold $\alpha$, data $y$

${\tilde{\phi}}\leftarrow\phi_{0}$ $\triangleright$ Initialize proposal
parameter

$\mathcal{Z}\leftarrow\\{z^{(i)}\sim q_{{\tilde{\phi}}}\\}_{i=1}^{N}$
$\triangleright$ Initialize samples

$\hat{L}_{\textsc{f}}(\phi\,;{\tilde{\phi}})=\sum_{i=1}^{N}\hat{w}^{(i)}_{{\tilde{\phi}}}\log\frac{p(y,z^{(i)})}{q_{\phi}(z^{(i)})}$
$\triangleright$ Initialize SAA

for $t=1,\ldots,T$ do

$\phi_{t}=\mathrm{optimizer}\text{-}\mathrm{step}(\hat{L}_{\textsc{f}},\phi_{t-1})$

if $\phi_{t}\notin S_{\mathcal{Z},\alpha}({\tilde{\phi}})$ then
$\triangleright$ Not inside trust region

${\tilde{\phi}}\leftarrow\phi_{t}$ $\triangleright$ Update proposal

$\mathcal{Z}\leftarrow\\{z^{(i)}\sim q_{{\tilde{\phi}}}\\}_{i=1}^{N}$
$\triangleright$ Refresh samples

$\hat{L}_{\textsc{f}}(\phi\,;{\tilde{\phi}})=\sum_{i=1}^{N}\hat{w}^{(i)}_{{\tilde{\phi}}}\log\frac{p(y,z^{(i)})}{q_{\phi}(z^{(i)})}$
$\triangleright$ Refresh SAA

end if

end for

#### Defining Trust regions.

To define a notion of a trust region, we compute an effective sample size
(ESS), which is a proxy measure for the variance of the importance weights.
Notably, the importance weights can be decomposed into two parts, (1) the
ratio of the variational density and the trust region density, and (2) the
ratio between posterior and variational density,

$\displaystyle w_{{\tilde{\phi}}}(z)=$ $\displaystyle
v_{\phi,{\tilde{\phi}}}(z)\>w_{\phi}(z),$ $\displaystyle
v_{\phi,{\tilde{\phi}}}(z)$
$\displaystyle=\frac{q_{\phi}(z)}{q_{{\tilde{\phi}}}(z)}.$

The variance of $w_{\phi}$ is independent of the fixed proposal parameters and
decreases during optimization. The variance of $v_{\phi,{\tilde{\phi}}}$
measures how similar $q_{\tilde{\phi}}$ is to $q_{\phi}$ and can be controlled
by updating the proposal parameters ${\tilde{\phi}}$ to the parameters of the
current variational approximation $\phi$ once the ESS drops below a certain
threshold.

We formalize this notion by defining a set-valued function which, for a given
threshold $\alpha$ and samples $\mathcal{Z}=\\{z^{(i)}\\}_{i=1}^{N}$, maps
each parameter to a corresponding trust region based on a scoring function,
which we choose to be the ESS,

$\displaystyle S_{\mathcal{Z},\alpha}({\tilde{\phi}})$
$\displaystyle=\\{\phi\in\Phi\ \mid
s_{\mathcal{Z}}({\tilde{\phi}},\phi)>\alpha\\}$ $\displaystyle
s_{\mathcal{Z}}({\tilde{\phi}},\phi)$
$\displaystyle=\frac{\left(\sum_{i=1}^{N}v_{{\tilde{\phi}},\phi}(z^{(i)})\right)^{2}}{N\sum_{i=1}^{N}\left(v_{{\tilde{\phi}},\phi}(z^{(i)})\right)^{2}}.$

In other words, for a given ESS threshold $\alpha$ we can verify $\phi\in
S_{\mathcal{Z},\alpha}({\tilde{\phi}})$ by checking
$s_{\mathcal{Z}}({\tilde{\phi}},\phi)>\alpha$. We visualize how new trust
regions, corresponding to different SAAs, are constructed sequentially during
optimization in Figure 1.

Figure 1: Visualization of parameter traces and trust regions corresponding to
different SAAs. If after an update $\phi\notin
S_{\mathcal{Z},\alpha({\tilde{\phi}})}$, we set ${\tilde{\phi}}\leftarrow\phi$
to construct a new SAA and corresponding trust region.

#### Effect of the Trust-Region on Convergence.

For a low enough threshold $\alpha$ the algorithm might converge to an optimal
parameter $\hat{\phi}$ within the trust region of the current SAA that does
not yet satisfy our global convergence criteria, i.e. sufficiently minimizes
the forward KL-divergence. In these cases, if the variational approximation is
not degenerate, we can try recover by decreasing $\alpha$ such that
$\hat{\phi}\notin S_{\alpha}({\tilde{\phi}})$ and continue optimization. As
$\alpha\to 1$ the frequency of sample acquisition increases and, for
$\alpha=1$ VISA reduces to standard IWFVI. In this work we chose $\alpha$ high
enough that the algorithm does not converge prematurely to an optima of an
intermediate SAA. In these cases we find that convergence of the training loss
can be used as indicator for convergence of the forward KL-divergence or
corresponding upper bound, which we verify in our experiments. We also
experimented with caching past sample sets to compute a secondary loss based
on the last $M$ SAAs. While we did not find it to add additional value in our
experiments, it be a useful tool to assess convergence in other settings.

#### Efficient Implementation.

To avoid recomputing density values for old sample locations, we cache both
sample location $z^{(i)}$ and the corresponding log-joint density of the model
$\log p(z^{(i)},y)$. If sampling from the proposal is cheap and memory is of
concern, e.g. for large samples set or if past sample sets are stored to
compute a validation loss, we can store the random seed instead of the sample
and rematerialize the sample when needed.

##  4 Related Work

#### VI with SAAs

Recent work that studies SAAs (Giordano et al., 2024; Burroni et al., 2023) in
the context of variational inference focuses on the reparameterized black-box
VI setting and optimizes a reverse KL-divergence. These methods rely on
reparameterization to move samples to areas of high posterior density while
keeping a fixed set of noise realizations from the base distribution, which
does not depend of the variational parameters. Optimizing a deterministic
objective allows the authors to use second-order optimization and linear
response methods (Giordano et al., 2015) to fit covariances. While these
methods allow to realize substantial gains in terms of inference quality and
efficiency, in contrast to VISA, they require differentiability of the model.

#### Stochastic second-order optimization.

There is also work outside of the context of SAAs that aims to incorporate
second order information to improve stochastic optimization and variational
inference. Byrd et al. (2016) propose batched-L-BFGS, which computes stable
curvature estimates by sub-sampled Hessian-vector products instead of
computing gradient differences at every iteration. This work has also been
been adopted to the variational inference setting by Liu & Owen (2021).
Pathfinder (Zhang et al., 2022) uses a quasi-Newton method to find the mode of
the a target density and construct normal approximations to the density along
the optimization path. The intermediate normal approximations are used to
define a variational approximation that minimizes an evidence lower bound.
Similar to SAA-based methods, pathfinder can reduce the number of model
evaluations by up to an order of magnitude compared to HMC, but requires a
differentiable model.

#### VI with forward KL-divergence.

VISA is also related to other methods that aim to optimize a forward KL-
divergence or its stochastic upper bound. This includes reweighted-wake sleep
(and wake-wake) methods (Bornschein & Bengio, 2014; Le et al., 2018) to which
we compare VISA in the experiment section, as well as their doubly-
reparameterized variants (Tucker et al., 2018; Finke & Thiery, 2019; Bauer &
Mnih, 2021), which are not directly comparable as they require a
differentiable model. While the methods above use a single importance sampling
step using the variational approximation as a proposal, other methods use more
complex proposal including MCMC proposals (Naesseth et al., 2020; Zhang et
al., 2023), approximate Gibbs kernels (Wu et al., 2020), or proposal defined
by probabilitic programs (Stites et al., 2021; Zimmermann et al., 2021). While
these methods do not necessarily require a differentiable proposal they are
not designed to be sample efficient but to approximate complex target
densities.

Figure 2: Symmetric KL-divergence as a function of the number of model
evaluations for a Gaussian target with diagonal covariance matrix (top row)
and dense covariance matrix (bottom row). For small learning rates (0.001,
0.005, 0.01) IWFVI and BBVI-SF, need a larger number of model evaluations to
converge. VISA converges much faster as it compensates for the small step size
by reusing samples. For a learning rate of 0.05 VISA fails to reliably
converge, while IWFVI still converges. Overall, VISA converges faster or at
the same rate as IWFVI and BBVI-SF with the same or higher learning rates.

##  5 Experiments

We compare VISA to standard importance-weighted VI (IWFVI) in terms of
inference quality and the number of model evaluations. We assess inference
quality by measuring or estimating the forward KL-divergence to the true
posterior, or when such an estimate is not available by assessing the log-
joint probability of samples under the model. In the Gaussian experiment we
also compare to black-box VI, which optimizes the reverse KL-divergence, and
therefore report inference quality by measuring the symmetric KL-divergence to
allow for a fair comparison.

###  5.1 Gaussians

To study the effect of different learning rates and ESS threshold parameters,
we first evaluate VISA on approximating medium- to high- dimensional Gaussians
and compare the inference performance over the number of model evaluations to
IWFVI, standard reparameterized variational inference (BBVI-RP) and
variational inference using a score-function gradient estimator (BBVI-SF).
Notably, we include BBVI-RP as a reference only, showcasing that faster
convergence can be achieved by leveraging the differentiability of the model,
and do not compare to it directly. To allow for a fair comparison between
methods that optimize a forward KL-divergence (VISA, baseline) and methods
that optimize a reverse KL-divergence (BBVI), we evaluate the inference
qualify in terms of the symmetric KL-divergence.

We study two different target densities, a (1) $D=128$ dimensional Gaussian
with a diagonal covariance matrix and (2) $D=32$ dimensional Gaussian with a
dense covariance matrix. For the diagonal covariance matrix, we scale the
variances of the individual dimension from $\sigma_{\min}=0.1$ to
$\sigma_{\max}=1$ such that the covariance matrix takes the form

$\displaystyle
C_{\mathrm{diag}}=\mathrm{diag}\left(\left[\sigma_{\min}+(i-1)*\frac{\sigma_{max}-\sigma_{min}}{D-1}\right]_{i=1}^{D}\right).$

To create the dense covariance matrix we first sample a random positive semi-
definite matrix $M=AA^{T}$, where $A_{ij}\sim\mathcal{U}(0,1)$ and then
construct the covariance matrix

$\displaystyle
C_{\mathrm{dense}}=\left(\frac{M}{||M||_{\mathrm{F}}}+0.1\mathbb{I}\right).$

Figure 2 shows the results for different learning rates
$\text{lr}\in\\{0.001,0.005,0.01,0.05\\}$ and ESS threshold parameters
$\alpha\in\\{0.9,0.95,0.99\\}$. We compute gradient estimates for VISA, IWFVI,
and BBVI-SF with $N=10$ samples, and gradient estimates for BBVI-RP using a
single sample.

We observe that VISA converges substantially faster than IWFVI and BBVI-SF at
lower learning rates (0.001, 0.005, 0.01) for both, targets with diagonal- and
dense covariance matrix. The difference in the convergence rate gets less
pronounced as the learning rate increases. For large learning rates (0.05)
VISA fail to converge reliably, while IWFVI still converges. For even higher
learning rates all methods but BBVI-RP fail to converge. If we compare
convergence across learning rates, we observe that VISA converges faster or at
the same rate as IWFVI with the same or higher learning rate. In the Gaussian
experiment, the overall effect of the threshold $\alpha$ on the convergence
rate is minor, but we observe slightly faster convergence for lower ESS
thresholds. We also find that, especially for lower ESS thresholds, VISA is
also more susceptible to underestimating posterior variance and might not
fully converge in the final phase of training. We hypothesize that this is due
to the deterministic optimization procedure, which overfits to high
probability samples if the samples are not refreshed frequently enough. A
possible remedy for this is to increase the ESS threshold during the final
phase of training, for $\alpha\to 1$ this recovers the behaviour of IWFVI.
However, since the overall effect of the threshold on the rate of convergence
is small, a viable strategy is to employ a fixed threshold value that is
somewhat conservative, such as $\alpha=0.95$ or $\alpha=0.99$.

Overall we find that VISA converges faster or at the same rate as IWFVI and
BBVI-SF for a range of different learning rates and ESS thresholds. Moreover
VISA is substantially more robust to the choice of learning rate. Using a
deterministic objective allows us to draw fresh samples less frequently and
consequently requires less evaluations of a potentially expensive to evaluate
model log-joint density.

###  5.2 Lotka-Volterra

The Lotka-Vorterra predator-prey population dynamics (Lotka, 1925; Volterra,
1927) are modeled by a pair of first-order ordinary differential equations
(ODEs),

$\displaystyle\frac{du}{dt}=(\alpha-\beta v)u,$
$\displaystyle\frac{dv}{dt}=(-\gamma+\delta u)v,$

where $v$ denotes the predator- and $u$ denotes the prey population. We will
in the following denote the pair of predator-prey populations at time $t$ with
$z_{t}=(u_{t},v_{t})$. The dynamics of the ODE are governed by its respective
population growth and shrinkage parameters
$\theta=(\alpha,\beta,\gamma,\delta)$, which we would like to infer together
with the initial conditions of the system given noisy observations
$y_{1:T}=(y_{1},\ldots,y_{T})$.

Figure 3: Results for Lotka-Volterra model with different learning rates. (Top
row) Training objective over number of model evaluations. (Middle row)
Approximate forward KL-divergence computed on reference samples obtained by
MCMC. For smaller step sizes (0.001, 0.005) VISA achieves comparable forward
KL-divergence to IWFVI while requiring significantly less model evaluations to
converge (see vertical lines). For larger step sizes (0.01) VISA only
converges with a high ess threshold (0.99) for which it requires approximately
the same number of evaluations as IWFVI. (Bottom row) Gradient steps over
number of batch evaluations of the model, each batch evaluation corresponds to
evaluating a batch of $N=100$ samples. VISA requires fewer evaluations per
gradient step compared to IWFVI.

Following Carpenter (2018), we place priors over the initial population size
$z_{0}$ and system parameters $\theta$

$\displaystyle z^{\mathrm{prey}}_{0},z^{\mathrm{pred}}_{0}$
$\displaystyle\sim\mathrm{LogNormal}(\log(10),1),$ (4)
$\displaystyle\alpha,\gamma$ $\displaystyle\sim\mathrm{Normal}(1,0.5),$ (5)
$\displaystyle\beta,\delta$ $\displaystyle\sim\mathrm{Normal}(0.05,0.05),$ (6)

and assume a fractional observation error,

$\displaystyle y^{\mathrm{prey}}_{t},y^{\mathrm{pred}}_{t}$
$\displaystyle\sim\mathrm{LogNormal}(\log z_{t},\sigma_{t})$
$\displaystyle\sigma_{t}$ $\displaystyle\sim\mathrm{LogNormal}(-1,1).$

Given an initial population $z_{0}$, system parameters $\theta$, and
observations $y_{1:T}$ we can solve the ODE numerically to obtain approximate
population sizes $z_{1:T}$ for time steps $1,\ldots,T$ which we use to compute
the likelihood of the observed predator-prey populations,

$\displaystyle p(y_{1:T}\mid z_{0},\theta)$
$\displaystyle=\prod_{t=0}^{T}p(y_{t}\mid z_{t}).$

Our goal is to learn an approximation to the posterior $p(\theta,z_{0}\mid y)$
by minimizing the evidence upper bound

$\displaystyle L_{\textsc{f}}(\phi):=\mathbb{E}_{(z_{0},\theta)\sim
p(\cdot,\cdot\mid
y)}\left[\log\frac{p(z_{0},y,\theta)}{q_{\phi}(z_{0},\theta)}\right].$ (7)

We model the variational approximation $q_{\phi}$ for the interaction
parameters $\theta$ and initial population sizes $z_{0}$ as jointly log-normal
and initialize $\phi$ such that the the marginal over $z_{0}$ matches the
prior (Equation 4) and the marginal over $\theta$ has similar coverage to the
prior. We approximate the objective and its gradient using with $N=100$
samples. To specify a common convergence criteria, we compute the highest
common test loss value for VISA with $\alpha=0.99$ and IWFVI that is not
exceeded by more than $1$ nat by all consecutive test loss values. The
convergence threshold computed this way is $-712.6$ nats.

To evaluate the inference performance, we first generate $N=10\,000$
approximate posterior samples using a No-U-Turn Sampler (NUTS) (Hoffman &
Gelman, 2014) with $10\,000$ burn-in steps and window-adaption, which
generally provides good performance out of the box111NUTS is an adaptive
Hamiltonian Monte Carlo sampler and uses the gradient information of the model
to guide the generation of proposals. As such it requires the log-joint
density model to be differentiable which is not required by VISA or IWFVI.. We
use the approximate posterior samples to approximate an “oracle” for the upper
bound in Equation 7,

$\displaystyle\hat{L}^{\textsc{nuts}}_{\textsc{f}}(\phi)=\frac{1}{N}\sum_{i=1}^{N}\log\frac{p(z^{(i)}_{0},y,\theta^{(i)})}{q_{\phi}(z_{0}^{(i)},\theta^{(i)})},$

which we evaluate along with the training objective during optimization to
assess convergence.

We find that VISA is able to obtain variational distribution of similar
quality to IWFVI while requiring significantly fewer model evaluations for
smaller learning rates (see Figure 3). Interestingly, VISA requires
significantly fewer model evaluations per gradient step during the early
stages of training, while requiring slightly more evaluations per gradient
step thereafter. We hypothesise that this is again a result of under
approximating posterior variance in the later stages of training. As a result,
even small changes in the variational distribution can lead to big changes in
the ESS, which triggers the drawing of fresh samples. For VISA, we also find a
more pronounced difference between different ESS thresholds and their
influence on convergence. Runs with a higher ESS thresholds converge more
stably and are able to achieve lower test loss in the final stages of
training.

###  5.3 Pickover Attractor

Figure 4: Results for Pickover attractor. (a) Approximate log-joint density
over number of batch-evaluations of model. (b) Log-joint approximation plotted
over domain of prior. The variational approximation capture the high density
area containing the data. (c) Visualization of pickover attractor with ground
truth parameters $\theta=[-2.3,1.25]$. (d) Visualization of attractor with
average system parameters computed over $10.000$ samples from the learned
variational approximation. Each evaluation in the plot corresponds to
evaluating a batch of $N=10$ samples.

Following Rainforth et al. (2016), we model a 3D Pickover attractor (Pickover,
1995) with parameters $\theta=(\beta,\eta)$,

$\displaystyle x_{t+1,1}$ $\displaystyle=\sin(\beta
x_{t,2})-\cos(2.5x_{t,1})x_{t,3}$ $\displaystyle x_{t+1,2}$
$\displaystyle=\sin(1.5x_{t,1})x_{t,3}-\cos(\eta x_{t,2})$ $\displaystyle
x_{t+1,3}$ $\displaystyle=\sin(x_{t,1}).$

Due to its chaotic nature the system is sensitive to small perturbations in
its initial state, i.e. even small variations in the initial state lead to
exponentially fast diverging trajectories. Therefore, to track the evolution
of the system, we employ a bootstrap particle filter (Gordon et al., 1993)
which assumes noisy observations $y:=y_{1:T}$ and introduces auxiliary
variables $z:=z_{1:T}$ to model the latent state of the system. We define the
prior over system parameters

$\displaystyle p(\theta)=\begin{cases}1/18&-3\leq\theta_{1}\leq
3,0\leq\eta\leq 3\\\ 0&\text{otherwise}\end{cases},$

and model the transition and observation as,

$\displaystyle z_{0}$
$\displaystyle\sim\mathcal{N}(\cdot\mid\mathbf{0}_{3},\mathbb{I}_{3}),$
$\displaystyle z_{t}$ $\displaystyle\sim\mathcal{N}(\cdot\mid
h(z_{t-1},\theta),\sigma_{z})\quad$ $\displaystyle\text{for}\ t>0,$
$\displaystyle y_{t}$ $\displaystyle\sim\mathcal{N}(\cdot\mid
z_{t},\sigma_{y})\ $ $\displaystyle\text{for}\ t\geq 0,$

where $\sigma_{z}=0.01$, $\sigma_{y}=0.2$, and $h$ evolves the system by one
time step using the equations of the Pickover attractor described above. The
particle filter is used to simulate $T=100$ time steps with $M=500$ particles,
which renders evaluating the model expensive.

To restrict the proposal to the same domain as the prior by first sampling
from a Normal, which we denote $\bar{q}_{\phi}$ and parameterized with a mean
and lower Cholesky factor, and then transform the sample by

$f(\theta)=(\tanh(3\cdot\theta_{1}),\tanh(1.5\cdot\theta_{2}+1.5)).$

The density of the transformed samples is
$q_{\phi}(\theta)=\bar{q}_{\phi}(\theta)\lvert\mathrm{det}\frac{df(\theta)}{d\theta}\rvert^{-1}$,
which is appropriately restricted to the domain of the prior.

We are interested in approximating the marginal posterior $p(\theta\mid y)$
over system parameters by optimizing the evidence upper bound

$\displaystyle\mathop{\mathbb{E}}_{p_{\theta}(\theta\mid
y)}\left[\log\frac{p_{\theta}(y,\theta)}{q_{\phi}(\theta)}\right]$
$\displaystyle=\mathop{\mathbb{E}}_{p_{\theta}(\theta\mid
y)}\left[\log\frac{\mathbb{E}_{p_{{}_{\mathrm{pf}}}}[\hat{p}(y\mid\theta)]}{q_{\phi}(\theta)}\right]$
$\displaystyle\leq\mathop{\mathbb{E}}_{p_{\theta}(\theta\mid
y)}\left[\mathbb{E}_{p_{{}_{\mathrm{pf}}}}\left[\log\frac{\hat{p}(y\mid\theta)}{q_{\phi}(\theta)}\right]\right].$

To obtain a tractable objective we replace the intractable marginal likelihood
$p(y\mid\theta)\approx\hat{p}_{\theta}(y\mid\theta)$ with the marginal
likelihood estimate obtained by running the particle filter (Naesseth et al.,
2019), similar to pseudo-marginal methods (Andrieu et al., 2010) and
approximate the gradient with $N=10$ samples. As the likelihood estimate is
non-differentiable due to the discrete ancestor choices made inside the
particle filter, we cannot run NUTS to obtain approximate posterior sample as
before, but instead report the log-joint density of the variational
distribution.

We observe that VISA converges more stably with fewer samples compared to
IWFVI and find that attractors corresponding to samples from the variational
approximation look qualitatively similar to those based on the true
parameters. We summarize the result in Figure. 4.

##  6 Discussion and Limitations

In this paper we developed VISA, a method for approximate inference for
expensive to evaluate models that optimizes the forward KL-divergence through
a sequence of SAAs. Each SAA is optimized deterministically and requires and
fixes a single set of samples, hereby requiring new model evaluations only
when the SAA is refreshed. To track the approximation quality of the current
SAA, VISA computes the ESS of the ratio between the current variational
distribution and the proposal distribution that was used to construct the SAA.
If the ESS falls below a predefined threshold, a new SAA approximation is
constructed based on fresh samples from the current variational distribution.
We observe gains of a factor 2 or more in terms of the number of required
model evaluations for conservatively chosen step sizes, while achieving
similar posterior approximation accuracy as IWFVI, the equivalent method that
does not employ the sequential sample-average approximation.

#### Underapproximation of posterior variance.

Both reparameterized VI, which optimizes the reverse KL-divergence, and IWFVI,
which optimizes the forward KL-divergence via importance sampling, are prone
to under approximating posterior variance. In the case of reparameterized VI,
this can often be attributed to the mode seeking behaviour of the reverse KL-
divergence, while in IWFVI the low effective samples sizes can lead to over
fitting to a small number of high-weight samples. We found that keeping the
samples fixed for too long, i.e. using an ESS threshold that is too low, can
exacerbate this problem, as the optimizer can take multiple steps towards the
same high-weight samples.

Giordano et al. (2024) and Burroni et al. (2023) showed that when applying SAA
to reparameterized VI, it is possible to make use of second-order methods. We
experimented with optimizing SAAs with L-BGFS, which is a quasi-Newton method
with line search. However, we found that in the setting of optimizing a
forward KL with relatively few samples, L-BGFS can amplify the problem of
overfitting, often leading to instabilities and collapsed variational
distributions.

#### Number of latent variables and parameters.

Because VISA employs a relatively small number of samples and does not refresh
samples at every iteration, we found that it is not well-suited to models with
a large number of latent variables or large number of parameters. This agrees
with theoretical findings by Giordano et al. (2024), who show that SAAs for a
full covariance Gaussian fail if the number of samples is not at least in the
same regime as the number of latent dimensions. Burroni et al. (2023) manage
to train full covariance normal approximation by using a sequence of SAAs
using increasingly large sample sizes, however, this is directly opposed to
our goal of reducing the number of model evaluations in expensive to evaluate
models.

## Acknowledgement

The authors would like to thank Tamara Broderick for helpful discussion about
sample-average approximations in the context of reparameterized variational
inference.

## References

  * Andrieu et al. (2010) Andrieu, C., Doucet, A., and Holenstein, R.  Particle Markov Chain Monte Carlo Methods.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 72(3):269–342, June 2010.
  * Bauer & Mnih (2021) Bauer, M. and Mnih, A.  Generalized Doubly Reparameterized Gradient Estimators.  In _International Conference on Machine Learning_ , pp. 738–747. PMLR, 2021.
  * Bornschein & Bengio (2014) Bornschein, J. and Bengio, Y.  Reweighted wake-sleep.  _arXiv preprint arXiv:1406.2751_ , 2014.
  * Burroni et al. (2023) Burroni, J., Domke, J., and Sheldon, D.  Sample Average Approximation for Black-Box VI.  _arXiv preprint arXiv:2304.06803_ , 2023.
  * Byrd et al. (2016) Byrd, R. H., Hansen, S. L., Nocedal, J., and Singer, Y.  A Stochastic Quasi-Newton Method for Large-Scale Optimization.  _SIAM Journal on Optimization_ , 26(2):1008–1031, January 2016.  ISSN 1052-6234.  doi: 10.1137/140954362.
  * Carpenter (2018) Carpenter, B.  Predator-Prey Population Dynamics: The Lotka-Volterra model in Stan.  https://mc-stan.org/users/documentation/case-studies/lotka-volterra-predator-prey.html#abstract, January 2018.
  * Finke & Thiery (2019) Finke, A. and Thiery, A. H.  On importance-weighted autoencoders.  _arXiv preprint arXiv:1907.10477_ , 2019.
  * Giordano et al. (2024) Giordano, R., Ingram, M., and Broderick, T.  Black box variational inference with a deterministic objective: Faster, more accurate, and even more black box.  _Journal of Machine Learning Research_ , 25(18):1–39, 2024.
  * Giordano et al. (2015) Giordano, R. J., Broderick, T., and Jordan, M. I.  Linear Response Methods for Accurate Covariance Estimates from Mean Field Variational Bayes.  In _Advances in Neural Information Processing Systems_ , volume 28. Curran Associates, Inc., 2015.
  * Glynn (1990) Glynn, P. W.  Likelihood ratio gradient estimation for stochastic systems.  _Communications of the ACM_ , 33(10):75–84, October 1990.  ISSN 0001-0782.  doi: 10.1145/84537.84552.
  * Gordon et al. (1993) Gordon, N. J., Salmond, D. J., and Smith, A. F. M.  Novel approach to nonlinear/non-Gaussian Bayesian state estimation.  _IEE Proceedings F (Radar and Signal Processing)_ , 140(2):107–113, April 1993.  ISSN 2053-9045.  doi: 10.1049/ip-f-2.1993.0015.
  * Hoffman & Gelman (2014) Hoffman, M. D. and Gelman, A.  The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo.  _Journal of Machine Learning Research_ , 15(47):1593–1623, 2014.  ISSN 1533-7928.
  * Kim et al. (2015) Kim, S., Pasupathy, R., and Henderson, S. G.  A Guide to Sample Average Approximation.  _Handbook of simulation optimization_ , pp. 207–243, 2015.
  * Kingma & Welling (2013) Kingma, D. P. and Welling, M.  Auto-Encoding Variational Bayes.  _ICLR 2014 conference submission_ , December 2013.
  * Le et al. (2018) Le, T. A., Kosiorek, A. R., Siddharth, N., Teh, Y. W., and Wood, F.  Revisiting Reweighted Wake-Sleep.  _International Conference on Learning Representations_ , September 2018.
  * Liu & Owen (2021) Liu, S. and Owen, A. B.  Quasi-Monte Carlo Quasi-Newton in Variational Bayes.  _Journal of Machine Learning Research_ , 22(243):1–23, 2021.  ISSN 1533-7928.
  * Lotka (1925) Lotka, A. J.  Principles of physical biology.  1925\.
  * Naesseth et al. (2020) Naesseth, C., Lindsten, F., and Blei, D.  Markovian Score Climbing: Variational Inference with KL(p || q).  In _Advances in Neural Information Processing Systems_ , volume 33, pp. 15499–15510. Curran Associates, Inc., 2020.
  * Naesseth et al. (2019) Naesseth, C. A., Lindsten, F., and Thomas B., S.  Elements of Sequential Monte Carlo.  _Foundations and Trends in Machine Learning_ , 2019.
  * Paisley et al. (2012) Paisley, J., Blei, D., and Jordan, M.  Variational bayesian inference with stochastic search.  In _Proceedings of the 37th International Conference on Machine Learning_ , 2012.
  * Pickover (1995) Pickover, C. A.  _The Pattern Book: Fractals, Art And Nature_.  World Scientific, June 1995.  ISBN 978-981-4504-03-4.
  * Rainforth et al. (2016) Rainforth, T., Le, T. A., van de Meent, J.-W., Osborne, M. A., and Wood, F.  Bayesian optimization for probabilistic programs.  _Advances in Neural Information Processing Systems_ , 29, 2016.
  * Ranganath et al. (2014) Ranganath, R., Gerrish, S., and Blei, D.  Black Box Variational Inference.  In _Proceedings of the Seventeenth International Conference on Artificial Intelligence and Statistics_ , pp. 814–822. PMLR, April 2014.
  * Rezende et al. (2014) Rezende, D. J., Mohamed, S., and Wierstra, D.  Stochastic Backpropagation and Approximate Inference in Deep Generative Models.  In _Proceedings of the 31st International Conference on Machine Learning_ , pp. 1278–1286. PMLR, June 2014.
  * Stites et al. (2021) Stites, S., Zimmermann, H., Wu, H., Sennesh, E., and van de Meent, J.-W.  Learning proposals for probabilistic programs with inference combinators.  In _Proceedings of the Thirty-Seventh Conference on Uncertainty in Artificial Intelligence_ , pp. 1056–1066. PMLR, December 2021.
  * Tucker et al. (2018) Tucker, G., Lawson, D., Gu, S., and Maddison, C. J.  Doubly Reparameterized Gradient Estimators for Monte Carlo Objectives.  _arXiv preprint arXiv:1810.04152_ , 2018.
  * van de Meent et al. (2018) van de Meent, J.-W., Paige, B., Yang, H., and Wood, F.  An Introduction to Probabilistic Programming.  _arXiv preprint arXiv:1809.10756_ , 2018.
  * Volterra (1927) Volterra, V.  Fluctuations in the Abundance of a Species considered Mathematically.  _Nature_ , 119(2983):12–13, January 1927.  ISSN 1476-4687.  doi: 10.1038/119012b0.
  * Wingate & Weber (2013) Wingate, D. and Weber, T.  Automated variational inference in probabilistic programming.  _arXiv preprint arXiv:1301.1299_ , 2013.
  * Wu et al. (2020) Wu, H., Zimmermann, H., Sennesh, E., Le, T. A., and Meent, J.-W. V. D.  Amortized Population Gibbs Samplers with Neural Sufficient Statistics.  In _Proceedings of the 37th International Conference on Machine Learning_ , pp. 10421–10431. PMLR, November 2020.
  * Zhang et al. (2022) Zhang, L., Carpenter, B., Gelman, A., and Vehtari, A.  Pathfinder: Parallel quasi-Newton variational inference.  _Journal of Machine Learning Research_ , 23(306):1–49, 2022.  ISSN 1533-7928.
  * Zhang et al. (2023) Zhang, L., Blei, D., and Naesseth, C. A.  Transport Score Climbing: Variational Inference Using Forward KL and Adaptive Neural Transport.  _Transactions on Machine Learning Research_ , May 2023.  ISSN 2835-8856.
  * Zimmermann et al. (2021) Zimmermann, H., Wu, H., Esmaeili, B., and van de Meent, J.-W.  Nested Variational Inference.  In _Advances in Neural Information Processing Systems_ , volume 34, pp. 20423–20435. Curran Associates, Inc., 2021.