# A Machine Learning Benchmark for Facies Classification

Yazeed Alaudah11footnotemark: 133footnotemark: 3 Patrycja
Michałowicz22footnotemark: 2 Motaz Alfarraj33footnotemark: 3 and Ghassan
AlRegib33footnotemark: 3

###### Abstract

The recent interest in using deep learning for seismic interpretation tasks,
such as facies classification, has been facing a significant obstacle, namely
the absence of large publicly available annotated datasets for training and
testing models. As a result, researchers have often resorted to annotating
their own training and testing data. However, different researchers may
annotate different classes, or use different train and test splits. In
addition, it is common for papers that apply machine learning for facies
classification to not contain quantitative results, and rather rely solely on
visual inspection of the results. All of these practices have lead to
subjective results and have greatly hindered the ability to compare different
machine learning models against each other and understand the advantages and
disadvantages of each approach.

To address these issues, we open-source a fully-annotated 3D geological model
of the Netherlands F3 Block. This model is based on the study of the 3D
seismic data in addition to 26 well logs, and is grounded on the careful study
of the geology of the region. Furthermore, we propose two baseline models for
facies classification based on a deconvolution network architecture and make
their codes publicly available. Finally, we propose a scheme for evaluating
different models on this dataset, and we share the results of our baseline
models. In addition to making the dataset and the code publicly available,
this work helps advance research in this area by creating an objective
benchmark for comparing the results of different machine learning approaches
for facies classification.

Alaudah et al.

##  1 Introduction

In recent years, there has been great interest in using fully-supervised deep
learning models for seismic interpretation tasks such as facies classification
Shi et al., (2018); Dramsch and Lüthje, (2018); Waldeland and Solberg, (2017);
Araya-Polo et al., (2017); Huang et al., (2017); Zhao, (2018); Di et al.,
(2018); Rutherford Ildstad and Bormann, (2017). Typically, deep learning
models — such as convolutional neural networks (CNNs) — have millions of free
parameters, and therefore require a large amount of annotated training data.
Unfortunately, and unlike other areas of research such as computer vision,
there is a lack of large publicly-available annotated datasets for seismic
interpretation that can be used to train and benchmark machine learning
models. To address this problem, some researchers resort to annotating their
own training and testing datasets. For example, in the Netherlands F3 block,
Zhao, (2018) annotated 40 inlines, Di et al., (2018) annotated 12 inlines,
while Rutherford Ildstad and Bormann, (2017) only annotated a single inline
for their model. The limited number of annotated sections is understandable
given that the annotation process is time-consuming, requires subject matter
expertise, and can be quite subjective. Nevertheless, such limited annotations
undermine the mass potential machine learning could have when deployed in such
a field.

Alternatively, there has been some research in attempting to avoid annotating
large amounts of data by using weakly-supervised learning approaches. Alaudah
and AlRegib, (2016) trained a facies classification model using seismic images
with image-level labels only. Later, Alaudah et al., (2019) proposed a method
for generating large amounts of training data using similarity-based retrieval
and a weakly-supervised label mapping algorithm. As few as one or two exemplar
images per class were enough to automatically generate a large amount of
training data. These automatically-generated training data were then used to
train a weakly-supervised deconvolution network Alaudah et al., (2018) for
facies classification. Other researchers avoid supervision all-together by
using traditional unsupervised machine learning techniques such as principal
component analysis or self-organizing maps. There is a very rich literature on
traditional supervised and unsupervised methods for facies classification,
e.g., Coléou et al., (2003); de Matos et al., (2006); Dubois et al., (2007).
Zhao et al., (2015) provides a review of some of the most commonly used
techniques. More recently, unsupervised techniques based on deep learning
models such as deep convolutional autoencoders have been explored Qian et al.,
(2018); Shafiq et al., (2018); Veillard et al., (2018).

Whether researchers annotate their own training data or use other techniques,
there still remains a lack of large publicly-available annotated datasets for
seismic interpretation that can be used for training different models and
comparing the performance of different approaches. Furthermore, it is common
for papers that apply machine learning for facies classification, or other
seismic interpretation tasks, to not contain quantitative results, but rather
rely solely on subjective visual inspection of the results. All of this leads
to highly subjective results and greatly hinders the ability of researchers to
compare different approaches against each other and understand the advantages
and disadvantages of each approach.

To address these issues, and to help make machine learning research in seismic
interpretation more reproducible, we open-source a fully-annotated 3D
geological model of the Netherlands F3 Block dGB Earth Sciences, (1987). This
model is grounded in the geology of the region and based on the study of both
the 3D seismic data and 26 well logs located within the F3 block or its
vicinity. The data also includes fault planes that we have extracted from the
F3 block. While we do not use the fault data in this work, we make the data
publicly available for those who are interested in exploring fault detection
within our model. Furthermore, we also present two baseline models for facies
classification based on a deconvolution network architecture. The first
baseline is a patch-based model that is trained using a large number of small
patches extracted from all the inlines and crosslines of the training set. The
second baseline is a section-based model that was trained directly on entire
inlines and crosslines of the training set. In addition, we have open-sourced
all the codes that were used to train and test our baseline models using the
PyTorch deep learning library***Both the code and the dataset can be viewed
at: www.github.com/olivesgatech/facies_classification_benchmark. Finally, we
propose a common procedure for evaluating different models on this dataset,
and we share the results for our baseline models. The next section provides a
quick overview of the geology of the Netherlands F3 block and introduces our
geological model.

##  2 A 3D Geological Model of the Netherlands F3 Block

Figure 1: The location of the F3 block. Adapted from Duin et al., (2006).

The North Sea is rich in hydrocarbon deposits, which is why this area is very
well studied in the literature Doornenbal, (2014). The North Sea continental
shelf, located off the shores of the Netherlands, is divided into geographical
zones described by different letters of the alphabet; within these zones are
smaller areas marked with numbers. One of these areas is a rectangle of
dimensions 16 km x 24 km known as the F3 block, see Figure 1. In 1987, the F3
block 3D seismic survey was conducted to identify the geological structures of
this area and to search for hydrocarbon reservoirs. In addition, many
boreholes were drilled within the F3 block throughout the years. The F3 block
became one of the most widely known and studied seismic surveys after dGB
Earth Sciences made the data obtained from the survey publicly available.

The aim of this section is to briefly describe the geology of the survey area
and introduce the 3D geological model that we have developed and how it was
obtained.

###  2.1 The geology of the F3 block

Within the continental shelf of the North Sea, ten groups of
lithostratigraphic units have been identified in the literature Van Adrichem
Bogaert and Kouwe, (1993); Mijnlieff, (2002); Scheck-Wenderoth and Lamarche,
(2005); Duin et al., (2006). These groups and their main lithostratigraphic
features are listed below from newest to oldest:

  1. 1.

Upper North Sea group: claystones and sandstones from Miocene to Quaternary.

  2. 2.

Lower and Middle North Sea groups: sands, sandstones, and claystones from
Paleocene to Miocene.

  3. 3.

Chalk group: carbonates of Upper Cretaceous and Paleocene.

  4. 4.

Rijnland group: clay formations with sandstones of Upper Cretaceous.

  5. 5.

Schieland, Scruff and Niedersachsen groups: claystones of Upper Jurassic and
Lower Cretaceous.

  6. 6.

Altena group: claystones and carbonates of Lower and Middle Jurassic.

  7. 7.

Lower and Upper Germanic Trias groups: sandstones and claystones of Triassic.

  8. 8.

Zechstein group: evaporites and carbonates of Zechstein.

  9. 9.

Upper and Lower Rotliegend groups: siliceous rocks and basalts of the Lower
Zechstein.

  10. 10.

Limburg group: Upper carboniferous siliceous rock, which are the bedrock for
hydrocarbons.

Figure 2: A) A geological cross-section of the North Sea continental shelf
along axis A-A’; B) A map of the location of the cross-section. Adapted from
Duin et al., (2006).

The F3 block is located on the border of two tectonic structures: the Step
Graben and the Dutch Central Graben (see Figure 2). These tectonic structures
are characterized by different lithostratigraphic units of varying thickness.
This varying thickness is a result of tectonic activity Ziegler, (1988, 1990),
which was started in the Variscan orogeny Schroot and De Haan, (2003). The
area within the Step Graben is strongly disturbed by salt diapirs, which were
active several times, from the Zechstein to the Paleogene period Remmelts,
(1996). On the other hand, and as a result of subsiding Jurassic rocks, the
Altena, Scruff, Schieland and Niedersachsen groups are observed only within
the Dutch Central Graben Duin et al., (2006).

###  2.2 The modeling process

To prepare our 3D geological model of the F3 block, we relied on both well
logs and 3D seismic data. The next two subsections describe this process.

####  2.2.1 3D model building using well logs data

The well log data were obtained from a website managed by the Geological
Survey of the Netherlands (www.nlog.nl). The data (including information
related to coordinates, true vertical depth, measured depth along the
curvature, inclinations, and individual horizons) were collected for 26
boreholes located within the F3 block or its vicinity. The exact locations of
these wells are visualized in Figure 3.

Figure 3: Locations of the boreholes that were used to create the geological
model.

Originally, the 26 wells contained 40 different horizons, so it was necessary
to assign these different horizons to the various lithostratigraphic units
that were adopted in literature and were presented in the previous subsection.
The next step was correlating wells with each other. After that, it was
possible to create a preliminary 3D model based on the well log data by using
Petrel’s make/edit surface tool. This process facilitated the preliminary
visualization of the range of individual horizons, which was very helpful in
the further interpretation of the 3D seismic data.

####  2.2.2 3D model building using seismic data

Figure 4: A 3D view of our geological model of the F3 block.

The F3 block data was migrated in time, not depth, so it was necessary to do
time-depth conversion since the structural model must be prepared in the depth
domain. OpendTect 5.0 was used to perform the time-depth conversion using a
velocity model that was provided with the F3 block data.

The next step in creating the model was faults-surface interpretation. Using
Petrel’s polygon editing tool, we interpreted the main fault surfaces and the
fault networks were created by using the fault framework modeling tool.
Horizons were interpreted in a similar fashion, but by using the seeded 3D
autotracking tool, which interpolated data automatically and took into account
the faults networks modeled previously.

Based on the interpreted horizons and faults, preliminary modeling was
conducted. This was done using the horizon modeling tool with volume-based
modeling, which is an advanced method of isochronous geological space
modeling. The preliminary model included several imperfections in the
interpretations of horizons and faults, so it was necessary to re-model
several faults and conduct small corrections in the interpreted horizons.

After this, it was possible to create the final three-dimensional model which
highlights the regions between individual horizons. Here, Petrel’s structural
modeling module in the horizon modeling tool was used in addition to the
create zone model function. The final 3D geological model is shown in Figure
4.

###  2.3 The 3D geological model

Within our 3D geological model of the F3 Block, we identified seven groups of
lithostratigraphic units (see Figure 4). These are (from newest to oldest):
the Upper North Sea group, the Middle North Sea group, the Lower North Sea
group, the Chalk group, the Rijnland group, the Scruff group, and the
Zechstein group.

