Solano et al

# Confirmation and characterization of neglected WDS systems using Gaia DR3
and the Virtual Observatory

E. Solano  I. Novalbos  A. J. Ros  M. Cortés-Contreras  C. Rodrigo
Departmento de Astrofísica, Centro de Astrobiología (CSIC-INTA), ESAC Campus,
Camino Bajo del Castillo s/n
E-28692 Villanueva de la Cañada, Madrid, Spain  O.A.N.L. BCN, E-08032
Barcelona, Spain  Asociación Astronómica de Cartagena. C/ Mayor 148, E-30394
Canteras, Cartagena, Spain  Departamento de Física de la Tierra y Astrofísica.
Facultad Ciencias Físicas Universidad Complutense de Madrid. Plaza Ciencias,
1, E-28040 Madrid, Spain<EMAIL_ADDRESS>
(26 April 2016; 6 June 2016; 6 June 2016)

###### Abstract

The aim of this paper is, making use of the Gaia DR3 catalogue and Virtual
Observatory tools, to confirm and characterize 428 binary and multiple stellar
systems classified as neglected (only one observation) in the Washington
Double Star Catalogue (WDS). The components of the stellar systems have the
same parallax and proper motion (within the errors) and are separated by less
than 50 000 AU, which minimizes the number of by-chance counterparts.
Effective temperatures calculated using VOSA were used to estimate stellar
masses. Binding energies were calculated for 42 binary systems confirming they
are physical pairs. Also we found 75 pairs with F/G- M spectral types which
are very interesting to improve the determination of the metallicity of the M
star from the higher-mass component.

###### keywords:

binaries:general, stars:fundamental parameters, astronomical databases:
miscellaneous

††articletype: Article Type

##  1 Introduction

It is well known that, while they are young, stars are not isolated but
grouped in clusters, associations of stars loosely bound by mutual
gravitational attraction. Cluster members share physical parameters like age
and metallicity as well as kinematic properties (distances, proper motions and
radial velocities). Typically, after a few hundred million years, open
clusters become disrupted by close encounters with other clusters and clouds
of gas, as they orbit the Galactic center. As a remnant of this process, a
significant fraction of main-sequence stars (the exact percentage depends on
the spectral type) are in binary and multiple systems (Duquennoy Mayor
(1991),Raghavan . (2010),Cortés-Contreras . (2017)). The fact that the
components of wide binary and multiple systems share physical properties and,
at the same time, evolve in an independent way due to their large separation,
makes them excellent testbenchs for stellar evolutionary models. In this work
we aim to increase the number of these systems through examination of
historical data.

The Washington Double Star Catalogue (WDS, Mason ., 2001) is an all-sky
survey, maintained by the US Naval Observatory (USNO), that represents one of
the most important databases of binary and multiple stellar systems. At the
time of writing, the catalogue contains 154 686
rows111https://vizier.cds.unistra.fr/viz-bin/VizieR-3?-source=B/wds/wds. Among
other information, each WDS row includes the WDS name, right ascension and
declination of the primary component and the position angle and separation
between the primary and secondary components. The WDS catalogue also includes
a category, called neglected, to flag primaries that have been observed only
once, either because the information on coordinates is wrong or simply because
they have not been re-observed yet.

