b=0. We are again using the RANDOM_INDEX to retrieve a random subsample of
Gaia point sources.

In the lower part, we plot a 2d-histogram, that is, a density plot for the
objects we have retrieved. As you can see, we access the list of all retrieved
galactic longitudes l by calling up `resultset[’l’]`, while we get the list of
latitudes b by calling up `resultset[’b’]`. This is the same as for other
tables, where we can retrieve a column by using the column name as an index.
The combination `(resultset[’l’]+180.0) % 360` is used to shift the galactic
longitude by 180 degrees. That way, the galactic center is not around l=0,
that is, at the left and right margin of the image, but in the center.

Figure 96: Plot of Gaia stars, retrieved with PyVO

The resulting image gaia-plot.pdf can be seen in Fig. 96. Once more, the Milky
Way, the Large Magellanic Cloud and the Small Magellanic Cloud are clearly
visible.

Just like the TOPCAT TAP access, there is a whole world of different services
to explore via this interface, from stars to exoplanets323232E.g. via the API
for the highly useful Extrasolar Planet Encyclopaedia website,
[http://exoplanet.eu/API/] to galaxies.

Another way of accessing various catalogues is the AstroQuery package,333333A
general description can be found in Ginsburg et al. 2019, “astroquery: An
Astronomical Web-Querying Package in Python,” https://arXiv:1901.04520 a
collection of modules to retrieve data from various astronomical catalogues
via the web. Information about the specific catalogues and the specifics of
accessing them can be found in the documentation at
[https://astroquery.readthedocs.io].

As a simple example, here is the same Gaia query we had already executed with
`pyvo`, but this time using `astroquery`:

⬇

from astroquery.gaia import Gaia

import matplotlib.pyplot as plt

job = Gaia.launch_job("""

SELECT TOP 1000000

l,b

FROM gaiadr2.gaia_source

ORDER BY RANDOM_INDEX

""")

resultset = job.get_results()

plt.clf()

plt.hist2d((resultset[’l’]+180.0) \%

360,resultset[’b’], bins=(200, 200),

cmap=plt.cm.jet)

plt.savefig(’gaia-plot2.pdf’,bbox_inches=’tight’)

The resulting figure looks exactly the same as Fig. 96.

##  10 Astronomical image manipulation with Python

All professional astronomical images are in the FITS format, file extension
.fits or .fit, which stands for “Flexible Image Transport System”. Where your
ordinary JPG file gives you 256 steps between darkest and brightest (8 bits,
separately for each color RGB), FITS gives you, by default, 65536, or 16 bits.
That’s a lot of contrast. In section 3, we learned how to use application
software to take a look at, and perform some measurements in, images in FITS
format. Let’s see how we can display and analyse FITS images in Python.

###  10.1 FITS files and python

In order to open a FITS image file in Python, we once more use the Astropy
sub-module `astropy.io`. Just as we did with FITS tables in section 9.1, we
first load the HDU list. For convenience, we use one of the Hubble Space
Telescope images we had already downloaded in section 3.1, namely the file
hst_05773_05_wfpc2_f502n_wf_drz.fits. We load the file’s HDU list like this:

⬇

from astropy.io import fits

fitsURL=’hst_05773_05_wfpc2_f502n_wf_drz.fits’

hdulist=fits.open(fitsURL)

If we call the function info on the hdulist, we are shown the different
header/data units of this fits file:

⬇

In: hdulist.info()

Out:

Filename: hst_05773_05_wfpc2_f502n_wf_drz.fits

No. Name Type Cards Dimensions Format

0 PRIMARY PrimaryHDU 509 ()

1 SCI ImageHDU 103 (2150, 2150) float32

2 WHT ImageHDU 124 (2150, 2150) float32

3 CTX ImageHDU 123 (2150, 2150) int32

We are interested in the science image SCI, the second HDU, hence the one with
index 1. The header component of the HDU variable will give you access to the
header, as shown in the following listing.

⬇

In: hdulist[1].header

Out: Out[279]:

XTENSION= ’IMAGE ’ / Image extension

BITPIX = -32 / array data type

NAXIS = 2 / number of array dimensions

NAXIS1 = 2150

NAXIS2 = 2150

PCOUNT = 0 / number of parameters

GCOUNT = 1 / number of groups

CRVAL1 = 274.7173822566667 / right ascension of reference pixel (deg)

CRVAL2 = -13.83106772194444 / declination of reference pixel (deg)

CRPIX1 = 1075.0 / x-coordinate of reference pixel

CRPIX2 = 1075.0 / y-coordinate of reference pixel

Once again, I am showing only a small part of the output here. You can access
specific header information by plugging in the respective keyword. This here,
for instance, will give you the number of pixels along the first image
dimension:

⬇

In: hdulist[1].header[’NAXIS1’]

Out: 2150

An interesting piece of information is the exposure time in seconds, which is
contained in the header of the primary HDU:

⬇

In: hdulist[0].header[’EXPTIME’]

Out: 2200.0

The data component of the science HDU will give you the image data. Let’s
define

⬇

imdata = hdulist[1].data

###  10.2 Displaying (showing) an image

Once we have the image data, we can display it, using matplotlib, using the
`imshow` function, as follows:

⬇

plt.clf()

plt.axes().set_aspect(’equal’)

plt.imshow(imdata,cmap=’gray’)

The `set_aspect(’equal’)` tells matplotlib that both x and y axes should have
the same scale, as behoves two spatial directions spanning a two-dimensional
plane. The cmap option tells imshow which colormap to use, in this case
grayscale; many others are possible. The result is at first rather dark, as
you can see in Fig. 97.

Figure 97: This is not the image of a black hole

Just like with DS9, we somehow need to map the high contrast of the FITS image
to our more modestly contrasted version. We can use the `clim` option to map a
more restricted range of values to our image. Let’s look at the 1st and 99th
percentile values of the image data (that is, the brightness value below which
the darkest 1% of the pixels fall, and the brightness value above which 1% of
the pixels fall). We can access those descriptive numbers by typing

⬇

In: np.percentile(imdata,1)

Out: -0.028674498219043016

In: np.percentile(imdata,99)

Out: 0.068247721269726738

With the results, we can plot the image setting more suitable brightness
limits (more generally, limits for our color map) like this:

⬇

plt.clf()

plt.axes().set_aspect(’equal’)

plt.imshow(imdata,cmap=’gray’,clim=(-0.03,0.088))

The result can be seen in Fig. 98.

Figure 98: HST image of M 16, with adapted colourmap

###  10.3 Pixelwise operations

The image data we have put on display is an array. Specifically, you can see
this as

⬇

In: type(imdata)

Out: numpy.ndarray

and

⬇

In: imdata.shape

Out: (2150, 2150)

which shows you that the image is a Numpy array $2150\times 2150$. Just like
with any other array, we can retrieve pixel-wise information. For instance,
the brightness value of the pixel $x=1200,y=1400$ can be retrieved by
accessing the respective element of the `imdata` array, namely as

⬇

In: imdata[1400][1200]

Out: 0.030607721

Note the counter-intuitive order — the first square parentheses contain the y
value, the second one the x value! We can change pixel values in the same way,
by assigning a new value to a specific `imdata[y][x]`. Being able to read out
and manipulate pixels individually gives us substantial power to analyse the
image. We can use all the tools the previous sections have put at our
disposal, comparing pixel values, or summing them up, or performing fits, or
trying to identify point sources, or doing much more complicated types of
analysis.

Let us perform at least some basic operations on the SDSS data file
frame-g-007923-5-0307.fits that we had downloaded and analysed in sections 3.6
and 3.7. It’s usually good idea to take an overall look at the image one
intends to analyse, which we have learned to do in this way:

⬇

hdulistS=fits.open(’frame-g-007923-5-0307.fits’)

imdataS=hdulistS[0].data

lowerOne = np.percentile(imdataS,1)

upperOne = np.percentile(imdataS,99)

plt.clf()

plt.imshow(imdataS,cmap=’gray’,

clim=(lowerOne,upperOne))

plt.savefig(’sdss-py.pdf’,bbox_inches=’tight’)

The result is shown in Fig. 99.

Figure 99: SDSS DR9 frame displayed with imshow

Next, let us see if we can re-do the simple aperture photometry measurements.
The stars we had compared in section 3.7 were located at $X=1819,\;Y=1215$ and
$X=1123,\;Y=32$, respectively. We can use `xlim` and `ylim` to zoom in onto
those locations; for instance, amending the plot with

⬇

centerX = 1819

centerY = 1215

plt.xlim(centerX-100,centerX+100)

plt.ylim(centerY-75,centerY+75)

we obtain a detailed view all around the first of the stars, 200 pixels in
width, 150 in height. Next, we will visualize both the position of the star we
are interested in, and the surrounding regions we will use in our aperture
photometry measurements. To this end, we will add a circle to our diagram,
centered on the star. The proper way of doing this is

⬇

thisCircle = plt.Circle((centerX, centerY), 10, color=’r’,fill=False,lw=2)

plt.gca().add_artist(thisCircle)

This has two parts. In the first part, we define a circle centered on
$(1819,1215)$ with radius 10, the colour red, which is not filled, just an
outline with linewidth 2. In matplotlib, that circle is a so-called “patch”, a
two-dimensional artist object. Artist objects are matplotlib’s way of drawing
specific shapes.

In the second line, `add_artist` adds this circle to the axes object of our
figure. We repeat those commands with radius 20; the outer circles mark the
area we will use for determining the background brightness. Note that in
repeating those commands, we need not even choose a different name for the
variable `thisCircle`.

Once we have pushed the object onto our diagram using `add_artist`, the
variable has done what it was meant to be, and we can use it to define, and
draw, another circle. The result can be see in Fig. 100: Our chosen location
is now surrounded by two circles, drawn in red.

Figure 100: Detailed view of the region around the first star, with an inner
circle and outer annulus marked

That is a cosmetic marking only, though. Next, we can do aperture photometry
as in section 3.7. First, we determine the sum of pixel brightness values in
the outer circle, as well as the area of that outer circle in pixels. We do so
by going over all pixels in a square with side length `2*radius`, and
including only the values of pixels whose distance is smaller than the value
of `radius`, in other words: those pixels that are in the circle we are
interested in:

⬇

radius=20

photCollector=np.array([])

for ii in range(centerX-radius, centerX+radius):

for jj in range(centerY-radius,centerY+radius):

distance = np.sqrt((ii-centerX)**2 + (jj-centerY)**2 )

if distance < radius:

photCollector= np.append(photCollector, imdataS[jj][ii])

C1 = np.sum(photCollector)

A1 = len(photCollector)

As a result, I obtain C1=107.65 and an area of A1 = 1245 pixels. For the inner
circle, with radius 10 pixel, I get C2 = 101.12 with an area A2 = 305 pixels.
Using the brightness formula (6),

$l=C_{2}-\frac{C_{1}-C_{2}}{A_{1}-A_{2}}\cdot A_{2},$ (18)

we subtract from the sum of brightness values $C_{2}$ the sum of background
brightness values (as estimated by the average of the brightness values in the
annulus bounded by the inner and the outer circle). For the first star, I get
the brightness $l_{A}=12.60$. I repeat the procedure for the second star, the
one at $X=1123,\;Y=32$, and obtain C1 = 107.47 with an area of A1 = 2809
pixels, and C2 = 107.65 for an area of A1= 1245 pixels. Using once more the
brightness formula (18), the star’s background-subtracted brightness is
$l_{B}=107.80$.

As I had already argued in section 3.7, given that the collecting area and the
exposure time are the same in both cases, the ratio of our values $l_{B}$ and
$l_{A}$ should be the same as the ratio of the intensities of the two stars;
inserting this into equation (1) for the astronomical magnitudes, we find

$m_{A}-m_{B}=-2.5\cdot\log\left(\frac{l_{A}}{l_{B}}\right)=-2.5\cdot\log\left(\frac{12.6}{107.8}\right)=2.33,$
(19)

again not far from the difference in the stars’ catalog g magnitudes (gmag) of
2.13.

This is a single example of how we can use our ability to access an
astronomical image pixel by pixel unlocks a treasure trove of analytical
possibilities.343434A more detailed guide for CCD image processing using
Python is the “CCD Data Reduction Guide” by Matt Craig and Lauren Chambers,
available at [https://mwcraig.github.io/ccd-as-book/00-00-Preface.html]. If
you really want to do aperture photometry on a larger scale, you would very
likely not want to code your analysis from scratch — you would use, for
instance, the Astropy-affiliated package Photutils
[https://photutils.readthedocs.io/] and its higher-level functions for
identifying point sources and measuring their brightness. But the main point
of this section was much more general: Once we have loaded FITS data, the full
force of coding can be brought into play. Our aperture photometry example was
comparatively simple. More complex operations are possible: We could fit a
brightness profile to an extended source, or try to reconstruct the point-
spread function for a point source. We could even try to have our script
identify objects such as stars automatically, by checking brightness
differences. The possibilities are virtually endless, definitely so in that
astronomers will continue to bring to bear new analysis tools as they are
being developed. Recent applications of machine learning in astronomy are a
pertinent example.

##  11 A simple simulation

So far, we have almost exclusively dealt directly with observational data,
either images/spectra or higher-level table data. This is not the only kind of
data that astronomers work with — there is a whole other branch of research
focused on simulated data. Such simulated data can play different roles. Some
are meant to provide a point of comparison for observations. For instance, if
you really want to understand all the details of the spectral line shape in a
stellar spectrum, you will need access to the results of simulations for the
genesis of those lines. The newest versions take into account complex effects
such as the absence of local thermodynamic equilibrium (a common
simplification). Compare the simulated spectra for different stellar parameter
values (such as chemical abundances) with your observations, and you can
deduce the properties of your star.353535At the time of this writing, a sample
set of tools from the group of Maria Bergemann at the Max Planck Institute for
Astronomy can be found on [http://nlte.mpia.de].

Figure 101: Details of a massive disk galaxy at redshift $z=1$, simulated with
the TNG50 simulation. Image: TNG collaboration

Other simulations are more ambitious in scope. The various runs of the
IllustrisTNG363636More information can be found at [http://www.tng-
project.org].  simulation, for instance, follow a cubic region within the
cosmos from shortly after the Big Bang to the present. The simulation runs are
of different degrees of coarseness, simulating either a very highly resolved
cube with a sidelength of 50 million parsecs, a less well resolved cube of 100
Mpc or a larger, but still less resolved cube 300 Mpc a side. The simulation
is mostly based on “particles” representing dark matter, stars, and gas, and
simulates numerous physical processes including the influence of magnetic
fields, the formation and evolution of stars in galaxies, the production of
heavier elements in stars, and interactions with the supermassive black holes
in the centres of galaxies. All in all, the IllustrisTNG simulations follow
the evolution of the universe over the full 13.8 billion years, including
diverse length scales from those on which the cosmos is, on average,
homogeneous, down to those of the sub-structure of galaxies, cf. Fig. 101.

As I have written before, the details of such simulations call for (at least!)
a whole set of lecture notes of their own. Most of that knowledge is far
beyond the scope of the present text. But I firmly believe that an overview of
working with astronomical data needs to include at least a brief introduction
to creating your own data based on models and on the laws of physics — that,
in my mind, is an integral part of working-with-astronomical-data literacy.
With this in mind, here is a simple example that illustrates at least some
basic techniques, and also some elementary pitfalls, of numerical simulations.
The physics in this case is simple, classical mechanics.

###  11.1 Step-by-step numerical integration: Euler method

We revert to what is probably the most basic (and arguably most useful!)
elementary system in all of physics: The harmonic oscillator. The basic set-up
is as shown in Fig. 102: A particle with mass $m$, which can move only in the
(horizontal) x direction, is fixed to the wall with a spring. If the particle
is displaced from its rest position at $x=0$, the spring exerts a force
following Hooke’s law, $F_{x}=-k\cdot x,$ with $k$ the spring constant.

Figure 102: Harmonic oscillator: particle on a spring

This system has the advantage that its equation of motion,

$m\ddot{x}=-k\cdot x,$ (20)

linking the x acceleration $\ddot{x}$ and the force using Newton’s second law,
is readily solved analytically, that is, in terms of a simple mathematical
function. The solution is

$x(t)=A\cdot\sin(\omega t),$ (21)

where the angular frequency $\omega$ is linked to the system’s oscillation
period $T$ by the standard definition

$\omega=\frac{2\pi}{T},$ (22)

and the equations of motion in this particular case demand

$\omega=\sqrt{\frac{k}{m}}.$ (23)

Differentiating the orbit equation (21) once with respect to time, we have

$\dot{x}=\omega\cdot A\cdot\cos(\omega t),$ (24)

and differentiating with respect to time once more,

$\ddot{x}=-\omega^{2}\cdot A\cdot\sin(\omega t)=-\frac{k}{m}\cdot x,$ (25)

which shows that our solution (21) indeed satisfies the equation of motion
(20). So far, so elementary. But now, pretend that we do not know of this
simple solution. How can we simulate the system, in other words: find a
solution not analytically, but numerically?

The basic idea is that, if we look at very small time interval, all of the
changes during such an interval will be approximately linear. This, is, after
all, the definition of a derivative: over an infinitesimally small interval
$\mathrm{d}t$, the change of the function $x(t)$, namely $\mathrm{d}x$, is
given by

$\mathrm{d}x=\dot{x}\cdot\mathrm{d}t.$ (26)

Replace the infinitesimally small interval $\mathrm{d}t$ by a finite small
interval $\Delta t$, and what was an equality in (26) becomes an
approximation, whose quality depends on the magnitude of $\Delta t$: the
smaller $\Delta t$, the better the approximation. Thus, if we know the
particle’s x position at one time $t$, we can estimate its position at a
slightly later time $t+\Delta t$ as

$x(t+\Delta t)=x(t)+v(t)\cdot\Delta t,$ (27)

where $v(t)$ is the particle’s velocity in x direction at the time $t$. What
this formula does not encode, of course, is how $v(t)$ changes over time. But
for the change of $v$, we can write down a similar equation. The rate of
change of the velocity, after all, is the acceleration, which by Newton’s
second law $F=m\ddot{x}$ is linked to the force acting on the particle. Thus,
to obtain the velocity at some time $t+\Delta t$, we can use the approximation

$v(t+\Delta t)=v(t)+\ddot{x}(t)\cdot\Delta t=v(t)+\frac{F(t)}{m}\cdot\Delta
t.$ (28)

Incidentally, in writing down the approximation equations (27) and (28), we
have applied a technique that can be used much more generally, when dealing
with higher-order differential equations: We have transformed a single second-
order differential equation (for $\ddot{x}$) into a system of (two) first
order equations, each of which we have solved approximately, in going from $t$
to $t+\Delta t$.

Now, we can discretize the whole problem: We consider time steps $t_{i}$, with
$i=1,\ldots,N$, and evaluate the position and the velocity of our particle at
each step. We choose the times $t_{i}$ equidistant, with $t_{i+1}-t_{i}=\Delta
t$ for all $i$, for some fixed, small $\Delta t$.

(What is small, and more specifically: what is sufficiently small? That
depends on the problem’s characteristic time scale, and requires physical
thought. With hindsight (or by looking at the analytical solution), we know
that we are dealing periodic, oscillatory motion, so whatever interval $\Delta
t$ we choose had better be much smaller than the system’s natural period $T$.
If we cannot find a physical time scale to ascribe to the system, we might
need to fall back on experimentation. If the solution changes significantly
when we repeat the simulation with ever smaller $\Delta t$, that is an
indication we have not reached the proper resolution yet. If, on the other
hand, the solution remains pretty much the same when we replace, say, $\Delta
t$ by $\Delta t/2$, that is an indication that we have reached the regime
where the finite size of $\Delta t$ does not exert significant influence on
our result any more.)

Differential equations do not, on their own, completely determine what is
happening. It is necessary to specify initial conditions in order to define a
unique solution. As an example, consider what happens when you throw a ball
vertically upward from the Earth’s surface. The differential equations tell
you how Earth’s gravity will accelerate the ball. But those equations alone
are not sufficient to tell you what happens. For a prediction, you will need
to specify both the ball’s initial position and it’s initial velocity. Both
those initial values are crucial in determining how the ball will move, and in
particular whether it will fall back to Earth or keep going forever (when its
velocity is larger than the position-dependent escape velocity).

In this case, let us choose an initial position $x_{0}$ and initial speed
$v_{0}$ for our particle. Let $x_{i}$ be the object’s position at time
$t_{i}$, $v_{i}$ its velocity in x direction at that time, $a_{i}$ its
acceleration and $F_{i}$ the force acting on it at the time. The simple, step-
wise evolution equations we have derived are

$\displaystyle v_{i+1}$ $\displaystyle=$ $\displaystyle v_{i}+a_{i}\cdot\Delta
t=v_{i}+\frac{1}{m}F_{i}\cdot\Delta t$ (29) $\displaystyle x_{i+1}$
$\displaystyle=$ $\displaystyle x_{i}+v_{i}\cdot\Delta t.$ (30)

The process of following the evolution step by step is called numerical
integration, and the simple algorithm we have given for going from one step to
the next is called Euler’s method. In our case, with the Hooke force (20), the
force depends only on the position, so the velocity equation reads

$v_{i+1}=v_{i}-\frac{k}{m}x_{i}\cdot\Delta t.$ (31)

This is readily implemented in Python, for instance in the following way,
using the loop function. (The array `tCollector` is only defined in
preparation of plotting positions against time later on, and not used in the
numerical integration itself.)

⬇

k=0.5

m=1.0

numberOfSteps = 30000

DeltaT = 0.001

finalT =numberOfSteps*DeltaT

tCollector=np.linspace(0,finalT,numberOfSteps+1)

# Initial conditions:

xCollector=[1.0]

vCollector=[0]

for ii in range(numberOfSteps):

xNew = xCollector[-1] + DeltaT*vCollector[-1]

vNew = vCollector[-1] + DeltaT*(-k/m*xCollector[-1])

xCollector.append(xNew)

vCollector.append(vNew)

Plotting the result, as in Fig. 103, shows that

Figure 103: Solution of the harmonic oscillator equation with initial position
1.0 and initial speed 0.0

we indeed obtain the proper sine shape, shifted so as to form a cosine
function. (Why cosine instead of sine? Because, by choosing our initial
condition to be $v=0$, we start our evolution at the maximum x value.)

###  11.2 Numerical errors

If I just were to plot the analytical solution and the numerical solution in
the same diagram, the curves would be overlaid so closely that no difference
would be visible upon direct inspection. As an alternative, Fig. 104 shows the
difference between the analytical and the numerical solution for our harmonic
oscillator solution, at each time.

Figure 104: Difference between our numerical solution and the analytical
solution

There is good news and there is bad news. The good news is that the
differences between the true solution and our simulation are very small,
namely 0.6% at most, and much better for most of the time. The bad news is
that the differences are getting larger over time. If that trend continues,
and we let our simulation run much further, its deviations from the true
solution will become so large as to be noticeable. Our numerical simulation is
unstable in this sense.

We can understand why that is. Imagine that, in truth, our particle reaches
our maximum x value at time $t_{i}$. From time $t_{i-1}$ to time $t_{i}$, we
change the x value by $\Delta x=+v_{i-1}\cdot\Delta t$. But $v_{i-1}$ is the
(non-zero) velocity at the beginning of the interval. Over the time $\Delta
t$, that velocity will change to zero (since velocity zero is what defines the
maximum value, the turning point). Thus, the average velocity over that
interval will be smaller than the initial velocity for that interval. We are
overestimating the amount $\Delta x$ that x increases during that time. We
pretend the object has flown with unchanged initial velocity $v_{i-1}$,
whereas, in reality, it has slowed down. That means our turning-point x will
be just a little further out, x a little larger, than in reality.

There is nothing in our simulation to compensate for that larger error. The
system knows nothing of its past. That slightly larger x, that slightly larger
amplitude at the turning point will be carried along for the rest of our
simulation. You can think about it in the following way: Because of the
slightly over-large x when the turning point comes, our system, from then on,
has a slightly larger potential energy than it should have. There is nothing
to compensate; once the extra energy has seeped into the system through the
numerical error, it will remain in the system.

Even worse: At every turning point, the same argument applies, so at every
turning point, the error of having a slightly larger amplitude will increase.
The errors systematically add up. That is what makes the system unstable.
Errors do not compensate each other; errors just add up over time, increasing
the overall deviation from the true time evolution.

Note that the systematic errors, and the instability, are a property of the
algorithm we have used to simulate the time evolution. There are several
different algorithms that, in the limit of infinitesimal $\Delta t$, all
amount to the same integration procedure, but which differ in their stability
properties and in how accurate they are for finite $\Delta t$. We will look at
one of them in the next section.

###  11.3 Velocity Verlet algorithm

One example for a better-behaved numerical integration scheme xsis the
velocity Verlet algorithm, which introduces “half-step” velocities, in order
to mitigate problems like those I described in section 11.2, as follows:

$\displaystyle v_{i+1/2}$ $\displaystyle=$ $\displaystyle
v_{i}+\frac{1}{2}\Delta t\cdot a_{i}$ (32) $\displaystyle x_{i+1}$
$\displaystyle=$ $\displaystyle x_{i}+v_{i+1/2}\cdot\Delta t$ (33)
$\displaystyle v_{i+1}$ $\displaystyle=$ $\displaystyle
v_{i+1/2}+\frac{1}{2}\Delta t\cdot a_{i+1},$ (34)

where the $a_{i}$ is again calculated from the force acting in time step $i$.
In this algorithm, we avoid always using the initial values of rates-of-change
at each time step, leading to systematic errors; instead, the calculation of
$x_{i+1}$ uses the intermediate velocity, while the transition from $v_{i}$ to
$v_{i+1}$ proceeds in two steps, one using the acceleration value at the
beginning, the second the one at the end of the interval $\Delta t$. The
implementation is, again, fairly straightforward:

⬇

k=0.5

m=1.0

numberOfSteps = 30000

DeltaT = 0.001

finalT =numberOfSteps*DeltaT

tCollector=np.linspace(0,finalT,numberOfSteps+1)

# Initial conditions:

xCollectorVV=[1.0]

vCollectorVV=[0]

for ii in range(numberOfSteps):

vHalf = vCollectorVV[-1] + 0.5*DeltaT*(-k/m*xCollectorVV[-1])

xNew = xCollectorVV[-1] + DeltaT*vHalf

vNew = vHalf + 0.5*DeltaT*(-k/m*xNew)

xCollectorVV.append(xNew)

vCollectorVV.append(vNew)

The comparison between the results can be seen in Fig. 105.

Figure 105: Difference between our numerical solution and the analytical
solution, once for the Euler algorithm, once for the velocity Verlet algorithm
(the latter here scaled up artificially by a factor $10^{4}$)

The velocity Verlet algorithm, too, appears to be unstable in the long-term,
with the error increasing over time. But the simple of expedient of adding the
half-step velocity has greatly improved the accuracy. After all, note that, in
this diagram, I have scaled up the difference for the velocity Verlet
algorithm by a whopping factor of $10^{4}$ to make it visible in comparison
with the Euler deviations!

This simple example shows the importance of implementing the evolution
algorithm carefully; finding the best ways of doing this is a science (and
possibly an art) of its own. If you decide to explore this further, you might
want to look at the Runge-Kutta family of iterative methods next.

###  11.4 A simple two-dimensional simulation

So far, our simple simulation was one-dimensional: motion in the x direction
over time. Let us choose a two-dimensional scenario next, and one that is very
important in astronomy: The motion of a test particle around a central mass
under the influence of the central mass’s (Newtonian) gravity, which provides
a good model for the orbit of a (not too massive) planet around a star.

Let us put the central mass into the origin of our coordinate system. We treat
the x and y components of the motion separately. Since we have two independent
directors, we will need to treat the force as a vector as well, separating its
x from its y component. Fig. 106 shows the geometry of the situation.

Figure 106: Planet with mass $m$ orbiting in the xy plane under the
gravitational influence of a central mass $M$ in the origin

Crucially, the triangle describing the x-y-position of the planet and the
triangle describing the decomposition of the gravitational force into x and y
component are similar, in the geometric sense: as the figure shows, they have
the same set of three angles. Similar triangles can only differ by an overall
length scale. In particular, ratios of the corresponding sides of such
triangles are the same. If we abbreviate the distance of the planet from the
origin (and thus from the central mass) as

$r=\sqrt{x^{2}+y^{2}}$ (35)

(having applied the Pythagorean theorem), then we have

$\frac{F_{x}}{F}=\frac{x}{r}$ (36)

and

$\frac{F_{y}}{F}=\frac{y}{r}.$ (37)

For the Newtonian gravitational force, we have

$F=-\frac{GMm}{r^{2}}.$ (38)

Using $F_{x}=ma_{x}$ and $F_{y}=ma_{y}$ for the link between the force and the
acceleration components, we have

$\displaystyle v_{x,i+1/2}$ $\displaystyle=$ $\displaystyle
v_{x,i}+\frac{1}{2}\Delta t\cdot a_{x,i}$ (39) $\displaystyle x_{i+1}$
$\displaystyle=$ $\displaystyle x_{i}+v_{x,i+1/2}\cdot\Delta t$ (40)
$\displaystyle v_{x,i+1}$ $\displaystyle=$ $\displaystyle
v_{x,i+1/2}+\frac{1}{2}\Delta t\cdot a_{x,i+1},$ (41)

and analogous equations linking the $v_{y}$ and $y$ values. Note that, via
$r$, the force at any given step, and thus the acceleration components, depend
on both coordinates.

Next, we need to choose suitable units. Let the central mass have one solar
mass, $M_{\odot}=2\cdot 10^{30}\;\mbox{kg}$. Going by Earth’s orbit, a
suitable unit of length is the astronomical unit (corresponding to the average
Earth-Sun distance), $1\;\mbox{au}=1.5\cdot 10^{11}\;\mbox{m}$. As our unit of
time, we choose the Julian year: 365.25 standard days, abbreviated as a for
the latin “annum” for year, related to the SI unit, the second, as
$1\;\mbox{a}=31~{}557~{}600\>\mbox{s}\sim\pi\cdot 10^{7}\>\mbox{s}$. The
corresponding unit for speed is related to the more usual one as

$1\>\frac{\mbox{km}}{\mbox{s}}=0.21\>\frac{\mbox{au}}{\mbox{a}}$ (42)

and the acceleration felt by the planet is given by

$F/m=-39.48\>\left(\frac{1\>\mbox{au}}{r}\right)^{2}\>\frac{\mbox{au}}{a^{2}}.$
(43)

The simulation code itself is listed here:

⬇

numberOfSteps = 30000

DeltaT = 0.0001

accFac = 39.48 # Corresponding to one solar mass, in au per square year

finalT=numberOfSteps*DeltaT

tCollector=np.linspace(0,finalT,numberOfSteps+1)

# Initial conditions:

xCollector=[1.5]

vxCollector=[0]

yCollector=[0.0]

vyCollector=[2.0]

for ii in range(numberOfSteps):

rNow = np.sqrt(xCollector[-1]**2+yCollector[-1]**2)

accNow = -accFac/rNow**2

accNowx= accNow*xCollector[-1]/rNow

accNowy= accNow*yCollector[-1]/rNow

vxHalf = vxCollector[-1] + 0.5*DeltaT*accNowx

vyHalf = vyCollector[-1] + 0.5*DeltaT*accNowy

xNew = xCollector[-1] + DeltaT*vxHalf

yNew = yCollector[-1] + DeltaT*vyHalf

rNew = np.sqrt(xNew**2 + yNew**2)

accNew = -accFac/rNew**2

accNewx = accNew*xNew/rNew

accNewy = accNew*yNew/rNew

vxNew = vxHalf + 0.5*DeltaT*accNewx

vyNew = vyHalf + 0.5*DeltaT*accNewy

xCollector.append(xNew)

yCollector.append(yNew)

vxCollector.append(vxNew)

vyCollector.append(vyNew)

For the half-step velocity, we calculate the accelerations in x and y
direction, starting with the magnitude of the acceleration, which follows
directly from Newton’s law. Then, we evolve the position one time step
further, re-calculate the acceleration for the new position, and use those to
update the x and y component of the velocity for the second half of the time
step.

To close this section, and put our simulation to the test, let us see if we
can recover Kepler’s three laws of planetary motion from our simulation. We
begin by plotting the shape of the orbit in Fig. 107.

Figure 107: Orbit of a planet around a central mass, simulated with the
velocity Verlet algorithm

That certainly looks like an elliptical orbit, with the central mass at the
origin serving as one of the focus points. But we should demonstrate this fact
in a more quantitative way. Taking the minimum and maximum x value on our
simulated orbit, we find that the left-most point of the orbit is at
$x=-0.12335823920305192$, while the rightmost point is, by construction, at
$x=1.5$ (since that is where we put the initial position of our particle, its
velocity pointing straight upwards).

Thus, our simulated orbit has a major half axis of $a=1.623\>\mbox{au}$. From
the minimum and maximum of the y coordinate value, we find that the minor half
axis is $b=0.860\>\mbox{au}$. For an ellipse, that would correspond to an
eccentricity of

$e=\sqrt{1-\frac{b^{2}}{a^{2}}}=0.848.$ (44)

The polar coordinate equation for an ellipse is

$r(\theta)=\frac{a(1-e^{2})}{1+e\cos\theta}.$ (45)

Inserting our parameter values $a$ and $e$, we can plot this reference ellipse
and compare with the simulated orbit. The resulting plot looks just like Fig.
107, with one ellipse directly on top of the other. At least qualitatively, we
have indeed confirmed Kepler’s first law: the orbit of a planet orbiting a
central mass is an ellipse, with the central mass in one of the focal points.
More quantitatively, we can compare analytical solution and simulation
directly. There are several possibilities for this; with the following piece
of code, I take each simulated point, calculate the position angle $\theta$
and distance $r$ from the focus point, and compute the absolute value of the
difference between the simulated value $r$ and the analytical value
$r(\theta)$ given by (45):

⬇

diffCollector=[]

for x,y in zip(xCollector,yCollector):

r=np.sqrt(x**2+y**2)

theta = np.arctan2(y,x)

anr = a*(1-e**2)/(1-e*np.cos(theta))

diffr = np.sqrt((r-anr)**2)

diffCollector.append(diffr)

The histogram of the values contained in `diffCollector` is shown in Fig. 108.

Figure 108: Histogram of the square difference between fitted ellipse and
simulated ellipse

It has a highly structured shape, and no doubt one could learn a lot about the
systematic errors involved in the simulation by understanding that cascade-
like structure. Such analysis is far beyond our current scope; for our
purposes, we note that this looks definitely non-random and, importantly, that
the largest deviation is 4/10 000 of an astronomical unit. Given that the
length scales of our orbit (half axis length, circumference) are on a scale of
1 astronomical unit, that deviation is fairly small, and our simulated planet
appears to have an elliptical orbit.

Next, for Kepler’s second law, which says that the connecting line between the
planet and the central mass sweeps out equal areas in equal time intervals.
Each of our time steps defines the same time interval, so if we calculate the
triangle swept out in each time step (whose three vertices are the planet’s
position at the beginning and at the end of the time step, and the location of
the central mass), we should always obtain the same area. For each such
triangle, we know the x and y coordinates of all three vertices, and thus can
calculate all the side lengths $a,b,c$ using the Pythagorean theorem. With
this information, we can use Heron’s formula to calculate the triangle’s area
as

$A=\sqrt{s(s-a)(s-b)(s-c)}$ (46)

where $s=(a+b+c)/2$ is the triangle’s semi-perimeter. (Alternatively, we can
use half of the cross product of the two position vectors, now viewed as
three-dimensional vectors, to obtain the same result.) The following bit of
code collects the relative deviation of each such triangle area from the mean
in an array `relativeDiff`:

⬇

areaCollector=np.array([])

for x1, x2, y1, y2 in zip(xCollector[1:], xCollector[:-1],yCollector[1:],
yCollector[:-1]):

a = np.sqrt(x1**2+y1**2)

b = np.sqrt(x2**2+y2**2)

c = np.sqrt((x1-x2)**2+(y1-y2)**2)

s = 0.5*(a+b+c)

A = np.sqrt(s*(s-a)*(s-b)*(s-c))

areaCollector = np.append(areaCollector,A)

averageArea=np.average(areaCollector)

relativeDiff = (areaCollector-averageArea)/averageArea

A histogram of the values in `relativeDiff` shows that all those areas, swept
out in the same time interval, are indeed very close to their average value,
as Fig. 109 shows.

Figure 109: Histogram of the relative difference from the average of the
triangle areas corresponding to each time step

The distribution shows that we do have a strong maximum at the average area
value, with small (a few parts in a trillion!) fluctuations to smaller and
larger values.

Next, to Kepler’s third law. In Kepler’s own version, this links the orbital
periods $T$ and major elliptical half axis $a$ of different planets orbiting
the same central mass, stating that the ratio $a^{3}/T^{2}$ is the same for
all of them. We did not simulate planets with different initial conditions
(although we could), and thus will check the more advanced form of Kepler’s
law found by Newton, which states that

$\frac{a^{3}}{T^{2}}=\frac{GM}{4\pi^{2}}$ (47)

(in the limit we have simulated, namely where the planetary mass $m$ is small
against the central mass $M$). We have already estimated $a$. Let us do the
same for $T$. To this end, I have plotted the y coordinate of our planet
against time in Fig. 110.

Figure 110: The y coordinate of our simulated planet over time

This is unsurprisingly periodical. In order to find out the period, we fold
the time evolution: we assume a value for the period $T$, and define the phase
$\phi$ in terms of this period as

$\phi=\frac{t}{T}\>\mbox{mod}\>1.$ (48)

All integer multiples of the period $T$ get mapped to 0, all times that can be
written as

$t=(n+f)\cdot T$ (49)

with integer $n$ and $0\leq f<1$ get mapped to $f$. The quick-and-dirty way of
finding the correct period $T$ is to vary the value by hand, and see whether
or not the resulting curves coincide. From Fig. 110, we can read off that the
period is somewhat less than one year. Fig. 111 shows the resulting plot for
$T=0.8\>\mbox{a}$.

Figure 111: Phase plot, with an assumed periodicity of $T=0.8\>\mbox{a}$, for
our simulated orbit

The fact that we still see separate, similar curves shows that we have not
quite hit on the correct period yet. By slowly decreasing $T$, and re-
plotting, I can bring those separate curves to coincidence. After a few dozen
tries, the final round with a microscopic look at the steepest curve region
via `plt.xlim(0.475,0.525)`, I arrive at $T=0.7313\;\mbox{a}$. Note that this
is very similar to how you determine the orbital of an exoplanet by folding
the light curve data (from the transit method) or Doppler shift data (in the
radial velocity method). We could think about automatising this, sorting the
phase values into different arrays indexed by the integer part of $t/T$, and
minimising the differences between those partial curves, but for our little
test, I consider the trial-and-error approach sufficient.

Thus, our final test is to see if this $T$ value indeed satisfies equation
(47). The right-hand side of that equation is

$1\;\frac{\mbox{au}^{3}}{\mbox{a}^{2}}$ (50)

by definition, since the Earth does have a semimajor axis of length
$1\;\mbox{au}$ (up to and including the sixth significant degree), and an
orbital period of 1 year. By our estimates, our simulation satisfies

$\left(\frac{a}{1\>\mbox{au}}\right)^{3}\cdot\left(\frac{1\>\mbox{a}}{T}\right)^{2}=0.99991.$
(51)

The two values coincide up to one part in 10 000. Our simulation reproduces
Kepler’s third law, as well.

This concludes our brief excursion into the realm of simulations. We have, as
ever, only scratched the surface, but along the way, you have encountered some
of the pitfalls and characteristics of numerical simulations, notably the
concept of numerical errors, the question of stability, and the importance of
choosing a suitable algorithm. What we did not encounter was the limitation
imposed by the available computing power. For our purposes, a simple script
which took a few seconds to run was sufficient. As simulations become more
complex, computing time increases, and can become a limiting factor. There are
various ways of addressing this problem, and pushing the limit. Powerful
hardware is one of them. Parallelising calculations, that is, having multiple
processors (or processor cores) tackle different parts of the problem
simultaneously, can be a powerful strategy for many (but not all) simulation
problems. Parallelisation is achieved by linking multiple processors in a way
that automatically distributes calculations among the available cores. GPUs,
graphics processing units, which commonly operate at lower frequencies but
consist of a greater number of cores, are frequently used for the purpose.
Projects like the IllustrisTNG simulation mentioned at the beginning of this
section push the envelope using an immense amount of computing power. The most
calculation-intensive run of that simulation, TNG300-1, involved 24 000 cores,
adding up to a total CPU time of nearly 34 million hours.373737The numbers
have been taken from the article by Dylan et al. 2017,
[http://adsabs.harvard.edu/abs/2018MNRAS.475..624N]

We have also only reproduced the simplest type of simulation, namely an N-body
simulation that follows the evolution of separate particles. Other
applications in astrophysics require simulations to deal with fluid mechanics,
dividing space into cells whose fluid content and properties (such as
temperature, entropy and flow properties) are then followed over time.383838A
step-by-step to writing a hydrodynamics code can be found in
[https://github.com/python-hydro/how_to_write_a_hydro_code], by Michael
Zingale.

As in the case of reducing and analysing observational data, there are
packages and libraries that will help you set up specific more advanced
simulations using Python.393939An example is GalPy by Jo Bovy, 2014,
[arXiv:1412.3451]. However, should you go further in this direction, you are
likely to reach a point where increased speed is so important that Python is
not fast enough for what you are about to do. More advanced simulations
typically use so-called compiled programming languages, such as C and its
different incarnations, where the code is first translated into a “compiled
code,” that is, into machine-readable instructions that can be carried out by
the computer directly. Compiled languages are typically faster than
interpreted languages like Python, where the translation into machine code
happens on the fly as the program is executed. An example is the Gadget-2 code
by Volker Springel et al.404040See [https://ascl.net/0003.001], as well as the
how-to for installing and running Gadget-2 by Nathan Goldbaum on astrobites,
at [https://astrobites.org/2011/04/02/installing-and-running-gadget-2/]. But
by the time you work at that level, you have moved far beyond the basic
introduction presented in this text.

##  12 Conclusion

Working with astronomical data requires a combination of skills. If you have
worked through this text, reproducing the analytical tasks on your own
computer, using DS9, TOPCAT and Python, you should have acquired basic
proficiency in a number of these skills, and you should now be familiar with
several key tools, both conceptual (histograms! diagrams!) and practical (how
to use different kinds of software to achieve specific purposes).

As we have seen on several occasions, we have only just scratched the surface.
But that is perfectly fine! If you dedicate your career to research, you will
continue to build on what you know. What you have learned here should allow
you to start that life-long process.

You will also have seen that the knowledge needed for data analysis falls into
different categories. It goes without saying that working with astronomical
data requires knowledge of physics, astronomy, and mathematics, specifically
statistics. When you begin working on a new topic, then likely as not in the
beginning, you will not fully understand what you are doing, and how the
different elements you are dealing with fit together properly. Your goal
should be to, eventually, reach the stage where you do understand what is
going on. Such knowledge is required if you want to understand your results,
but crucially also if you want to understand the limitations of your data, and
possible ways of improvement.

As you learn to use new tools, and more advanced tools, your newly acquired
knowledge might also open up new opportunities for exciting science. Applying
a tool that others did not think to apply could be the key step towards
finding a new result. In this respect, it pays to have an eye out for
neighbouring fields (or at least sub-fields). Could what they are doing in
their field help you achieve something new in yours?

Some of the knowledge you need for data analysis is a matter of convention.
How to perform a certain operation in TOPCAT, or plot a certain kind of
diagram, is nothing you can deduce logically from previous fundamental
knowledge (although previous experience will help you find the right answer).
If the software in question is well documented, it makes sense to familiarise
yourself with the basics; the more common strategy is to google what you are
looking for. Everybody does it. Reminding yourself how to “matplotlib equal
axis ratio” is just one search field away. Also, if you are working in an
institute, chances are there will be experienced people you can ask for help
with specific problems. Last but not least, if you need to accomplish a
complex task, and you happen to have a script which does something similar, it
makes sense to get to understand that other script, to try out variations on
what it does, and eventually to adapt it to your purposes.

Finally, there is meta-knowledge about working with astronomical data. One
piece of good advice is to always look at your data in simple form before
attempting complex operations on it. Make a few basic histograms and diagrams
to get a feeling for your data, look at an image, or make a quick plot of a
spectrum — incidentally, you might find out that your data is somehow
completely different from what you expected, and that is a good thing to know
early on!

A fairly universal truth you are likely to learn early on is how easy it is to
make mistakes. One reason is that, even while each of the elements of your
analysis might be straightforward, code can become fairly voluminous fairly
quickly once you combine all the different necessary steps. So what do you do
if, at the end of your analysis, your result is surprising, possibly wrong, or
even obviously wrong? Conversely, if things turn out as expected, how can you
be sure that this is not the result of several mistakes cancelling each other
out, or almost cancelling? Such questions would not be important if your goal
was just to create a visually pleasing astronomical image, for instance. But
when you are doing scientific research, you had better understand, and check,
every step of what you are doing; otherwise, you cannot be sure of your
result. Your analysis should contain as many cross-checks as you can come up
with, and have time to implement. After every step, think about what that step
is meant to do, and what you can do to check that the desired result has
indeed be achieved. If you write up your research, such cross-checks and
safeguards are likely to make up an important part of your description of what
you did.

Also, make sure your code is comprehensible. The most important element of
this is adding descriptive comments. A person reading your code should be able
to follow what you are doing step by step, guided by your code, but also by
your description in the comments. Remember that the person in question could
be you in a few years — it is amazing how incomprehensible ill-documented code
can become once you return to it a few months or years after having written
it! Meaningful variable names can help with comprehensibility, too.

Last but not least documenting your code is a matter of scientific
accountability. Science should be reproducible. Your research publications
should tell your colleagues what you have done, so they can build on your
results, but also critically examine what you have done. In an ideal world,
every scientific article would be accompanied by data files and script files;
running the script on the data, you should be able to reproduce the article’s
results on your own computer.414141Cf. Weiner et al. 2009,
[https://arxiv.org/abs/0903.3971]. (And, in digging into the script, you would
have a complete, unambiguous representation of what the article’s authors have
done with their data!) We’re not there yet, but why not introduce those good
practices right now? Comment and document your code. Preserve the definite
version of your analysis scripts together with your article. That way, you
will always be able to understand what you did earlier. If colleagues ask, you
can give them your script and your data, and they can check for themselves
what you have done. (And yes, I know you might be worried about giving away
your trade secrets. If that is the case, you could still make your scripts
available after some time. The competitive edge a specific script gives you is
likely to become less important over time in any case.)

If you use someone else’s data, or Python module, or software, you should make
sure to give them appropriate credit. Published astronomical data typically
comes with accompanying publications. Citing those publications when using the
data is the appropriate way of giving credit to the astronomers involved. As
an example, consider the THINGS survey data used in Fig. 22. Using this image
on a slide, you would add the information “Walter 2008,” or “Walter et al.
2008,” and every astronomer would be able to use the Astrophysics Data Service
ADS at https://ui.adsabs.harvard.edu/ to find the article. If you are writing
in a medium where you can add a link, add a link. If you are writing in a more
formal setting, put an entry into your reference section, and refer to that
entry in your text. For specific software or larger data sets, there is often
a standard paragraph of acknowledgement (“boilerplate”) you are asked to
include in your text. You can find examples in the Acknowledgements section,
below. Science is a networked activity, and giving credit where credit is due
is an important part of good scientific practice.

Astronomical data has never been as accessible as it is now, and computing
power never as cheap. Observatory and telescope archives provide images and
spectra, catalogs higher-level data like never before. This is likely to get
even better as new facilities come online, and new tools become available.
It’s an exciting time to work with astronomical data!

## Acknowledgements

I would like to thank the anonymous referee and also Wolfgang Brandner, Roland
Gredel and Carolin Liefke for helpful comments, Fabian Walter for helping with
the THINGS data and Catharina Hock for helpful discussions.

Working with astronomical data would not be as (comparatively) easy without
the countless individuals who have invested time and effort in programming the
software, packages and libraries used throughout this text. If you are using
the results, please acknowledge them appropriately. A number of scientific
software products have associated articles that you should cite when using
those products in your scientific work (e.g. see the footnotes for the table
on p. 5.3).

SAOImageDS9 development has been made possible by funding from the Chandra
X-ray Science Center (CXC) and the High Energy Astrophysics Science Archive
Center (HEASARC) with additional funding from the JWST Mission office at Space
Telescope Science Institute, cf. the article Joye & Mandel 2003,
[http://adsabs.harvard.edu/abs/2003ASPC..295..489J].

TOPCAT was developed, and keeps being developed further, by Mark Taylor
(University of Bristol). A description of the application can be found in
Taylor 2005, [http://adsabs.harvard.edu/abs/2005ASPC..347…29T].

This work makes use of observations from the LCOGT network, which is described
in Brown et al. 2013, [http://adsabs.harvard.edu/abs/2013PASP..125.1031B].

Funding for SDSS-III, from which I have used some data, has been provided by
the Alfred P. Sloan Foundation, the Participating Institutions, the National
Science Foundation, and the U.S. Department of Energy Office of Science. The
SDSS-III web site is [http://www.sdss3.org/]. SDSS-III is managed by the
Astrophysical Research Consortium for the Participating Institutions of the
SDSS-III Collaboration including the University of Arizona, the Brazilian
Participation Group, Brookhaven National Laboratory, Carnegie Mellon
University, University of Florida, the French Participation Group, the German
Participation Group, Harvard University, the Instituto de Astrofisica de
Canarias, the Michigan State/Notre Dame/JINA Participation Group, Johns
Hopkins University, Lawrence Berkeley National Laboratory, Max Planck
Institute for Astrophysics, Max Planck Institute for Extraterrestrial Physics,
New Mexico State University, New York University, Ohio State University,
Pennsylvania State University, University of Portsmouth, Princeton University,
the Spanish Participation Group, University of Tokyo, University of Utah,
Vanderbilt University, University of Virginia, University of Washington, and
Yale University.

Some of the sample files used here are based on observations made with the
NASA/ESA Hubble Space Telescope, and obtained from the Hubble Legacy Archive,
which is a collaboration between the Space Telescope Science Institute
(STScI/NASA), the Space Telescope European Coordinating Facility (ST-ECF/ESA)
and the Canadian Astronomy Data Centre (CADC/NRC/CSA).

The astronomical data query language ADQL was developed by the Virtual
Observatory Query Language Working Group, cf. the description in Ortiz et al.
2008, [http://www.ivoa.net/documents/latest/ADQL.html].

This work has made use of data from the European Space Agency (ESA) astrometry
mission Gaia ([https://www.cosmos.esa.int/gaia]), as processed by the Gaia
Data Processing and Analysis Consortium (DPAC, see more at
[https://www.cosmos.esa.int/web/gaia/dpac/consortium]). Funding for the DPAC
has been provided by national institutions, in particular the institutions
participating in the Gaia Multilateral Agreement; cf. Gaia Collaboration 2016,
[http://adsabs.harvard.edu/abs/2016A&A…595A…1G], and the description in Gaia
Collaboration 2018, [http://adsabs.harvard.edu/abs/2018arXiv180409365G].


# Fisher-Weighted Merge of Contrastive Learning Models
in Sequential Recommendation

Jung Hyun Ryu  Jaeheyoung Jeon  Jewoong Cho  Myungjoo Kang

###### Abstract

Along with the exponential growth of online platforms and services,
recommendation systems have become essential for identifying relevant items
based on user preferences. The domain of sequential recommendation aims to
capture evolving user preferences over time. To address dynamic preference,
various contrastive learning methods have been proposed to target data
sparsity, a challenge in recommendation systems due to the limited user-item
interactions. In this paper, we are the first to apply the Fisher-Merging
method to Sequential Recommendation, addressing and resolving practical
challenges associated with it. This approach ensures robust fine-tuning by
merging the parameters of multiple models, resulting in improved overall
performance. Through extensive experiments, we demonstrate the effectiveness
of our proposed methods, highlighting their potential to advance the state-of-
the-art in sequential learning and recommendation systems.

Machine Learning, ICML

Figure 1: Overview of Parameter Merging. Sharing model architecture, each
models differ in how the contrastive loss is constructed. We merge models by
using a weighted sum, where the weights are determined based on the posterior
distribution of each model’s parameters, assuming a Gaussian distribution.

##  1 Introduction

With the exponential growth of online platforms and services, a significant
amount of data is being generated daily. Recommendation systems have become
crucial in utilizing this data effectively. The system aim to identify
relevant items based on user preferences and interests. As user preferences
evolve over time, sequential recommendation has gained attention as a subfield
in this area. We address the problem of sequential recommendation as follows.

Let $\mathcal{U}$ be the set of users
$\mathcal{U}=\\{u_{1},u_{2},\cdots,u_{|\mathcal{U}|}\\}$, and $\mathcal{V}$ be
the set of items as $\mathcal{V}=\\{v_{1},v_{2},\cdots,v_{|\mathcal{V}|}\\}$.
The sequence of user-item interaction for $u_{i}$ is a list with chronological
order,
$s_{i}=[v_{1}^{u_{i}},v_{2}^{u_{i}},\cdots,v_{t}^{u_{i}},\cdots,v_{n_{u_{i}}}^{u_{i}}]$.
Here user $u_{i}\in\mathcal{U}$, $v_{t}^{u_{i}}\in\mathcal{V}$, and user
$u_{i}$ interact item $v_{t}^{u_{i}}$ in time step $t$. The length of sequence
for user $u_{i}$ is $n_{u_{i}}$, and our object is to build a model predicting
the item with which user is interact in the next time step, i.e.

$p\left(v_{n_{u_{i}}+1}^{u_{i}}=v\,\middle|\,s_{i}\right).$ (1)

The previous methodologies typically employ similar model structures but
utilize various learning frameworks (Xie et al., 2022; Qiu et al., 2022).
Prior research has shown that ensemble methods yield significant benefits when
multiple learning frameworks are employed (Gontijo-Lopes et al., 2021).

We propose a practical and feasible method to ensemble the parameters of
models trained with different contrastive learning techniques in a sequential
recommendation. The purpose of this study is to effectively aggregate the
obtained parameters $\theta$ in various learning frameworks and hyperparameter
settings, building on previous research and experiments.

Assuming the posterior distribution of parameters $\theta_{m}$ for each $m$-th
model, Sec 3.1, we achieved more effective ensemble results. This approach
allowed us to capture the uncertainty associated with each model’s parameter
estimates and leverage this information to enhance the ensemble process. By
considering the posterior distributions, we were able to account for the
variability in parameter values across different models and obtain a more
robust and comprehensive ensemble outcome.

##  2 Related Works

Researchers have explored various ensemble methods, including bootstrapping,
bagging, and boosting, to improve model performance (Ganaie et al., 2022;
Breiman, 1996, 2001; Natekin & Knoll, 2013; Liu et al., 2014). Ovadia et al.
(Ovadia et al., 2019) demonstrated the accuracy of ensembles even in the
presence of distribution shift, while Mustafa et al. (Mustafa et al., 2020)
proposed a method that combines fine-tuned subsets of pre-trained models to
achieve high accuracy and resilience to distribution shift. Parameter merging
is another technique to reduce model size and computational requirements
(Houlsby et al., 2019). However, ensemble methods often require additional
training, which can be computationally expensive and time-consuming.

###  2.1 Diverse Learning Framework

Wenzel et al. (2020) and Zaidi et al. (2021) investigated the role of random
hyperparamters and architectures in ensemble. Gontijo et al. (2021)
demonstrated the ensemble effect across various training methodologies;
initialization, hyperparameter, architecture, framework, and dataset levels.
Diverse training methodologies exhibit different generalization capabilities,
ultimately lead to uncorrelated errors. Models tend to specialize in
subdomains within the data and highlights the crucial role of ensemble
techniques in enhancing overall performance.

###  2.2 Merging Methods

##### Model Soup

Model Soup (Wortsman et al., 2022) presents an effective approach for
combining parameters without additional training. It demonstrates research
findings that improve the performance of trained models by constructing a
”recipe” composed of diverse models and averaging their parameters. The study
introduces three methods for creating the recipe: the uniform soup, which
averages the parameter values of all models; the greedy soup, which
sequentially adds models based on their performance ranking; and the learned
soup, which identifies the optimal model interpolation through training. These
approaches contribute to enhancing the overall performance of the model
without the need for additional training.

##### Fisher Merging

Within the scope of related works, parameter merging is interpreted as a
process that maximizes the joint likelihood of model parameters’ posteriors
(Matena & Raffel, 2022). Previous study (Wortsman et al., 2022) consider
averaging as a scenario where the posteriors of these models are assumed to
follow an isotropic Gaussian distribution, and the joint likelihood is
maximized accordingly. To refine this approach, efforts have been made to
approximate the posterior of the model using Laplace approximation (Daxberger
et al., 2021). In this case, the distribution of each model is modeled by
assuming the mean as the observed, which can be interpreted as trained
parameter and the variance as the Gaussian distribution’s Fisher matrix. By
employing this formulation, the joint likelihood is calculated.

###  2.3 Sequential Recommendation System

SASRec (Kang & McAuley, 2018) employ Transformer layers to dynamically assign
weights to previous items. BERT4Rec (Sun et al., 2019) demonstrate an
improvement by incorporating user behavior information from both directions
using a bidirectional Transformer. CL4SRec (Xie et al., 2022) employed three
data augmentation techniques, namely item cropping, item masking, and item
reordering, to create pairs for contrastive learning. DuoRec (Qiu et al.,
2022) integrated two types of contrastive loss. Firstly, it incorporated
unsupervised augmentation using dropout-based model-level augmentation to
generate positive pairs. Secondly, it incorporated supervised positive
sampling, which involves creating pairs by considering sequences with the same
target item as positive samples.

##  3 Methodology

We perform model ensemble based on different types of loss functions. BERT4Rec
(Sun et al., 2019), CL4SRec (Xie et al., 2022), and DuoRec (Qiu et al., 2022)
share the basic structure of BERT4Rec (Sun et al., 2019). However, they differ
in the sense of constructing positive pairs, a key component of their learning
framework of constrastive learning.

Figure 1 represents the overview of parameter merging process. By sharing the
structure of the model, which is parameterized with diverse learning
frameworks, we can leverage ensemble methods to our advantage. Furthermore,
inspired by previous studies demonstrating the effectiveness of ensemble
models trained using various learning methods, we apply parameter merging
techniques, namely Parameter Averaging and Fisher-weighted Parameter Merging,
described in Section 3.1, to combine these models.

###  3.1 Understanding Model Ensemble

We follow the work of Matena et al (2022). Consider a scenario where we have
models with the same structure,
$\text{model}_{1},\text{model}_{2},\cdots,\text{model}_{M}$, with
corresponding parameters $\theta_{1},\theta_{2},\cdots,\theta_{M}$. Our
objective is to find the parameter $\theta^{*}$ that maximizes the joint
likelihood of the posteriors of these parameters.

The posterior of $\theta_{m}$ can be represented as
$p\left(\theta|\theta_{m}\right)$. Since obtaining this posterior directly is
generally challenging, it can employ approximation methods such as Laplace
approximation to make assumptions and seek the parameter $\theta^{*}$ (MacKay,
1992; Daxberger et al., 2021). Let us interpret the process of finding
$\theta^{*}$ as maximizing the joint likelihood, $\sum_{m}\log
p\left(\theta\mid\theta_{m}\right)$. Assuming that
$p\left(\theta\mid\theta_{m}\right)$ follows a Gaussian distribution, we set
the mean of this Gaussian distribution as the observed $\theta_{m}$ and
examine the procedure for averaging parameters and Fisher merging separately,
depending on the method used to assume the variance.

##### Averaging Parameters

Assume that the posterior $p\left(\theta\mid\theta_{m}\right)$ follows a
Gaussian distribution $\mathcal{N}\left(\hat{\theta}_{m},I\right)$. Here,
$\hat{\theta}_{m}$ represents the parameters of the trained $m$-th model, and
$I$ denotes the identity matrix. In this case, the desired solution
$\theta^{*}$ can be obtained as the average of the parameters of the candidate
models, as shown in eq.2:

$\theta^{*}=\underset{\theta}{\arg\max}\sum_{m}\log
p\left(\theta\mid\theta_{m},I\right)=\frac{1}{M}\sum_{m}\theta_{m}.$ (2)

##### Fisher Merging

Let us consider the posterior $p\left(\theta\mid\theta_{m}\right)$ following a
Gaussian distribution $\mathcal{N}\left(\hat{\theta}_{m},H^{-1}\right)$. Here,
$\hat{\theta}_{m}$ represents the parameters of the trained $m$-th model, and
$H$ corresponds to the Hessian matrix of $\theta_{m}$ obtained through the
second-order Taylor expansion at the mode of the posterior. It has been
established that the Hessian matrix in this distribution coincides with the
Fisher information, but for computational efficiency, we only utilize the
diagonal elements of the Fisher matrix (Kirkpatrick et al., 2017).

The desired solution $\theta^{*}$ can be expressed as eq.3, capturing the
essence of the Fisher likelihood :

$\theta^{*}=\underset{\theta}{\arg\max}\sum_{m}\lambda_{m}\log
p\left(\theta\mid\theta_{m},F_{m}\right),$ (3)

where $F_{m}=\mathbb{E}_{x}\mathbb{E}_{y\sim
p_{\theta}\left(y|x\right)}\nabla_{\theta}\log{p_{\theta}\left(y|x\right)}\nabla_{\theta}\log{p_{\theta}\left(y|x\right)}^{T}$.
The closed-form solution for $\theta^{*}$ can be obtained as shown in eq.4,
which directly incorporates the Fisher matrix. In practice, we utilize an
empirical estimate of the Fisher matrix, denoted as $\hat{F}$, as shown eq.4
(Kirkpatrick et al., 2017).

$\theta^{*\left(j\right)}=\frac{\sum_{m}\lambda_{m}F_{m}^{\left(j\right)}\theta_{m}^{\left(j\right)}}{\sum_{m}\lambda_{m}F_{m}^{\left(j\right)}},$
(4)

where $F_{m}=\frac{1}{N}\mathbb{E}_{y\sim p_{\theta}\left(y\mid
x\right)}\left(\nabla_{\theta}\log{p_{\theta}\left(y\mid x\right)}\right)^{2}$
and $j=1,\cdots,|\theta|$, considering as element-wise multiplication.

###  3.2 Applying Model Ensemble

By expressing the Fisher matrix we intend to compute in eq.4 in terms of
recommendation factors, we can decompose it into the following components:

$\displaystyle\mathbb{E}_{x_{i}}\mathbb{E}_{y\sim p_{\theta}\left(y\mid
x_{i}\right)}\left(\nabla_{\theta}\log{p_{\theta}\left(y\mid
x_{i}\right)}\right)^{2}$ (5)
$\displaystyle=\frac{1}{N}\sum_{i}\sum_{j}p_{\theta}\left(y_{j}\mid
x_{i}\right)\left(\nabla_{\theta}\log{p_{\theta}\left(y_{j}\mid
x_{i}\right)}\right)^{2}$
$\displaystyle=\frac{1}{|\mathcal{U}|}\sum_{i}^{|\mathcal{U}|}\sum_{j}^{|\mathcal{V}|}p_{\theta}\left(v_{j}\mid
s_{i}\right)\left(\nabla_{\theta}\log{p_{\theta}\left(v_{j}\mid
s_{i}\right)}\right)^{2}.$

There are two computational challenges associated with the above equation.
First, calculations need to be performed for each individual sample $s_{i}$.
Second, calculations need to be performed for each item $v_{j}$ within a
single sample. The reason why these points acts as a drawback in
recommendation systems is due to the large number of users and items in the
data. For instance, in the case of MovieLens-1M dataset (Harper & Konstan,
2015), there are about 6000 users and 3500 items. However, performing Fisher
matrix calculations that require differentiation with respect to $\theta$ for
each user and item becomes a computational burden.

####  3.2.1 Sampling sequences

##### Batch-wise Computation

To address the first challenge of performing computations on individual
samples, we reinterpret the equation and carry out the calculations on a batch
basis. It should be noted that $p_{\theta}\left(v_{j}|s_{i}\right)$ can vary
for each sample $s_{i}$. Therefore, we perform the sorting of
$p_{\theta}\left(v_{j}|s\right)$ to address this variation, where BS indicates
batch size:

$\displaystyle\sum_{i}^{|\mathcal{U}|}\sum_{j}^{|\mathcal{V}|}p_{\theta}\left(v_{j}\mid
s_{i}\right)\left(\nabla_{\theta}\log{p_{\theta}\left(v_{j}\mid
s_{i}\right)}\right)^{2}$ (6)
$\displaystyle=\sum_{\text{BS}_{k}}\sum_{j}^{|\mathcal{V}|}\left(\sum_{i}^{\text{BS}_{k}}p_{\theta}\left(v_{j}|s_{i}\right)\right)\left(\nabla_{\theta}\sum_{i}^{\text{BS}_{k}}\log{p_{\theta}\left(v_{j}|s_{i}\right)}\right)^{2}.$

####  3.2.2 Sampling items

To alleviate the computational burden associated with iterating over all $j$
values, which scales with $|\mathcal{V}|$, we employ a sampling-based approach
within the methodology. This sampling strategy aims to reduce the
computational cost while maintaining the representativeness of the
calculations.

##### Random Sampling

We compute the eq.5 by randomly sampling $j$ from the total number of items.
This process was performed to calculate the Fisher matrix without any specific
assumptions or prior knowledge.

##### Top-$k$ Sampling

The probability which is output by the model can be interpreted as the
preference or likelihood of the recommended items for a given sample. Based on
this interpretation, we select a set of $n$ items that are most likely to be
of interest to the corresponding user, i.e. $p_{\theta}\left(v_{j}\mid
s_{i}\right)$. Subsequently, we compute the Fisher matrix with these selected
items as the focal points. By focusing on this subset of items that are
expected to be of highest interest, we aim to capture the relevant information
for optimizing the model’s performance effectively.

$\displaystyle\sum_{j}^{|\mathcal{V}|}p_{\theta}\left(v_{j}\mid
s\right)\left(\nabla_{\theta}\log{p_{\theta}\left(v_{j}\mid
s\right)}\right)^{2}$ (7)
$\displaystyle\approx\sum_{j}^{\text{top-}k}p_{\theta}\left(v_{j}\mid
s\right)\left(\nabla_{\theta}\log{p_{\theta}\left(v_{j}\mid
s\right)}\right)^{2}.$

##### Model-based Sampling

To select a subset of items for further analysis, we randomly sampled items
based on their conditional probability $p_{\theta}(v_{j}|s_{i})$ using a
weighted random selection process. The selection probability of each item was
determined by its associated probability stored in the model’s output. By
selecting items with higher probabilities, we focused on a specific number of
items that were more likely to align with the user’s preferences or interests.
This allowed us to analyze and evaluate the subset of items based on their
associated probabilities obtained from the model’s output. With $N$ denotes
the sample size, his approximation can be represented as:

$\displaystyle\mathbb{E}_{y\sim p_{\theta}\left(y\mid
x\right)}\left(\nabla_{\theta}\log{p_{\theta}\left(y\mid x\right)}\right)^{2}$
(8) $\displaystyle\approx\frac{1}{N}\sum_{v_{j}\sim p_{\theta}\left(v_{j}\mid
s\right)}^{N}\left(\nabla_{\theta}\log{p_{\theta}\left(v_{j}\mid
s\right)}\right)^{2}.$

##### Calculate with target item

We compute the Fisher matrix based on the target item, disregarding other
items with limited direct relevance. By employing this approach, we focus
solely on the target item and its associated information to calculate the
Fisher matrix. Our rationale behind this decision is to prioritize the target
item’s impact on the model’s optimization process, as it is directly linked to
the specific objective or task at hand. Consequently, we exclude items with
minimal direct relevance to ensure a more targeted and meaningful computation
of the Fisher matrix.

$p_{\theta}\left(v_{j}^{*}\mid
s\right)\left(\nabla_{\theta}\log{p_{\theta}\left(v_{j}^{*}\mid
s\right)}\right)^{2},$ (9)

where $v_{j}^{*}$ is the target item.

##  4 Experiments

We use MovieLens-1M dataset (Harper & Konstan, 2015) for experiments. For each
user, we have sequential data consisting of movies purchased in chronological
order. We adopt next-item prediction task (i.e. leave-one-out evaluation),
following previous works (Sun et al., 2019; Xie et al., 2022; Qiu et al.,
2022). The last movie is considered as the test set, and the validation data
is used to predict the preceding movies. During training, we adopt a masked
language modeling approach similar to BERT (Devlin et al., 2018), where we
mask certain movies in the sequentially ordered list and task the model with
predicting them.

The evaluation method used in this study is the Normalized Discounted
Cumulative Gain at 10 (NDCG@10), which is a ranking-based evaluation approach
(He et al., 2017). It ranks the top 10 items predicted by the model based on
their perceived preference and considers the actual ranking of the preferred
items. A higher NDCG value, closer to 1, indicates better performance.
Different NDCG values can be obtained depending on the selection of items,
such as from the full item pool, a random set of 100 items, or the top 100
most popular items.

###  4.1 Results of Model Merging

Examine the results through Table 1 and Table 2. Table 1 presents the results
obtained by training models, namely BET4Rec (Sun et al., 2019), CL4SRec (Xie
et al., 2022), and DuoRec (Qiu et al., 2022). We merge these models using
Fisher methods. While Table 1 demonstrates the results of models trained
solely from scratch. Table 2 represents the results of fine-tune setting. We
train the baseline model without contrastive loss for 20 epochs, which is the
convergence point of the baseline experiment without any additional
contrastive loss, similar to BERT4Rec. Following this, each model; BERT4Rec
(Sun et al., 2019), CL4SRec (Xie et al., 2022), DuoRec (Qiu et al., 2022),
underwent fine-tuning according to their respective methods, and the results
were merged using Fisher methods. In both conditions, we fine-tuned addtional
epoch after merging process.

Fisher merge fails to improve the performance of individual models in baseline
setting. When Fisher merge is applied during the fine-tuning setting, it leads
to improved performance compared to individual models. This finding aligns
with previously reported phenomena (Ganaie et al., 2022) where individual
models tend to achieve higher performance than merged model in the baseline
setting. However, the results of Fisher merge in the fine-tuning setting show
comparable performance with the individual models in baseline setting, while
the individual recipe models of fine-tuning setting do not exceed. Also, the
results indicate that even for models that have not been sufficiently trained
such as CL4SRec in our setting, merging parameters resulted in comparable
performance to other models, demonstrating robustness.

Table 1: Results of parameter merging; Fisher-merge on Baseline Settings. The
recipes used for merging were trained for the same number of epochs. ‘POS.’
refers to the method of constructing positive pair, ‘sup’ to supervised
augmentation and ‘-’ to supervised and unsupervised augmentation in subsection
2.3
model | pos | full | random | popular
---|---|---|---|---
baseline |  | 0.1398 | 0.5651 | 0.0482
cl4srec |  | 0.0955 | 0.043 | 0.0429
duorec | sup | 0.1348 | 0.5575 | 0.044
| unsup | 0.1301 | 0.5592 | 0.0438
| - | 0.1382 | 0.5588 | 0.0464
fisher |  | 0.1289 | 0.5495 | 0.0472
Table 2: Results of parameter merging; Fisher-merge on fine-tune Settings. The
recipes used for merging were trained on baseline model (without contrastive
loss) and fine-tuned on each model.
model | pos | full | random | popular
---|---|---|---|---
baseline |  | 0.135 | 0.5573 | 0.0426
cl4srec |  | 0.0585 | 0.0513 | 0.0466
duorec | sup | 0.1346 | 0.5547 | 0.0454
| unsup | 0.1358 | 0.5594 | 0.0445
| - | 0.1351 | 0.554 | 0.0423
fisher |  | 0.1386 | 0.5618 | 0.0428
Table 3: Effect of Sampling Method and Sampling Size. We merge models in
settings of Table 2, the fine-tune setting. We merged models with 4 sampling
methods; random sampling, top-k sampling, model-based sampling, and calculate
on target item, on 3 different sampling size; n=10, n=30 and n=50. Bold
represents the best variant in each evaluation setting, and underlines
indicates the second best variation.
| sample size | FULL | RANDOM | POPULAR
---|---|---|---|---
|  | NDCG10 | NDCG20 | NDCG10 | NDCG20 | NDCG10 | NDCG20
baseline |  | 0.135 | 0.1601 | 0.5573 | 0.5786 | 0.0426 | 0.0706
CL4SRec |  | 0.0585 | 0.0751 | 0.0513 | 0.043 | 0.0466 | 0.0701
DuoRec (sup.) |  | 0.1346 | 0.1591 | 0.5547 | 0.58 | 0.0454 | 0.068
DuoRec (unsup.) |  | 0.1358 | 0.1609 | 0.5594 | 0.5782 | 0.0445 | 0.0742
DuoRec (sup.&unsup.) |  | 0.1351 | 0.1599 | 0.554 | 0.5732 | 0.0423 | 0.0724
random sampling | 10 | 0.1379 | 0.1638 | 0.5606 | 0.5825 | 0.0457 | 0.0691
| 30 | 0.1366 | 0.1624 | 0.5584 | 0.58 | 0.0477 | 0.0726
| 50 | 0.1386 | 0.1636 | 0.5598 | 0.5813 | 0.0419 | 0.0419
top-k sampling | 10 | 0.1364 | 0.1624 | 0.5602 | 0.5817 | 0.0446 | 0.0689
| 30 | 0.1373 | 0.1616 | 0.5637 | 0.5835 | 0.0457 | 0.0708
| 50 | 0.1387 | 0.1635 | 0.5592 | 0.5807 | 0.0424 | 0.0672
model-based sampling | 10 | 0.1358 | 0.1619 | 0.5564 | 0.5782 | 0.044 | 0.0696
| 30 | 0.1385 | 0.1646 | 0.5579 | 0.5784 | 0.0446 | 0.0689
| 50 | 0.138 | 0.1632 | 0.5605 | 0.5814 | 0.0465 | 0.0719
calculate on target item | 1 | 0.1386 | 0.1628 | 0.5618 | 0.5806 | 0.0428 | 0.0725

###  4.2 The Validity of Batch-wise Computation

We performed batch-wise computations with the aim of implementing an efficient
Fisher matrix calculation. Compared to computing on individual samples,
grouping samples into batches allowed us to achieve computational efficiency.

The following Figure 4 in Appendix A.3 illustrates the method for minimizing
errors when performing calculations on a batch basis. The figure demonstrates
that within a batch containing 10 samples, denoted as $s_{i}$, there is a
phenomenon where the probabilities of item $v_{j}$ decrease in a similar
manner. By sorting the samples $s_{i}$ based on the probability of $v_{j}$,
even when grouping them into batches, it is possible to minimize the error
described by the eq.6. Furthermore, the figure illustrates the rationale
behind top-k sampling. For the top-k items, the probabilities hold meaningful
information, whereas for the remaining items, the probabilities are nearly
zero or close to it.

###  4.3 Effect of Sampling Methods and Size

To investigate the effect of sampling methods, we conduct experiments by
varying the number of sampled samples and the sampling techniques employed.
Specifically, we consider three sample sizes: $n=10,n=30\text{ and }n=50$, and
four different sampling methods: random sampling, top-k sampling, model-based
sampling, and calculate with target item. The results of these experiments can
be observed in Table 3. The table provides insights into the performance of
each sampling method under different sample sizes, allowing us to analyze
their respective effects on the task at hand. Note that this result is
calculated on batched data. To examine the results of parameter merging, we
conducted experiments in fine-tuning setting, explained in 2.

The experiments revealed effective ensemble results, particularly showcasing
the robust performance of CL4SRec (Xie et al., 2022). Despite having
significantly lower performance compared to other models during the parameter
merging process, the model with poor performance exhibited robust performance
in the Fisher merge results. Regarding the sampling methods, top-k sampling
demonstrated the best performance. This can be attributed to the concentration
of probabilities assigned to specific items by the model, effectively
approximating the Fisher criterion sought in the evaluation. Also, the model-
based sampling method exhibits a more pronounced improvement in performance as
the sampling size increases compared to other models. We interpret these
results as being rooted in the direct interpretation of the equation defined
for Fisher merging. Interestingly, despite the fact that calculating Fisher
matrix on target item has a single sample, the method demonstrated sampling
effectiveness by achieving good performance even with a small sample size.
These findings shed light on the interpretation of experimental results in the
context of deep learning research.

Figure 2: Measured Time Consumed for Each Sampling Method and Size. We sample
items from MovieLens-1M dataset. Time is measured in batch-wise setting, where
batch size is 256.

###  4.4 Computational Cost

Figure 2 demonstrates computational cost in terms of time consumed during
calculating Fisher matrix for single model. The concept of parameter merging
involves additional computation on the existing parameters. Therefore, it is
important to ensure efficiency in this process. To achieve efficiency,
considerations such as calculating the Fisher matrix in batch units and
performing sampling are necessary. It is observed that, except for the
calculation on the target item, the computational complexity increases
linearly with the sampling size. As for the calculation on the target item,
the sampling size remains fixed at 1 since each sequence has a single target
item. Thus, our research is significant as it approximates the Fisher matrix
calculation with a much smaller number of items (around 3000) compared to
calculating it on the entire item set.

Figure 3: Visualization for Weights of Merged Models. Based on the plane
containing 64-dimension parameters of three model, we visualised its weight,
100 samples from each posteriors and merged parameters

###  4.5 Visualization of Merged Weights

We present a visual illustration to aid in the intuitive understanding of the
merged weights. Figure 3 represents the fine-tuning setting of 2, where the
three centroids correspond to the weights of individual models. The plane
visualized in 3 encompasses these three weights. The scattered points,
projected onto the plane, depict 100 samples drawn from
$\mathcal{N}(\theta_{m},F_{m})$. It is observed that the baseline weight
exhibits the largest variance. This can be attributed to the experimental
setup where the baseline is pre-trained and then fine-tuned with CL4SRec (Xie
et al., 2022) and DuoRec (Qiu et al., 2022). The weights obtained through
uniform merging are represented as the average of the three centroid points,
while the weights obtained through Fisher merging take into account the
variances of these recipe weights. It can be seen that the weights obtained
through Fisher merging considered posterior and variance with Laplace
approximation and provides nice initial point for fine-tune.

##  5 Conclusion

We apply ensemble technique, Fisher merging, for sequential models, enabling
robust fine-tuning through parameter merging. Our experimental results
demonstrate the effectiveness of these proposed methods in improving
recommendation performance. These contributions have the potential to advance
the field of sequential learning and recommendation systems, offering valuable
insights for future research and practical applications.

## Acknowledgements

This work was support by the NRF grant [2012R1A2C3010887] and the MSIT/IITP
[1711117093, 2021-0-00077, 2021-0-01343, Artificial Intelligence Graduate
School Program(SNU)].

## References

  * Breiman (1996) Breiman, L.  Bagging predictors.  _Machine learning_ , 24:123–140, 1996.
  * Breiman (2001) Breiman, L.  Random forests.  _Machine learning_ , 45:5–32, 2001.
  * Daxberger et al. (2021) Daxberger, E., Kristiadi, A., Immer, A., Eschenhagen, R., Bauer, M., and Hennig, P.  Laplace redux-effortless bayesian deep learning.  _Advances in Neural Information Processing Systems_ , 34:20089–20103, 2021.
  * Devlin et al. (2018) Devlin, J., Chang, M.-W., Lee, K., and Toutanova, K.  Bert: Pre-training of deep bidirectional transformers for language understanding.  _arXiv preprint arXiv:1810.04805_ , 2018.
  * Ganaie et al. (2022) Ganaie, M. A., Hu, M., Malik, A., Tanveer, M., and Suganthan, P.  Ensemble deep learning: A review.  _Engineering Applications of Artificial Intelligence_ , 115:105151, 2022.
  * Geirhos et al. (2020) Geirhos, R., Meding, K., and Wichmann, F. A.  Beyond accuracy: quantifying trial-by-trial behaviour of cnns and humans by measuring error consistency.  _Advances in Neural Information Processing Systems_ , 33:13890–13902, 2020.
  * Gontijo-Lopes et al. (2021) Gontijo-Lopes, R., Dauphin, Y., and Cubuk, E. D.  No one representation to rule them all: Overlapping features of training methods.  _arXiv preprint arXiv:2110.12899_ , 2021.
  * Harper & Konstan (2015) Harper, F. M. and Konstan, J. A.  The movielens datasets: History and context.  _Acm transactions on interactive intelligent systems (tiis)_ , 5(4):1–19, 2015.
  * He et al. (2017) He, X., Liao, L., Zhang, H., Nie, L., Hu, X., and Chua, T.-S.  Neural collaborative filtering.  In _Proceedings of the 26th international conference on world wide web_ , pp. 173–182, 2017.
  * Houlsby et al. (2019) Houlsby, N., Giurgiu, A., Jastrzebski, S., Morrone, B., De Laroussilhe, Q., Gesmundo, A., Attariyan, M., and Gelly, S.  Parameter-efficient transfer learning for nlp.  In _International Conference on Machine Learning_ , pp. 2790–2799. PMLR, 2019.
  * Kang & McAuley (2018) Kang, W.-C. and McAuley, J.  Self-attentive sequential recommendation.  In _2018 IEEE international conference on data mining (ICDM)_ , pp. 197–206. IEEE, 2018.
  * Kirkpatrick et al. (2017) Kirkpatrick, J., Pascanu, R., Rabinowitz, N., Veness, J., Desjardins, G., Rusu, A. A., Milan, K., Quan, J., Ramalho, T., Grabska-Barwinska, A., et al.  Overcoming catastrophic forgetting in neural networks.  _Proceedings of the national academy of sciences_ , 114(13):3521–3526, 2017.
  * Liu et al. (2014) Liu, P., Han, S., Meng, Z., and Tong, Y.  Facial expression recognition via a boosted deep belief network.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pp. 1805–1812, 2014.
  * MacKay (1992) MacKay, D. J.  A practical bayesian framework for backpropagation networks.  _Neural computation_ , 4(3):448–472, 1992.
  * Matena & Raffel (2022) Matena, M. S. and Raffel, C. A.  Merging models with fisher-weighted averaging.  _Advances in Neural Information Processing Systems_ , 35:17703–17716, 2022.
  * Mustafa et al. (2020) Mustafa, B., Riquelme, C., Puigcerver, J., Pinto, A. S., Keysers, D., and Houlsby, N.  Deep ensembles for low-data transfer learning.  _arXiv preprint arXiv:2010.06866_ , 2020.
  * Natekin & Knoll (2013) Natekin, A. and Knoll, A.  Gradient boosting machines, a tutorial.  _Frontiers in neurorobotics_ , 7:21, 2013.
  * Ovadia et al. (2019) Ovadia, Y., Fertig, E., Ren, J., Nado, Z., Sculley, D., Nowozin, S., Dillon, J., Lakshminarayanan, B., and Snoek, J.  Can you trust your model’s uncertainty? evaluating predictive uncertainty under dataset shift.  _Advances in neural information processing systems_ , 32, 2019.
  * Qiu et al. (2022) Qiu, R., Huang, Z., Yin, H., and Wang, Z.  Contrastive learning for representation degeneration problem in sequential recommendation.  In _Proceedings of the fifteenth ACM international conference on web search and data mining_ , pp. 813–823, 2022.
  * Sun et al. (2019) Sun, F., Liu, J., Wu, J., Pei, C., Lin, X., Ou, W., and Jiang, P.  Bert4rec: Sequential recommendation with bidirectional encoder representations from transformer.  In _Proceedings of the 28th ACM international conference on information and knowledge management_ , pp. 1441–1450, 2019.
  * Wang & Isola (2020) Wang, T. and Isola, P.  Understanding contrastive representation learning through alignment and uniformity on the hypersphere.  In _International Conference on Machine Learning_ , pp. 9929–9939. PMLR, 2020.
  * Wenzel et al. (2020) Wenzel, F., Snoek, J., Tran, D., and Jenatton, R.  Hyperparameter ensembles for robustness and uncertainty quantification.  _Advances in Neural Information Processing Systems_ , 33:6514–6527, 2020.
  * Wortsman et al. (2022) Wortsman, M., Ilharco, G., Gadre, S. Y., Roelofs, R., Gontijo-Lopes, R., Morcos, A. S., Namkoong, H., Farhadi, A., Carmon, Y., Kornblith, S., et al.  Model soups: averaging weights of multiple fine-tuned models improves accuracy without increasing inference time.  In _International Conference on Machine Learning_ , pp. 23965–23998. PMLR, 2022.
  * Xie et al. (2022) Xie, X., Sun, F., Liu, Z., Wu, S., Gao, J., Zhang, J., Ding, B., and Cui, B.  Contrastive learning for sequential recommendation.  In _2022 IEEE 38th international conference on data engineering (ICDE)_ , pp. 1259–1273. IEEE, 2022.
  * Yosinski et al. (2015) Yosinski, J., Clune, J., Nguyen, A., Fuchs, T., and Lipson, H.  Understanding neural networks through deep visualization.  _arXiv preprint arXiv:1506.06579_ , 2015.
  * Zaidi et al. (2021) Zaidi, S., Zela, A., Elsken, T., Holmes, C. C., Hutter, F., and Teh, Y.  Neural ensemble search for uncertainty estimation and dataset shift.  _Advances in Neural Information Processing Systems_ , 34:7898–7911, 2021.

##  Appendix A Appendix

###  A.1 Motivation : Error Inconsistency

Table 4: Effect of Constructing Pair of Contrastive Loss (%). We observe that
models with divergent training methodologies exhibit distinct generalization
behavior, resulting in highly uncorrelated errors.
| similar |  | disssimilar
---|---|---|---
CL4SRec (Xie et al., 2022) | 8.05 | $<$ | 11.41
DuoRec (Qiu et al., 2022) | 8.67 | $<$ | 11.18

Previous research (Gontijo-Lopes et al., 2021; Yosinski et al., 2015)
demonstrated the increased effectiveness of ensemble methods as error
inconsistency grows. Building upon the existing research discourse, we
conducted the current experiment. In this study, we analyze the impact of
Fisher merging in the context of sequential recommendation systems,
attributing its effectiveness to the selection of recipe models trained using
different frameworks .

In our experiments, we employ a model based on the BERT4Rec (Sun et al., 2019)
architecture as our baseline. To enhance the performance of the model, we
apply various data augmentation techniques to enable contrastive learning.

To analyze the effects of contrastive loss, we divide the training frameworks
into two categories: similar and dissimilar. The similar learning frameworks
are trained using the same loss function but with slight variations such as
different seeds and hyperparameters, indicating the relationship among models
trained with small changes. On the other hand, the dissimilar learning
frameworks involve different data augmentation techniques, resulting in
variations in the construction of positive and negative pairs for contrastive
loss (Wang & Isola, 2020).

Error inconsistency (Geirhos et al., 2020) refers to the percentage of data
where two models have different classification results, with one model making
correct predictions while the other model makes incorrect predictions. Since
we are not dealing with classification, we considered a model to have made a
correct prediction if the value of NDCG@10 is above 0.5.

By comparing the error inconsistency between similar framework and dissimilar
framework, we observe the effectiveness of contrastive loss. An observation
that can be inferred from Table 4 is that the constructing positive pair for
contrastive loss significantly affects the similarity of the samples that the
models predict accurately. As the method for constructing positive pair
varies, the models demonstrate a considerable difference in their ability to
predict samples correctly. This finding highlights the sensitivity of the
models to the specific construction of the contrastive loss, which in turn
impacts their predictive performance.

###  A.2 Robustness of Fisher Merging; Recipe Selection

We compared two different recipe selection in Table 5; Fisher merged
parameters with least performance model and Fisher merged parameters without
the model. In our experimental setup, CL4SRec did not exhibit superior
performance compared to other models, considering the chosen hyperparameter
settings and other factors. Therefore, we aim to leverage the elements of the
recipe to demonstrate the robustness effect of Fisher merge. Our findings
confirm that by removing underperforming models as individual components and
applying Fisher merge, the resulting ensemble demonstrates robustness.

Table 5: Ablation Study of Results of parameter merging; Fisher-merge on fine-
tune Settings. The recipes used for merging were trained on baseline model
(without contrastive loss) and fine-tuned on each model. Ablation Study
illustrates the situation of recipe without the model with least performance.
model | pos | full | random | popular
---|---|---|---|---
baseline |  | 0.135 | 0.5573 | 0.0426
(cl4srec) |  | (0.0585) | (0.0513) |  (0.0466)
duorec | sup | 0.1346 | 0.5547 | 0.0454
| unsup | 0.1358 | 0.5594 | 0.0445
| - | 0.1351 | 0.554 | 0.0423
fisher (with) |  | 0.1386 | 0.5618 | 0.0428
fisher (w.o.) |  | 0.1373 | 0.5603 | 0.0487

###  A.3 Visualization of Sorted Probability

The figure displays the sorted probabilities of the top 50 items for 10
sequences, where single line represents single sequence. The cumulative
probability values for sample sizes of 10, 30, and 50 are 0.381, 0.569, and
0.658, respectively. With the exception of a few largest ones, the majority of
probabilities approximate $0$.

Figure 4: Sorted Probability $p_{\theta}\left(v_{j}|s_{i}\right)$.