These groups can be divided into three structural levels: Cenozoic (Lower,
Middle, and Upper North Sea groups), Mesozoic (Scruff, Rijnland, and Chalk
groups), and Permian (Zechstein group).

As is evident in Figure 4, the F3 Block is characterized by highly variable
geological structures, both in the horizontal and vertical range, which is
manifested by the differential thicknesses of individual units and by the
expanded faults network related to salt tectonics. The area of the F3 Block
can be divided into two regions: Eastern and Western. The Eastern region is
disturbed by the occurrence of Zechstein diapirs and irregular faults network.
The Western region is characterized by regular fault networks and a more
uniform thickness of lithostratigraphic units.

The Upper North Sea group is the youngest and the flattest lithostratigraphic
unit within our model. The top of the Upper North Sea group is the bottom of
the North Sea at the same time, which is about -40 meters above sea level (m
a.s.l). Differences in the depth of the ocean floor are small, and they are
maximally 6 meters within the whole F3 Block. It can be noted that the depth
of this top decreases from SW to NE. The thickness of the Upper North Sea
group varies from about 1000 m (in places deformed by Permian diapirs) to
about 1320 m in the northern part of the research area (see Figure 4).

Below the Upper North Sea group lays the Middle North Sea group. The depth of
the top of this unit ranges from -1000 m a.s.l. within the diapir in NE part
of the F3 Block to about -1360 m a.s.l. in the northern part of this area,
between diapirs. The thickness of the Middle North Sea group is from 20 to 150
m. As in the case of the Upper North Sea group, there is a clear relationship
between the occurrence of Zechstein salts and the depth and thickness of this
unit. Differences in the thickness of this unit between both sites of faults
are also visible.

The next unit is the Lower North Sea group. This unit contains similar
lithostratigraphic units to the Middle North Sea group, but is visually
distinct in the seismic data. The top is at a depth from -1100 m a.s.l., while
the thickness is from about 180 to 750 m.

The top of the Chalk group is at depth from -1300 m a.s.l. (above the diapirs
in the NE part of the survey) to -2100 m a.s.l. (in the Eastern part of the
survey, which is undisturbed by diapirs). The minimum thickness of this unit
is 25 m, while above the salt diapirs in NE part of the F3 Block, this
substantially increases to 525 m.

The Rijnland group is submerged in the NNE direction, while it is the
shallowest in the SW part of the F3 Block and above some Zechstein diapirs at
the center of the survey (see Figure 6). The maximum thickness of the Rijnland
group is about 200 m (above some diapirs), while in the other parts of the F3
block it can be less than 20 m or does not occur at all.

The Scruff group, similar to the Rijnland group, is thinned out in NNE
direction, more or less in the middle of the F3 Block, where the top of this
layer has a depth of -2180 m a.s.l. This layer is shallowest (-1500 m a.s.l.)
in the SW part of the F3 block and above the Zechstein diapirs in the Southern
part of the survey. The thickness of the Scruff group within our model
boundaries ranges from 100 m to almost 700 m, but is much larger in reality
and can reach several kilometers Duin et al., (2006).

The Zechstein group occurs only in the eastern part of the survey, as
irregularly-shaped salt diapirs. The shallowest part of the Zechstein group is
at a depth of -1500 m a.s.l. while the maximum thickness of the Zechstein
group within the research area is about 700 m. However, as in the case of the
Scruff group, the depth is much bigger. According to the literature, it can
reach several kilometers Duin et al., (2006).

Figure 5: An overhead view of 3D fault planes from three different generations
of faults that we have identified in the F3 block.

In addition to the identified groups of lithostratigraphic units mentioned
above, we have also identified three generations of faults. The first
generation are reverse, oblique-slip, sinistral faults with an SSW-NNE
orientation. This direction is connected with the course of the tectonic axis
of the Dutch Central Graben, which (similar to the whole Graben) has an SSW-
NNE orientation. The second generation of faults are normal, oblique-slip,
dextral faults with a W-E orientation. Finally, the third generation are
faults that are genetically linked with faults from the first and second
generations, but were disturbed by the Permian halokinesis. Figure 5 shows an
overhead view of the three generations of faults that we have identified.
Also, Figure 6 shows two diagonal cross sections along the SW-NE and NW-SE
axis in our 3D model shown in Figure 4.

Figure 6: Two diagonal cross sections of our 3D geological model in Figure 4.

##  3 Deconvolution Network Baseline

In addition to the geological model we have introduced in the last section, we
propose two baseline models for facies classification based on a deconvolution
network architecture. In this section, we introduce deconvolution networks and
describe our two baseline models.

###  3.1 Deconvolution networks

Figure 7: The architecture of the deconvolution network used in this work.
White layers are convolution or deconvolution layers. Red layers are max-
pooling layers, while green layers are unpooling layers.

Convolutional neural networks (CNNs) have seen great success in a wide number
of visual applications, from image classification to semantic labeling. In the
previous few years, it was not well established that end-to-end CNN
architectures can perform very well in semantic labeling tasks (such as facies
classification). A major hurdle to the success of end-to-end CNN architectures
in these tasks was what seemed like a trade-off between classification and
localization accuracy. Deeper networks that have many convolution and pooling
layers have proven to be the most successful models in image classification
tasks. However, their large receptive fields and increased spatial invariance
(due to pooling and convolutional layers) make it difficult to infer the
locations of various objects within the image. In other words, the deeper we
go into a network, the more it seemed we lose the location information of
various objects within the image. Some researchers have attempted to overcome
this hurdle by using various pre- or post-processing techniques. However, the
introduction of fully convolutional network architectures, such as FCN Long et
al., (2015) and DeconvNet Noh et al., (2015) have shown that it is possible to
achieve good semantic labeling results using a convolutional network only,
with no pre- or post-processing steps required. FCN accomplish this by
replacing the fully-connected layers of the CNN with 1D convolutional layers
that produce coarse feature maps. These coarse feature maps are then
upsampled, and concatenated with the scores from intermediate feature maps in
the network to generate the output. These upsampling steps, however, result in
a blurred output that loses some of the resolution of the original image.

Deconvolution networks, on the other hand, overcome this problem by using a
symmetric encoder-decoder style architecture composed of stacks of convolution
and pooling layers in the encoder, and stacks of deconvolution and unpooling
layers in the decoder that mirror the encoders architecture. The role of the
encoder can be seen as doing object detection and classification, while the
decoder is used for accurate localization of these objects within the image.
This architecture can achieve finer and more accurate results than those of
the FCN, and therefore is adopted in our work.

A few recent papers have illustrated the successful application of
deconvolution networks for seismic interpretation applications Alaudah et al.,
(2018); Di et al., (2018).

Figure 7 illustrates the architecture of the deconvolution network used for
both of our baseline models. Every convolution or deconvolution layer (in
white) is followed by a rectified linear unit (ReLU) non-linearity. The layers
in red perform $2\times 2$ max pooling to select the maximum filter response
within small windows. The indices of the maximum responses for every pooling
layer are then shared with their respective unpooling layers (in green) to
undo this pooling operation and get a higher resolution image.

###  3.2 Baseline Models

In this work, we use two baseline models; a patch-based and a section-based
model. These two models use the exact same architecture, optimizer, and
hyperparameters but differ in the way they are trained and the way they are
used to label the seismic volume.

####  3.2.1 Patch-based model:

The patch-based model is trained on small patches extracted from the inlines
and crosslines of the training data. For very large seismic volumes, this
approach can be more feasible than using entire sections for training. At
training time, the patches of seismic data and their associated labels are
sampled randomly from the inlines and crosslines of the training set. During
test time, the model samples overlapping patches in the inline and crossline
direction and averages the results to generate a 2D labeled version of the
test inline or crossline. This is done for all inlines and crosslines in the
test sets.

####  3.2.2 Section-based model:

The section-based model is trained on entire inline and crossline sections.
The advantage of this approach is two-fold. First, since the network is fed an
entire section, it can easily learn the relationships between different
lithostratigraphic units and can take the depth information into account when
labeling the section. The second advantage is more practical. Training and
testing entire sections at once means the network can be trained or tested
very quickly since there are only a relatively small number of seismic inlines
and crosslines†††This is assuming the GPU memory is large enough to handle the
size of the seismic sections. On our Nvidia Titan X GPU, we trained the
baseline section-based network – eight sections at a time – in about 70
minutes.. One advantage of using a fully convolutional architecture (such as
the one we are using) is that the size of the network input does not have to
be fixed. The size of the output of the network changes as the size of its
input change. Therefore, the different size of the inline and crossline
sections does not pose any problem to the training of this network‡‡‡While the
sizes of the inlines and crosslines do not need to match, their resolutions
(in terms of meters/pixel) should. In our case, pixels in the inline and
crossline directions are both 25m$\times$25m..

####  3.2.3 Other variations:

In addition to the baseline patch- and section-based models, we have trained
other variations of these models to test how they can be improved. We have
tested the following variations:

  * •

Baseline + data augmentation: data augmentation applies different label-
preserving transformations to the training data such as small rotations,
random horizontal flipping, and the addition of Gaussian noise. This can help
increase the training sample size, and help the network generalize better to
the test data.

  * •

Baseline + data augmentation + skip connections: we further improve on the
previous model by adding skip connections. In a deep neural network, the
output of a layer is typically passed on as the input to the next layer in the
network. Skip connections allow the output of a layer to be also passed as an
input to a layer farther up the network, skipping intermediate layers in the
process. These connections are implemented by directly adding the outputs of
various layers in the encoder part of the deconvolution network to the outputs
of the corresponding layers in the decoder. Skip connections help networks
overcome the vanishing gradient problem Hochreiter et al., (2001) by providing
“shortcuts” for the computed gradients to propagate to the lower layers of the
network.

##  4 Experimental Setup

Figure 8: A 3D view of the F3 block from above with the  Zechstein  Group
shown in red, while the  Rijnland/Chalk  Group is shown in a semi-transparent
beige color. Inline 300 and crossline 1000 divide the survey into four
regions. The NW region of the survey is used for training, while the SW region
constitutes test set #1. The remaining region East of crossline 1000
constitutes test set #2.

In this section, we will introduce the main elements of the experimental
setup, including how the final geological model was produced, how the model is
split into training and testing sets, and what metrics are used to objectively
evaluate the performance.

###  4.1 The geological model

The final geological model that we use to train and test our models is not the
entire volume shown in Figure 4. The time-depth conversion process of the
seismic data resulted in some artifacts. These artifacts were concentrated
along the sigmoidal structure in the Upper North Sea group. Due to these
artifacts, and missing data on the sides of the survey, we only use the data
between inlines 100 and 701, crosslines 300 and 1201, and depth between 1005
and 1877 meters. Furthermore, we combine the Rijnland and the Chalk groups in
our final model to a single class due to various issues with processing the
Rijnland/Chalk boundary when generating the final model. Table 1 shows the
percentage of different classes in our training set.