The Virtual Observatory (VO222http://www.ivoa.net) is an international
initiative aiming at optimizing the usage of the scientific information hosted
in astronomical archives. VO has developed tools and services which enormously
facilitate the access and analysis of astronomical data. In particular, Simbad
(Wenger ., 2000), Vizier (Ochsenbein ., 2000), TOPCAT (Taylor, 2005), Aladin
(Bonnarel ., 2000; Boch Fernique, 2014), and VOSA (Bayo ., 2008) have been
intensively used in this paper.

The paper is structured as follows. In Section 2, we describe the methodology
used to obtain the sample of objects studied in this paper together with the
results of their visual inspection. Physical parameters are estimated in
Section 3. The Virtual Observatory SED Analyzer tool
(VOSA333http://svo2.cab.inta-csic.es/theory/vosa/) was used to compute the
effective temperatures of our objects as well as to identify unresolved
binaries among them. Making use of the Gaia DR3 information on colours and
distances, the objects were placed on a colour - absolute magnitude diagram
(CMD) allowing to separate them into main sequence objects, subgiants/giants
or white dwarfs. Also, masses were estimated from effective temperatures and
used to calculate binding energies. Finally, in Section 4 we summarize the
main results of the paper. A brief description of the Virtual Observatory
compliant archive that contains detailed information on the candidates is
given in the Appendix.

##  2 Sample selection

The sample of objects analyzed in this paper was obtained after applying a
workflow which consists of the following steps:

  * •

Filtering on the number of observations: We kept only those primaries flagged
as neglected in the WDS catalogue, that is, with just one observation
(Nobs=1). This search reduced the number of rows from 154 686 to 18 314. Some
of these primaries have not been observed for many years (for instance, over
two hundred were observed more than a century ago).

  * •

Cross-matching: The 18 314 primary components obtained in the previous step
were cross-matched with Gaia DR3. Primary components not having counterparts
in Gaia at less than 5 arcsec were rejected. If there is more than one Gaia
counterpart at less than 5 arcsec, only the nearest one was considered. The 5
arcsec search radius was adopted as a compromise solution to avoid an
unmanageable number of false positives. We also forced the association to be
symmetric in the sense that the nearest object to the Gaia counterpart must
coincide with our original primary component. After cross-matching the number
of rows was reduced to 17 598.

  * •

Filtering on parallaxes and proper motions: We used the Gaia DR3 information
on parallaxes and proper motions to keep only primary components with relative
errors of less than 10 per cent in $\rm{PMRA}$ and $\rm{PMDEC}$ and less than
20 per cent in parallax. The condition in parallax is necessary to have a
reliable estimation of distance (Luri ., 2018). After this filtering, the
number of rows reduced to 14 808. For each one of these 14 808 entries, we
cross-matched the primary components with Gaia DR3 in a 180 arcmin radius,
keeping all the counterparts that also fulfil the conditions in the relative
errors in parallax and proper motion previously mentioned. We adopted this
large value for the search radius to keep all the physically bound pairs.
After this, we obtained 865 974 pairs.

  * •

Comoving pairs: From the 865 974 pairs, we kept only those for which the
differences in parallax and in proper motion, both in right ascension and
declination, were less than 3 times the corresponding errors. After applying
this condition, 2 735 pairs were kept. Each pair is formed by the WDS source
(primary component) and a Gaia counterpart (secondary component).

  * •

Physical separation: According to Jiménez-Esteban . (2019), the great majority
of chance alignment counterparts occurs at physical separations between
components larger than 50 000 AU. Using this number as an upper limit, 603
pairs were left. Physical separations were estimated by using the formulae
$s=\rho\times d$ , where $\rho$ is the angular separation between components
and $d$ is the distance to the pair (estimated from parallaxes). Moreover,
according to Gaia Collaboration . (20211), the minimum angular separation
above which a pair can be considered by Gaia as resolved is 180 mas.
Therefore, we did not impose any condition on the minimum separation between
components since this limit is much lower than what can be achieved from
ground-based observations.

  * •

RUWE. The Gaia Renormalised Unit Weight Error (RUWE) helps to identify objects
with problematic astrometric solution (Arenou ., 2018; Lindegren ., 2018,
2021). We adopted a conservative value of RUWE $<$ 1.4 to keep stars with good
astrometry. We found 429 pairs with RUWE $<$ 1.4 in both the primary and
secondary components.

  * •

Radial velocities: If a pair is physically bound, then, the primary and
secondary components should present similar radial velocities (RVs) within the
errors. However, as stated in Jiménez-Esteban . (2019), the radial velocity
dispersion in field stars is large and the probability of chance alignment
(same values within errors) is not negligible. Therefore, similar RVs cannot
be used as a proof to fully confirm that the pair is physically bound but, on
the contrary, it is a good approach to discard pairs with different RV values.

The Gaia DR3 catalog provides the RVs of both components for only 13 pairs.
One of them, WDS07222-2558, showed clearly discrepant RVs (35.38$\pm$1.65 km
s-1 and 0.16$\pm$3.91 km s-1 for the primary and secondary, respectively) and
was, thus, discarded.

After all this process, we ended up with 428 pairs comprising 354 primaries
and 372 secondaries (the same primary can be associated to more than one
secondary). The sky distribution of primary components is shown in Fig. 1
while Fig. 2 provides information on the distribution of the separation
between the primary and secondary components according to their distance.
Detailed information on these pairs can be found at the SVO archive of
neglected systems (see appendix).

Figure 1: Sky distribution in Galactic coordinates (Aitoff projection) of the
primary components of the 428 pairs identified in Section 2 (filled red
circles). A 2MASS coloured image is displayed in the background.

Figure 2: Separation (both physical and angular) vs distance of the 428 pairs
identified in Section 2.

###  2.1 Visual inspection

The 428 pairs were visually inspected taking advantage of the scripting
capabilities of Aladin444https://aladin.u-strasbg.fr/, an interactive software
sky atlas. For each pair we conducted the following steps:

  * •

Upload of an 5 arcmin image of the Second Palomar Sky Survey (Reid ., 1991,
POSS II) centered on the position of the primary component. The POSS II image
is used as background source.

  * •

Identification of the primary component using the WDS coordinates.

  * •

Identification of the secondary component using the information on separation
and position angle available in the WDS catalogue.

  * •

Upload of the Gaia DR3 sources lying in the region of the sky covered by the
POSS II image.

  * •

Graphical representation of the proper motion of the Gaia sources using
arrows.

An example of the output of the script is given in Fig. 3. During the visual
inspection of the 428 pairs the following cases were identified:

  * •

The WDS information on coordinates of the primary component, position angle
and separation is correct but it was superseded by the astrometric information
provided by Gaia due to its superior accuracy (Fig. 3). We found 32 pairs
lying in this category.

  * •

The WDS information is wrong. There is no Gaia DR3 source at the
separation/position angle given in WDS. The secondary component is found at a
different separation and/or position angle (Fig. 4). 189 pairs belong to this
category. For 127 pairs, the primary component is flagged in Simbad as close
binary itself. However, in all cases, they have associated a single entry in
Gaia DR3 with RUWE < 1.4, indicating that they are, most likely, single.

  * •

For eight pairs there is a Gaia DR3 source at the expected position of the
secondary according to the separation/position angle given in WDS.
Nevertheless, the parallax/proper motion of the Gaia DR3 source lying at the
expected position of the secondary are different by more than 3$\sigma$ from
the parallax/proper motion of the primary component. The secondary component
is, actually, found at a different separation/position angle.

  * •

35 of our pairs belong to a multiple system according to WDS, that is, each
one of the 35 primaries has associated more than one secondary. For three of
the primaries we found that none of the secondaries reported in WDS have
similar parallaxes and proper motions and, thus, are not physically
associated. For 29 primaries only one of the secondaries reported in WDS can
be considered as physical pair (Fig. 5), while we confirm that three primaries
belong to a multiple system. Finally, we found three primaries, reported as
double in WDS but that, actually, are triple systems according to Gaia
parallaxes and proper motions.

  * •

The primary is part of a larger structure like an open cluster or a stellar
association (Fig. 6). 161 primaries belong to this category.

The previously listed categories are properly flagged in the archive (see
Appendix). Also, the results obtained after the visual inspection will be
reported to USNO for their ingestion in the WDS catalogue.

Figure 3: Example of a primary (WDS 00243+5753) for which the information
provided by the WDS on separation and position angle is correct (the small red
arc indicates the expected position of the secondary according to the WDS
information), information that is, anyway, superseded by the Gaia parameters
due to its superior performance. Blue arrows represent the Gaia DR3 proper
motions.

Figure 4: Example of a primary (WDS 02199+5236) for which the information
provided by the WDS on separation and/or position angle is wrong. The
secondary is found at a different separation/position angle (at the center of
the circle that appears furthest to the right). The expected position of the
secondary according to the WDS information is marked by a small red arc at the
North-East of the crosshair. Blue arrows represent the Gaia DR3 proper
motions.

Figure 5: Example of a triple system in WDS (WDS 04241+2413) for which one of
the secondaries (small red arc below the crosshair) is not physically bound
based on its Gaia proper motion information (blue arrows).

Figure 6: Example of a primary component (WDS 03447+3206, marked with a cross
and a red circle) member of a stellar cluster (IC 348).

##  3 Physical parameters

###  3.1 Effective temperatures

We used VOSA to estimate effective temperatures for our pairs. VOSA is a tool
developed by the Spanish Virtual Observatory designed to build the Spectral
Energy Distributions (SEDs) of thousands of objects at a time from a large
number of photometric catalogues, ranging from the ultraviolet to the
infrared. VOSA compares catalogue photometry with different collections of
theoretical models and determines which model best reproduces the observed
data following different statistical approaches. Physical parameters can then
be estimated for each object from the model that best fits the data.

Using VOSA we queried the GALEX (Bianchi ., 2017), SDSS DR12 (Alam ., 2015),
APASS DR9 (Henden ., 2015), Gaia DR3 (Gaia Collaboration ., 20211), 2MASS-PSC
(Skrutskie ., 2006), and ALLWISE (Wright ., 2010) catalogs to build the SEDs
from the ultraviolet to the infrared. Observational SEDs were then compared to
the grid of BT-Settl model atmospheres (Allard ., 2012). We assumed 1 200 K
$\leq$ $T_{\rm eff}$ $\leq$ 12 000 K; 3 $\leq$ logg $\leq$ 4.5 and solar
metallicity.

Extinction can play an important role in shaping the SED in particular at
short wavelengths. If extinction is underestimated, the slope of the SED will
appear flattened at short wavelengths and the derived effective temperature
will be lower. To account for this effect and to minimize the extinction -
effective temperature degeneracy in the SED fitting, we decided to leave
extinction as a free parameter in the fitting process taking values in the
range 0 $\leq$ Av $\leq$ 1 mag and keep only objects at distances $<$ 1000 pc
(1 Kpc roughly corresponds to an extinction of 1 mag in the optical regime.
See, for instance, Fig. 8 in Lallement . (2022)). Therefore, effective
temperatures were estimated only for this subset of objects.

The goodness of fit of the SED in VOSA can be assessed with the vgfb
parameter, a pseudo-reduced $\chi^{2}$ internally used by VOSA that is
calculated by forcing $\sigma(F_{\rm obs})>0.1\times F_{\rm obs}$, where
$\sigma(F_{\rm obs})$ is the error in the observed flux ($F_{\rm obs}$). This
parameter is useful to avoid the risk of overweighting photometric points with
under-estimated photometric errors. Only sources with vgfb < 15 (which is an
indicator of good fit) were kept.

VOSA also allows the identification of flux excess in a SED, excess that could
be ascribed to the presence of a disc or the existence of a non-resolved
companion. This way, we can check whether each of the primaries and
secondaries are, themselves, single or double objects.A detailed description
of how VOSA manages the flux excess can be found in the VOSA
documentation555https://bit.ly/2KRCv9x. Out of 354 primaries, VOSA did not
find excess for 289 of them (216 at less than 1000 pc), while 34 primaries
show flux excess, 24 of which at less than 1 Kpc. The rest of primaries (31)
showed a bad SED fitting due to, for instance, poor quality photometry or lack
of enough photometric points. Similarly, for the 372 secondaries, VOSA did not
find excess for 256 of them (189 at less than 1000 pc), while 43 secondaries
show flux excess (32 at less than 1 Kpc). 73 secondaries were discarded
because of their poor SED fitting. Physical parameters (effective temperature,
luminosity, stellar radius) of the 216 primaries and 189 secondaries at less
than 1000 pc and not showing flux excess can be found at the SVO archive of
neglected systems (see Appendix). Examples of the VOSA SED fitting are shown
in Fig. 7.

Figure 7: Top: Example of a VOSA SED fitting showing no flux excess
(WDS16219-6729). The red dots represent the observed photometry while the blue
spectrum indicates the theoretical model that fits best. Bottom: Example of a
VOSA SED fitting showing flux excess (WDS06290+3717). As before, red dots
represent the observed photometry, purple and green dots the synthetic
photometry of the models that best fits and blue dots the composite synthetic
photometry. The yellow inverted triangles indicate photometric values
corresponding to upper limits. They are not considered in the fitting process.

Fig. 8 shows the distribution in effective temperature of the primaries and
secondaries classified by VOSA as single stars (i.e., not showing flux excess
in their SEDs). The primaries reach the maximum of the distribution at $\sim$
6 500 K while the secondaries have it at $\sim$ 3 500 K with $\sim$ 50% of
them in the range 3 000 $\leq$ $T_{\rm eff}$ $\leq$ 4 000 K. Of special
interest are the 75 pairs composed of a primary with F-G spectral types (5 300
K $<$ $T_{\mathrm{eff}}$$<$ 7 300 K, according to the updated version of Table
4 in Pecaut Mamajek
(2013)666https://www.pas.rochester.edu/~emamajek/EEM_dwarf_UBVIJHK_colors_Teff.txt)
and an M-dwarf secondary ($T_{\mathrm{eff}}$$<$ 3 900 K) as the metallicity of
the M-dwarf can be inferred from the hotter component.

Figure 8: Distribution of the effective temperatures of the primaries (blue)
and secondaries (pink) classified as single objects by VOSA.

###  3.2 HR diagram

The absolute Gaia magnitude in the $G$ band was estimated using

$M_{G}=G+5\log{\varpi}+5,$ (1)

where $\varpi$ is the parallax in arcseconds and $G$ the apparent magnitude.
With the absolute magnitude and the BP-RP colour, we built a colour - absolute
magnitude diagram (CMD). Fig. 9 shows the position of our pairs in the CMD.

Among our pairs we identified a white dwarf (Gaia DR3 150264795265471616)
already reported in Gentile Fusillo . (2019) and two sources (Gaia DR3
4009115472737321600 and Gaia DR3 6409528409064496000) lying in the locus
occupied by the white dwarf - main sequence binaries Rebassa-Mansergas .
(2021).

Figure 9: Colour-magnitude diagram built using Gaia DR3 sources with
parallaxes larger than 10 mas, photometric errors in $G$, $G_{\mathrm{BP}}$
and $G_{\mathrm{RP}}$ less than 10 per cent and RUWE $<$ 1 (black dots). On
top of it we have plotted the primary components classified by VOSA as singles
(green bullets) and unresolved binaries (red bullets). Likewise, the secondary
components classified as singles and unresolved binaries are plotted in pink
and blue crosses, respectively . The correction by the Gaia colour excess was
applied following the approach described in Riello . (2021).

###  3.3 Binding energies

Stellar masses of objects lying on the main sequence were derived from
effective temperatures by interpolating in Table 4 in Pecaut Mamajek (2013).
With these values and the projected physical separations we computed reduced
binding energies as in Caballero (2009)

$U_{g}=-GM_{1}M_{2}/s$ (2)

where G is the gravitational constant, M1 and M2 the masses of the primary and
the secondary and s the projected angular separation. Fig. 10 shows the
binding enery - total mass distribution for the 42 pairs with mass
determinations for both components. We consider a physically bound pair when
binding energy is over 1033 J (Caballero, 2009). All the pairs show binding
energies well above this threshold.

Figure 10: Binding energy-total mass diagram. The red line indicates the
minimum energy above which the pair is supposed to be bound.

##  4 Conclusions

Starting with Washington Double Stars catalogue (154 686 rows), we selected
those primary stars having just one observation, and look for counterparts
(secondaries) sharing the same Gaia parallaxes and proper motions (and radial
velocities, whenever available) within the errors, and a RUWE value typical of
objects with good astrometrical solution. This returned 428 pairs that were
visually inspected using Aladin to check the WDS information on separation and
position angle.

Effective temperatures, luminosities and radii of both the primaries and the
secondaries were estimated using VOSA. In order to minimize the impact of the
effective temperature - extinction degeneracy, physical parameters were
derived only for objects at less than 1 Kpc. VOSA also allows to identify
unresolved binaries by identifying the flux excess in the SED distribution.
This way, out of 354 primaries, VOSA classified 289 as singles and 34 as
unresolved binaries according to their SED. Similarly, for the 372
secondaries, VOSA identified 256 as singles and 43 as unresolved binaries. We
were also able to identify 75 M dwarf + F/G pairs for the subsample for which
the primary and secondary are both single stars. These pairs are very helpful
to accurately estimate the metallicity of the M dwarf from the hotter
companion.

Using Gaia DR3 parallaxes and magnitudes, we plotted the 428 pairs on a
Hertzsprung-Russell diagram (HRD). According to their position in the HRD we
found a white dwarf and two sources lying in the region of the parameter space
typically occupied by the white dwarf - main sequence binaries. Finally, we
computed the binding energies for 42 pairs finding that all of them are
consistent with being gravitationally bound.

Detailed information of the pairs can be found in the Virtual Observatory
compliant archive described in the Appendix.

## Acknowledgments

This work has been funded by MCIN/AEI/10.13039/501100011033 through grant
PID2020-112949GB-I00 at Centro de Astrobiología (CSIC-INTA). This research
used the Washington Double Star Catalog maintained at the U.S. Naval
Observatory. This publication makes use of VOSA, developed under the Spanish
Virtual Observatory project. This research has made use of Aladin sky atlas
developed at CDS, Strasbourg Observatory, France. Vizier, Simbad, and TOPCAT
have also been widely used in this paper. This work has made use of data from
the European Space Agency (ESA) mission Gaia
(https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and
Analysis Consortium (DPAC,
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

##  Appendix A Virtual Observatory compliant, online catalogue

In order to help the astronomical community with using our catalogue of
neglected WDS objects, we developed an archive system that can be accessed
from a webpage777http://svocats.cab.inta-csic.es/wds_neglected/ or through a
Virtual Observatory ConeSearch888e.g.http://svocats.cab.inta-
csic.es/wds_list/cs.php?RA=31.825&DEC=-7.905&SR=0.1&VERB=2.

The archive system implements a very simple search interface that allows
queries by coordinates and radius as well as by other parameters of interest.
The user can also select the maximum number of sources (with values from 10 to
unlimited). The result of the query is a HTML table with all the sources found
in the archive fulfilling the search criteria. The result can also be
downloaded as a VOTable or a CSV file. Detailed information on the output
fields can be obtained placing the mouse over the question mark located close
to the name of the column. The archive also implements the
SAMP999http://www.ivoa.net/documents/SAMP (Simple Application Messaging)
Virtual Observatory protocol. SAMP allows Virtual Observatory applications to
communicate with each other in a seamless and transparent manner for the user.
This way, the results of a query can be easily transferred to other VO
applications, such as, for instance, Topcat.

The query syntaxes to recover the different subsets identified in this work
are the following:

  * •

32 pairs for which the WDS information has been superseded by the Gaia
astrometry: http://svocats.cab.inta-csic.es/wds_list/index.php?action=search
and flag=0

  * •

189 pairs for which the WDS information on separation and/or position angle is
wrong: http://svocats.cab.inta-csic.es/wds_list/index.php?action=search and
flag=11

  * •

8 pairs for which there is a Gaia DR3 source at the separation/position angle
given in WDS but with a different parallax/proper motion:
http://svocats.cab.inta-csic.es/wds_list/index.php?action=search and flag=22

  * •

3 primaries in WDS for which none of the secondaries are physical pairs
according to Gaia DR3 parallaxes and proper motions: http://svocats.cab.inta-
csic.es/wds_list/index.php?action=search and flag=550

  * •

29 primaries for which only one of the WDS secondaries can be considered as a
physical pair according to Gaia DR3 parallaxes and proper motions:
http://svocats.cab.inta-csic.es/wds_list/index.php?action=search and flag=551

  * •

3 primaries, reported as part of double systems in WDS, but that actually
belong to triple systems according to Gaia: http://svocats.cab.inta-
csic.es/wds_list/index.php?action=search and flag=552

  * •

Physical parameters (effective temperature, luminosity, stellar radius) of the
216 primaries at less than 1000 pc: http://svocats.cab.inta-
csic.es/wds_primary

  * •

Physical parameters (effective temperature, luminosity, stellar radius) of the
189 secondaries at less than 1000 pc: http://svocats.cab.inta-
csic.es/wds_secondary

  * •

Effective temperatures of the 75 FG + M pairs: http://svocats.cab.inta-
csic.es/wds_fgkm

  * •

Binding energies for the 42 pairs with mass determinations for both
components: http://svocats.cab.inta-csic.es/wds_binding/

## References

  * Alam . (2015) Alam15Alam, S., Albareti, FD., Allende Prieto, C. et al.  201507,  ApJS219112.  10.1088/0067-0049/219/1/12
  * Allard . (2012) Allard12Allard, F., Homeier, D. Freytag, B.  201206,  Philosophical Transactions of the Royal Society of London Series A37019682765-2777.  10.1098/rsta.2011.0269
  * Arenou . (2018) Arenou2018Arenou, F., Luri, X., Babusiaux, C. et al.  201808,  A&A616A17.  10.1051/0004-6361/201833234
  * Bayo . (2008) Bayo08Bayo, A., Rodrigo, C., Barrado Y Navascués, D., Solano, E., Gutiérrez, R., Morales-Calderón, M. Allard, F.  200812,  A&A492277-287.  10.1051/0004-6361:200810395
  * Bianchi . (2017) Bianchi17Bianchi, L., Shiao, B. Thilker, D.  201706,  ApJS230224.  10.3847/1538-4365/aa7053
  * Boch Fernique (2014) Boch14Boch, T. Fernique, P.  201405,  Aladin Lite: Embed your Sky in the Browser Aladin Lite: Embed your Sky in the Browser.  N. Manset P. Forshay (), Astronomical Data Analysis Software and Systems XXIII Astronomical Data Analysis Software and Systems XXIII 485, 277.
  * Bonnarel . (2000) Bonnarel00Bonnarel, F., Fernique, P., Bienaymé, O. et al.  200004,  A&AS14333-40.  10.1051/aas:2000331
  * Caballero (2009) Caballero09Caballero, JA.  200911,  A&A5071251-259.  10.1051/0004-6361/200912596
  * Cortés-Contreras . (2017) CC17Cortés-Contreras, M., Béjar, VJS., Caballero, JA. et al.  201701,  A&A597A47.  10.1051/0004-6361/201629056
  * Duquennoy Mayor (1991) Duq91Duquennoy, A. Mayor, M.  199108,  A&A248485.
  * Gaia Collaboration . (20211) Gaiacollab21Gaia Collaboration, Brown, AGA., Vallenari, A. et al.  2021105,  A&A649A1.  10.1051/0004-6361/202039657
  * Gaia Collaboration . (20212) Gaia21Gaia Collaboration, Brown, AGA., Vallenari, A. et al.  2021205,  A&A649A1.  10.1051/0004-6361/202039657
  * Gentile Fusillo . (2019) Fusillo19Gentile Fusillo, NP., Tremblay, PE., Gänsicke, BT. et al.  201902,  MNRAS48244570-4591.  10.1093/mnras/sty3016
  * Henden . (2015) Henden15Henden, AA., Levine, S., Terrell, D. Welch, DL.  201501,  APASS - The Latest Data Release APASS - The Latest Data Release.  American Astronomical Society Meeting Abstracts #225 American Astronomical Society Meeting Abstracts #225 225, 336.16.
  * Jiménez-Esteban . (2019) Fran19Jiménez-Esteban, FM., Solano, E. Rodrigo, C.  201902,  AJ157278.  10.3847/1538-3881/aafacc
  * Lallement . (2022) Lallement22Lallement, R., Vergely, JL., Babusiaux, C. Cox, NLJ.  202205,  A&A661A147.  10.1051/0004-6361/202142846
  * Lindegren . (2018) Lindegren2018Lindegren, L., Hernández, J., Bombrun, A. et al.  201808,  A&A616A2.  10.1051/0004-6361/201832727
  * Lindegren . (2021) Lindegren2021Lindegren, L., Klioner, SA., Hernández, J. et al.  202105,  A&A649A2.  10.1051/0004-6361/202039709
  * Luri . (2018) Luri18Luri, X., Brown, AGA., Sarro, LM. et al.  201808,  A&A616A9.  10.1051/0004-6361/201832964
  * Makarov . (2008) Makarov08Makarov, VV., Zacharias, N. Hennessy, GS.  200811,  ApJ6871566-578.  10.1086/591638
  * Mason . (2001) Mason01Mason, BD., Wycoff, GL., Hartkopf, WI., Douglass, GG. Worley, CE.  200112,  AJ12263466-3471.  10.1086/323920
  * Ochsenbein . (2000) Ochsenbein00Ochsenbein, F., Bauer, P. Marcout, J.  200004,  A&AS14323-32.  10.1051/aas:2000169
  * Pecaut Mamajek (2013) pecaut13Pecaut, MJ. Mamajek, EE.  201309,  ApJS20819.  10.1088/0067-0049/208/1/9
  * Raghavan . (2010) Rag10Raghavan, D., McAlister, HA., Henry, TJ. et al.  201009,  ApJS19011-42.  10.1088/0067-0049/190/1/1
  * Rebassa-Mansergas . (2021) alberto21Rebassa-Mansergas, A., Solano, E., Jiménez-Esteban, FM. et al.  202110,  MNRAS50645201-5211.  10.1093/mnras/stab2039
  * Reid . (1991) Reid91Reid, IN., Brewer, C., Brucato, RJ. et al.  199107,  PASP103661.  10.1086/132866
  * Riello . (2021) Riello21Riello, M., De Angeli, F., Evans, DW. et al.  202105,  A&A649A3.  10.1051/0004-6361/202039587
  * Skrutskie . (2006) Skrutskie06Skrutskie, MF., Cutri, RM., Stiening, R. et al.  200602,  AJ13121163-1183.  10.1086/498708
  * Taylor (2005) Taylor05Taylor, MB.  200512,  TOPCAT & STIL: Starlink Table/VOTable Processing Software TOPCAT & STIL: Starlink Table/VOTable Processing Software.  P. Shopbell, M. Britton R. Ebert (), Astronomical Data Analysis Software and Systems XIV Astronomical Data Analysis Software and Systems XIV 347, 29.
  * Tokovinin (2014) Tokovinin14Tokovinin, A.  201404,  AJ147487.  10.1088/0004-6256/147/4/87
  * Wenger . (2000) Wenger00Wenger, M., Ochsenbein, F., Egret, D. et al.  200004,  A&AS1439-22.  10.1051/aas:2000332
  * Wright . (2010) Wright10Wright, EL., Eisenhardt, PRM., Mainzer, AK. et al.  201012,  AJ14061868-1881.  10.1088/0004-6256/140/6/1868


# INTRODUCTION TO DATA ASSIMILATION
ALGORITHMS FOR PARAMETER ESTIMATION

Loc Luong  New Mexico Tech

## Abstract

In this study, two classes of methods including statistical and variational
data assimilation algorithms will be described. In statistical methods, the
model state is updated sequentially based on the previous estimate.
Variational methods, on the other hand, seek an estimation in space and time
by minimizing a cost function. Both of these methods require estimates of
background state which is the prior information of the system and its error
covariances. In terms of linear and Gaussian problems, they have the same
solution. In the family of Kalman Filter algorithms, the conventional Kalman
Filter (KF) and Ensemble Kalman Filter (EnKF) will be implemented. A three-
dimension variational method (3D-Var) will be employed to illustrate the
variational approaches. A simple case of an ordinary differential equation
(ODE) is coupled to highlight the difference between these algorithms. Namely,
a mass-spring system governed by a second order differential equation will be
examined. We also look at the situation where a periodic external force is
applied to the system. Then different data assimilation algorithms will be
applied to this system. Results from the experiments will be analyzed to
showcase the advantages and disadvantages of each method.

##  1 Introduction

The process of integrating model predictions with actual system measurements
to provide an optimal approximation of the current state of the system or
predictions of the future of the system is known as data assimilation. In
practice, the observations or measurements usually contain noise. It is
important to incorporate noisy data in efficient way to provide more accurate
estimates of system state. And data assimilation is by far one of the most
efficient way to tackle such tasks. This technique combines different types of
information which are forecasting data and observation of state to minimize
prediction errors. The method initially developed in weather forecasting, and
later gain a great deal of attention to other fields of science.
Over time, data assimilation became very popular in many fields of science
such as geoscience, petroleum engineering, seismology, medicine, control
theory, etc. Data assimilation can be classified in two classes: statistical
DA and variational DA. In statistical aspect, we seek an optimal solution with
minimum variance, whereas in variational DA we seek an optimal solution that
minimizes a cost function. Kalman filter is one of the simplest method in data
assimilation, and is first developed by R. E. Kalman.

##  2 Methodology

Consider a dynamic system governed by the following equation:

$x_{k}=Ax_{k-1}+Bu_{k-1}+w_{k-1}$ (1)

where $x_{k}$ is the state vector of the process at time k, A is the state
transition model which get the previous state vector $x_{k-1}$ as input, B is
the control input model which get the previous control variable $u_{k-1}$ as
input, $w_{k-1}$ is the process noise that is assumed to have a multivariate
normal distribution with mean 0 and a covariance matrix Q.
The measurements or observations $z_{k}$ at time k is modeled by:

$z_{k}=Hx_{k}+v_{k}$ (2)

where H is the observation model, and is assumed stationary over time, $v_{k}$
is the associated observation error that is assumed to be normally distributed
with mean 0 and covariance matrix R.
We want to estimate the state of the system in the sense of minimizing the
norm of $\mathinner{\\!\left\lVert z_{k}-Hx_{k}\right\rVert}$ as well as the
uncertainty. The Kalman filter helps to model the prior values (background
information) of the state and measured values (observations) into the
following form:

$\hat{x}_{k}=\hat{x}^{\prime}_{k}+K_{k}(z_{k}-H\hat{x}^{\prime}_{k})$ (3)

where $\hat{x}^{\prime}_{k}$ is the prior estimate of $x_{k}$, $K_{k}$ is the
Kalman gain which will be derived later. The second term
$z_{k}-H\hat{x}^{\prime}_{k}$ is known as the innovation.
The covariance of the new update is:

$P_{k}=cov(\hat{x}_{k})$ (4)
$\begin{split}P_{k}&=cov((I-K_{k}H)\hat{x}^{\prime}_{k}+K_{k}z_{k})\\\
&=(I-K_{k}H)cov(\hat{x}^{\prime}_{k})(I-K_{k}H)^{T}+K_{k}cov(z_{k})K_{k}^{T}\\\
&=(I-K_{k}H)P^{\prime}_{k}(I-K_{k}H)^{T}+K_{k}cov(z_{k})K_{k}^{T}\\\
&=(I-K_{k}H)P^{\prime}_{k}(I-K_{k}H)^{T}+K_{k}RK_{k}^{T}\\\ \end{split}$ (5)

where $P^{\prime}_{k}$ is the prior estimate of $P_{k}$.
Expand the covariance $P_{k}$, we can have:

$\begin{split}P_{k}&=P^{\prime}_{k}-K_{k}HP^{\prime}_{k}-P^{\prime}_{k}H^{T}K_{k}^{T}+K_{k}HP^{\prime}_{k}H^{T}K_{k}^{T}+K_{k}RK_{k}^{T}\\\
&=P^{\prime}_{k}-K_{k}HP^{\prime}_{k}-P^{\prime}_{k}H^{T}K_{k}^{T}+K_{k}(HP^{\prime}_{k}H^{T}+R)K_{k}^{T}\end{split}$
(6)

The sum of diagonal elements of the covariance matrix $P_{k}$ is total
variation of the system. Therefore, we want to minimize the trace of the
covariance matrix $P_{k}$.
Taking the trace operator of both sides, we have:

$\begin{split}Tr(P_{k})&=Tr(P^{\prime}_{k})-2Tr(K_{k}HP^{\prime}_{k})+Tr(K_{k}(HP^{\prime}_{k}H^{T}+R)K_{k}^{T})\\\
\end{split}$ (7)

Using vector calculus, the derivative with respect to $K_{k}$ can be expressed
as:

$\begin{split}\frac{\partial T(P_{k})}{\partial
K_{k}}&=-2(HP^{\prime}_{k})^{T}+2K_{k}(HP^{\prime}_{k}H^{T}+R)\end{split}$ (8)

Setting the derivative to be 0, we have:

$\begin{split}-2(HP^{\prime}_{k})^{T}+2K_{k}(HP^{\prime}_{k}H^{T}+R)=0\end{split}$
(9)

This leads to:

$\begin{split}K_{k}&=(HP^{\prime}_{k})^{T}(HP^{\prime}_{k}H^{T}+R)^{-1}\\\
K_{k}&=P^{\prime}_{k}H^{T}(HP^{\prime}_{k}H^{T}+R)^{-1}\\\ \end{split}$ (10)

Then, we can substitute the Kalman gain $K_{k}$ to simplify $P_{k}$ as
follows:

$\begin{split}P_{k}&=P^{\prime}_{k}-K_{k}HP^{\prime}_{k}-P^{\prime}_{k}H^{T}K_{k}^{T}+K_{k}(HP^{\prime}_{k}H^{T}+R)K_{k}^{T}\\\
&=P^{\prime}_{k}-P^{\prime}_{k}H^{T}(HP^{\prime}_{k}H^{T}+R)^{-1}HP^{\prime}_{k}-P^{\prime}_{k}H^{T}(P^{\prime}_{k}H^{T}(HP^{\prime}_{k}H^{T}+R)^{-1})^{T}+\\\
&+P^{\prime}_{k}H^{T}(HP^{\prime}_{k}H^{T}+R)^{-1}(HP^{\prime}_{k}H^{T}+R)(P^{\prime}_{k}H^{T}(HP^{\prime}_{k}H^{T}+R)^{-1})^{T}\\\
&=P^{\prime}_{k}-P^{\prime}_{k}H^{T}(HP^{\prime}_{k}H^{T}+R)^{-1}HP^{\prime}_{k}-P^{\prime}_{k}H^{T}(P^{\prime}_{k}H^{T}(HP^{\prime}_{k}H^{T}+R)^{-1})^{T}+\\\
&+P^{\prime}_{k}H^{T}(P^{\prime}_{k}H^{T}(HP^{\prime}_{k}H^{T}+R)^{-1})^{T}\\\
&=P^{\prime}_{k}-P^{\prime}_{k}H^{T}(HP^{\prime}_{k}H^{T}+R)^{-1}HP^{\prime}_{k}\\\
&=P^{\prime}_{k}-K_{k}HP^{\prime}_{k}\\\
&=(I-K_{k}H)P^{\prime}_{k}\end{split}$ (11)

Finally, the Kalman Filter algorithm can be summarized as follows:

Algorithm 1 Kalman Filter algorithm

Input: A, B, H, Q, R, $P_{0}$, $x_{0}$

1:For k = 1:k

2: $x^{\prime}_{k}=Ax_{k-1}+Bu_{k}$

3: $P^{\prime}_{k}=AP_{k-1}A^{T}+Q$

4: $K_{k}=P^{\prime}_{k}H^{T}(HP^{\prime}_{k}H^{T}+R)^{-1}$

5: $x_{k}=x^{\prime}_{k}+K_{k}(z_{k}-Hx^{\prime}_{k}$.

6: $P_{k}=(I-K_{k}H)P^{\prime}_{k}$

7:End For

Output: $x_{k}$, $P_{k}$

###  2.1 Ensemble Kalman Filter

In the Ensemble Kalman Filter, we use Monte Carlo simulation to generate a
collection of the state variable which represent the uncertainty of the
system. Firstly, using multivariate normal distribution we can initialize the
ensemble for each random vector and perform the estimate update in a way
similar to Kalman Filter:

$x_{k}^{i}=\hat{x}_{k}^{i}+K(y-H\hat{x}_{k}^{i})$ (12)

where $\hat{x}_{k}^{i}$ is the forecast state vector, and is estimated as:

$\hat{x}_{k}^{i}=M(\hat{x}_{k-1}^{i})$ (13)

where M is the system operators that can be linear or non-linear.
In Kalman Filter, the Kalman gain is:

$\begin{split}K_{k}&=P^{\prime}_{k}H^{T}(HP^{\prime}_{k}H^{T}+R)^{-1}\\\
\end{split}$ (14)

We can estimate the Kalman gain using the statistic data from ensemble
vectors. Firstly, the forecast covariance $P^{\prime}_{k}$ can be estimated
as:

$P^{\prime}_{k}=\frac{1}{m-1}\sum_{i=1}^{m}(\hat{x}_{k-1}^{i}-\bar{x})(\hat{x}_{k-1}^{i}-\bar{x})^{T}$
(15)

where $\bar{x}$ is the mean of ensemble vectors:

$\bar{x}=\frac{1}{m}\sum_{i=1}^{m}\hat{x}_{k-1}^{i}$ (16)

Then, the Kalman gain can be calculated similarly to KF using equation 14.
The Ensemble Kalman Filter algorithm can be summarized as follows:

Algorithm 2 Ensemble Kalman Filter algorithm

Input: Observation model H, Observation covariance matrix R, system operator H

Initialize ensemble vectors $\hat{x}_{k}^{i}$, i=1,2..,m.

2:For k = 1:k

Perturbation the observation with noise $z_{i,k}=z_{i}+u_{i}$,
$u_{i}\sim\mathcal{N}(0,R)$

4: Compute ensemble mean and covariance

Compute Kalman gain

6: Update ensemble vectors

Compute ensemble forecast

8:End For

Output: $x_{k}$, $P_{k}$

###  2.2 3D-Variational method

The key concept behind variational methods is to strike a balane between the
background (historical state) and the measured values (observations). In terms
of variational approach, we seek for solutions that are not far away from the
observation by minimizing a cost function. The cost function is:

$J(x)=\frac{1}{2}(x-x^{b})^{T}D^{-1}(x-x^{b})+\frac{1}{2}(y-Hx)^{T}R^{-1}(y-Hx)$
(17)

where R and D are the observation and background error covariance matrices,
respectively. H is the observation operator.
The gradient of J is:

$\nabla J=D^{-1}(x-x^{b})-H^{T}R^{-1}(y-Hx)$ (18)

Solving equation 18, we can have:

$\begin{split}D^{-1}(x-x^{b})&=H^{T}R^{-1}(y-Hx)\\\
D^{-1}x+H^{T}R^{-1}Hx&=H^{T}R^{-1}y+D^{-1}x^{b}\\\
x&=(D^{-1}+H^{T}R^{-1}H)^{-1}(H^{T}R^{-1}y+D^{-1}x^{b})\\\
x&=(D^{-1}+H^{T}R^{-1}H)^{-1}(D^{-1}x^{b}+H^{T}R^{-1}Hx^{b}+\\\
&+H^{T}R^{-1}y-H^{T}R^{-1}Hx^{b})\\\
x&=(D^{-1}+H^{T}R^{-1}H)^{-1}((D^{-1}+H^{T}R^{-1}H)x^{b}+\\\
&+H^{T}R^{-1}(y-Hx^{b})\\\ x&=x^{b}+K(y-Hx^{b})\\\ \end{split}$

Finally, we have:

$\begin{split}x&=x^{b}+K(y-Hx^{b})\end{split}$ (19)

where the gain matrix K is given by:

$\begin{split}K=(D^{-1}+H^{T}R^{-1}H)^{-1}H^{T}R^{-1}\end{split}$

Equation 19 express the solution of 3D-Variational method. The solution is
very similar to the Kalman filter form (Equation 3) where we model the prior
values of the state and measured values. The gain matrix K is weighted average
of the background (or prior) and the observations of the state.
The background covariance matrix D describe the errors covariances of model
variables. The matrix D in particular must be symmetric positive definite, and
it must be reasonable in expressing observation parameters. In practice, the
matrices involved in equation 19 are often in large dimensions, thus the
direct calculation is unfeasible. Therefore, to minimize this cost function of
3D-Variation method, gradient descent (GD) is often implemented. The 3D-Var
algorithm is given as:

Algorithm 3 3D-Var algorithm

Initialization: $x_{0}=x^{b}$, n =0

1:While $\mathinner{\\!\left\lVert\nabla J\right\rVert}>\epsilon$ or $n\leq
n_{max}$

2: Compute J

3: Compute $\nabla J$

4: Update $x_{j+1}$

5: j = j+1

6:End While

Output: $x_{k}$

By using a change of variable, we can reformulate the variational cost
function into Tikhonov regularization problem.
Let $D=\sigma_{D}^{2}C_{D}$, $R=\sigma_{R}^{2}C_{R}$, and ignore the constant
term $\frac{1}{2}$ in the cost function, we have:

$\begin{split}J(x)&=(x-x^{b})^{T}D^{-1}(x-x^{b})+(y-Hx)^{T}R^{-1}(y-Hx)\\\
&=\frac{1}{\sigma_{D}^{2}}(x-x^{b})^{T}C_{D}^{-1/2}C_{D}^{-1/2}(x-x^{b})+\frac{1}{\sigma_{R}^{2}}(y-Hx)^{T}C_{R}^{-1/2}C_{R}^{-1/2}(y-Hx)\end{split}$
(20)

Since $D$ and $R$ are symmetric and positive definite, we obtain:

$\begin{split}J(x)&=\frac{1}{\sigma_{D}^{2}}\mathinner{\\!\left\lVert
C_{D}^{-1/2}(x-x^{b})\right\rVert}_{2}^{2}+\frac{1}{\sigma_{R}^{2}}\mathinner{\\!\left\lVert
C_{R}^{-1/2}(y-Hx)\right\rVert}_{2}^{2}\\\
&=\frac{1}{\sigma_{D}^{2}}\mathinner{\\!\left\lVert
C_{D}^{-1/2}(x-x^{b})\right\rVert}_{2}^{2}+\frac{1}{\sigma_{R}^{2}}\mathinner{\\!\left\lVert
C_{R}^{-1/2}(y-Hx^{b})-C_{R}^{-1/2}(Hx-Hx^{b})\right\rVert}_{2}^{2}\\\
&=\frac{1}{\sigma_{D}^{2}}\mathinner{\\!\left\lVert
C_{D}^{-1/2}(x-x^{b})\right\rVert}_{2}^{2}+\frac{1}{\sigma_{R}^{2}}\mathinner{\\!\left\lVert
C_{R}^{-1/2}(y-Hx^{b})-C_{R}^{-1/2}HC_{D}^{1/2}C_{D}^{-1/2}(x-x^{b})\right\rVert}_{2}^{2}\\\
\end{split}$ (21)

By letting:

$\begin{split}A&=R^{-1/2}HD^{1/2}\\\ b&=R^{-1/2}(y-Hx^{b})\\\
z&=D^{-1/2}(x-x^{b})\\\
\mu^{2}&=\frac{\sigma_{R}^{2}}{\sigma_{D}^{2}}\end{split}$

the cost function $J(x)$ became:

$\begin{split}J(z)&=\mathinner{\\!\left\lVert
b-Az\right\rVert}_{2}^{2}+\mu^{2}\mathinner{\\!\left\lVert
z\right\rVert}_{2}^{2}\end{split}$ (22)

where $\mu$ can be interpreted as regularization parameter that determined the
relative weight between the background state and observations. It is evident
that the classical variational problem can be formulated into the Tikhonov
regularization problem as shown above.

##  3 Results and discussions

To run experiments, we will consider the a spring-mass system that governed by
the following equation:

$y^{\prime\prime}(t)+0.45y^{\prime}(t)+y(t)=3cos(2t)$ (23)

with initial condition $y(0)=2$, $y^{\prime}(0)=0$.
By letting:

$\begin{split}x_{1}(t)&=y(t)\\\ x_{2}(t)&=y^{\prime}(t)\end{split}$ (24)

we can convert this equation into a system of two first order differential
equations. We have:

$\begin{split}x_{1}^{\prime}(t)&=x_{2}(t)\\\
x_{2}^{\prime}(t)&=-x_{1}(t)-0.45x_{2}(t)+3cost(2t)\end{split}$ (25)

This system can be rewritten in:

$\begin{split}x^{\prime}(t)=Ax(t)+Bu(t)\end{split}$ (26)

where

$A=\begin{bmatrix}0&1\\\ -1&-0.45\end{bmatrix}$ (27) $B=\begin{bmatrix}1&0\\\
0&1\end{bmatrix}$ (28)

and

$u(t)=\begin{bmatrix}0\\\ 3cost(2t)\end{bmatrix}$ (29)

Using the finite difference approximation method, we can discrete the
governing equation and solve for the true state of the system. We assume that
the system has process noise with covariance matrix:

$Q=\begin{bmatrix}0.0005&0\\\ 0&0.0005\end{bmatrix}$ (30)

, and

$H=\begin{bmatrix}1&0\end{bmatrix}.$ (31)

The initial guess is: $x_{0}=\begin{bmatrix}0\\\ 0\end{bmatrix}$, and
$P_{0}=\begin{bmatrix}0.05&0\\\ 0&0.05\end{bmatrix}$.

Figure 1: The system states x1 and x2 Figure 2: Kalman filter prediction
Figure 3: Ensemble Kalman filter prediction Figure 4: 3D-Variational
prediction

Figure 1 depicts the true states of the system. Figure 2 shows the estimate of
$x_{1}$ using Kalman Filter (KF) algorithm. Figure 3 illustrates the estimate
of $x_{1}$ using Ensemble Kalman Filter (EnKF) algorithm with ensemble size of
100. It is observed that the Ensemble Kalman Filter performed better than
Kalman Filter. In both cases, although the initial guesses are far away from
observations, after several iterations, the algorithms can be able to quickly
get close to the true state of the system.
Figure 4 demonstrates the results from 3D-Variational method for predicting
$x_{1}$. It can be seen that there are some discrepancies between actual
states and predictions from 3D-Var method. This method also takes more
iterations to get close to the observations.
On the advantage of the Kalman Filter, when new data or observations are
available, it is easy to incorporate the previous optimal estimates without
recomputing anything. However, the Kalman Filter is not applicable for non-
linear system. In contrast, the Ensemble Kalman Filter has been proven in many
high-dimensional, non-linear data assimilation problems. By simulating a
collection of state vectors this method avoid computing large covariance
matrices, and then update the ensemble through time will improve the
computational overhead. One of the significant drawbacks of the Ensemble
Kalman Filter is that the full error covariance matrix of a high-dimensional
system can not be represented with only few ensemble m especially in a rank
deficiency problem where m is significant less than n.
The 3D-Variational method allows us to solve an optimization problem in one
run. A major drawback of this method is the convergence to the global minimum
of the cost function to be minimized. Another significant difficulty of 3D-Var
method is the need to design an appropriate background covariance matrix D
that properly defines the covariances of the model parameters.

##  4 Conclusion

In this study, we have implemented different data assimilation algorithms
including the Kalman Filter (KF), Ensemble Kalman Filter (EnKF), and
3D-Variational (3D-Var) to solve the spring-mass system governed by a second
order differential equations. Results from experiments show that the Ensemble
Kalman Filter algorithm has better performance compared to the Kalman Filter
algorithm. Although there are some discrepancies in 3D-Var, the method can
reasonably predict the states of the system. Results from 3D-Var method are
highly depend on the accuracy of the background error covariance matrix D. If
D is not define properly, the method may fail to find an optimal solution of
the cost function. In further study, we can compare these algorithms with
other variants of variational methods as well as filtering methods.

## References

  * [1] Welch, Greg, and Gary Bishop.  "An introduction to the Kalman filter." (1995): 127-132.
  * [2] Evensen, Geir. "The ensemble Kalman filter: Theoretical formulation and practical implementation." Ocean dynamics 53, no. 4 (2003): 343-367.
  * [3] Barker, Dale M., Wei Huang, Yong-Run Guo, A. J. Bourgeois, and Q. N. Xiao. "A three-dimensional variational data assimilation system for MM5: Implementation and initial results." Monthly Weather Review 132, no. 4 (2004): 897-914.
  * [4] Evensen, Geir. Data assimilation: the ensemble Kalman filter. Springer Science & Business Media, 2009.
  * [5] Asch, Mark, Marc Bocquet, and Maëlle Nodet. Data assimilation: methods, algorithms, and applications. Society for Industrial and Applied Mathematics, 2016.


# SFVInt: Simple, Fast and Generic Variable-Length Integer Decoding using Bit
Manipulation Instructions

Gang Liao Ye Liu Yonghua Ding Le Cai Jianjun Chen  ByteDance Infrastructure
System Lab<EMAIL_ADDRESS>{ye.liu, yonghua.ding, le.cai,
<EMAIL_ADDRESS>
(2018)

###### Abstract.

The ubiquity of variable-length integers in data storage and communication
necessitates efficient decoding techniques. In this paper, we present SFVInt,
a simple and fast approach to decode the prevalent Little Endian Base-128
(LEB128) varints. Our approach, distilled into a mere 500 lines of code,
effectively utilizes the Bit Manipulation Instruction Set 2 (BMI2) in modern
Intel and AMD processors, achieving significant performance improvement while
maintaining simplicity and avoiding overengineering. SFVInt, with its generic
design, effectively processes both 32-bit and 64-bit unsigned integers using a
unified code template, marking a significant leap forward in varint decoding
efficiency. We thoroughly evaluate SFVInt’s performance across various
datasets and scenarios, demonstrating that it achieves up to a 2x increase in
decoding speed when compared to varint decoding methods used in established
frameworks like Facebook Folly and Google Protobuf.

† Work performed at ByteDance Infrastructure System Lab

††copyright: acmcopyright††journalyear: 2018††doi:
XXXXXXX.XXXXXXX††conference: Make sure to enter the correct conference title
from your rights confirmation emai; June 03–05, 2018; Woodstock, NY††price:
15.00††isbn: 978-1-4503-XXXX-X/18/06

##  1\. Introduction

In data processing and storage, handling variable-length integers efficiently
is of paramount importance. These integers feature prominently across numerous
systems such as search engines like Apache Lucene (Foundation, 2011),
databases including IBM DB2 (Bhattacharjee et al., 2009) and Apache Kudu
(Lipcon et al., 2015), and popular file formats such as ORC (Foundation,
2013a) and Parquet (Foundation, 2013b), acting as a key component for compact
data representation. Their significance extends to serialization frameworks,
notably Google Protobuf (Google, 2001), and the WebAssembly binary encoding
(Consortium, 2022); they are also an essential part of the Go programming
language’s default API (Griesemer, 2011), highlighting its widespread adoption
across various critical technological domains (Liao et al., 2023;
Bhattacherjee et al., 2021; Liao and Abadi, 2023; Liao, 2022). The Little
Endian Base-128 (LEB128) encoding (Wikipedia, 2023), used for these variable-
length integers (varints), is widely adopted for its optimal balance of space
efficiency and implementation simplicity. However, a critical challenge arises
in varint decoding, particularly due to the unpredictable lengths of encoded
integers, which often lead to frequent branch mispredictions and thus become a
bottleneck in performance.

Recognizing this challenge, we turn our attention to the Bit Manipulation
Instruction Set 2 (BMI2) (bmi, 2013), an advancement in processor instruction
sets available in modern Intel (Intel Corporation, 2016) and AMD CPUs (AMD,
2023). BMI2 offers a suite of instructions that are particularly adept at
handling bitwise operations, a core component of varint decoding processes.
Leveraging these capabilities, we propose SFVInt, a generic, simple, and fast
approach to accelerate the decoding of LEB128 varints. Condensed into around
500 lines of code, our approach significantly outperforms traditional methods,
offering up to a 2x increase in decoding speed over techniques used in well-
established systems and libraries. Our approach is both simple and highly
effective, capable of handling 32-bit and 64-bit unsigned integers with a
unified C++ template. We provide a detailed analysis, highlighting how the use
of BMI2 instructions contributes to this substantial performance enhancement.

To the best of our knowledge, this is the first time the BMI2 technique has
been applied in the context of varint decoding, marking a significant
advancement in this field. This research presents a critical enhancement in
the efficiency of varint decoding, opening pathways for more refined data
processing optimizations. The ramifications of this discovery are extensive,
potentially influencing the design and performance optimization of databases,
network systems, and big data applications, crucial in a data-centric era.

##  2\. Background

###  2.1. Variable-Length Integers

Varints, or variable-length integers, are a method of encoding integers using
one or more bytes. They are widely adopted in computing for their ability to
optimize space usage. The principal idea behind varints is to use a dynamic
number of bytes for representation, depending on the integer’s magnitude.
Smaller numbers consume fewer bytes, making this method highly space-efficient
for data that predominantly consists of smaller integers.

Varint encoding, essential in data format optimization, encompasses several
schemes, including Group Varint (Dean, 2009), Prefix Varint (Dean, 2009), and
Stream Varint (Lemire et al., 2018), each with its unique characteristics.
Among these, LEB128 stands out for its simplicity and widespread adoption.
LEB128 operates on a simple yet effective principle: it encodes an integer
into a sequence of bytes, where each byte holds 7 bits of the integer plus a
continuation bit. This continuation bit indicates the presence of subsequent
bytes. In this scheme, the integer’s least significant 7-bit group is
positioned in the first byte, adhering to the little-endian convention.
Despite the existence of more complex varint formats, LEB128’s simplicity
makes it a prevalent choice in various applications, underlining its
effectiveness in efficiently handling integer encoding.

The decoding process of LEB128 involves reading each byte, extracting the 7
bits of the integer, and assembling them to reconstruct the original integer.
This process continues byte by byte until a byte with the continuation bit set
to 0 is encountered, indicating the end of the encoded integer.

LEB128 decoding, a bottleneck in high-volume data processing, faces challenges
due to encoded integers’ unpredictable lengths. This variability frequently
leads to branch mispredictions, exacerbating CPU overhead. The need for per-
byte decision-making during decoding amplifies this issue, especially in
large-scale or real-time processing environments. Addressing these
inefficiencies is crucial for optimizing LEB128’s decoding performance.

###  2.2. Bit Manipulation Instructions (BMI)

Bit Manipulation Instructions (BMI2) (bmi, 2013) are advanced CPU instruction
sets designed for efficient bit-level operations, introduced with Intel’s
Haswell architecture and adopted by AMD. These instructions, including PDEP
and PEXT, enable faster and more effective manipulation of bits in registers
or memory. BMI2’s inclusion in modern processors optimizes critical operations
in high-performance computing, data processing, and cryptographic
applications, reflecting a broader trend towards enhancing computational
efficiency through specialized hardware capabilities.

PDEP (Parallel Bit Deposit): This instruction streamlines the process of
depositing bits from a source operand into a destination operand based on a
mask. It allows for selectively placing bits in specific positions, enabling
precise control over bit placement. PDEP is particularly useful in scenarios
where bits from a variable need to be reorganized or masked in a specific
pattern, as illustrated in Figure 1.

Figure 1. PDEP Example.

PEXT (Parallel Bit Extract): PEXT complements PDEP by enabling the extraction
of bits from a source operand based on a mask. This instruction allows for the
selective extraction of bits, effectively ’filtering’ the source operand
through the mask to produce the desired bit pattern in the destination
operand. PEXT is invaluable for operations that require isolating specific
bits from a larger set, such as decoding varints where bits need to be
extracted from a sequence of bytes. An example of PEXT in action is depicted
in Figure 2.

Figure 2. PEXT Example.

Our motivation for utilizing BMI2 in LEB128 varint decoding lies in its
ability to simplify the bit extraction and assembly process, crucial for
efficient decoding. Traditional LEB128 decoding, characterized by sequential
bitwise operations, benefits from the parallelism and efficiency offered by
BMI2 instructions.

##  3\. SIMD Implementation

This section outlines traditional varint encoding/decoding methods, then
focuses on optimizing operations like varint sequence skipping and storage
size determination. We explore BMI2-enhanced implementation, detailing how it
improves decoding efficiency and performance, demonstrating a progression from
basic to advanced techniques in varint processing.

###  3.1. Basic Varint Operations

Varint Encoding: LEB128 encodes integers into a compact, variable-length
format by dividing each integer into 7-bit blocks, each represented by a
single byte. All bytes, barring the final one, mark their most significant bit
as 1 to denote the sequence’s continuation. The final byte’s most significant
bit is set to 0, indicating the end of the encoding. This method is efficient
for storing integers, particularly smaller values. The pseudocode below
details this encoding process:

    Algorithm: LEB128 Integer Encoding
    Input: Integer val, BufferPointer buf, Output: None
    1. while val >= 0x80 do
    2.   *buf = 0x80 | (val & 0x7f);
    3.   val >>= 7; buf++;
    5. end while
    6. *buf = val; buf++;

If val is 128 (0x80) or greater, the function stores the lowest 7 bits of val
into the buf, setting the 8th bit of this byte (denoted by 0x80) to signal
that more bytes follow. The buffer pointer buf is then advanced, and val is
right-shifted by 7 bits to prepare the next 7 bits for encoding. This process
repeats until val is less than 128. At this point, val is small enough to fit
into the final byte, indicating the end of the variable-length integer.

Varint Decoding: LEB128 extracts the original integer from its encoded
representation. It involves iteratively processing each byte, reconstructing
the integer by aggregating 7-bit segments.

    Algorithm: LEB128 Integer Decoding (Basic Version)
    Input: BufferPointer buf
    Output: T res (either uint32_t or uint64_t)
    1. res = 0; max_shift = std::is_same_v<T, uint32_t> ? 28 : 63;
    2. for shift = 0 to max_shift step 7 do
    3.   res |= ((*buf & 0x7f) << shift);
    4.   if (likely(!(*buf++ & 0x80))) break;
    5. end for

The decoding algorithm reads each byte, combining its lower 7 bits into the
result. The most significant bit indicates if subsequent bytes are part of the
integer. The decoding concludes when a byte with the most significant bit
unset is encountered, indicating the final byte of the varint. In line 1, the
maximum shift values for LEB128 decoding – 28 for 32-bit integers and 63 for
64-bit integers – are determined by the integer size and the 7-bit grouping
used in this encoding. For example, for a 32-bit unsigned integer, which has a
maximum value of 4,294,967,295 ($2^{32}$ \- 1), it would require up to 5 bytes
for representation, as 32 divided by 7 is approximately 4.57.

Varint Skipping: In columnar storage systems, efficient skip operations are
vital for effective filtering and scanning of varint-encoded data streams. We
employ a technique that processes 64-bit data blocks, enabling simultaneous
handling of multiple varints.

    Algorithm: LEB128 Integer Skipping
    Input: uint64_t n (# integers to skip), BufferPointer buf
    Output: Updated BufferPointer
    1. res = 0; w = (const uint64_t*)(buf);
    2. while n >= 8 do
    3.   n -= __builtin_popcountll(~(*w++) & 0x8080808080808080);
    4. end while
    5. buf = (const char*)(w);
    6. while n-- do
    7.   while *buf++ & 0x80 do {;} end while // fallback
    8. end while

In line 3, applying the bitwise NOT operation ($\sim$) inverts all bits of a
64-bit integer, turning continuation bits from ’1’s to ’0’s and vice versa.
With the mask 0x8080808080808080, highlighting MSBs in each 8-bit segment, the
AND operation (&) then isolates these MSBs. Utilizing SIMD instructions,
specifically POPCNT, this step efficiently identifies varint terminations
within the 64-bit block.

Varint Sizing: In varint-based systems, accurately estimating storage needs is
crucial. We utilize precomputed lookup tables (LUTs) to quickly calculate the
exact byte size needed for varint encoding, based on the position of the most
significant set bit in integers. This approach optimizes computational
efficiency with low overhead in caching, facilitating rapid estimations for
dynamic storage allocation.

    Algorithm: LEB128 Integer Sizing
    Input: vec<T> values (T is either uint32_t or uint64_t)
    Output: uint64_t size
    1. kLUT32 = {5, 5, 5, 5,
          4, 4, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3,
          2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1};
    2. kLUT64 = {10,
          9, 9, 9, 9, 9, 9, 9, 8, 8, 8, 8, 8, 8, 8,
          7, 7, 7, 7, 7, 7, 7, 6, 6, 6, 6, 6, 6, 6,
          5, 5, 5, 5, 5, 5, 5, 4, 4, 4, 4, 4, 4, 4,
          3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2,
          1, 1, 1, 1, 1, 1, 1};
    3. size = 0;
    4. for each val in values do
    5.   if constexpr (std::is_same_v<val, uint32_t>) then
    6.       size += kLUT32[__builtin_clz  (val | 1)]
    7.   else
    8.       size += kLUT64[__builtin_clzll(val | 1)]
    9. end for

This algorithm computes the total storage requirement for varint-encoded
integers by iteratively applying the LZCNT/BSR instruction to count leading
zeros, adjusted for zero values with (val|1). The outcome indexes into
predefined lookup tables (LUTs) for 32-bit and 64-bit integers (kLUT32,
kLUT64) to retrieve byte sizes.

###  3.2. BMI2-Enhanced Bulk Varint Decoding

This section enhances varint decoding with BMI2 instructions, advancing beyond
the basic methods in Section 3.1.

Mask Configuration: The first step involves employing a 64-bit mask with the
_pext_u64 BMI2 instruction for efficient extraction of most significant bits
(MSBs) from byte sequences. Configurations like 0x0000808080808080 for six
bytes, 0x0080808080808080 for seven bytes, and 0x8080808080808080 for eight,
cater to different ranges of extracted values (0-63, 0-127 and 0-255,
respectively), influencing the number of processed bytes. Our implementation
adopts the six-byte mask 0x0000808080808080, based on findings that longer
masks risk L1 instruction cache spillage, compromising performance. This
choice ensures an optimized balance between decoding efficiency and CPU
resource utilization, in line with the operational constraints of modern CPU
architectures.

Utilizing _pext_u64(word, mask), we decipher the count and positions of
integers within a 64-bit block. For example, using an 8-byte mask
(0x8080808080808080) that yields zero indicates eight complete integers in the
segment. Initially setting aside carryover or boundary effects for simplicity,
we employ _pext_u64 instructions strategically to isolate individual integers:

             int1 = _pext_u64(word, 0x000000000000007f);
             int2 = _pext_u64(word, 0x0000000000007f00);
             int3 = _pext_u64(word, 0x00000000007f0000);
             int4 = _pext_u64(word, 0x000000007f000000);
             int5 = _pext_u64(word, 0x0000007f00000000);
             int6 = _pext_u64(word, 0x00007f0000000000);
             int7 = _pext_u64(word, 0x007f000000000000);
             int8 = _pext_u64(word, 0x7f00000000000000);

For a masked result of 63, indicating two integers spanning six and two bytes
respectively in the original varint sequence. To extract them efficiently, we
apply the _pext_u64 instruction with specific masks tailored to these lengths:

             int1 = _pext_u64(word, 0x00007f7f7f7f7f7f);
             int2 = _pext_u64(word, 0x7f7f000000000000);

This method showcases adaptability to varint structures, identified by mask
results, reducing the need for conditional byte-by-byte checks that can cause
branch mispredictions and pipeline stalls.

Cross-Boundary Cases: In addressing cross-boundary cases where integers span
two or even three consecutive blocks, we adopt a refined approach crucial for
accurate and complete decoding. Our approach employs two key variables:
shift_bits and partial_value. shift_bits tracks the necessary bit displacement
for integers crossing block boundaries, while partial_value retains the
already decoded segment from a previous block. This system ensures proper
alignment and integration of integer segments across blocks, merging
partial_value with new segments to reconstruct complete integers.

Continuing from our previous example where a masked result of 63 suggests two
integers with the first spanning six bytes and the second two, we encounter
the complexity of int1 extending across 64-bit block boundaries. Here,
partial_value retains int1’s initial segment from the preceding block, while
shift_bits tracks the decoded bit count. The subsequent processing extracts
int1’s remainder, aligns it using shift_bits, and merges it with
partial_value, forming the complete integer, as demonstrated below:

            int1 = (_pext_u64(word, 0x00007f7f7f7f7f7f)
                    << shift_bits) | partial_value;
            int2 =  _pext_u64(word, 0x7f7f000000000000);

This method effectively handles integers crossing 64-bit block boundaries,
ensuring their precise and complete reconstruction.

Algorithm: The following presents an advanced varint decoding algorithm
optimized for bulk processing with a 6-byte MASK, integrating cross-boundary
case handling and fully leveraging BMI2 instructions for improved efficiency
and speed. Line 7 employs _pext_u64 to extract MSBs from a 64-bit word using a
defined mask, revealing the varint structure, including the count and
distribution of integers. The algorithm then uses this output (mval) in a
switch-case to customize the decoding for each integer’s unique encoding
pattern.

    Algorithm: Optimized Bulk Varint Decoding with BMI2
    Input: uint64_t n, BufferPointer buf, T* res (uint32/uint64)
    Output: Updated BufferPointer in the data stream
    1:  mask_length = 6; mask = 0x0000808080808080;
    2:  shift_bits = 0; pt_val = partial_value = 0;
    3:  buf -= mask_length;
    4:  while n >= 8 do
    5:    buf += mask_length;
    6:    word = *(const uint64_t*)(buf);
    7:    mval = _pext_u64(word, mask);
    8:    switch (mval)
    9:      case 0:
    10:         cu_val = _pext_u64(word, 0x000000000000007f);
    11:         *res++ = (cu_val << shift_bits) | pt_val;
    12:         *res++ = _pext_u64(word, 0x0000000000007f00);
    13:         *res++ = _pext_u64(word, 0x00000000007f0000);
    14:         *res++ = _pext_u64(word, 0x000000007f000000);
    15:         *res++ = _pext_u64(word, 0x0000007f00000000);
    16:         *res++ = _pext_u64(word, 0x00007f0000000000);
    17:         shift_bits = 0; pt_val = 0; n -= 6;
    18:      case 15:
    19:         cu_val = _pext_u64(word, 0x0000007f7f7f7f7f);
    20:         *res++ = (cu_val << shift_bits) | pt_val;
    21:         *res++ = _pext_u64(word, 0x00007f0000000000);
    22:         shift_bits = 0; pt_val = 0; n -= 2;
    23:      case 32:
    24:         cu_val = _pext_u64(word, 0x000000000000007f);
    25:         *res++ = (cu_val << shift_bits) | pt_val;
    26:         *res++ = _pext_u64(word, 0x0000000000007f00);
    27:         *res++ = _pext_u64(word, 0x00000000007f0000);
    28:         *res++ = _pext_u64(word, 0x000000007f000000);
    29:         *res++ = _pext_u64(word, 0x0000007f00000000);
    30:         pt_val = _pext_u64(word, 0x00007f0000000000);
    31:         shift_bits = 7; n -= 5;
    32:      case 45:
    33:         cu_val = _pext_u64(word, 0x0000000000007f7f);
    34:         *res++ = (cu_val << shift_bits) | pt_val;
    35:         *res++ = _pext_u64(word, 0x0000007f7f7f0000);
    36:         pt_val = _pext_u64(word, 0x00007f0000000000);
    37:         shift_bits = 7; n -= 2;
    38:      case 62:
    39:         cu_val = _pext_u64(word, 0x000000000000007f);
    40:         *res++ = (cu_val << shift_bits) | pt_val;
    41:         pt_val = _pext_u64(word, 0x00007f7f7f7f7f00);
    42:         shift_bits = 35; n -= 1;
    43:      case 63:
    44:         pt_val |= _pext_u64(word, 0x00007f7f7f7f7f7f)
    45:                << shift_bits;
    46:         shift_bits += 42;
    47:      case ... // Other cases omitted for brevity
    48: end while // fallback later

Figure 3. BMI2-enhanced Bulk Decoding.

In an example from lines 38-42, the algorithm manages a case where a
significant part of an integer starts in a preceding 64-bit block and
continues into the current one. The binary pattern 62 (represented as 111110)
under a 6-byte MASK indicates an integer’s end within the current block. Using
_pext_u64, the algorithm merges this block’s decoded segment with the
preceding block’s pt_val, completing the integer. Additionally, the presence
of five leading 1s (11111) signals an incomplete integer extending into the
next block, which is extracted and stored in pt_val for subsequent decoding,
with shift_bits adjusted to 35 for alignment, i.e., the five 7-byte sequences.
This decrement in n by 1 ensures the algorithm’s seamless management of cross-
boundary integers, maintaining decoding continuity and precision.

Case 63, depicted in lines 44-45, addresses a binary representation of 111111,
signaling an integer spanning three data blocks, with the current block as its
mid-segment. This specific case requires extracting six bytes (42 bits) via
_pext_u64, appending this data to the partial value to maintain integer
continuity across blocks.

##  4\. Performance Evaluation

Figure 4. Workload 1 (W1): uniform distribution - generated by
folly::Random::secureRand32() or secureRand64().

Figure 5. W2: varint byte distribution — 1 byte: 90.08%, 2 bytes: 4.63%, 3
bytes: 3.22%, 4 bytes: 1.20%, 5 bytes: 0.88%.

Figure 6. W3: varint byte distribution — 1 byte: 81.22%, 2 bytes: 7.31%, 3
bytes: 6.16%, 4 bytes: 4.20%, 5 bytes: 1.10%.

Figure 7. W4: varint byte distribution — 1 byte: 72.13%, 2 bytes: 12.31%, 3
bytes: 8.53%, 4 bytes: 5.31%, 5 bytes: 1.72%.

In our evaluation of the BMI2-enhanced varint decoding technique, we leveraged
various AWS EC2 instances with diverse Intel and AMD CPU architectures (refer
to Table 1). We created datasets encoded as varints, following uniform and
skewed distributions. The uniform distribution served as a stable benchmark,
while the skewed distribution, biased towards smaller integers (1-byte or
2-bytes LEBs), mirrored real-world varint patterns.

We employed four distinct workloads for evaluation: W1 featured uniformly
distributed 32-bit integers as a benchmark; W2 used real-world LEB lengths
from the WebAssembly build suite (web, [n. d.]), providing practical insights;
and Workloads W3 and W4 derived from our proprietary systems’ data
distributions, offering a realistic view of their application. Notably, W2 to
W4 exhibited a skew towards 1-byte integers, aligning with typical data trends
observed in practical scenarios. The specific distributions for these
workloads are detailed in the figure captions of Figures 4 to 7.

Our evaluation compares our BMI2-based decoding technique against established
varint encoding solutions from Apache Kudu (a columnar storage system)(Lipcon
et al., 2015), Google Protobuf (a serde library)(Google, 2001), Apache ORC (a
columnar format)(Foundation, 2013a), and Facebook Folly (a C++ std
lib)(Facebook, 2012), which all support 32-bit and 64-bit integer encoding.
This benchmark situates our approach within the current industry
practices111Our assessment directly interacted with the varint decoding
interfaces of these libraries, ensuring no additional overhead is incurred.,
focusing on in-memory decoding performance given the widespread use of the
LEB128 algorithm. Performance gains were assessed in a single-threaded
environment to attribute improvements directly to our method.

EC2 | CPU Architecture
---|---
m6i | Intel Xeon Ice Lake 8375C
m5 | Intel Xeon Skylake 8175M / Cascade Lake 8259CL
m4 | Intel Xeon Broadwell E5-2686 / Haswell E5-2676
c6i | Intel Xeon Ice Lake 8375C
c5 | Intel Xeon Cascade Lake
c4 | Intel Xeon Haswell E5-2666 v3
m6a | 3rd Gen AMD EPYC (Milan)
c6a | 3rd Gen AMD EPYC (Milan)
m5a | 2nd Gen AMD EPYC 7571
c5a | 2nd Gen AMD EPYC 7R32
Table 1. AWS EC2 Types and CPU Architectures

###  4.1. Intel CPU Microarchitectures

In Figures 4-7, the left panel showcases the varint decoding performance
across six EC2 instance types featuring Intel CPUs, detailed in Table 1. Each
figure within the panel includes benchmarks for each EC2 instance type, with
the y-axis representing the time required per iteration (in milliseconds),
where each iteration involves decoding one million integers generated
according to a specific distribution. Here, SFVInt consistently excels,
surpassing all other systems in decoding speed across every instance type and
under various distributions. For example, in W1 (Figure 4), SFVInt
demonstrates a significant performance edge, nearly doubling the speed of
Protobuf. Similarly, in W4 (Figure 7), SFVInt outperforms all other systems by
nearly a factor of two.

The varint byte distribution evolves from W2 to W4, indicating a transition to
a more varied spread of integer sizes. In W2 (Figure 5), the distribution
heavily favors 1-byte integers, constituting 90.08% of the data, with a minor
representation of 2-byte at 4.63%, and even less for 3-byte (3.22%), 4-byte
(1.20%), and 5-byte integers (0.88%). Progressing to W4, there’s a noticeable
decrease in 1-byte integers to 72.13%, with corresponding increases across the
board—2-byte integers rise to 12.31%, 3-byte to 8.53%, 4-byte to 5.31%, and
5-byte to 1.72%. This shift implies an increase in the number of bytes subject
to the decoding process, logically resulting in prolonged iteration times as
the byte distribution broadens. Specifically, on the m6i.4xlarge instance, the
iteration time for SFVInt extends from 1.6ms in W2, to 2.2ms in W3, ultimately
reaching 3.1ms in W4. Such an increment reflects the additional computational
overhead necessitated by the decoding of a higher proportion of multi-byte
integers.

On the other hand, the performance gap between SFVInt and other systems widens
most notably in W4. Taking the c6i.4xlarge instance as an example, SFVInt
achieves a 2x speed increase compared to Protobuf in W4, while in W3 and W2,
the gains are 45% and 19%, respectively. This indicates that SFVInt’s
BMI2-fortified decoding mechanism gains a significant edge when dealing with
varints exceeding the single-byte range. This is because with shorter byte-
length encodings, our approach requires a sequence of serial instructions
similar to traditional byte-by-byte processing. This is exemplified by case 0
(lines 10 - 17) in Figure 3, which involves a greater number of sequential
executions than case 63 (lines 44 - 46), hence, SFVInt’s relatively subdued
enhancement for 1-byte integers reflects this similarity to conventional
decoding approaches.

###  4.2. AMD CPU Microarchitectures

In our analysis detailed in Figures 4-7, we observe that SFVInt demonstrates
up to 40% faster varint decoding on 3rd generation AMD EPYC processors (m6a
and c6a), in contrast to the slower performance on 2nd generation CPUs (m5a
and c5a). This performance disparity is consistent with the known emulation
overhead of PEXT/PDEP instructions on AMD platforms, which can significantly
lag behind Intel’s native execution (Kullberg, 2017; red, 2017). Studies
indicate that on AMD hardware, these instructions’ latency can vary widely,
from 18 to approximately 300 cycles (InstLatX64, 2019). The improved results
with newer AMD generations suggest advancements in BMI2 instruction handling.
For enhanced consistency across CPU architectures, we can incorporate a
dynamic selection mechanism between BMI2-accelerated and standard decoding or
employs ZP7 (Wegner, 2019), a branchless alternative for AMD processors.

##  5\. Related Works

Varint Formats. In byte-oriented integer compression, various techniques aim
to improve decoding efficiency. Traditional LEB128, using a 7-bit structure,
lacks SIMD optimization. In contrast, Google’s VARINT-GB (Dean, 2009) offers
an alternative by encoding a constant block of integers with a unified control
byte, streamlining decoding for data with variable lengths. The patented
VARINT-G8IU (Stepanov et al., 2011) further advances this innovation,
employing SIMD to process variable counts of integers within a fixed byte
structure, controlled by an initial byte. Stream VByte (Lemire et al., 2018)
segregates control and data bytes into separate streams to fully leverage
SIMD’s parallel processing strengths. However, these techniques, despite their
SIMD-enabled efficiency gains, diverge from the prevalent LEB128 standard,
limiting their seamless integration into existing infrastructures. Conversely,
the Masked VByte method (Plaisance et al., 2017) retains the LEB128 form while
harnessing SIMD for improved performance, yet it is not without its
drawbacks—namely, a non-generic and maintenance-intensive codebase exceeding
2300 lines, tailored specifically for 32-bit unsigned integers. Our proposed
SFVInt is a simple, elegant, fast, and generic solution that effectively
addresses the limitations present in earlier methods.

BMI in DBMS. The pursuit of SIMD optimizations for data scanning operations
and operator vectorization has been prominent in data management research (Li
and Patel, 2013; Polychroniou and Ross, 2015; Sitaridi et al., 2016;
Polychroniou and Ross, 2014). Kudu was seminal in applying BMI2 for data scan
filtering (Lipcon et al., 2015; Lipcon, 2020), a similar idea extended by
Parquet-Select (Li et al., 2023), which utilizes BMI2 for optimized selection
pushdown in columnar storage.

##  6\. Conclusions

In this paper, we introduced SFVInt, a simple, fast and generic varint
decoding approach leveraging BMI2 instructions, significantly enhancing the
performance of varint decoding. Through rigorous evaluation across diverse CPU
architectures, SFVInt consistently demonstrated up to a 2x increase in
decoding speed over traditional methods on modern processors. This performance
boost is crucial for applications dealing with large-scale data processing.

## References

  * (1)
  * web ([n. d.]) [n. d.].  _WebAssembly Build Suite_.  Retrieved Jan 2, 2024 from https://github.com/WebAssembly/build-suite
  * bmi (2013) 2013\.  _x86 Bit manipulation instruction set_.  Retrieved Jan 6, 2023 from https://en.wikipedia.org/wiki/X86_Bit_manipulation_instruction_set
  * red (2017) 2017\.  _Ryzen and BMI2: Strange behavior and high latencies_.  Retrieved Jan 2, 2024 from https://www.reddit.com/r/Amd/comments/60i6er/ryzen_and_bmi2_strange_behavior_and_high_latencies/
  * AMD (2023) AMD. 2023.  _AMD64 Architecture Programmer’s Manual_.  Retrieved Jan 6, 2023 from https://www.amd.com/content/dam/amd/en/documents/processor-tech-docs/programmer-references/24594.pdf
  * Bhattacharjee et al. (2009) Bishwaranjan Bhattacharjee, Lipyeow Lim, Timothy Malkemus, George Mihaila, Kenneth Ross, Sherman Lau, Cathy McArthur, Zoltan Toth, and Reza Sherkat. 2009.  Efficient Index Compression in DB2 LUW.  _Proc. VLDB Endow._ 2, 2 (aug 2009), 1462–1473.  https://doi.org/10.14778/1687553.1687573
  * Bhattacherjee et al. (2021) Souvik Bhattacherjee, Gang Liao, Michael Hicks, and Daniel J. Abadi. 2021.  BullFrog: Online Schema Evolution via Lazy Evaluation. In _Proceedings of the 2021 International Conference on Management of Data_ (Virtual Event, China) _(SIGMOD ’21)_. Association for Computing Machinery, New York, NY, USA, 194–206.  https://doi.org/10.1145/3448016.3452842
  * Consortium (2022) World Wide Web Consortium. 2022.  _WebAssembly Binary Format_.  Retrieved Jan 6, 2023 from https://webassembly.github.io/spec/core/binary/values.html
  * Dean (2009) Jeffrey Dean. 2009.  Challenges in building large-scale information retrieval systems: invited talk. In _WSDM ’09: Proceedings of the Second ACM International Conference on Web Search and Data Mining_. New York, NY, USA, 1–1.  http://doi.acm.org/10.1145/1498759.1498761
  * Facebook (2012) Facebook. 2012.  _Folly: An open-source C++ library developed and used at Facebook_.  Retrieved Dec 30, 2023 from https://github.com/facebook/folly
  * Foundation (2011) Apache Software Foundation. 2011.  _Apache Lucene_.  Retrieved Jan 4, 2023 from https://lucene.apache.org/
  * Foundation (2013a) Apache Software Foundation. 2013a.  _Apache ORC: the smallest, fastest columnar storage for Hadoop workloads_.  Retrieved Dec 30, 2023 from https://orc.apache.org
  * Foundation (2013b) Apache Software Foundation. 2013b.  _Apache Parquet: an open source, column-oriented data file format_.  Retrieved Jan 4, 2023 from https://parquet.apache.org/
  * Google (2001) Google. 2001.  _Protocol Buffers: Google’s data interchange format_.  Retrieved Dec 30, 2023 from https://github.com/protocolbuffers/protobuf
  * Griesemer (2011) Robert Griesemer. 2011.  _Support for varint encoding in Go_.  Retrieved Jan 6, 2023 from https://github.com/golang/go/commit/f30719dc89c2a41502fa584b790943170ad2d1ce
  * InstLatX64 (2019) InstLatX64. 2019.  _Achilles heel of AMD Zens: data dependency of PDEP/PEXT instructions_.  Retrieved Jan 6, 2023 from https://mobile.twitter.com/InstLatX64/status/1209095219087585281
  * Intel Corporation (2016) Intel Corporation. 2016.  _Intel® 64 and IA-32 Architectures Software Developer’s Manual_.  Retrieved Jan 6, 2023 from https://www.intel.com/content/dam/www/public/us/en/documents/manuals/64-ia-32-architectures-software-developer-vol-2b-manual.pdf
  * Kullberg (2017) Kullberg. 2017.  _Ryzen Schach Performance - BMI2 Problem_.  Retrieved Jan 2, 2024 from https://www.hardwareluxx.de/community/threads/ryzen-schach-performance-bmi2-problem.1156117/
  * Lemire et al. (2018) Daniel Lemire, Nathan Kurz, and Christoph Rupp. 2018.  Stream VByte: Faster byte-oriented integer compression.  _Inform. Process. Lett._ 130 (Feb. 2018), 1–6.  https://doi.org/10.1016/j.ipl.2017.09.011
  * Li et al. (2023) Yinan Li, Jianan Lu, and Badrish Chandramouli. 2023.  Selection Pushdown in Column Stores Using Bit Manipulation Instructions.  _Proc. ACM Manag. Data_ 1, 2, Article 178 (jun 2023), 26 pages.  https://doi.org/10.1145/3589323
  * Li and Patel (2013) Yinan Li and Jignesh M. Patel. 2013.  BitWeaving: Fast Scans for Main Memory Data Processing. In _Proceedings of the 2013 ACM SIGMOD International Conference on Management of Data_ (New York, New York, USA) _(SIGMOD ’13)_. Association for Computing Machinery, New York, NY, USA, 289–300.  https://doi.org/10.1145/2463676.2465322
  * Liao (2022) Gang Liao. 2022.  The Evolution of Cloud Data Architectures: Storage, Compute, and Migration.  https://drum.lib.umd.edu/items/e591f36a-a240-42db-8252-196ed4facee9.
  * Liao and Abadi (2023) Gang Liao and Daniel J. Abadi. 2023.  FileScale: Fast and Elastic Metadata Management for Distributed File Systems. In _Proceedings of the 2023 ACM Symposium on Cloud Computing_ (, Santa Cruz, CA, USA,) _(SoCC ’23)_. Association for Computing Machinery, New York, NY, USA, 459–474.  https://doi.org/10.1145/3620678.3624784
  * Liao et al. (2023) Gang Liao, Amol Deshpande, and Daniel J. Abadi. 2023.  Flock: A Low-Cost Streaming Query Engine on FaaS Platforms.  arXiv:2312.16735 [cs.DB] https://github.com/flock-lab/flock
  * Lipcon (2020) Todd Lipcon. 2020.  _Core algorithms for columnar serialization_.  Retrieved Jan 6, 2023 from https://github.com/apache/kudu/commit/0ba6cb8d6b38a992786e5027528349a43802fd31
  * Lipcon et al. (2015) Todd Lipcon, David Alves, Dan Burkert, Jean-Daniel Cryans, Adar Dembo, Mike Percy, Silvius Rus, Dave Wang, Matteo Bertozzi, Colin Patrick McCabe, et al. 2015\.  Kudu: Storage for fast analytics on fast data.  _Cloudera, inc_ 28 (2015), 36–77.
  * Plaisance et al. (2017) Jeff Plaisance, Nathan Kurz, and Daniel Lemire. 2017.  Vectorized VByte Decoding.  arXiv:1503.07387 [cs.IR]
  * Polychroniou and Ross (2014) Orestis Polychroniou and Kenneth A. Ross. 2014.  Vectorized Bloom Filters for Advanced SIMD Processors. In _Proceedings of the Tenth International Workshop on Data Management on New Hardware_ (Snowbird, Utah) _(DaMoN ’14)_. Association for Computing Machinery, New York, NY, USA, Article 6, 6 pages.  https://doi.org/10.1145/2619228.2619234
  * Polychroniou and Ross (2015) Orestis Polychroniou and Kenneth A. Ross. 2015.  Efficient Lightweight Compression Alongside Fast Scans. In _Proceedings of the 11th International Workshop on Data Management on New Hardware_ (Melbourne, VIC, Australia) _(DaMoN’15)_. Association for Computing Machinery, New York, NY, USA, Article 9, 6 pages.  https://doi.org/10.1145/2771937.2771943
  * Sitaridi et al. (2016) Evangelia Sitaridi, Orestis Polychroniou, and Kenneth A. Ross. 2016.  SIMD-Accelerated Regular Expression Matching. In _Proceedings of the 12th International Workshop on Data Management on New Hardware_ (San Francisco, California) _(DaMoN ’16)_. Association for Computing Machinery, New York, NY, USA, Article 8, 7 pages.  https://doi.org/10.1145/2933349.2933357
  * Stepanov et al. (2011) Alexander A. Stepanov, Anil R. Gangolli, Daniel E. Rose, Ryan J. Ernst, and Paramjit S. Oberoi. 2011.  SIMD-Based Decoding of Posting Lists. In _Proceedings of the 20th ACM International Conference on Information and Knowledge Management_ (Glasgow, Scotland, UK) _(CIKM ’11)_. Association for Computing Machinery, New York, NY, USA, 317–326.  https://doi.org/10.1145/2063576.2063627
  * Wegner (2019) Zach Wegner. 2019.  _ZP7: Zach’s Peppy Parallel-Prefix-Popcountin’ PEXT/PDEP Polyfill_.  Retrieved Jan 6, 2023 from https://github.com/zwegner/zp7
  * Wikipedia (2023) Wikipedia. 2023.  _LEB128 (Little Endian Base 128)_.  Retrieved Jan 6, 2023 from https://en.wikipedia.org/wiki/LEB128