In addition to the final model labels and seismic data, we also release the
original horizons for all the lithostratigraphic units, in addition to the
extracted fault planes from all three generations.

Table 1: The percentage of pixels from different classes in the training set. Zechstein | Scruff | Rijnland/Chalk | Lower N. S. | Middle N. S. | Upper N. S.
---|---|---|---|---|---
1.48% | 3.17% | 6.53% | 48.44% | 11.89% | 28.49%

###  4.2 The train/test split

Careful selection of the training and testing sets is crucial in any machine
learning application. This is especially true in seismic data, where
neighboring sections are highly correlated. Selecting the training and testing
sections randomly will lead to artificially good test results, that are not
representative of the actual generalization performance of the tested models.
Therefore, it is important to minimize the correlation between the training
and testing sets as much as possible. It is also important to ensure that both
the training and testing sets have adequate representation of all the classes
in the dataset.

Therefore, we decide to split the data as shown in Figure 8. Namely, the data
is split into the following three sets:

  1. 1.

Training set: This includes all the data in the range of inlines [300,700] and
crosslines [300,1000].

  2. 2.

Test set #1: This set includes all the data in the range of inlines [100,299]
and crosslines [300,1000].

  3. 3.

Test set #2: This sets includes all the data in the ranges of inlines
[100,700] and crosslines [1001,1200]. This set includes a large Zechstein
diapir in the NE of the survey that is never seen in the training set.

For a fair comparison with others who might use this benchmark in the future,
it is important to note that the test sets should not be used more than once.
Testing a model on the test set, then retraining that model with different
parameters means that the test set has been used for validation, which defeats
its purpose.

###  4.3 Evaluation metrics

To objectively evaluate the performance of different models on our two test
sets, we use the following metrics: pixel accuracy (PA), class accuracy (CA)
for each individual class, mean class accuracy (MCA) for all classes, and
frequency-weighted intersection over union (FWIU). These metrics are detailed
in the appendix.

##  5 Results

Table 2: Results of our two baseline models, with other variations, when
tested on both test splits of our dataset. All metrics are in the range
$[0,1]$, with larger values being better. The best performing model for every
metric is highlighted in bold.

Model Metric | PA | Class Accuracy | MCA | FWIU
---|---|---|---|---
Zechstein | Scruff | Rijnland/Chalk | Lower N. S. | Middle N. S. | Upper N. S.
Patch-based model | 0.788 | 0.264 | 0.074 | 0.499 | 0.992 | 0.804 | 0.754 | 0.565 | 0.640
Patch-based + aug. | 0.852 | 0.434 | 0.221 | 0.707 | 0.974 | 0.884 | 0.916 | 0.689 | 0.743
Patch-based + aug + skip | 0.862 | 0.458 | 0.286 | 0.673 | 0.974 | 0.912 | 0.926 | 0.705 | 0.757
Section-based model | 0.879 | 0.219 | 0.539 | 0.744 | 0.951 | 0.872 | 0.973 | 0.716 | 0.789
Section-based + aug. | 0.901 | 0.714 | 0.423 | 0.812 | 0.979 | 0.940 | 0.956 | 0.804 | 0.844
Section-based + aug + skip | 0.905 | 0.602 | 0.674 | 0.772 | 0.941 | 0.938 | 0.974 | 0.817 | 0.832

(a) Seismic data

(b) Ground truth labels

(c) Patch-based baseline

(d) Section-based baseline

(e) Patch-based + aug

(f) Section-based + aug

(g) Patch-based + aug + skip

(h) Section-based + aug + skip

Figure 9: The results of the different models on inline 200 from test set #1.
The color map is shown in Table 2.

After the final geological model was created, we train each of the models
described earlier on the training set until its training loss converges. We
note that on our Nvidia Titan X GPU, the baseline patch-based model, and the
augmented version, converged after 16 hours of training. The patch-based model
with skip connections required less than 5 hours to converge. The section-
based models required significantly less time, all of them converging in less
than 90 minutes. We test these models by using them to label all inlines and
crosslines in both test sets, and computing the performance metrics on the
final result. Table 2 summarizes the objective results for all the models that
we have tested on both test sets§§§The results for test set #1 and test set #2
are available from https://bit.ly/2D4QTbH, while Figure 9 shows inline 200 of
test set #1 labeled using the six different models we have tested. In the
remainder of this section, we will discuss these results, and suggest various
methods to improve upon them.

####  5.0.1 Patch-based vs section-based models:

Since the patch-based models are trained on patches from different depths in
the data, they can easily confuse various classes that typically exist at
different depths. For example, the patch-based models in Figure 9 often
confuse the Scruff group in the bottom left of the image with the Lower North
Sea group, while the section-based models do not suffer from these problems as
often. Figure 10 shows the confusion matrices for the baseline patch and
section models. It shows how the patch-based model confuses many classes in
our test sets with the Lower North Sea group. The baseline section-based model
is better at classifying the other classes as Figure 10(b) shows.

Table 2 shows that both patch and section-based models perform fairly well on
the North Sea groups, with the section-based models performing better.
However, for smaller classes such as the Scruff and Zechstein groups, the
section-based models show a clear advantage. The MCA score shows a 15%
improvement of the section-based baseline model vs. the patch-based model.
Overall, section-based models are superior to patch-based models due to their
ability to incorporate spatial and contextual information within each seismic
section. They also have the advantage of being faster to train and test.

(a) Patch-based baseline model

(b) Section-based baseline model

Figure 10: Confusion matrices for our two baseline models on both test set #1
and #2. Each row shows the distribution of the model output for each class.

####  5.0.2 Imbalanced classes:

As Table 1 shows, our dataset is highly imbalanced. The Zechstein and Scruff
groups are far smaller than the Lower or Upper North Sea groups. This means
while training, the network is trained on far more examples of Lower or Upper
North Sea groups than of Zechstein or Scruff groups for example. This leads to
the networks being biased towards classifying pixels as Lower or Upper North
Sea groups, and therefore artificially achieve high class accuracy scores for
those classes. However, this is at the expense of the very poor performance
for the smaller classes. In addition, the larger classes have a more diverse
visual appearance compared to smaller ones, and therefore, it would be easier
for the network to confuse features from smaller classes with those learned
from larger ones.

We do not make any changes to the baseline models to overcome this imbalance.
However, using various techniques to overcome this class imbalance can
significantly improve the results, especially for the smaller classes, such as
Zechstein and Scruff groups.

####  5.0.3 Data augmentation:

Data augmentation is a technique to increase the size of the training set
artificially. This is quite useful when training a large network with a
limited amount of training data. We use simple augmentation operations
including randomly rotating the patch or the section by up to $\pm
15^{\circ}$, adding random Gaussian noise, and randomly flipping the patches
or the sections horizontally. Using data augmentation significantly improved
the results for both baseline models, but especially the patch-based model.
The FWIU and MCA scores increased by more than 10% in the patch-based model,
and significantly improved the results for smaller classes such as the
Zechstein and Scruff groups. The results of the section-based model were also
improved by using data augmentation, although to a lesser degree.

####  5.0.4 Skip connections:

For both patch- and section-based models, adding skip connections can improve
the results, and speed up the training. This is especially noticeable in the
patch-based model where adding skip connections to the baseline + aug model
improved the results by about 1% in the PA metric and about 1.5% in the MCA
and FWIU metrics. The improvement in the results in the section-based models
are more subtle, as adding skip connections only improved the PA result by
0.1%. Interestingly, the Scruff group which worst is the performing class in
both the patch and section-based model seemed to benefit the most from the
addition of skip connections. The class accuracy score for the Scruff group
increased by 6.5% and 25% in the patch and section based models respectively.
Overall, adding skip connections seems to help improve the results. It also
can speed up to the training process. In the case of the patch-based model,
the skip connection model converged four times faster than the baseline.

##  6 Conclusions

In conclusion, we have introduced and made publicly available a new annotated
dataset for facies classification. This dataset includes six different
lithostratigraphic classes based on the underlying geology of the Netherlands
F3 block. The dataset also includes fault planes from three different
generations that we have identified in the F3 block. In addition, we present
two baseline deep learning models for facies classification, a patch- and a
section-based model, both based on a deconvolution network architecture. We
train these models using our dataset, and we analyze their performance.
Furthermore, we make the code for training and testing these models publicly
available for others to use. It is our hope that this dataset, and the code
that we have released, will help facilitate more research in this area, and
help create an objective benchmark for comparing the results of different
machine learning approaches for facies classification.

##  7 Acknowledgements

We want to acknowledge the support of the Center for Energy and Geo Processing
(CeGP) at the Georgia Institute of Technology and King Fahd University of
Petroleum and Minerals (KFUPM). We would also like to thank dGB Earth Sciences
for making the Netherlands F3 data publicly available, and Schlumberger for
providing educational licenses for the Petrel software.

Evaluation Metrics To objectively evaluate the performance of our models on
this dataset, we use a set of evaluation metrics that are commonly used in the
computer vision literature. If we denote the set of pixels that belong to
class $i$ as $G_{i}$, and the set of pixels classified as class $i$ as
$F_{i}$. Then, the set of correctly classified pixels is $G_{i}\cap F_{i}$. We
use $|\cdot|$ to denote the number of elements in a set. Now, we can define
the following metrics:

  * •

Pixel Accuracy (PA) is the percentage of pixels over all classes that are
correctly classified,

$\text{PA}=\frac{\sum_{i}|F_{i}\cap G_{i}|}{\sum_{i}|G_{i}|}.$ (1)

  * •

Class Accuracy for class $i$ ($\text{CA}_{i}$) is the percentage of pixels
that are correctly classified in a class $i$.

$\text{CA}_{i}=\frac{|F_{i}\cap G_{i}|}{|G_{i}|}.$ (2)

We will also define the Mean Class Accuracy (MCA) as the average of CA over
all classes,

$\text{MCA}=\frac{1}{n_{c}}\sum_{i}\text{CA}_{i}=\frac{1}{n_{c}}\sum_{i}\frac{|F_{i}\cap
G_{i}|}{|G_{i}|},$ (3)

where $n_{c}$ is the number of classes.

  * •

Intersection over Union ($\text{IU}_{i}$) is defined as the number of elements
of the intersection of $G_{i}$ and $F_{i}$ over the number of elements of
their union set,

$\text{IU}_{i}=\frac{|F_{i}\cap G_{i}|}{|F_{i}\cup G_{i}|}.$ (4)

This metric measures the overlap between the two sets and it should be $1$ if
and only if all pixels were correctly classified. Further, when we average IU
over all classes, we arrive at the Mean Intersection over Union (Mean IU),

$\text{Mean
IU}=\frac{1}{n_{c}}\sum_{i}\text{IU}_{i}=\frac{1}{n_{c}}\sum_{i}\frac{|F_{i}\cap
G_{i}|}{|F_{i}\cup G_{i}|}.$ (5)

To prevent this metric from being overly sensitive to small classes, it is
common to weigh each class by its size. The resulting metric is known as
Frequency-Weighted Intersection over Union (FWIU),

$\text{FWIU}=\frac{1}{\sum_{i}|\mathcal{G}_{i}|}\cdot\sum_{i}|\mathcal{G}_{i}|\cdot\frac{|\mathcal{F}_{i}\cup\mathcal{G}_{i}|}{|\mathcal{F}_{i}\cap\mathcal{G}_{i}|}.$
(6)

## References

  * Alaudah et al., (2019) Alaudah, Y., M. Alfarraj, and G. AlRegib, 2019, Structure label prediction using similarity-based retrieval and weakly supervised label mapping: GEOPHYSICS, 84, V67–V79.
  * Alaudah and AlRegib, (2016) Alaudah, Y., and G. AlRegib, 2016, Weakly-supervised labeling of seismic volumes using reference exemplars: 2016 IEEE International Conference on Image Processing (ICIP), 4373–4377.
  * Alaudah et al., (2018) Alaudah, Y., S. Gao, and G. AlRegib, 2018, Learning to label seismic structures with deconvolution networks and weak labels, in SEG Technical Program Expanded Abstracts 2018: Society of Exploration Geophysicists, 2121–2125.
  * Araya-Polo et al., (2017) Araya-Polo, M., T. Dahlke, C. Frogner, C. Zhang, T. Poggio, and D. Hohl, 2017, Automated fault detection without seismic processing: The Leading Edge, 36, 208–214.
  * Coléou et al., (2003) Coléou, T., M. Poupon, and K. Azbel, 2003, Unsupervised seismic facies classification: A review and comparison of techniques and implementation: The Leading Edge, 22, 942–953.
  * de Matos et al., (2006) de Matos, M. C., P. L. Osorio, and P. R. Johann, 2006, Unsupervised seismic facies analysis using wavelet transform and self-organizing maps: Geophysics, 72, P9–P21.
  * dGB Earth Sciences, (1987) dGB Earth Sciences, 1987, The Netherlands Offshore, The North Sea, F3 Block - Complete: https://bit.ly/2B44GO0.
  * Di et al., (2018) Di, H., Z. Wang, and G. AlRegib, 2018, Real-time seismic-image interpretation via deconvolutional neural network, in SEG Technical Program Expanded Abstracts 2018: Society of Exploration Geophysicists, 2051–2055.
  * Doornenbal, (2014) Doornenbal, J. C., 2014, Kilka uwag o kartografii wglebnej i modelowaniu geologicznym w holandii: Przeglad Geologiczny, 62, 806.
  * Dramsch and Lüthje, (2018) Dramsch, J. S., and M. Lüthje, 2018, Deep-learning seismic facies on state-of-the-art cnn architectures, in SEG Technical Program Expanded Abstracts 2018: Society of Exploration Geophysicists, 2036–2040.
  * Dubois et al., (2007) Dubois, M. K., G. C. Bohling, and S. Chakrabarti, 2007, Comparison of four approaches to a rock facies classification problem: Computers & Geosciences, 33, 599 – 617.
  * Duin et al., (2006) Duin, E., J. Doornenbal, R. Rijkers, J. Verbeek, and T. E. Wong, 2006, Subsurface structure of the netherlands-results of recent onshore and offshore mapping: Netherlands Journal of Geosciences, 85, 245.
  * Hochreiter et al., (2001) Hochreiter, S., Y. Bengio, P. Frasconi, J. Schmidhuber, et al., 2001, Gradient flow in recurrent nets: the difficulty of learning long-term dependencies, in A Field Guide to Dynamical Recurrent Networks: IEEE.
  * Huang et al., (2017) Huang, L., X. Dong, and T. E. Clee, 2017, A scalable deep learning platform for identifying geologic features from seismic attributes: The Leading Edge, 36, 249–256.
  * Long et al., (2015) Long, J., E. Shelhamer, and T. Darrell, 2015, Fully convolutional networks for semantic segmentation: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 3431–3440.
  * Mijnlieff, (2002) Mijnlieff, H., 2002, Top pre-permian distribution map and some thematic regional geologic maps of the netherlands: ICCP.
  * Noh et al., (2015) Noh, H., S. Hong, and B. Han, 2015, Learning deconvolution network for semantic segmentation: Proceedings of the IEEE International Conference on Computer Vision, 1520–1528.
  * Qian et al., (2018) Qian, F., M. Yin, X.-Y. Liu, Y.-J. Wang, C. Lu, and G.-M. Hu, 2018, Unsupervised seismic facies analysis via deep convolutional autoencoders: GEOPHYSICS, 83, A39–A43.
  * Remmelts, (1996) Remmelts, G., 1996, Salt tectonics in the southern north sea, the netherlands, in Geology of Gas and Oil under the Netherlands: Springer, 143–158.
  * Rutherford Ildstad and Bormann, (2017) Rutherford Ildstad, C., and P. Bormann, 2017, MalenoV: Tool for training and classifying SEGY seismic facies using deep neural networks: https://github.com/bolgebrygg/MalenoV.
  * Scheck-Wenderoth and Lamarche, (2005) Scheck-Wenderoth, M., and J. Lamarche, 2005, Crustal memory and basin evolution in the central european basin system—new insights from a 3d structural model: Tectonophysics, 397, 143–165.
  * Schroot and De Haan, (2003) Schroot, B., and H. De Haan, 2003, An improved regional structural model of the upper carboniferous of the cleaver bank high based on 3d seismic interpretation: Geological Society, London, Special Publications, 212, 23–37.
  * Shafiq et al., (2018) Shafiq, M. A., M. Prabhushankar, H. Di, and G. AlRegib, 2018, Towards understanding common features between natural and seismic images, in SEG Technical Program Expanded Abstracts 2018: Society of Exploration Geophysicists, 2076–2080.
  * Shi et al., (2018) Shi, Y., X. Wu, and S. Fomel, 2018, Automatic salt-body classification using a deep convolutional neural network, in SEG Technical Program Expanded Abstracts 2018: Society of Exploration Geophysicists, 1971–1975.
  * Van Adrichem Bogaert and Kouwe, (1993) Van Adrichem Bogaert, H., and W. Kouwe, 1993, 1997: Stratigraphic nomenclature of the Netherlands, revision and update. Mededelingen Rijks Geologische Dienst, 50.
  * Veillard et al., (2018) Veillard, A., O. Morére, M. Grout, and J. Gruffeille, 2018, Fast 3d seismic interpretation with unsupervised deep learning: Application to a potash network in the north sea: Presented at the 80th EAGE Conference and Exhibition 2018.
  * Waldeland and Solberg, (2017) Waldeland, A., and A. Solberg, 2017, Salt classification using deep learning: Presented at the 79th EAGE Conference and Exhibition 2017.
  * Zhao, (2018) Zhao, T., 2018, Seismic facies classification using different deep convolutional neural networks, in SEG Technical Program Expanded Abstracts 2018: Society of Exploration Geophysicists, 2046–2050.
  * Zhao et al., (2015) Zhao, T., V. Jayaram, A. Roy, and K. J. Marfurt, 2015, A comparison of classification techniques for seismic facies recognition: Interpretation, 3, SAE29–SAE58.
  * Ziegler, (1988) Ziegler, P. A., 1988, Evolution of the arctic-north atlantic and the western tethys: A visual presentation of a series of paleogeographic-paleotectonic maps: AAPG memoir, 43, 164–196.
  * Ziegler, (1990) ——–, 1990, Geological atlas of western and central europe: Presented at the , Geological Society of London, Shell Internat. Petroleum Maatschappij.


# Faster Population Counts Using AVX2 Instructions

Wojciech Muła<EMAIL_ADDRESS>Nathan Kurz<EMAIL_ADDRESS>Daniel Lemire  ⋆Université du Québec (TELUQ), Canada<EMAIL_ADDRESS>5800
Saint-Denis, Montreal (Quebec) H2S 3L5 Canada

###### Abstract

Counting the number of ones in a binary stream is a common operation in
database, information-retrieval, cryptographic and machine-learning
applications. Most processors have dedicated instructions to count the number
of ones in a word (e.g., popcnt on x64 processors). Maybe surprisingly, we
show that a vectorized approach using SIMD instructions can be twice as fast
as using the dedicated instructions on recent Intel processors. The benefits
can be even greater for applications such as similarity measures (e.g., the
Jaccard index) that require additional Boolean operations. Our approach has
been adopted by LLVM: it is used by its popular C compiler (Clang).

###### keywords:

Software Performance; SIMD Instructions; Vectorization; Bitset; Jaccard Index

##  1 Introduction

We can represent all sets of integers in $\\{0,1,\ldots,63\\}$ using a single
64-bit word. For example, the word 0xAA (0b10101010) represents the set
$\\{1,3,5,7\\}$. Intersections and unions between such sets can be computed
using a single bitwise logical operation on each pair of words (AND, OR). We
can generalize this idea to sets of integers in $\\{0,1,\ldots,n-1\\}$ using
$\lceil n/64\rceil$ 64-bit words. We call such data structures _bitsets_ ;
they are also known as a bit vectors, bit arrays or bitmaps. Bitsets are
ubiquitous in software, found in databases [1], version control systems [2],
search engines [SPE:SPE2325, 3, 4], and so forth. Languages such as Java and
C++ come with their own bitset classes (java.util.BitSet and std::bitset
respectively).

The cardinality of a bitset (the number of one bits, each representing an
element in the set) is commonly called a population count, a popcount, a
Hamming weight, a sideways addition, or sideways sum. For example, the
population counts of the words 0xFFFF, 0xAA and 0x00 are 16, 4 and 0
respectively. A frequent purpose for the population count is to determine the
size of the intersection or union between two bitsets. In such cases, we must
first apply a logical operation on pairs of words (AND, OR) and then compute
the population count of the resulting words. For example, the cardinality of
the intersection of the sets $A=\\{4,5,6,7\\}$ and $B=\\{1,3,5,7\\}$
represented by the words 0xF0 and 0xAA can be computed as $|A\cap
B|=\operatorname{popcount}(\texttt{0xF0}\operatorname{AND}\texttt{0xAA})=\operatorname{popcount}(\texttt{0xA0})=2$.

Population-count functions are used in cryptography [5], e.g., as part of
randomness tests [6] or to generate pseudo-random permutations [7]. They can
help find duplicated web pages [8]. They are frequently used in bioinformatics
[9, 10, 11], ecology [12], chemistry [13], and so forth. Gueron and Krasnov
use population-count instructions as part of a fast sorting algorithm [14].

The computation of the population count is so important that commodity
processors have dedicated instructions: popcnt for x64 processors and cnt for
the 64-bit ARM architecture.111The x64 popcnt instruction was first available
in the Nehalem microarchitecture, announced in 2007 and released in November
2008. The ARM cnt instruction was released as part of the Cortex-A8
microarchitecture, published in 2006 [15]. The x64 popcnt instruction is fast:
on recent Intel processors, it has a throughput of one instruction per cycle
[16] (it can execute once per cycle) and a latency of 3 cycles (meaning that
the result is available for use on the third cycle after execution). It is
available in common C and C++ compilers as the intrinsic mmpopcntu64. In Java,
it is available as the Long.bitCount intrinsic.

Commodity PC processors also support Single-Instruction-Multiple-Data (SIMD)
instructions. Starting with the Haswell microarchitecture (2013), Intel
processors support the AVX2 instruction set which offers rich support for
256-bit vector registers. The contest between a dedicated instruction
operating on 64-bits at a time (popcnt) and a series of vector instructions
operating on 256-bits at a time (AVX2) turns out to be interesting. In fact,
we show that we can achieve twice the speed of the an optimized popcnt-based
function using AVX2: 0.52 versus 1.02 cycles per 8 bytes on large arrays. Our
claim has been thoroughly validated: at least one major C compiler (LLVM’s
Clang) uses our technique [17].

Thus, in several instances, SIMD instructions might be preferable to dedicated
non-SIMD instructions if we are just interested in the population count of a
bitset. But what if we seek the cardinality of the intersection or union, or
simply the Jaccard index between two bitsets? Again, the AVX2 instructions
prove useful, more than doubling the speed ($2.4\times$) of the computation
against an optimized function using the popcnt instruction.

##  2 Existing Algorithms and Related Work

Conceptually, one could compute the population count by checking the value of
each bit individually by calling count += (word >> i) & 1 for i ranging from 0
to 63, given that word is a 64-bit word. While this approach scales linearly
in the number of input words, we expect it to be slow since it requires
multiple operations for each bit in each word. It is $O(n)$—the best that we
can do—but with a high constant factor.

population count of 4 bitssum of bits 1 & 21${}^{\text{st}}$
bit2${}^{\text{nd}}$ bitsum of bits 3 & 43${}^{\text{rd}}$
bit4${}^{\text{th}}$ bit Figure 1:  A tree of adders to compute the population
count of four bits in two steps.

Instead, we should prefer approaches with fewer operations per word. We can
achieve the desired result with a tree of adders and bit-level parallelism. In
Fig. 1, we illustrate the idea over words of 4 bits (for simplicity). We
implement this approach with two lines of code.

  1. 1.

We can sum the individual bits to 2-bit subwords with the line of C code: ( x
& 0b0101 ) + ( ( x >> 1 ) & 0b0101 ). This takes us from the bottom of the
tree to the second level. We say to this step exhibits _bit-level parallelism_
since two sums are executed at once, within the same 4-bit word.

  2. 2.

We can then sum the values stored in the 2-bit subwords into a single 4-bit
subword with another line of C code: ( x & 0b0011 ) + ( ( x >> 2 ) & 0b0011 ).

Fig. 2 illustrates a non-optimized (naive) function that computes the
population count of a 64-bit word in this manner.

⬇ uint64_t c1 = UINT64_C(0x5555555555555555); uint64_t c2 =
UINT64_C(0x3333333333333333); uint64_t c4 = UINT64_C(0x0F0F0F0F0F0F0F0F);
uint64_t c8 = UINT64_C(0x00FF00FF00FF00FF); uint64_t c16 =
UINT64_C(0x0000FFFF0000FFFF); uint64_t c32 = UINT64_C(0x00000000FFFFFFFF);
uint64_t count(uint64_t x) { x = (x & c1) + ((x >> 1) & c1); x = (x & c2) +
((x >> 2) & c2); x = (x & c4) + ((x >> 4) & c4); x = (x & c8) + ((x >> 8) &
c8); x = (x & c16)+ ((x >> 16)) & c16); return (x & c32) + ((x >> 32) & c32);
}
---
Figure 2: A naive tree-of-adders function in C ⬇ uint64_t c1 =
UINT64_C(0x5555555555555555); uint64_t c2 = UINT64_C(0x3333333333333333);
uint64_t c4 = UINT64_C(0x0F0F0F0F0F0F0F0F); uint64_t count(uint64_t x) { x -=
(x >> 1) & c1; x = (( x >> 2) & c2) + (x & c2); x = ( x + (x >> 4) ) & c4; x
*= UINT64_C(0x0101010101010101); return x >> 56; }
---
Figure 3: The Wilkes-Wheeler-Gill function in C

A fast and widely used tree-of-adder function to compute the population count
has been attributed by Knuth [18] to a 1957 textbook by Wilkes, Wheeler and
Gill [19]: see Fig. 3. It involves far fewer than 64 instructions and we
expect it to be several times faster than a naive function checking the values
of each bit and faster than the naive tree-of-adder approach on processor with
a sufficiently fast 64-bit integer multiplication (which includes all x64
processors).

  * •

The first two lines in the count function correspond to the first two levels
of our simplified tree-of-adders count function from Fig. 1. The first line
has been optimized. We can verify the optimization by checking that for each
possible 2-bit word, we get the sum of the bit values:

    * –

0b11 - 0b01 = 0b10 = 2,

    * –

0b10 - 0b01 = 0b01 = 1,

    * –

0b01 - 0b00 = 0b01 = 1,

    * –

0b00 - 0b00 = 0b00 = 0.

  * •

After the first two lines, we have 4-bit population counts (in
$\\{0b0000,0b0001,0b0010,0b0011,0b0100\\}$) stored in 4-bit subwords. The next
line sums consecutive 4-bit subwords to bytes. We use the fact that the most
significant bit of each 4-bit subword is zero.

  * •

The multiplication and final shift sum all bytes in an efficient way.
Multiplying x by 0x0101010101010101 is equivalent to summing up x, x << 8, x
<< 16, …, x << 56. The total population count is less than 64, so that the sum
of all bytes from x fits in a single byte value (in $[0,256)$). In that case,
the most significant 8 bits from the product is the sum of all eight byte
values.

Knuth also attributes another common technique to Wegner [20] (see Fig. 4)
that could be competitive when the population count is relatively low (e.g.,
less than 4 one bit per 64-bit word). When the population count is expected to
be high (e.g., more than 60 one bit per 64-bit words), one could simply negate
the words prior to using the function so as to count the number of zeros
instead. The core insight behind the Wegner function is that the line of C
code x &= x - 1 sets to zero the least significant bit of x, as one can
readily check. On an x64 processor, the expression x &= x - 1 might be
compiled to the blsr (reset lowest set bit) instruction. On current generation
processors, this instruction achieves a throughput of two instructions per
cycle with a latency of one cycle [16]. The downside of the Wegner approach
for modern processors is that the unpredictable loop termination adds a
mispredicted branch penalty of at least 10 cycles [21], which for short loops
can be more expensive than the operations performed by the loop.

⬇ int count(uint64_t x) { int v = 0; while(x != 0) { x &= x - 1; v++; } return
v; }
---
Figure 4: The Wegner function in C.

Another simple and common technique is based on tabulation. For example, one
might create a table that contains the corresponding population count for each
possible byte value, and then look up and sum the count for each byte. Such a
table would require only 256 bytes. A population count for a 64-bit word would
require only eight table look-ups and seven additions. On more powerful
processor, with more cache, it might be beneficial to create a larger table,
such as one that has a population count for each possible short value (2
bytes) using 64 KB. Each doubling of the bit-width covered by the table halves
the number of table lookups, but squares the memory required for the table.

We can improve the efficiency of tree-of-adders techniques by _merging_ the
trees across words [22]. To gain an intuition for this approach, consider that
in the Wilkes-Wheeler-Gill approach, we use 4-bit subwords to store the
population count of four consecutive bits. Such a population count takes a
value in $\\{0,1,2,3,4\\}$, yet a 4-bit integer can represent all integers in
$[0,16)$. Thus, as a simple optimization, we could accumulate the 4-bit counts
across three different words instead of a single one. Next consider that if
you sum two 4-bit subwords (representing integers in $[0,16)$) the result is
in $[0,32)$ whereas an 8-bit subword (a byte) can represent all integers in
$[0,256)$, a range that is four times larger. Hence, we can accumulate the
counts over four triple of words. These two optimizations combined lead to a
function to compute the population count of twelve words at once (see Fig. 5)
faster than would be possible if we processed each word individually.

⬇ uint64_t count(uint64_t *input) { uint64_t m1 =
UINT64_C(0x5555555555555555); uint64_t m2 = UINT64_C(0x3333333333333333);
uint64_t m4 = UINT64_C(0x0F0F0F0F0F0F0F0F); uint64_t m8 =
UINT64_C(0x00FF00FF00FF00FF); uint64_t m16= UINT64_C(0x0000FFFF0000FFFF);
uint64_t acc = 0; for (int j = 0; j < 12; j += 3) { uint64_t count1 = input[j
+ 0]; uint64_t count2 = input[j + 1]; uint64_t half1 = input[j + 2]; uint64_t
half2 = input[j + 2]; half1 &= m1; half2 = (half2 >> 1) & m1; count1 -=
(count1 >> 1) & m1; count2 -= (count2 >> 1) & m1; count1 += half1; count2 +=
half2; count1 = (count1 & m2) + ((count1 >> 2) & m2); count1 += (count2 & m2)
+ ((count2 >> 2) & m2); acc += (count1 & m4) + ((count1 >> 4) & m4); } acc =
(acc & m8) + ((acc >> 8) & m8); acc = (acc + (acc >> 16)) & m16; acc = acc +
(acc >> 32); return acc; }
---
Figure 5: The Lauradoux population count in C for sets of 12 words.

However, even before Lauradoux proposed this improved function, Warren [23]
had presented a superior alternative attributed to a newsgroup posting from
1997 by Seal, inspired from earlier work by Harley. This approach, henceforth
called Harley-Seal, is based on a carry-save adder (CSA). Suppose you are
given three bit values ($a,b,c\in\\{0,1\\}$) and you want to compute their sum
($a+b+c\in\\{0,1,2,3\\}$). Such a sum fits in a 2-bit word. The value of the
least significant bit is given by $(a\operatorname{XOR}b)\operatorname{XOR}c$
whereas the most significant bit is given by
$(a\operatorname{AND}b)\operatorname{OR}((a\operatorname{XOR}b)\operatorname{AND}c)$.
Table 1 illustrates these expressions: the least significant bit
($(a\operatorname{XOR}b)\operatorname{XOR}c$) takes value 1 only when $a+b+c$
is odd and the most significant bit takes value 1 only when two or three of
the input bits ($a,b,c$) are set to 1. There are many possible expressions to
compute the most significant bit, but the chosen expression is convenient
because it reuses the $a\operatorname{XOR}b$ expression from the computation
of the least significant bit. Thus, we can sum three bit values to a 2-bit
counter using 5 logical operations. We can generalize this approach to work on
all 64-bits in parallel. Starting with three 64-bit input words, we can
generate two new output words: $h$, which holds the 64 most significant bits,
and $l$, which contains the corresponding 64 least significant bits. We
effectively compute 64 sums in parallel using bit-level parallelism. Fig. 6
presents an efficient implementation in C of this idea. The function uses 5
bitwise logical operations (two XORs, two ANDs and one OR): it is optimal with
respect to the number of such operations [24, 7.1.2]. However, it requires at
least three cycles to complete due to data dependencies.

⬇ void CSA(uint64_t* h, uint64_t* l, uint64_t a, uint64_t b, uint64_t c) { uint64_t u = a ^ b; *h = (a & b) | (u & c); *l = u ^ c; }
---
Figure 6: A C function implementing a bitwise parallel carry-save adder (CSA). Given three input words $a,b,c$, it generates two new words $h,l$ in which each bit represents the high and low bits in the bitwise sum of the bits from $a$, $b$, and $c$.  Table 1: Sum of three bits $a+b+c$. We use $\oplus$ for XOR, $\land$ for AND and $\lor$ for OR. $a$ | $b$ | $c$ | $a+b+c$ |  $(a\oplus b)\oplus c$ |  $(a\land b)\lor\allowbreak((a\oplus b)\land c)$
---|---|---|---|---|---
0 | 0 | 0 | 0 |  0 |  0
0 | 0 | 1 | 1 |  1 |  0
0 | 1 | 0 | 1 |  1 |  0
1 | 0 | 0 | 1 |  1 |  0
0 | 1 | 1 | 2 |  0 |  1
1 | 0 | 1 | 2 |  0 |  1
1 | 1 | 0 | 2 |  0 |  1
1 | 1 | 1 | 3 |  1 |  1
… twos
(input) ones
(input) $d_{\\_}i$ (input) $d_{\\_}{i+1}$ (input) CSA $d_{\\_}{i+2}$ (input)
$d_{\\_}{i+3}$ (input) CSA…CSA fours
(output) twos
(output) ones
(output) Figure 7: Harley-Seal algorithm aggregating four new inputs
($d_{\\_}i,d_{\\_}{i+1},d_{\\_}{i+2}$, $d_{\\_}{i+3}$) to inputs ones and
twos, producing new values of ones, twos and fours.

From such a CSA function, we can derive an efficient population count. Suppose
we start with three words serving as counters (initialized at zero): one for
the least significant bits (henceforth ones), another one for the second least
significant bits (twos, so named because each bit set represents 2 input
bits), and another for the third least significant bits (fours, representing 4
input bits). We can proceed as follows; the first few steps are illustrated in
Fig. 7. We start with a word serving as a population counter $c$ (initialized
at zero). Assume with we have a number of words $d_{\\_}1,d_{\\_}2,\ldots$
divisible by 8. Start with $i=0$.

  * •

Load two new words ($d_{\\_}i,d_{\\_}{i+1}$). Use the CSA function to sum
ones, $d_{\\_}i$ and $d_{\\_}{i+1}$, write the least significant bit of the
sum to ones and store the carry bits in a temporary register (noted twosA). We
repeat with the next two input words. Load $d_{\\_}{i+2},d_{\\_}{i+3}$, use
the CSA function to sum ones, $d_{\\_}i$ and $d_{\\_}{a+i}$, write the least
significant bit of the sum to ones and store the carry bits in a temporary
register (noted twosB).

  * •

At this point, we have three words containing second least significant bits
(twos, twosA, twosB). We sum them up using a CSA, writing back the result to
twos and the carry bits to a temporary register foursA.

  * •

We do with $d_{\\_}{i+4},d_{\\_}{i+5}$ and $d_{\\_}{i+6},d_{\\_}{i+7}$ as we
did with $d_{\\_}i,d_{\\_}{i+1}$ and $d_{\\_}{i+2},d_{\\_}{i+3}$. Again we
have three words containing second least significant bits (twos, twosA,
twosB). We sum them up with CSA, writing the result to twos and to a carry-bit
temporary register foursB.

  * •

At this point, we have three words containing third least significant bits
(fours, foursA, foursB). We can sum them up with a CSA, write the result back
to fours, storing the carry bits in a temporary register eights.

  * •

We compute the population count of the word eights (e.g, using the Wilkes-
Wheeler-Gill population count) and increment the counter $c$ by the population
count.

  * •

Increment $i$ by 8 and continue for as long as we have new words.

When the algorithm terminates, multiply $c$ by 8. Compute the population count
of fours, multiply the result by 4 and add to $c$. Do similarly with twos and
ones. The counter $c$ contains the population count. If the number of input
words is not divisible by 8, adjust accordingly with the leftover words (e.g,
using the Wilkes-Wheeler-Gill population count).

In that particular implementation of this idea, we used blocks of eight words.
More generally, the Harley-Seal approach works with blocks of $2^{n}$ words
for $n=3,4,5,\ldots$ (8, 16, 32, …). We need $2^{n}-1$ CSA function calls when
using $2^{n}$ words, and one call to an auxiliary function (e.g., Wilkes-
Wheeler-Gill). If we expect the auxiliary function to be significantly more
expensive than the CSA function calls, then larger blocks should lead to
higher performance, as long as we have enough input data and many available
registers. In practice, we found that using blocks of sixteen words works well
on current processors (see Fig. 8). This approach is only worthwhile if we
have at least 16 input words (64-bits/word $\times$ 16 words $=128$ bytes).

⬇ uint64_t harley_seal(uint64_t * d, size_t size) { uint64_t total = 0, ones =
0, twos = 0, fours = 0, eights = 0, sixteens = 0; uint64_t twosA, twosB,
foursA, foursB, eightsA, eightsB; for(size_t i = 0; i < size - size % 16; i +=
16) { CSA(&twosA, &ones, ones, d[i+0], d[i+1]); CSA(&twosB, &ones, ones,
d[i+2], d[i+3]); CSA(&foursA, &twos, twos, twosA, twosB); CSA(&twosA, &ones,
ones, d[i+4], d[i+5]); CSA(&twosB, &ones, ones, d[i+6], d[i+7]); CSA(&foursB,
&twos, twos, twosA, twosB); CSA(&eightsA, &fours, fours, foursA, foursB);
CSA(&twosA, &ones, ones, d[i+8], d[i+9]); CSA(&twosB, &ones, ones,
d[i+10],d[i+11]); CSA(&foursA, &twos, twos, twosA, twosB); CSA(&twosA, &ones,
ones, d[i+12],d[i+13]); CSA(&twosB, &ones, ones, d[i+14],d[i+15]);
CSA(&foursB, &twos, twos, twosA, twosB); CSA(&eightsB, &fours, fours, foursA,
foursB); CSA(&sixteens, &eights, eights, eightsA, eightsB); total +=
count(sixteens); } total = 16 * total + 8 * count(eights) + 4 * count(fours) +
2 * count(twos) + count(ones); for(size_t i = size - size % 16 ; i < size;
i++) total += count(d[i]); return total; }
---
Figure 8: A C function implementing the Harley-Seal population count over an
array of 64-bit words. The count function could be the Wilkes-Wheeler-Gill
function.

The functions we presented thus far still have their uses when programming
with high-level languages without convenient access to dedicated functions
(e.g., JavaScript, Go) or on limited hardware. However, they are otherwise
obsolete when a sufficiently fast instruction is available, as is the case on
recent x64 processors with popcnt. The popcnt instruction has a reciprocal
throughput222The reciprocal throughput is the number of processor clocks it
takes for an instruction to execute. of one instruction per cycle. With a
properly constructed loop, the load-popcnt-add sequence can be executed in a
single cycle, allowing for a population count function that processes 64-bits
per cycle.

###  2.1 Existing Vectorized Algorithms

To our knowledge, the first published vectorized population count on Intel
processor was proposed by Muła in 2008 [25]. It is a vectorized form of
tabulation on 4-bit subwords. Its key ingredient is the SSSE3 vector
instruction pshufb (see Table 2). The pshufb instruction shuffles the input
bytes into a new vector containing the same byte values in a (potentially)
different order. It takes an input register $v$ and a control mask $m$,
treating both as vectors of sixteen bytes. Starting from
$v_{\\_}0,v_{\\_}1,\ldots,v_{\\_}{16}$, it outputs a new vector
$(v_{\\_}{m_{\\_}0},v_{\\_}{m_{\\_}1},v_{\\_}{m_{\\_}2},v_{\\_}{m_{\\_}3},\ldots,v_{\\_}{m_{\\_}{15}})$
(assuming that $0\leq m_{\\_}i<16$ for $i=0,1,\ldots,15$). Thus, for example,
if the mask $m$ is $0,1,2,\ldots,15$, then we have the identify function. If
the mask $m$ is $15,14,\ldots,0$, then the byte order is reversed. Bytes are
allowed to be repeated in the output vector, thus the mask $0,0,\ldots,0$
would produce a vector containing only the first input byte, repeated sixteen
times. It is a fast instruction with a reciprocal throughput and latency of
one cycle on current Intel processors, yet it effectively “looks up” 16 values
at once. In our case, we use a fixed input register made of the input bytes
$0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4$ corresponding to the population counts of
all possible 4-bit integers $0,1,2,3,\ldots,15$. Given an array of sixteen
bytes, we can call pshufb once, after selecting the least significant 4 bits
of each byte (using a bitwise AND) to gather sixteen population counts on
sixteen 4-bit subwords. Next, we right shift by four bits each byte value, and
call pshufb again to gather sixteen counts of the most significant 4 bits of
each byte. We can sum the two results to obtain sixteen population counts,
each corresponding to one of the sixteen initial byte values. See Fig. 9 for a
C implementation. If we ignore loads and stores as well as control
instructions, Muła’s approach requires two pshufb, two pand, one paddb, and
one psrlw instruction, so six inexpensive instructions to compute the
population counts of sixteen bytes. The Muła algorithm requires fewer
instructions than the part of Wilkes-Wheel-Gill that does the same work (see
Fig. 3), but works on twice as many input bytes per iteration.

⬇ __m128i count_bytes(__m128i v) { __m128i lookup =
_mm_setr_epi8(0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4); __m128i low_mask =
_mm_set1_epi8(0x0f); __m128i lo = _mm_and_si128(v, low_mask); __m128i hi =
_mm_and_si128( _mm_srli_epi16(v, 4), low_mask); __m128i cnt1 =
_mm_shuffle_epi8(lookup, lo); __m128i cnt2 = _mm_shuffle_epi8(lookup, hi);
return _mm_add_epi8(cnt1, cnt2); }
---
Figure 9: A C function using SSE intrinsics implementing Muła’s algorithm to
compute sixteen population counts, corresponding to sixteen input bytes.

The countbytes function from Fig. 9 separately computes the population count
for each of the sixteen input bytes, storing each in a separate byte of the
result. As each of these bytes will be in $[0,8]$, we can sum the result
vectors from up to 31 calls to countbytes using the mmaddepi8 intrinsic
without risk of overflow before using the psadbw instruction (using the
mmsadepu8 intrinsic) to horizontally sum the individual bytes into two 64-bit
counters. In our implementation, we found it adequate to call the countbytes
function eight times between each call to psadbw.

Morancho observed that we can use both a vector approach, like Muła’s, and the
popcnt in a hybrid approach [26]. Morancho proposed a family of hybrid schemes
that could be up to 22% faster than an implementation based on popcnt for
sufficiently large input arrays.

Table 2: Relevant SSE instructions with latencies and reciprocal throughput in CPU cycles on recent (Haswell) Intel processors . instruction | C intrinsic |  description | latency |  rec. throughput
---|---|---|---|---
paddb | mmaddepi8 |  add sixteen pairs of 8-bit integers | 1 | 0.5
pshufb | mmshuffleepi8 |  _shuffle_ sixteen bytes | 1 | 1
psrlw | mmsrliepi16 |  shift right eight 16-bit integers | 1 | 1
pand | mmandsi128 |  128-bit AND | 1 | 0.33
psadbw | mmsadepu8 |  sum of the absolute differences of the byte values to the low 16 bits of each 64-bit word | 5 | 1

##  3 Novel Vectorized Algorithms

Starting with the Haswell microarchiture released in 2013, Intel processors
support the AVX2 instruction set with 256-bit vectors, instead of the shorter
128-bit vectors. It supports instructions and intrinsics that are analogous to
the SSE intrinsics (see Table 2).333To our knowledge, Muła was first to
document the benefits of AVX2 for the population count problem in March 2016
[25].

The Muła function provides an effective approach to compute population counts
at a speed close to an x64 processor’s popcnt instruction when using 128-bit
vectors, but after upgrading to AVX2’s 256-bit vectors, it becomes faster than
functions using the popcnt instruction. We present the basis of such a
function in Fig. 10 using AVX2 intrinsics; the AVX2 intrinsics are analogous
to the SSE intrinsics (see Fig. 2). It returns a 256-bit word that can be
interpreted as four 64-bit counts (each having value in $[0,64]$). We can then
add the result of repeated calls with the mm256addepi64 intrinsic to sum
64-bit counts.

⬇ __m256i count(__m256i v) { __m256i lookup = _mm256_setr_epi8(0, 1, 1, 2, 1,
2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4, 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3,
4); __m256i low_mask = _mm256_set1_epi8(0x0f); __m256i lo = =
_mm256_and_si256(v, low_mask); __m256i hi =
_mm256_and_si256(_mm256_srli_epi32(v, 4), low_mask); __m256i popcnt1 =
_mm256_shuffle_epi8(lookup, lo); __m256i popcnt2 = _mm256_shuffle_epi8(lookup,
hi); __m256i total = _mm256_add_epi8(popcnt1,popcnt2); return
_mm256_sad_epu8(total, _mm256_setzero_si256()); }
---
Figure 10: A C function using AVX2 intrinsics implementing Muła’s algorithm to
compute the four population counts of the four 64-bit words in a 256-bit
vector. The 32 B output vector should be interpreted as four separate 64-bit
counts that need to be summed to obtain the final population count.

For a slight gain in performance, we can call the Muła function several times
while skipping the call to mm256sadepu8, adding the byte values with
mm256addepi8 before calling mm256sadepu8 once. Each time we call the Muła
function, we process 32 input bytes and get 32 byte values in $[0,8]$. We can
add sixteen totals before calling mm256sadepu8 to sum the results into four
64-bit words (since $8\times 16=128<2^{8}$), thus processing a block of 512
bytes per call.444We could call the Muła function up to 31 times, since
$8\times 31=248<2^{8}$.

Of all the non-vectorized (or _scalar_) functions, Harley-Seal approaches are
fastest. Thus we were motivated to port the approach to AVX2. The carry-save
adder that worked on 64-bit words (see Fig. 6) can be adapted in a straight-
forward manner to work with AVX2 intrinsics (see Fig. 11).

⬇ void CSA(__m256i* h, __m256i* l, __m256i a, __m256i b, __m256i c) { __m256i
u = _mm256_xor_si256(a , b); *h = _mm256_or_si256(_mm256_and_si256(a , b) ,
_mm256_and_si256(u , c) ); *l = _mm256_xor_si256(u , c); }
---
Figure 11: A C function using AVX2 intrinsics implementing a bitwise parallel
carry-save adder (CSA).

Fig. 12 presents an efficient Harley-Seal function using an AVX2 carry-save
adder. The function processes the data in blocks of sixteen 256-bit vectors
(512 B). It calls Muła’s AVX2 function (see Fig. 10).

⬇ uint64_t avx_hs(__m256i* d, uint64_t size) { __m256i total =
_mm256_setzero_si256(); __m256i ones = _mm256_setzero_si256(); __m256i twos =
_mm256_setzero_si256(); __m256i fours = _mm256_setzero_si256(); __m256i eights
= _mm256_setzero_si256(); __m256i sixteens = _mm256_setzero_si256(); __m256i
twosA, twosB, foursA, foursB, eightsA, eightsB; for(uint64_t i = 0; i < size;
i += 16) { CSA(&twosA, &ones, ones, d[i], d[i+1]); CSA(&twosB, &ones, ones,
d[i+2], d[i+3]); CSA(&foursA, &twos, twos, twosA, twosB); CSA(&twosA, &ones,
ones, d[i+4], d[i+5]); CSA(&twosB, &ones, ones, d[i+6], d[i+7]); CSA(&foursB,&
twos, twos, twosA, twosB); CSA(&eightsA,&fours, fours, foursA,foursB);
CSA(&twosA, &ones, ones, d[i+8], d[i+9]); CSA(&twosB, &ones, ones,
d[i+10],d[i+11]); CSA(&foursA, &twos, twos, twosA, twosB); CSA(&twosA, &ones,
ones, d[i+12],d[i+13]); CSA(&twosB, &ones, ones, d[i+14],d[i+15]);
CSA(&foursB, &twos, twos, twosA, twosB); CSA(&eightsB, &fours, fours, foursA,
foursB); CSA(&sixteens, &eights, eights, eightsA, eightsB);
total=_mm256_add_epi64(total, count(sixteens)); } total =
_mm256_slli_epi64(total, 4); total = _mm256_add_epi64(total,
_mm256_slli_epi64(count(eights), 3)); total = _mm256_add_epi64(total,
_mm256_slli_epi64(count(fours), 2)); total = _mm256_add_epi64(total,
_mm256_slli_epi64(count(twos), 1)); total
=_mm256_add_epi64(total,count(ones)); return _mm256_extract_epi64(total, 0) +
_mm256_extract_epi64(total, 1) + _mm256_extract_epi64(total, 2) +
_mm256_extract_epi64(total, 3); }
---
Figure 12: A C function using AVX2 intrinsics implementing Harley-Seal’s
algorithm. It assumes, for simplicity, that the input size in 256-bit vectors
is divisible by 16. See Fig. 10 for the count function.

Processors execute complex machine instructions using low-level instructions
called $\mu$ops.

  * •

Using the dedicated popcnt instruction for the population of an array of words
requires loading the word (movq), counting the bits (popcnt), and then adding
the result to the total (addq). The load and the popcnt can be combined into a
single assembly instruction, but internally they are executed as separate
µops, and thus each 64-bit word requires three $\mu$ops. Apart from minimal
loop overhead, these three operations can be executed in a single cycle on a
modern x64 superscalar processor, for a throughput of just over one cycle per
8 B word.

  * •

The AVX2 Harley-Seal function processes sixteen 256-bit vectors (512 B) with
98 $\mu$ops: 16 loads (vpmov), 32 bitwise ANDs (vpand), 15 bitwise ORs (vpor),
and 30 bitwise XORs (vpxor). Each 64-bit word (8 B) thus takes just over 1.5
$\mu$ops—about half as many as required to use the builtin popcnt instruction
on the same input.

While fewer $\mu$ops does does not guarantee faster execution, for
computationally intensive tasks such as this it often proves to be a
significant advantage. In this case, we find that it does in fact result in
approximately twice the speed.

##  4 Beyond population counts

In practice, we often want to compute population counts on the result of some
operations. For example, given two bitsets, we might want to determine the
cardinality of their intersection (computed as the bit-wise logical AND) or
the cardinality of their union (computed as the bit-wise logical OR). In such
instances, we need to load input bytes from the two bitsets, generate a
temporary word, process it to determine its population count, and so forth.
When computing the Jaccard index, given that we have no prior knowledge of the
population counts, we need to compute both the intersection and the union, and
then we need to compute the two corresponding population counts (see Fig. 13).

⬇ void popcnt_jaccard_index(uint64_t* A, uint64_t* B, size_t n) { double s = 0; double i = 0; for(size_t k = 0; k < n; k++) { s += _mm_popcnt_u64(A[k] | B[k]); i += _mm_popcnt_u64(A[k] & B[k]); } return i / s; }
---
Figure 13: A C function using the mmpopcntu64 intrinsic to compute the Jaccard
index of a pair of 64-bit inputs.

Both loads and logical operations benefit greatly from vectorization, and
hybrid scalar/vector approaches can be difficult because inserting and
extracting elements into and from vectors adds overhead. With AVX2, in one
operation, we can load a 256-bit register or compute the logical AND between
two 256-bit registers. This is four times the performance of the corresponding
64-bit operations. Thus we can expect good results from fast population count
functions based on AVX2 adapted for the computation of the Jaccard index, the
cardinality of the intersection or union, or similar operations.

##  5 Experimental Results

We implemented our software in C. We use a Linux server with an Intel i7-4770
processor running at $3.4\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$. This
Haswell processor has $32\text{\,}\mathrm{k}\mathrm{B}$ of L1 cache and
$256\text{\,}\mathrm{k}\mathrm{B}$ of L2 cache per core with
$8\text{\,}\mathrm{M}\mathrm{B}$ of L3 cache. The machine has
$32\text{\,}\mathrm{G}\mathrm{B}$ of RAM (DDR3-1600 with double-channel). We
disabled Turbo Boost and set the processor to run at its highest clock speed.
Our software is freely available (https://github.com/CountOnes/hamming_weight)
and was compiled using the GNU GCC 5.3 compiler with the “-O3 -march=native”
flags.

For all experiments, we use randomized input bits. However, we find that the
performance results are insensitive to the data values.

Table 3 presents our results in number of cycles per word, for single-threaded
execution. To make sure our results are reliable, we repeat each test 500
times and check that the minimum and the average cycle counts are within 1% of
each other. We report the minimum cycle count divided by the number of words
in the input. All the scalar methods (WWG, Laradoux, and HS) are significantly
slower than the native popcnt-based function. We omit tabulation-based
approaches from Table 3 because they are not competitive: 16-bit tabulation
uses over 5 cycles even for large arrays. We can see that for inputs larger
than 4 kB, the AVX2-based Harley-Seal approach is twice as fast as our
optimized popcnt-based function, while for small arrays (fewer 64 words) the
popcnt-based function is fastest.

Table 3: Number of cycles per 64-bit input word required to compute the population of arrays of various sizes. array size | WWG | Lauradoux | HS | popcnt | AVX2 Muła | AVX2 HS
---|---|---|---|---|---|---
256 B | 6.00 | 4.50 | 3.25 | 1.12 | 1.38 | —
512 B | 5.56 | 2.88 | 2.88 | 1.06 | 0.94 | —
1 kB | 5.38 | 3.62 | 2.66 | 1.03 | 0.81 | 0.69
2 kB | 5.30 | 3.45 | 2.55 | 1.01 | 0.73 | 0.61
4 kB | 5.24 | 3.41 | 2.53 | 1.01 | 0.70 | 0.54
8 kB | 5.24 | 3.36 | 2.42 | 1.01 | 0.69 | 0.52
16 kB | 5.22 | 3.36 | 2.40 | 1.01 | 0.69 | 0.52
32 kB | 5.23 | 3.34 | 2.40 | 1.01 | 0.69 | 0.52
64 kB | 5.22 | 3.34 | 2.40 | 1.01 | 0.69 | 0.52

We present the results for Jaccard index computations in Table 4. Contrary to
straight population counts, the Jaccard-index AVX2 Muła remains faster than
the popcnt-based function even for small blocks (256 B). AVX2 HS provides the
best speed, requiring only 1.15 cycles to calculate the Jaccard similarity
between each pair of 64-bit inputs. This is more than twice as fast
($2.4\times$) as the popcnt-based function. Since the population count is done
for each input of the pair, the speed of the similarity is only slightly
greater than the speed of calculating the two population counts individually.
That is, using AVX2 for both the Boolean operation and both population counts
gives us the Boolean operation almost for free.

Table 4: Number of cycles per pair of 64-bit input words required to compute the Jaccard index of arrays of various sizes. array size | popcnt | AVX2 Muła | AVX2 HS
---|---|---|---
256 B | 3.00 | 2.50 | —
512 B | 2.88 | 2.00 | —
1 kB | 2.94 | 2.00 | 1.53
2 kB | 2.83 | 1.84 | 1.33
4 kB | 2.80 | 1.76 | 1.22
8 kB | 2.78 | 1.75 | 1.16
16 kB | 2.77 | 1.75 | 1.15
32 kB | 2.76 | 1.75 | 1.15
64 kB | 2.76 | 1.74 | 1.15

##  6 Conclusion

On recent Intel processors, the fastest approach to compute the population
count on moderately large arrays (e.g., 4 kB) relies on a vectorized version
of the Harley-Seal function. It is twice as fast as functions based on the
dedicated instruction (popcnt). For the computation of similarity functions
between two bitsets, a vectorized approach based on the Harley-Seal function
is more than twice as fast ($2.4\times$) as an optimized approach based on the
popcnt instruction.

Future work should consider other computer architectures. Building on our work
[25], Sun and del Mundo tested various population-count functions on an Nvidia
GPU and found that its popcount intrinsic gave the best results [27].

We are grateful to S. Utcke for pointing out a typographical error in one of
the figures.

## Funding

This work was supported by the Natural Sciences and Engineering Research
Council of Canada [261437].

## References

  * [1] Lemire, D., Ssi-Yan-Kai, G., and Kaser, O. (2016) Consistently faster and smaller compressed bitmaps with roaring.  Software: Practice and Experience, 46, 1547–1569.
  * [2] Martí, V. (2015).  Counting objects.  GitHub Engineering Blog, http://githubengineering.com/counting-objects/ [last checked December 2016].
  * [3] Lemire, D., Boytsov, L., and Kurz, N. (2016) SIMD compression and the intersection of sorted integers.  Software: Practice and Experience, 46, 723–749.
  * [4] Grand, A. (2015).  Frame of Reference and Roaring Bitmaps.  https://www.elastic.co/blog/frame-of-reference-and-roaring-bitmaps [last checked December 2016].
  * [5] Hilewitz, Y., Shi, Z. J., and Lee, R. B. (2004) Comparing fast implementations of bit permutation instructions.  Signals, Systems and Computers, 2004. Conference Record of the Thirty-Eighth Asilomar Conference on, Washington, DC, USA, Nov, pp. 1856–1863. IEEE Computer Society.
  * [6] Suciu, A., Cobarzan, P., and Marton, K. (2011) The never ending problem of counting bits efficiently.  2011 RoEduNet International Conference 10th Edition, Washington, DC, USA, Jun, pp. 1–4. IEEE Computer Society.
  * [7] Stefanov, E. and Shi, E. (2012) FastPRP: Fast Pseudo-Random Permutations for Small Domains. Technical Report 338. IACR Cryptology ePrint Archive.
  * [8] Manku, G. S., Jain, A., and Das Sarma, A. (2007) Detecting near-duplicates for web crawling.  Proceedings of the 16th International Conference on World Wide Web, New York, NY, USA, May WWW ’07, pp. 141–150. ACM.
  * [9] Prokopenko, D., Hecker, J., Silverman, E. K., Pagano, M., Nöthen, M. M., Dina, C., Lange, C., and Fier, H. L. (2016) Utilizing the Jaccard index to reveal population stratification in sequencing data: a simulation study and an application to the 1000 Genomes Project.  Bioinformatics, 32, 1366–1372.
  * [10] Lacour, A., Schüller, V., Drichel, D., Herold, C., Jessen, F., Leber, M., Maier, W., Noethen, M. M., Ramirez, A., Vaitsiakhovich, T., and Becker, T. (2015) Novel genetic matching methods for handling population stratification in genome-wide association studies.  BMC Bioinformatics, 16, 84.
  * [11] Li, Y., Patel, J. M., and Terrell, A. (2012) WHAM: a high-throughput sequence alignment method.  ACM T. Database Syst., 37, 28.
  * [12] Dambros, C. S., Cáceres, N. C., Magnus, L., and Gotelli, N. J. (2015) Effects of neutrality, geometric constraints, climate, and habitat quality on species richness and composition of Atlantic forest small-mammals.  Global Ecol. Biogeogr., 24, 1084–1093.
  * [13] Zhang, B., Vogt, M., Maggiora, G. M., and Bajorath, J. (2015) Design of chemical space networks using a Tanimoto similarity variant based upon maximum common substructures.  J. Comput. Aided Mol., 29, 937–950.
  * [14] Gueron, S. and Krasnov, V. (2016) Fast quicksort implementation using AVX instructions.  The Computer Journal, 59, 83–90.
  * [15] Hill, S. (2006) Design of a reusable 1GHz, superscalar ARM processor.  2006 IEEE Hot Chips 18 Symposium (HCS), Washington, DC, USA, Aug, pp. 1–18. IEEE Computer Society.
  * [16] Fog, A. (2016) Instruction tables: Lists of instruction latencies, throughputs and micro-operation breakdowns for Intel, AMD and VIA CPUs. Technical report. Copenhagen University College of Engineering, Copenhagen, Denmark.  http://www.agner.org/optimize/instruction_tables.pdf [last checked December 2016].
  * [17] Carruth, C. (2015).  Implement a faster vector population count based on the PSHUFB.  https://reviews.llvm.org/rL238636 [last checked December 2016].
  * [18] Knuth, D. E. (2009) Bitwise Tricks & Techniques, The Art of Computer Programming, 4. Addison-Wesley, Boston, Massachusetts.
  * [19] Wilkes, M. V., Wheeler, D. J., and Gill, S. (1957) The Preparation of Programs for an Electronic Digital Computer, second edition. Addison-Wesley Publishing, Boston, USA.
  * [20] Wegner, P. (1960) A technique for counting ones in a binary computer.  Commun. ACM, 3, 322–.
  * [21] Rohou, E., Swamy, B. N., and Seznec, A. (2015) Branch prediction and the performance of interpreters: Don’t trust folklore.  Proc. 13th IEEE/ACM International Symposium on Code Generation and Optimization, Washington, DC, USA, February CGO ’15, pp. 103–114. IEEE Computer Society.
  * [22] Lauradoux, C. (2009).  Hamming weight.  (It was posted on the author’s website, but is apparently no longer available online.).
  * [23] Warren, H. S., Jr. (2007) The quest for an accelerated population count.  In Wilson, G. and Oram, A. (eds.), Beautiful Code: Leading Programmers Explain How They Think, chapter 10, pp. 147–160. O’Reilly Media, Sebastopol, California.
  * [24] Knuth, D. E. (2011) Combinatorial Algorithms, Part 1, The Art of Computer Programming, 4A. Addison-Wesley, Boston, Massachusetts.
  * [25] Muła, W.  SSSE3: fast popcount.  http://0x80.pl/articles/sse-popcount.html [last checked December 2016].
  * [26] Morancho, E. (2014) A hybrid implementation of Hamming weight.  2014 22nd Euromicro International Conference on Parallel, Distributed, and Network-Based Processing, Sankt Augustin, Germany, Feb, pp. 84–92. Euromicro.
  * [27] Sun, C. and del Mundo, C. C. (2016) Revisiting POPCOUNT operations in CPUs/GPUs.  ACM Student Research Competition Posters at the ACM/IEEE International Conference for High Performance Computing, Networking, Storage, and Analysis.
  * [28] Intel Corporation (2016).  Intel Architecture Instruction Set Extensions Programming Reference.  https://intel.ly/2hurf64 [last checked December 2016].

##  Appendix A Population Counts in AVX-512

Though Intel processors through the current Kaby Lake generation do not yet
support the AVX-512 instruction set, it is straight-forward to generalize our
vectorized algorithms to 512-bit vectors. However, even beyond the increase in
vector width, it should be possible to implement the carry-save adder more
efficiently with AVX-512, which also adds the vpternlogd instruction.
Available through C intrinsics as mm512ternarylogicepi64, this instruction
allows us to compute arbitrary three-input binary functions in a single
operation. Utilizing this, we can replace the 5 logical instructions we needed
for AVX2 with just two instructions. The vpternlogd instruction relies on an
integer parameter $i$ that serves as a look-up table. Given the input bits
$x,y,z$, the value given that the $(x+2y+4z)^{\mathrm{th}}$ bit of the
parameter $i$ is returned. For example, to compute XOR of the inputs, the $i$
parameter needs to have a 1-bit at indexes 1, 2, 4 and 7 (i.e.,
$i=2^{1}+2^{2}+2^{4}+2^{7}=150$ or 0x96 in hexadecimal). Similarly, to compute
the most significant bit of the carry-save adder, the $i$ parameter needs to
have a 1-bit at indexes 3, 5, 6 and 7 (0xe8 in hexadecimal). Fig. 14 presents
a C function implementing a carry-save adder (CSA) using AVX-512 intrinsics.

⬇ void CSA(__m512i* h, __m512i* l, __m512i a, __m512i b, __m512i c) { *l =
_mm512_ternarylogic_epi32(c, b, a, 0x96); *h = _mm512_ternarylogic_epi32(c, b,
a, 0xe8); }
---
Figure 14: A C function using AVX-512 intrinsics implementing a bitwise
parallel carry-save adder (CSA).

Intel has also announced that future processors might support the AVX-512
vpopcnt instruction [28] which computes the population count of each word in a
vector. To our knowledge, no available processor supports vpopcnt at this
time.