These are not the k-mers you are looking for: efficient online k-mer counting
using a probabilistic data structure

Qingpeng Zhang1, Jason Pell1, Rosangela Canino-Koning1, Adina Chuang Howe2,3,
C. Titus Brown1,2∗

1 Department of Computer Science and Engineering, Michigan State University,
East Lansing, MI, USA
2 Department of Microbiology and Molecular Genetics, Michigan State
University, East Lansing, MI, USA
3 Department of Plant, Soil, and Microbial Sciences, Michigan State
University, East Lansing, MI, USA
$\ast$ E-mail<EMAIL_ADDRESS>
## Abstract

K-mer abundance analysis is widely used for many purposes in nucleotide
sequence analysis, including data preprocessing for de novo assembly, repeat
detection, and sequencing coverage estimation. We present the khmer software
package for fast and memory efficient online counting of k-mers in sequencing
data sets. Unlike previous methods based on data structures such as hash
tables, suffix arrays, and trie structures, khmer relies entirely on a simple
probabilistic data structure, a Count-Min Sketch. The Count-Min Sketch permits
online updating and retrieval of k-mer counts in memory which is necessary to
support online k-mer analysis algorithms. On sparse data sets this data
structure is considerably more memory efficient than any exact data structure.
In exchange, the use of a Count-Min Sketch introduces a systematic overcount
for k-mers; moreover, only the counts, and not the k-mers, are stored. Here we
analyze the speed, the memory usage, and the miscount rate of khmer for
generating k-mer frequency distributions and retrieving k-mer counts for
individual k-mers. We also compare the performance of khmer to several other
k-mer counting packages, including Tallymer, Jellyfish, BFCounter, DSK, KMC,
Turtle and KAnalyze. Finally, we examine the effectiveness of profiling
sequencing error, k-mer abundance trimming, and digital normalization of reads
in the context of high khmer false positive rates. khmer is implemented in C++
wrapped in a Python interface, offers a tested and robust API, and is freely
available under the BSD license at github.com/ged-lab/khmer.

## Introduction

The goal of k-mer counting is to determine the number of occurrences for each
fixed-length word of length k in a DNA data set [1]. Efficient k-mer counting
plays an important role in many bioinformatics approaches, including data
preprocessing for de novo assembly, repeat detection, and sequencing coverage
estimation [2].

Short-read shotgun sequencing data is both relatively sparse in k-mers and
contains many erroneous k-mers. For typical values of k such as 32 these data
sets are sparse, as only a small fraction of the total possible number of
k-mers ($4^{32}$) are actually present in any genome or read data sets derived
from the genome. The high error rate (e.g. Illumina has a 0.1-1% per-base
error rate [3]) generates many unique k-mers. As the total number of generated
reads increases, the total number of errors grows with it linearly. This leads
to data sets where the erroneous k-mers vastly outnumber the true k-mers [4].
Tracking and counting the resulting large number of k-mers, most of which are
erroneous, has become an unavoidable and challenging task in sequence analysis
[5].

A variety of k-mer counting approaches, and standalone software packages
implementing them, have emerged in recent years; this includes Tallymer,
Jellyfish, BFCounter, DSK, KMC, Turtle and KAnalyze [2, 1, 6, 7, 8, 9, 10].

These approaches and implementations each offer different algorithmic trade-
offs and enable a non-overlapping set of functionality. Tallymer uses a suffix
tree to store k-mer counts in memory and on disk [2]. Jellyfish stores k-mer
counts in in-memory hash tables, and makes use of disk storage to scale to
larger data sets [1]. BFCounter uses a Bloom filter as a pre-filter to avoid
counting unique k-mers, and is the first published probabilistic approach to
k-mer counting [6]. DSK adopts an approach to k-mer counting that enables
time- and memory-efficient k-mer counting with an explicit trade-off between
disk and memory usage [7]. KMC and KAnalyze rely primarily on fast and
inexpensive disk access to count k-mers in low memory [8, 10]. Turtle provides
several different containers that offer different false positive and false
negative tradeoffs when counting k-mers [9].

Our motivation for exploring efficient k-mer counting comes from our work with
metagenomic data, where we routinely encounter data sets that contain
$300\times 10^{9}$ bases of DNA and over 50 billion distinct k-mers [11]. To
efficiently filter, partition, and assemble these data, we need to store
counts for each of these k-mers in main memory, and query and update them in
realtime — a set of functionality not readily offered by current packages.
Moreover, we wish to enable the use of cloud and desktop computers, which may
have poor I/O performance or limited memory. These needs have dictated our
exploration of efficient in-memory k-mer counting techniques.

Below, we describe an implementation of a simple probabilistic data structure
for k-mer counting. This data structure is based on a Count-Min Sketch [12], a
generalized probabilistic data structure for storing the frequency
distributions of distinct elements. Our implementation extends an earlier
implementation of a Bloom filter [13], which has been previously used in
bioinformatics applications, such as sequence matching [14], k-mer counting
[6], and de Bruijn graph storage and traversal [15, 16]. Many other variations
of Bloom filters have been proposed [17], including counting Bloom filters
[18], multistage filters [19], and spectral Bloom filters [20], which are
related to the Count-Min Sketch and our khmer implementation.

Probabilistic approaches can be particularly memory efficient for certain
problems, with memory usage significantly lower than any exact data structure
[15]. However, their use introduces set membership or counting false
positives, which have effects that must be analyzed in the context of specific
problems. Moreover, unlike existing techniques, the Count-Min Sketch stores
only counts; k-mers must be retrieved from the original data set. In exchange,
the low memory footprint enabled by this probabilistic approach enables online
updating and retrieval of k-mer counts entirely in memory, which in turn
supports streaming applications such as digital normalization [21].

We use the Amazon cloud to compare time, memory, and disk usage of our k-mer
counting implementation with that of other k-mer counting software packages,
for two problems. First, we generate a k-mer abundance distribution for large
data sets; and second, we query many individual k-mer counts at random from a
previously constructed k-mer count database. We show that khmer is competitive
in speed, memory, and disk usage for these problems. We also analyze the
effects of counting error on calculations of the k-mer count in sequencing
data sets, and in particular on metagenomic data sets. Finally, we discuss
khmer’s miscount performance in the context of two specific applications: low-
abundance k-mer trimming of reads, and digital normalization.

The khmer software [22] is implemented in C++ in a Python wrapper, enabling
flexible use and reuse by users with a wide range of computational expertise.
The software package is freely available for academic and commercial use and
redistribution under the BSD license at github.com/ged-lab/khmer/. khmer comes
with substantial documentation and many tutorials, and contains extensive unit
tests. Moreover, we have built several applications on top of khmer, including
memory-efficient de Bruijn graph partitioning [15] and lossy compression of
short-read data sets for assembly [21].

## Results

### Implementing a Count-Min Sketch for k-mers

The two basic operations supported by khmer are increment_count(kmer)  and c =
get_count(kmer).  Both operate on the data structure in memory, such that
neither incrementing a count nor retrieving a count involves disk access.

The implementation details are similar to those of the Bloom filter in [15],
but with the use of 8 bit counters instead of 1 bit counters. Briefly, Z hash
tables are allocated, each with a different size of approximately H bytes
($H_{1},H_{2},...,H_{Z}$); the sum of these hash table sizes must fit within
available main memory. To increment the count for a particular k-mer, a single
hash is computed for the k-mer, and the modulus of that hash with each hash
table’s size H gives the location for each hash table; the associated count in
each hash table is then incremented by 1. We use different sizes for each hash
table so as to vary the hash function. Even if two k-mers have the same
modulus in one hash table (a collision), they are unlikely to collide in the
other hash tables. To retrieve the count for a k-mer, the same hash is
computed and the minimum count across all hash tables is computed. While
different in implementation detail from the standard Bloom filter, which uses
a single hash table with many hash functions, the performance details are
identical [15]. One particularly important feature of the Count-Min Sketch is
that the counting error is one-sided [12]. Because counts are only
incremented, collisions result in inflated miscounts; if there is no collision
for a particular k-mer, the count is correct.

An additional benefit of the Count-Min Sketch is that it is extremely easy to
implement correctly, needing only about 3 dozen lines of C++ code for a simple
threadsafe implementation. (We have described how khmer scales with multiple
threads in [23].)

To determine the expected false positive rate — the average frequency with
which a given k-mer count will be incorrect when retrieved — we can look at
the hash table load. Suppose N distinct k-mers have been counted using Z hash
tables, each with size H. The probability that no collisions happened in a
specific entry in one hash table is $(1-1/H)^{N}$, or approximately
$e^{-N/H}$. The individual collision rate in one hash table is then $\approx
1-e^{-N/H}$. The total collision rate, which is the probability that a
collision occurred in each entry where a k-mer maps across all Z hash tables,
is $\approx(1-e^{-N/H})^{Z}$, which is also the expected false positive rate.

While the false positive rate can easily be calculated from the hash table
load, the average miscount — the degree to which the measured count differs
from the true count — depends on the k-mer frequency distribution, which must
be determined empirically. We analyze the effects of this below.

### Choosing number and size of hash tables used for k-mer counting

The false positive rate depends on the number of distinct k-mers $N$, the
number of hash tables $Z$, and the size of the hash tables $H$:
$f\approx(1-e^{-N/H})^{Z}$, with an associated memory usage of $M=HZ$. We face
two common scenarios: one in which we have a fixed number of k-mers $N$ and
fixed memory $M$ and we want to calculate the optimal number of hash tables
$Z$; and one in which we have a desired maximum false positive rate $f$ and a
fixed number of k-mers $N$, and we want to calculate the minimum memory usage
required to achieve $f$.

For fixed memory $M$ and number of distinct k-mers $N$, the optimal number of
hash tables can be found by minimizing $f$; taking the derivative, $df/dZ$,
with $f\approx\exp(Z\log(1-e^{-ZN/M}))$ and solving for 0, we find that $f$ is
minimized when $Z=\log(2)*(M/N)$ (see [24] for details).

Given a desired false positive rate $f$ and a fixed number of k-mers $N$, the
optimal memory usage can be calculated as follows. First, the optimal number
of hash tables is determined by the expected false positive rate alone:
$Z=\log_{0.5}f$. Using this $Z$, the minimum average hash table size $H$
necessary to achieve $f$ can be calculated as $H=(\log_{0.6185}(f)\times N)/Z$
(see [24] for details).

A remaining problem is that the number of distinct k-mers $N$ is typically not
known. However, memory- and time-efficient algorithms for calculating $N$ do
exist and we plan to implement this in khmer in the future [25].

### khmer efficiently calculates k-mer abundance histograms

We measured time and memory required to calculate k-mer abundance histograms
in five soil metagenomic read data sets using khmer, Tallymer, Jellyfish, DSK,
KMC, Turtle, and KAnalyze (Table 1; Figures 1 and 2). We chose to benchmark
abundance histograms because this functionality is common to all the software
packages, and is a common analysis approach for determining assembly
parameters [26]. We applied each package to increasingly large subsets of a
50m read soil metagenome data set [11]. For the BFCounter, KMC, Turtle and
KAnalyze packages, which do not generate k-mer abundance distribution
directly, we output the frequency of each k-mer to a file but do no further
analysis.

Figure 1 shows that the time usage of the khmer approach is comparable to DSK
and BFCounter, and, as expected, increases linearly with data set size.
Tallymer is the slowest of the four tools in this testing, while KMC, Turtle,
and Jellyfish are the fastest.

From Figure 2, we see that the memory usage of Jellyfish, Tallymer, BFCounter,
and Turtle increases linearly with data set size. Tallymer uses more memory
than Jellyfish generally, while BFCounter and Turtle have considerably lower
memory usage. DSK, KMC, and KAnalyze use constant memory across the data sets,
but at the cost of more limited functionality (discussed below).

The memory usage of khmer also increases linearly with data set size as long
as we hold the false positive rate constant. However, the memory usage of
khmer varies substantially with the desired false positive rate: we can
decrease the memory usage by increasing the false positive rate as shown in
Figure 2. We also see that with a low false positive of 1%, the memory usage
is competitive with Tallymer and Jellyfish; with a higher 5% false positive
rate, the memory usage is lower than all but the disk-based DSK; with an false
positive rate as high as 20%, the memory usage is further lower, close to DSK,
KAnalyze, and KMC.

We also measured disk usage during counting. Figure 3 shows that the disk
usage also increases linearly with the number of k-mers in the data set. For a
high-diversity metagenomic data set of 5 GB, the disk usage of both Jellyfish
and Tallymer is around 30 GB. khmer counts k-mers entirely in working memory
and does not rely on any on-disk storage to store or retrieve k-mer counts,
although for practicality the hash tables can be saved for later reuse; the
uncompressed disk usage for khmer in Figure 3 is the same as its memory. At
the expense of more time, khmer supports saving and loading gzip-compressed
hash tables, which are competitive in size to DSK’s on-disk database (Figure
3, dashed line).

### khmer accesses k-mer counts efficiently

We measured the time it took to access 9.7m 22-mers across five different data
sets after the initial databases had been built (Figure 4). Note that
Tallymer, Jellyfish, and khmer all support random access to k-mer counts,
while BFCounter, DSK, KMC, Turtle and KAnalyze do not. Here, khmer performed
well, dramatically outperforming Jellyfish and Tallymer. In all three cases,
system time dominated the overall time required to retrieve k-mers, suggesting
that the primary reason for the increase in retrieval time was due to the
increased size of the database on the disk (data not shown). In particular,
khmer is independent of the size of the database in retrieval time once the
hash tables are loaded into memory.

### The measured counting error is low on short-read data

Due to the use of Count-Min Sketch and its lack of collision tracking, khmer
will report some incorrect counts for k-mers; these counts are always higher
than the true counts, up to the bound of 255 (a limit imposed by our use of
8-bit counters). The frequency with which incorrect counts are reported can be
estimated from the hash table load. However, the expected miscount — the
difference between the true k-mer frequency and the reported k-mer frequency —
cannot be calculated without knowing the distribution of k-mer abundances; in
general, the average miscount will be small if the data is left-skewed. As
noted by Melsted and Pritchard, a large number of k-mers in short-read data
are low-abundance, leading to precisely the skew that would yield low
miscounts [6]. Here we use both real and simulated data sets to evaluate the
counting performance in practice.

Figure 5 shows the relationship between average miscount and counting false
positive rate for five different test data sets with similar numbers of
distinct k-mers: one metagenome data set; a simulated set of random k-mers; a
simulated set of reads, chosen with 3x coverage and 1% error; a simulated set
of reads (3x) with no error; and a set of E. coli reads (Table 2). Even when
the counting false positive rate is as high as 0.9 — where 90% of k-mers have
an incorrect count — the average miscount is still below 4.

We separately analyzed the average percentage miscount between true and false
k-mers; e.g. an miscount of 4 for a k-mer whose true count is 1 would be 400%.
Figure 6 shows the relationship between average miscount and counting false
positive rate for the same five data sets as in Figure 5. For a false positive
rate of 0.1 (10% of k-mer counts are incorrect), the average percentage
miscount is less than 10% for all five data sets; this will of course
generally be true, because the average miscount is bounded by the product of
the false positive rate with k-mer abundance.

We see here that for a fixed false positive rate, the simulated reads without
error have the highest average miscount, and the randomly generated k-mers
have the lowest average miscount. This is because these two abundance
distributions have the least and most left-skew, respectively: the simulated
reads without error have no abundance-1 k-mers, while the randomly generated
k-mers are entirely low abundance.

### Sequencing error profiles can be measured with k-mer abundance profiles

One specific use for khmer is detecting random sequencing errors by looking at
the k-mer abundance distribution within reads [27]. This approach, known also
as “k-mer spectral analysis”, was first proposed in by [28] and further
developed in [29]. The essential idea is that low-abundance k-mers contained
in a high-coverage data set typically represent random sequencing errors.

A variety of read trimming and error correcting tools use k-mer counting to
reduce the error content of the read data set, independent of quality scores
or reference genomes [30]. This is an application where the counting error of
the Count-Min Sketch approach used by khmer may be particularly tolerable: it
will never falsely call a high-abundance k-mer as low-abundance because khmer
never underestimates counts.

In Figure 7, we use khmer to examine the sequencing error pattern of a 5m-read
subset of an Illumina reads data set from single-colony sequencing of E. coli
[31]. The high rate of occurrence of unique k-mers close to the 3’ end of
reads is due to the increased sequencing error rate at the 3’ end of reads.

### khmer can be applied iteratively to read trimming

We next evaluated the effect of false-positive induced miscounts on read
trimming, in which reads are truncated at the first low-abundance k-mer.
Because the Count-Min Sketch never undercounts k-mers, reads will never be
erroneously trimmed at truly high-abundance k-mers; however, reads may not be
trimmed correctly when miscounts inflate the count of low-abundance k-mers. In
cases where many errors remain, read trimming can potentially be applied
multiple times, with each round reducing the total number of k-mers and hence
resulting in lower false positive rates for the same memory usage.

We performed six iterations of unique k-mer trimming on 5 million Illumina
reads from sequencing of E. coli, with memory usage less than 30 MB. For each
iteration we measured empirical false positive rate compared with number of
bases trimmed as well as the total number of k-mers (Table 3). In the first
round, the estimated false positive rate was 80.0%, and 13.5% of the total
bases were removed by trimming reads at low-abundance k-mers; the second
iteration had a false positive rate of 37.7%, and removed only 1.5% additional
data; and by the fourth iteration the false positive rate was down to 23.2%
with 0.0% of the data removed.

The elimination of so many unique k-mers (column 5) in the first pass was
unexpected: the high false positive rate should have resulted in fewer k-mers
being identified as unique, were the erroneous k-mers independent of each
other. Upon examination, we realized that in Illumina data erroneous k-mers
typically come from substitution errors that yield runs of up to k erroneous
k-mers in a row [30]. When trimming reads with high false positive rates,
these runs are typically trimmed after the first few unique k-mers, leaving
unique k-mers at the 3’ end. Because of this we hypothesized that high-FP rate
trimming would result in the retention of many unique k-mers at the 3’ end of
the read, and this was confirmed upon measurement (Table 3, column 6, pass 1
vs pass 2).

In comparison to quality-based trimming software such as seqtk and FASTX,
trimming at unique k-mers performed very well: in this data set, all unique
k-mers represent errors, and even with an initial false positive rate of 80%,
khmer outperformed all but the most stringent seqtk run (Table 3). With a
lower false positive rate or multiple passes, khmer eliminates more erroneous
k-mers than seqtk or FASTX. The tradeoff here is in memory usage: for larger
data sets, seqtk and FASTX will consume the same amount of memory as on
smaller data sets, while khmer’s memory usage will need to grow with the data
set size.

### Using khmer for digital normalization, a streaming algorithm

Digital normalization is a lossy compression algorithm that discards short
reads based on saturating coverage of a de Bruijn graph [21]. While several
non-streaming implementations exist, including Trinity’s in silico
normalization [32, 33], digital normalization can be efficiently implemented
as a streaming algorithm. In the streaming implementation, if a read is not
kept, it is not loaded into the Count-Min Sketch structure, and the false
positive rate does not increase. For high coverage data sets, the digital
normalization algorithm is sublinear in memory because it does not collect the
majority of k-mers in those data sets [21]. This has the advantage of enabling
low-memory preprocessing of both high-coverage genomic data sets, as well as
mRNAseq or metagenomic data sets with high-coverage components [21, 11].

While digital normalization is already implemented inside khmer, previous work
did not explore the lower bound on memory usage for effective digital
normalization. In particular, the effects of high false positive rates have
not been examined in any prior work.

We applied digital normalization to the E. coli data set used above, and chose
seven different Count-Min Sketch sizes to yield seven different false positive
rates 4. The data set was normalized to a k-mer coverage of 20 and the
resulting data were evaluated for retention of true and erroneous k-mers, as
in [21] (Table 4). The results show that digital normalization retains the
same set of underlying “true” k-mers until the highest false positive rate of
100% (Table 4, column 5), while discarding only about 2% additional reads
(Table 4, column 6).

To evaluate the effect of digital normalization with high false positive rates
on actual genome assembly, we next performed normalization to a coverage of 20
with the same range of false positive rates as above. We then assembled this
data with Velvet [34] and compared the resulting assemblies to the known E.
coli MG1655 genome using QUAST (Table 5). To our surprise, we found that even
after executing digital normalization with a false positive rate of 83.2%, a
nearly complete assembly was generated. No progressive increase in
misassemblies (measured against the real genome with QUAST) was seen across
the different false positive rates (data not shown). This suggests that below
83.2% FP rate, the false positive rate of digital normalization has little to
no effect on assembly quality with Velvet. (Note that the Velvet assembler
itself used considerably more memory than digital normalization.)

While these results are specific to Velvet and the coverage parameters used in
digital normalization, they do suggest that no significant information loss
occurs due to false positive rates below 80%. Further evaluation of assembly
quality in response to different normalization parameters and assemblers is
beyond the scope of of this paper.

## Discussion

### khmer enables fast, memory-efficient online counting

khmer enables memory- and time-efficient online counting (Figures 1, 2, and
4). This is particularly important for the streaming approaches to data
analysis needed as data set sizes increase. Because query and updating of
k-mer counts can be done directly as data is being loaded, with no need for
disk access or an indexing step, khmer can also perform well in situations
with poor disk I/O performance. (Note that BFCounter also supports online
k-mer counting [6].)

### khmer is a generally useful k-mer counting approach

In addition to online counting, khmer offers a general range of useful
performance tradeoffs for disk I/O, time and memory. From the performance
comparison between khmer and other k-mer counting packages in calculating
k-mer abundance distributions, khmer is comparable with existing packages. In
time, khmer performs competitively with DSK and BFCounter (Figure 1); khmer
also provides a way to systematically trade memory for miscounts across a wide
range of parameters (Figure 2). khmer’s uncompressed disk storage is
competitive with Jellyfish, and, in situations where disk space is at a
premium, khmer can take advantage of gzip compression to provide storage
similar to that of DSK (Figure 3, purple line with boxes).

KMC, DSK, and KAnalyze perform especially well in memory usage for calculating
the abundance distribution of k-mers. However, in exchange for this
efficiency, retrieving specific k-mer counts at random is likely to be quite
slow, as DSK is optimized for iterating across partition sets of k-mers rather
than randomly accessing k-mer counts.

For retrieving the counts of individual k-mers, khmer is significantly faster
than both Tallymer and Jellyfish. This is not surprising, since this was a
primary motivation for the development of khmer.

### khmer memory usage is fixed and low

The memory usage of the basic Count-Min Sketch approach is fixed: khmer’s
memory usage does not increase as data is loaded. While this means that khmer
will never crash due to memory limitations, and all operations can be
performed in main memory without recourse to disk storage, the false positive
rate may grow too high. Therefore the memory size must be chosen in light of
the false positive rate and miscount acceptable for a given application. In
practice, we recommend choosing the maximum available memory, because the
false positive rate decreases with increasing memory and there are no negative
effects to minimizing the false positive rate.

For any given data set, the size and number of hash tables will determine the
accuracy of k-mer counting with khmer. Thus, the user can control the memory
usage based on the desired level of accuracy (Figure 2). The time usage for
the first step of k-mer counting, consuming the reads, depends on the total
amount of data, since we must traverse every k-mer in every read. The second
step, k-mer retrieval, is algorithmically constant for fixed k; however, for
practicality, the hash tables are usually saved to and loaded from disk,
meaning that k-mer retrieval time depends directly on the size of the database
being queried.

The memory usage of khmer is particularly low for sparse data sets, especially
since only main memory is used and no disk space is necessary beyond that
required for the read data sets. This is no surprise: the information
theoretic comparison in [15] shows that, for sparse sequencing data sets,
Bloom filters require considerably less memory than any possible exact
information storage for a wide range of false positive rates and data set
sparseness.

In our implementation we use 1 byte to store the count of each k-mer in the
data structure. Thus the maximum count for a k-mer will be 255. In cases where
tracking bigger counts is required, khmer also provides an option to use an
STL map data structure to store counts above 255, with the trade-off of
significantly higher memory usage. In the future, we may extend khmer to
counters of arbitrary bit sizes.

### False positive rates in k-mer counting are low and predictable

The Count-Min Sketch is a probabilistic data structure with a one-sided error
that results in random overestimates of k-mer frequency, but does not generate
underestimates.

In the Count-Min Sketch, the total memory usage is fixed; the memory usage,
the hash functions, and the total number of distinct objects counted all
influence the accuracy of the count. While the probability of an inaccurate
count can easily be estimated based on the hash table load, the miscount size
is dependent on details of the frequency distribution of k-mers [12].

More specifically, in the analysis of the Count-Min Sketch, the difference
between the incorrect count and actual count is related to the total number of
k-mers in a data set and the size of each hash table [12]. Further study has
shown that the behavior of Count-Min Sketch depends on specific
characteristics of the data set under consideration, especially left-skewness
[35, 36]. These probabilistic properties suit short reads from next generation
sequencing data sets: the miscounts are low because of the highly left-skewed
abundance distribution of k-mers in these data sets.

Figures 5 and 6 demonstrate these properties well. We see more correct
counting for error-prone reads from a genome than for error-free reads from a
genome, with a normal distribution of k-mer abundance. Thus, this counting
approach is especially suitable for high diversity data sets, such as
metagenomic data, in which a larger proportion of k-mers are low abundance or
unique due to sequencing errors.

### Real-world applications for khmer

For many applications, an approximate k-mer count is sufficient. For example,
when eliminating reads with low abundance k-mers, we can tolerate a certain
number of low-frequency k-mers remaining in the resulting data set falsely. If
RAM-limited we can do the filtering iteratively so that at each step we are
making more effective use of the available memory.

In practice, we have found that a false positive rate of between 1% and 10%
offers acceptable miscount performance for a wide range of tasks, including
error profiling, digital normalization and low-abundance read-trimming.
Somewhat surprisingly, false positive rates of up to 80% can still be used for
both read trimming and digital normalization in memory-limited circumstances,
although multiple passes across the data may be needed.

For many applications, the fact that khmer does not break an imposed memory
bound is extremely useful, since for many data sets — especially metagenomic
data sets — high memory demands constrain analysis [11, 37]. Moreover, because
the false positive rate is straightforward to measure, the user can be warned
that the results should be invalidated when too little memory is used. When
combined with the graceful degradation of performance for both error trimming
and digital normalization, khmer readily enables analysis of extremely large
and diverse data sets [38]. In an experiment to assemble the reads of a soil
metagenomic sample collected from Iowa prairie, the number of reads to
assemble drops from 3.3 million to 2.2 million and the size of the data set
drops from 245GB to 145GB accordingly after digital normalization [11]. 240GB
memory was used in the process. This also shows that khmer works well to
analyze large, real-world metagenomic data sets.

### Conclusion

K-mer counting is widely used in bioinformatics, and as sequencing data set
sizes increase, graceful degradation of data structures in the face of large
amounts of data has become important. This is especially true when the
theoretical and practical effects of the degradation can be predicted (see
e.g. [6, 15, 9]). This is a key property of the Count-Min Sketch approach, and
its implementation in khmer.

The khmer software implementation offers good performance, a robust and well-
tested Python API, and a number of useful and well-documented scripts. While
Jellyfish, DSK, KMC, and Turtle also offer good performance, khmer is
competitive, and, because it provides a Python API for online counting, is
flexible. In memory-limited situations with poor I/O performance, khmer is
particularly useful, because it will not break an imposed memory bound and
does not require disk access to store or retrieve k-mer counts. However, in
exchange for this memory guarantee, counting becomes increasingly incorrect as
less memory is used or as the data set size grows large; in many situations
this may be an acceptable tradeoff.

### Future considerations

Applying khmer to extremely large data sets with many distinct k-mers requires
a large amount of memory: approximately 446 GB of memory is required to
achieve an false positive rate of 1% for $50\times 10^{9}$ k-mers. It is
possible to reduce the required memory by dividing k-mer space into multiple
partitions and counting k-mers separately for each partition. Partitioning
k-mer space into $M$ partitions results in a linear decrease in the number of
k-mers under consideration, thus reducing the occupancy by a constant factor
$M$ and correspondingly reducing the collision rate. Partitioning k-mer space
is a generalization of the systematic prefix filtering approach, where one
might first count all k-mers starting with AA, then AC, then AG, AT, CA, etc.,
which is equivalent to partitioning k-mer space into 16 equal-sized
partitions. These partitions can be calculated independently, either across
multiple machines or iteratively on a single machine, and the results stored
for later comparison or analysis. This is similar to the approach taken by DSK
[7], and could easily be implemented in khmer.

Further optimization of khmer on single machines, e.g. for multi-core
architectures, is unlikely to achieve significantly greater speed. Past a
certain point k-mer counting is fundamentally I/O bound [23].

Perhaps the most interesting future direction for probabilistic k-mer counting
is that taken by Turtle [9], in which several data structures are provided,
each with different tradeoffs, but with a common API. We hope to pursue this
direction in the future by integrating such approaches into khmer.

## Methods

### Code and data set availability

The version of khmer used to generate the results below is available at
http://github.com/ged-lab/khmer.git, tag ’2013-khmer-counting’. Scripts
specific to this paper are available in the paper repository at
https://github.com/ged-lab/2013-khmer-counting. The IPython[39] notebook file
and data analysis to generate the figures are also available in that github
repository. Complete instructions to reproduce all of the results in this
paper are available in the khmer-counting repository; see README.rst.

### Sequence Data

One human gut metagenome reads data set (MH0001) from the MetaHIT
(Metagenomics of the Human Intestinal Tract) project [40] was used. It
contains approximately 59 million reads, each 44bp long; it was trimmed to
remove low quality sequences.

Five soil metagenomics reads data sets with different size were taken from the
GPGC project for benchmark purpose (see Table 1). These reads are from soil in
Iowa region and they are filtered to make sure there are less than 30% Ns in
the read and each read is longer than 30 bp. The exact data sets used for the
paper are available on Amazon S3 and the instructions to acquire these data
sets are available in the paper repository on github.com.

We also generated four short-read data sets to assess the false positive rate
and miscount distribution. One is a subset of a real metagenomics data set
from the MH0001 data set, above. The second consists of randomly generated
reads. The third and fourth contain reads simulated from a random, 1 Mbp long
genome. The third has a substitution error rate of 3%, and the fourth contains
no errors. The four data sets were chosen to contain identical numbers of
distinct 22-mers. The scripts necessary to regenerate these data are available
in the paper repository on github.com.

### Count-Min Sketch implementation

We implemented the Count-Min Sketch data structure, a simple probabilistic
data structure for counting distinct elements [12]. Our implementation uses
$Z$ independent hash tables, each containing a prime number of counters
$H_{i}$. The hashing function for each hash table is fixed, and reversibly
converts each DNA k-mer (for $k\leq 32$) into a 64-bit number to which the
modulus of the hash table size is applied. This provides $Z$ distinct hash
functions.

To increment the count associated with a k-mer, the counter associated with
the hashed k-mer in each of the $N$ hash tables is incremented. To retrieve
the count associated with a k-mer, the minimum count across all $N$ hash
tables is chosen.

In this scheme, collisions are explicitly not handled, so the count associated
with a k-mer may not be accurate. Because collisions only falsely increment
counts, however, the retrieved count for any given k-mer is guaranteed to be
no less than the correct count. Thus the counting error is one-sided.

### Hash function and khmer implementation

The current khmer hash function works only for $k\leq 32$ and converts DNA
strings exactly into 64-bit numbers. However, any hash function would work.
For example, a cyclic hash would enable khmer to count k-mers larger in size
than 32; this would not change the scaling behavior of khmer at all, and is a
planned extension.

By default khmer counts k-mers in DNA, i.e. strandedness is disregarded by
having the hash function choose the lower numerical value for the exact hash
of both a k-mer and its reverse complement. This behavior is configurable via
a compile-time option.

### Comparing with other k-mer counting programs

We generated k-mer abundance distribution from five soil metagenomic reads
data sets of different sizes using khmer, Tallymer, Jellyfish, DSK, BFCounter,
KMC, Turtle and KAnalyze. If the software does not include function to
generate k-mer abundance distribution directly, we output the frequency of
each k-mer in an output file. We fixed k at 22 unless otherwise noted.

#### khmer:

For khmer, we set hash table sizes to fix the false positive rate at either
1%, 5% or 20%, and used 8 threads in loading the data.

We did the khmer random-access k-mer counting benchmark with a custom-written
Python script khmer-count-kmers which loaded the database file and then used
the Python API to query each k-mer individually.

#### Tallymer:

Tallymer is from the genometools package version 1.3.4. For the suffixerator
subroutine we used: -dna -pl -tis -suf -lcp.

We use the mkindex subroutine to generate k-mer abundance distribution, we
used: -mersize 22.

The Tallymer random access k-mer counting benchmark was done using the
’tallymer search’ routine against both strands; see the script tallymer-
search.sh.

#### Jellyfish:

The Jellyfish version used was 1.1.10 and the multithreading option is set to
8 threads.

Jellyfish uses a hash table to store the k-mers and the size of the hash table
can be modified by the user. When the specified hash table fills up, Jellyfish
writes it to the hard disk and initializes a new hash table. Here we use a
similar strategy as in [6] and chose the minimum size of the hash tables for
Jellyfish so that all k-mers were stored in memory.

We ran Jellyfish with the options as below:

jellyfish count -m 22 -c 2 -C for k=22.

The Jellyfish random access k-mer counting benchmark was performed using the
’query’ routine and querying against both strands; see the script jelly-
search.sh.

#### DSK:

We ran DSK with default parameters with -histo option to generate k-mer
abundance distribution. The DSK version used was 1.5031.

#### BFCounter:

The BFcounter version used was 1.0 and the multithreading option is set to 8
threads.

We ran BFCounter count subroutine with the options as below:

BFCounter count -k 22 -t 8 -c 100000. -n option representing the estimated
number of k-mers is adjusted to the different test data sets.

This subroutine produces the actual count of k-mers in input files.

We ran BFCounter dump subroutine with the options as below: BFCounter dump -k
22.

This subroutine can write k-mer occurrences into a tab-separated text file.

#### KMC:

The KMC version used was 0.3. We ran both kmc and kmc_dump subroutines with
default parameters.

#### Turtle:

The Turtle version used was 0.3. We ran scTurtle32 with the multithreading
option set to 8 threads and -n option representing expected number of frequent
k-mers is adjusted to different test data sets.

#### KAnalyze:

The KAnalyze version used was 0.9.3. We ran count subroutine with default
parameters.

## Acknowledgments

We thank Eric McDonald for technical assistance with optimizing the khmer
codebase.

## References

  * 1. Marçais G, Kingsford C (2011) A fast, lock-free approach for efficient parallel counting of occurrences of k-mers.  Bioinformatics 27: 764–770.
  * 2. Kurtz S, Narechania A, Stein JC, Ware D (2008) A new method to compute K-mer frequencies and its application to annotate large repetitive plant genomes.  BMC Genomics 9: 517.
  * 3. Metzker M (2010) Sequencing technologies - the next generation.  Nat Rev Genet 11: 31-46.
  * 4. Conway TC, Bromage AJ (2011) Succinct data structures for assembling large genomes.  Bioinformatics 27: 479-86.
  * 5. Minoche AE, Dohm JC, Himmelbauer H (2011) Evaluation of genomic high-throughput sequencing data generated on illumina hiseq and genome analyzer systems.  Genome Biol 12: R112.
  * 6. Melsted P, Pritchard JK (2011) Efficient counting of k-mers in DNA sequences using a bloom filter.  BMC bioinformatics 12: 333.
  * 7. Rizk G, Lavenier D, Chikhi R (2013) Dsk: k-mer counting with very low memory usage.  Bioinformatics 29: 652-3.
  * 8. Deorowicz S, Debudaj-Grabysz A, Grabowski S (2013) Disk-based k-mer counting on a pc.  BMC Bioinformatics 14: 160.
  * 9. Roy RS, Bhattacharya D, Schliep A (2014) Turtle: Identifying frequent k-mers with cache-efficient algorithms.  Bioinformatics: Advance Access published March 10, 2014 : doi: 10.1093/bioinformatics/btu132.
  * 10. Audano P, Vannberg F (2014) Kanalyze: A fast versatile pipelined k-mer toolkit.  Bioinformatics: Advance Access published March 18, 2014 : doi: 10.1093/bioinformatics/btu152.
  * 11. Howe AC, Jansson JK, Malfatti SA, Tringe SG, Tiedje JM, et al. (2014) Tackling soil diversity with the assembly of large, complex metagenomes.  Proc Natl Acad Sci U S A 111: 4904-9.
  * 12. Cormode G, Muthukrishnan S (2005) An improved data stream summary: the count-min sketch and its applications.  Journal of Algorithms 55: 58–75.
  * 13. Bloom BH (1970) Space/time trade-offs in hash coding with allowable errors.  Commun ACM 13: 422-426.
  * 14. Malde K, O’Sullivan B (2009) Using bloom filters for large scale gene sequence analysis in haskell.  In: Gill A, Swift T, editors, PADL. Springer, volume 5418 of _Lecture Notes in Computer Science_ , pp. 183-194.
  * 15. Pell J, Hintze A, Canino-Koning R, Howe A, Tiedje JM, et al. (2012) Scaling metagenome sequence assembly with probabilistic de bruijn graphs.  Proc Natl Acad Sci U S A 109: 13272-7.
  * 16. Jones DC, Ruzzo WL, Peng X, Katze MG (2012) Compression of next-generation sequencing reads aided by highly efficient de novo assembly.  Nucleic Acids Res 40: e171.
  * 17. Broder AZ, Mitzenmacher M (2003) Survey: Network applications of bloom filters: A survey.  Internet Mathematics 1: 485-509.
  * 18. Fan L, Cao P, Almeida J, Broder AZ (2000) Summary cache: A scalable wide-area web cache sharing protocol.  IEEE/ACM Trans Netw 8: 281–293.
  * 19. Estan C, Varghese G (2002) New directions in traffic measurement and accounting.  In: SIGCOMM. ACM, pp. 323-336.
  * 20. Cohen S, Matias Y (2003) Spectral bloom filters.  In: Halevy AY, Ives ZG, Doan A, editors, SIGMOD Conference. ACM, pp. 241-252.
  * 21. Brown CT, Howe A, Zhang Q, Pyrkosz AB, Brom TH (2012) A reference-free algorithm for computational normalization of shotgun sequencing data.  arXiv : 1203.4802.
  * 22. Crusoe M, Edvenson G, Fish J, Howe A, McDonald E, et al. (2014) The khmer software package: enabling efficient sequence analysis.  URL http://dx.doi.org/10.6084/m9.figshare.979190.
  * 23. McDonald E, Brown CT (2013) Working with big data in bioinformatics.  In: Armstrong T, editor, The Performance of Open Source Applications, lulu.com, chapter 12. p. 151.
  * 24. Broder A, Mitzenmacher M (2004) Network applications of bloom filters: A survey.  Internet mathematics 1: 485–509.
  * 25. Flajolet P, Fusy É, Gandouet O, Meunier F (2008) Hyperloglog: the analysis of a near-optimal cardinality estimation algorithm.  DMTCS Proceedings .
  * 26. Chikhi R, Medvedev P (2014) Informed and automated k-mer size selection for genome assembly.  Bioinformatics 30: 31-7.
  * 27. Medvedev P, Scott E, Kakaradov B, Pevzner P (2011) Error correction of high-throughput sequencing datasets with non-uniform coverage.  Bioinformatics 27: i137-41.
  * 28. Pevzner PA, Tang H, Waterman MS (2001) An eulerian path approach to dna fragment assembly.  Proc Natl Acad Sci U S A 98: 9748-53.
  * 29. Li X, Waterman MS (2003) Estimating the repeat structure and length of dna sequences using l-tuples.  Genome Res 13: 1916-22.
  * 30. Kelley DR, Schatz MC, Salzberg SL (2010) Quake: quality-aware detection and correction of sequencing errors.  Genome Biol 11: R116.
  * 31. Chitsaz H, Yee-Greenbaum J, Tesler G, Lombardo M, Dupont C, et al. (2011) Efficient de novo assembly of single-cell bacterial genomes from short-read data sets.  Nat Biotechnol 29: 915-21.
  * 32. Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, et al. (2013) De novo transcript sequence reconstruction from rna-seq using the trinity platform for reference generation and analysis.  Nat Protoc 8: 1494-512.
  * 33. Brown CT (2012) What does trinity’s in silico normalization do?  URL http://dx.doi.org/10.6084/m9.figshare.98198.
  * 34. Zerbino DR, Birney E (2008) Velvet: algorithms for de novo short read assembly using de bruijn graphs.  Genome Res 18: 821-9.
  * 35. Rusu F, Dobra A (2008) Sketches for size of join estimation.  ACM Transactions on Database Systems 33: 1–46.
  * 36. Cormode G, Muthukrishnan S (2005) Summarizing and mining skewed data streams.  In: Kargupta H, Srivastava J, Kamath C, Goodman A, editors, SDM. SIAM, pp. 44-55.
  * 37. Luo W, Friedman MS, Shedden K, Hankenson KD, Woolf PJ (2009) Gage: generally applicable gene set enrichment for pathway analysis.  BMC Bioinformatics 10: 161.
  * 38. Howe AC, Pell J, Canino-Koning R, Mackelprang R, Tringe S, et al. (2012) Illumina sequencing artifacts revealed by connectivity analysis of metagenomic datasets.  arXiv : 1212.0159.
  * 39. Pérez F, Granger B (2007) Ipython: A system for interactive scientific computing.  Computing in Science Engineering 9: 21-29.
  * 40. Qin J, Li R, Raes J, Arumugam M, Burgdorf KS, et al. (2010) A human gut microbial gene catalogue established by metagenomic sequencing.  Nature 464: 59–65.

## Figure Legends

Figure 1: Comparison of the time it takes for k-mer counting tools to
calculate k-mer abundance histograms, with time (y axis, in seconds) against
data set size (in number of reads, x axis). All programs executed in time
approximately linear with the number of input reads.

Figure 2: Memory usage of k-mer counting tools when calculating k-mer
abundance histograms, with maximum resident program size (y axis, in GB)
plotted against the total number of distinct k-mers in the data set (x axis,
billions of k-mers).

Figure 3: Disk storage usage of different k-mer counting tools to calculate
k-mer abundance histograms in GB (y axis), plotted against the number of
distinct k-mers in the data set (x axis). ∗Note that khmer does not use the
disk during counting or retrieval, although its hash tables can be saved for
reuse.

Figure 4: Time for several k-mer counting tools to retrieve the counts of 9.7m
randomly chosen k-mers (y axis), plotted against the number of distinct k-mers
in the data set being queried (x axis). BFCounter, DSK, Turtle, KAnalyze, and
KMC do not support this functionality.

Figure 5: Relation between average miscount — amount by which the count for
k-mers is incorrect — on the y axis, plotted against false positive rate (x
axis), for five data sets. The five data sets were chosen to have the same
total number of distinct k-mers: one metagenome data set; a set of randomly
generated k-mers; a set of reads, chosen with 3x coverage and 1% error, from a
randomly generated genome; a simulated set of error-free reads (3x) chosen
from a randomly generated genome and a set of E. coli reads.

Figure 6: Relation between percent miscount — amount by which the count for
k-mers is incorrect relative to its true count — on the y axis, plotted
against false positive rate (x axis), for five data sets. The five data sets
are the same as in Figure 5.

Figure 7: Number of unique k-mers (y axis) by starting position within read (x
axis) in an untrimmed E. coli 100-bp Illumina shotgun data set, for k=17 and
k=32. The increasing numbers of unique k-mers are a sign of the increasing
sequencing error towards the 3’ end of reads. Note that there are only 69
starting positions for 32-mers in a 100 base read.

## Tables

Table 1:  Benchmark soil metagenome data sets for k-mer counting performance, taken from [11]. Data set | size of file (GB) | number of reads | number of distinct k-mers | total number of k-mers
---|---|---|---|---
subset 1 | 1.90 | 9,744,399 | 561,178,082 | 630,207,985
subset 2 | 2.17 | 19,488,798 | 1,060,354,144 | 1,259,079,821
subset 3 | 3.14 | 29,233,197 | 1,445,923,389 | 1,771,614,378
subset 4 | 4.05 | 38,977,596 | 1,770,589,216 | 2,227,756,662
entire data set | 5.00 | 48,721,995 | 2,121,474,237 | 2,743,130,683
Table 2:  Data sets used for analyzing miscounts. Data set | Size of data set file | Number of total k-mers | Number of distinct k-mers
---|---|---|---
Real metagenomics reads | 7.01M | 2,917,200 | 1,944,996
Totally random reads with randomly generated k-mers | 3.53M | 2,250,006 | 1,973,059
Simulated reads from simulated genome with error | 5.92M | 3,757,479 | 2,133,592
Simulated reads from simulated genome without error | 9.07M | 5,714,973 | 1,989,644
Real E. coli reads | 4.85M | 4,004,911 | 2,079,302
Table 3:  Iterative low-memory k-mer trimming. The results of trimming reads at unique (erroneous) k-mers from a 5m read E. coli data set (1.4 GB) in under 30 MB of RAM. After each iteration, we measured the total number of distinct k-mers in the data set, the total number of unique (and likely erroneous) k-mers remaining, and the number of unique k-mers present at the 3’ end of reads. | FP rate | bases trimmed | distinct k-mers | unique k-mers | unique k-mers at 3’ end
---|---|---|---|---|---
untrimmed | - | - | 41.6m | 34.1m | 30.4%
khmer iteration 1 | 80.0% | 13.5% | 13.3m | 6.5m | 29.8%
khmer iteration 2 | 40.2% | 1.7% | 7.6m | 909.9k | 12.3%
khmer iteration 3 | 25.4% | 0.3% | 6.8m | 168.1k | 3.1%
khmer iteration 4 | 23.2% | 0.1% | 6.7m | 35.8k | 0.7%
khmer iteration 5 | 22.8% | 0.0% | 6.6m | 7.9k | 0.2%
khmer iteration 6 | 22.7% | 0.0% | 6.6m | 1.9k | 0.0%
filter by FASTX | - | 9.1% | 26.6m | 20.3m | 26.3%
filter by seqtk(default) | - | 8.9% | 17.7m | 12.1m | 12.3%
filter by seqtk(-q 0.01) | - | 15.4% | 9.9m | 5.1m | 5.2%
filter by seqtk(-b 3 -e 5) | - | 8.0% | 34.5m | 27.7m | 25.3%
Table 4:  Low-memory digital normalization. The results of digitally normalizing a 5m read E. coli data set (1.4 GB) to C=20 with k=20 under several memory usage/false positive rates. The false positive rate (column 1) is empirically determined. We measured reads remaining, number of “true” k-mers missing from the data at each step, and the number of total k-mers remaining. Note: at high false positive rates, reads are erroneously removed due to inflation of k-mer counts. memory | FP rate | retained reads | retained reads % | true k-mers missing | total k-mers
---|---|---|---|---|---
before diginorm | - | 5,000,000 | 100.0% | 170 | 41.6m
2400 MB | 0.0% | 1,656,518 | 33.0% | 172 | 28.1m
240 MB | 2.8% | 1,655,988 | 33.0% | 172 | 28.1m
120 MB | 18.0% | 1,652,273 | 33.0% | 172 | 28.1m
60 MB | 59.1% | 1,633,182 | 32.0% | 172 | 27.9m
40 MB | 83.2% | 1,602,437 | 32.0% | 172 | 27.6m
20 MB | 98.8% | 1,460,936 | 29.0% | 172 | 25.7m
10 MB | 100.0% | 1,076,958 | 21.0% | 185 | 20.9m
Table 5:  E. coli genome assembly after low-memory digital normalization. A comparison of assembling reads digitally normalized with low memory/high false positive rates. The reads were digitally normalized to C=20 (see [21] for more information) and were assembled using Velvet. We measured total length of assembly, as well as percent of true MG1655 genome covered by the assembly using QUAST. memory | FP rate | N contigs | total length(bases) | % of true genome covered
---|---|---|---|---
before diginorm | - | 106 | 4,546,051 | 97.84%
2400 MB | 0.0% | 617 | 4,549,235 | 98.05%
240 MB | 2.8% | 87 | 4,549,253 | 98.04%
120 MB | 18.0% | 86 | 4,549,335 | 98.04%
60 MB | 59.1% | 90 | 4,548,619 | 98.03%
40 MB | 83.2% | 89 | 4,550,599 | 98.11%
20 MB | 98.8% | 85 | 4,550,014 | 98.04%
10 MB | 100.0% | 97 | 4,545,871 | 97.97%


# Parallel sparse interpolation
using small primes

Mohamed Khochtali  Computer Science Department
United States Naval Academy
Annapolis, Maryland, USA Daniel S. Roche  Computer Science Department
United States Naval Academy
Annapolis, Maryland, USA Xisen Tian  Computer Science Department
United States Naval Academy
Annapolis, Maryland, USA

###### Abstract

To interpolate a supersparse polynomial with integer coefficients, two
alternative approaches are the Prony-based “big prime” technique, which acts
over a single large finite field, or the more recently-proposed “small primes”
technique, which reduces the unknown sparse polynomial to many low-degree
dense polynomials. While the latter technique has not yet reached the same
theoretical efficiency as Prony-based methods, it has an obvious potential for
parallelization. We present a heuristic “small primes” interpolation algorithm
and report on a low-level C implementation using FLINT and MPI.

##  1 Introduction

Given a way to evaluate or _sample_ an unknown function or procedure,
interpolation is the fundamental and important problem of recovering a formula
which accurately and completely describes that unknown function. As
discovering an arbitrary unknown function from a finite set of evaluations
with any reliability would be impossible, some constraints on the size and
form of the output are inevitably required.

Here we consider the problem of _sparse polynomial interpolation_ , in which
we are guaranteed that the unknown function is a multivariate polynomial with
bounded degree. Sparse interpolation algorithms date to the 18th century, but
have been the focus of considerable recent work in numeric and symbolic
computation, with applications ranging from power consumption in medical
devices, to reducing intermediate expression swell in mathematical
computations Cuyt and Lee (2008); Kaltofen (2010); Kaltofen et al. (2011);
Boyer et al. (2012).

Specifically, this work focuses on algorithms to interpolate an unknown
_supersparse_ polynomial with integer coefficients, which make efficient use
of modern _parallel_ computing hardware. Focusing on the “supersparse” (a.k.a.
“lacunary”) case means that our running time will be in terms of the number of
variables, number of nonzero terms, and the logarithm of the output degree.

The first algorithms to solve this problem in polynomial-time were based on
the exponential sums technique of Prony, and can efficiently solve the integer
problem by working in a large finite field modulo a single “big prime.” A
number of theoretical improvements and practical implementations have been
made in this vein, including work on fast parallel implementations Javadi and
Monagan (2010); van der Hoeven and Lecerf (2015).

We consider another type of approach, whereby the unknown sparse integer
polynomial is reduced in degree modulo many small primes. This technique,
first used by Garg and Schost to avoid the need for discrete logarithm
computations in arbitrary finite fields, can also be applied to the integer
polynomial case. While it cannot yet match the theoretical efficiency of the
big primes algorithms, we will show that the “small primes” method is very
effectively parallelized. Furthermore, we develop a practical heuristic
version of this method which reduces further the size and number of primes
required based on experimental results rather than on the theoretical worst-
case bounds.

###  1.1 Related work

While it would be impossible to list all of the related work on sparse
interpolation, we will mention some of the most recent results which are most
closely connected to the current study, and which may provide the reader with
useful background.

The now-classical approach to sparse interpolation is variously credited to
Prony, Blahut Blahut (1979), or Ben-Or and Tiwari Ben-Or and Tiwari (1988).
Only the last of these considered explicitly the case of integer coefficients,
but all share the key property of requiring the minimal number $O(T)$ of
evaluations in order to recover a polynomial with $T$ nonzero terms.

We refer to these approaches as “big prime” techniques, as the more modern
variants Kaltofen and Yagati (1989); Kaltofen (2010) adapt to the case of
integer coefficients by choosing carefully a single large modulus, then
perform the interpolation over a finite field in order to avoid exponential
growth in the bit-length of evaluations.

The approach which we take in this work is based on that of Garg and Schost
Garg and Schost (2009), who developed the first polynomial-time supersparse
interpolation algorithm over an arbitrary finite field. By reducing the
unknown polynomial modulo $(z^{p}-1)$, the full coefficients are discovered
immediately, but the exponents are only discovered after repeating for
multiple values of $p$. This “small primes” approach, described in more
details below, has the considerable advantage of relying only on low-degree,
dense polynomial arithmetic.

There are, of course, other sparse interpolation methods which do not fit
nicely into this big/small prime characterization. Most notable are Zippel’s
algorithm and hybrid variants of it Zippel (1990); Kaltofen and Lee (2003),
the symbolic-numeric method of Mansour (1995), and the Newton-Hensel lifting
approach of Avendaño et al. (2006).

We point out some recent work on efficient implementations which are of
particular interest to the current study. In Javadi and Monagan (2010), a new
variant of the big prime approach is developed which can be performed variable
by variable, in parallel. More recently, van der Hoeven and Lecerf (2015)
investigated a number of tricks and techniques towards practical, efficient
sparse interpolation, and posed some new benchmark problems. Their methods are
also based on the big prime approach; to our knowledge there has been no
reported implementation work on the small primes technique other than the
numerical interpolation code reported in Giesbrecht and Roche (2011).

###  1.2 Our contributions

Suppose $f\in\mathbb{Z}[x_{1},x_{2},\ldots x_{n}]$ is an unknown polynomial in
$n$ variables, with partial degrees less than $D$ in each variable, at most
$T$ nonzero terms, and coefficients less than $H$ in absolute value. Then $f$
can be written as

$f=\sum_{i=1}^{T}c_{i}x_{1}^{e_{i1}}x_{2}^{e_{i2}}\cdots{}x_{n}^{e_{in}},$

where each exponent $e_{ij}<D$ and each $|c_{i}|<H$.

Given a way to evaluate $f(\theta_{1},\ldots,\theta_{n})\bmod q$, for any
modulus $q\in\mathbb{Z}$ and $n$-tuple of evaluation points
$(\theta_{1},\ldots,\theta_{n})\in(\mathbb{Z}/q\mathbb{Z})^{n}$, the _sparse
integer polynomial interpolation problem_ is to determine the coefficients
$c_{i}$ and exponent tuples $(e_{i1},\ldots,e_{in})$, for each $1\leq i\leq
T$.

We will actually consider a slight relaxation of this problem, wherein
evaluations are of the form

$(q,p,d_{1},\ldots,d_{n})\mapsto
f(z^{d_{1}},\ldots,z^{d_{n}})\bmod(z^{p}-1)\in(\mathbb{Z}/q\mathbb{Z})[z].$

That is, the $n$ variables are replaced by a single indeterminate $z$, all
coefficients are reduced modulo $q$, and the exponents are reduced modulo $p$.
The reduction modulo $(z^{p}-1)$ is possible without affecting the overall
complexity whenever the unknown polynomial $f$ is given as a straight-line
program or algebraic circuit, or if the prime $q$ is chosen so that
$\mathbb{Z}/q\mathbb{Z}$ has a $p$th root of unity.

An algorithm for this problem is said to handle the _supersparse_ case if it
requires a number of evaluations and running time which are polynomial in $n$,
$T$, $\log D$, and $\log H$. This corresponds to the size of the sparse
representation of $f$ as a list of coefficient-exponent tuples, which requires
$O(nT\log D+T\log H)$ bits in memory.

We present a randomized algorithm for the sparse integer polynomial
interpolation problem, derived from the existing literature, whose running
time*** Here and throughout we use the _soft-oh_ notation in order to simplify
the stated running times: a running time is said to be $\widetilde{O}(\phi)$
for some function $\phi$ if and only if it is $O(\phi(\log\phi)^{O(1)})$.  for
provable correctness is

$\widetilde{O}\left(\left(1+\tfrac{1}{m}n\log D\right)nT\log D(\log D+\log
H)\right),$ (1)

where $m\geq 1$ is the number of parallel processors available for the task.
Furthermore, we demonstrate a heuristic variant on this algorithm, which works
well in our experimental testing, and reduces the running time further to

$\widetilde{O}\left(n\log D\log H+\left(1+\tfrac{1}{m}n\log D\right)T\log
H\right)$ (2)

in the typical case that $n$ and $\log D$ are both $O(T\log H)$.

We have implemented this heuristic approach using the C library FLINT for
dense polynomial arithmetic and MPI for parallelization. Our experiments
demonstrate the smallest effective settings for the parameters in our
heuristic approach. With those parameters, we show that the heuristic method
is competitive with the state of the art in the single-processor setting, and
that its running time scales well with increasing numbers of parallel
processors.

Specifically, our contributions are:

  1. 1.

A sparse interpolation algorithm whose potential parallel speedup is $O(n\log
D)$, compared with the $O(n)$ parallel speedup that has been shown in previous
work for other sparse interpolation algorithms.

  2. 2.

A heuristic variant of our algorithm, which is demonstrated to be effective on
our (limited) random experiments, that brings the running time complexity of
the small primes approach to be competitive with that of the big prime
approach.

  3. 3.

An efficient C implementation of our interpolation algorithm which
demonstrates its competitiveness on a standard benchmark problem.

###  1.3 Organization of the paper

We first outline the two main classes of existing approaches to supersparse
integer polynomial interpolation, which we call the “big prime” and the “small
primes” methods, in Section 2. We then present our own heuristic method in
Section 3, which is based on previously-known “small primes” algorithms, but
goes beyond theoretical worst-case bounds on the sizes needed in order to
further improve the efficiency.

Section 4 reports the details of our parallel implementation of this heuristic
method, and Section 5 presents the preliminary experimental results which
demonstrate the efficacy of this approach. Finally, we state some conclusions
and directions for further investigation in Section 6.

##  2 Existing algorithms for supersparse interpolation

###  2.1 “Big prime” methods

The original sparse interpolation algorithm of Ben-Or and Tiwari (1988) used
evaluations at powers of the first $n$ prime numbers along with Prony’s method
to deterministically recover the nonzero integer coefficients and exponents of
an unknown polynomial. This cannot be considered a “supersparse” algorithm, as
it performs the evaluations over the integers directly and requires working
with integers with more than $D$ bits.

However, it was soon recognized that, by choosing a single large prime $p\geq
D^{n}$, with $p-1$ smooth††† An integer is said to be _smooth_ if it has only
small prime factors.  so as to facilitate discrete logarithms, a supersparse
integer polynomial could be interpolated efficiently modulo $p$ Kaltofen et
al. (1990); Kaltofen (2010). The basic steps of this approach are as follows:

  1. 1.

Choose prime $p$ so that $(p-1)$ has a large “smooth” factor greater than
$D^{n}$ and let $\theta$ be a primitive element modulo $p$.

  2. 2.

For $i=0,1,2,\ldots,2T-1$, evaluate

$v_{i}=f(\theta^{i},\theta^{Di},\ldots,\theta^{D^{n-1}i})\bmod p.$

  3. 3.

Compute the minimum polynomial $\Gamma\in\mathbb{Z}_{p}[z]$ of the sequence
$(v_{i})_{i\geq 0}$ with the Berlekamp-Massey algorithm.

  4. 4.

Factor $\Gamma$ over $\mathbb{Z}_{p}[z]$; each root of $\Gamma$ can be written
as $\theta^{e_{1}+e_{2}D+\cdots+e_{n}D^{n-1}}$, corresponding to a single term
$cx_{1}^{e_{i1}}\cdots x_{n}^{e_{in}}$ of $f$.

  5. 5.

Compute $T$ discrete logarithms of the roots, and then the $D$-adic expansion
of each one, to discover the actual exponents $(e_{i1},\ldots,e_{in})$, for
$1\leq i\leq T$.

  6. 6.

Once the exponents are known, the coefficients can be computed from the
evaluations $v_{i}$ by solving a transposed Vandermonde system of dimension
$T$.

The two steps which can be trivially parallelized are the evaluations and
discrete log computations on steps 2 and 5. However, the dominating cost in
the complexity is factoring a polynomial over $\mathbb{Z}_{p}[x]$ on step 4\.
This is also confirmed to be the dominating cost in practice by van der Hoeven
and Lecerf (2015), and it is not clear how to efficiently parallelize the
factorization. Using the fastest known algorithms, the running time of this
step is $\widetilde{O}(T\log^{2}p)$ Grenet et al. (2015), which is at least
$\widetilde{O}(n^{2}T\log^{2}D)$ bit operations from the condition $p>D^{n}$.

In the description above, the evaluations at powers of $\theta$ on step 2
amount to a Kronecker substitution from multivariate to univariate. Of note is
the algorithm of Javadi and Monagan (2010), which uses a different approach
than Kronecker substitution in order to work one variable at a time, gaining a
potential $n$-fold parallel speedup.

###  2.2 “Small primes” methods

The algorithm described above works the same over an arbitrary finite field,
except that the discrete logarithms required in step 5 cannot be performed in
polynomial-time in general. This difficulty was first overcome in Garg and
Schost (2009), where the idea is as follows:

  1. 1.

Choose a series of small primes $p_{1},p_{2},\ldots$

  2. 2.

Apply the Kronecker substitution and for each $i=1,2,\ldots$ compute
$f_{i}=f(z,z^{D},\ldots,z^{D^{n-1}})\bmod(x^{p_{i}}-1)$.

  3. 3.

Each $f_{i}$ of maximal sparsity contains all the coefficients of $f$, and all
the exponents modulo $p_{i}$. Collect sufficiently many modular images of the
exponents in order to recover the full exponents over $\mathbb{Z}$.

  4. 4.

Recover the multivariate exponents by $D$-adic expansion of each univariate
exponent, and use any $f_{i}$ of maximal sparsity to discover the coefficient
of each term.

There are two significant challenges of this approach. The first is that, in
reducing the polynomial modulo $x^{p_{i}}-1$, it is possible that some
exponents are equivalent modulo $p_{i}$, and then multiple terms in the
original polynomial will _collide_ to form a single term in $f_{i}$. By
choosing random primes whose values are roughly $\widetilde{O}(T^{2}\log D)$,
the probability of encountering any collisions can be made arbitrarily small.

The second challenge is how to correlate the exponents from different
$f_{i}$’s in step 3, in order to recover the full exponents via Chinese
remaindering. The approach of Garg and Schost (2009) was to compute an
auxiliary polynomial whose roots are the unknown exponents; however this
increases the overall running time to $\widetilde{O}(n^{2}T^{3}\log^{2}D\log
H)$.

The technique of _diversification_ in Giesbrecht and Roche (2011) is another
randomization which chooses a random element $\alpha$ and interpolates
$f(\alpha z)$ instead of $f(z)$ itself. With high probability, the
“diversified” polynomial $f(\alpha z)$ has distinct coefficients, which can
then be used to correlate the exponents in different $f_{i}$’s. This avoids
the factoring step and reduces the complexity by a factor of $T$.

Subsequently, and separately, Arnold et al. (2013) showed how to allow the
magnitude of each $p_{i}$ to decrease by a factor of $T$, by allowing some
constant fraction of the terms in each $f_{i}$ to collide. These separate
approaches are combined and further improved in Arnold et al. (2015), which in
the typical case that $n\ll\log D\ll\log H$, brings the theoretical complexity
down to $\widetilde{O}(nT\log^{2}D\log H)$. Observe that this is competitive
with the “big primes” algorithm, but could be slower by up to a factor of
$\log H$.

The algorithm we describe next first shows how to effectively parallelize this
small primes approach, and then further reduces the complexity through a
heuristic argument. After the gains from parallelization or from the heuristic
improvement, the complexity of our algorithm will be less than the “big
primes” approach as well.

##  3 Our parallel small primes algorithm

[tb] SparseInterp($f,n,T,D,H,k,\ell$) Input: Bounds $T,D,H$ for an $n$-variate
sparse polynomial $f$ in $n$ variables, with $T$ nonzero terms, partial
degrees less than $D$, and integer coefficients less than $H$ in absolute
value; and parameters $k,\ell,Q\in\mathbb{Z}_{>0}$. Output: A list of $t\leq
T$ coefficients $c_{i}$ and exponent tuples $(e_{i1},\ldots,e_{in})$. If
$k,\ell$ are sufficiently large, these coefficients and exponents comprise the
sparse representation of $f$ with high probability. $q\leftarrow$ random prime
in the range $[Q,2Q]$

2 $\alpha\leftarrow$ random element of $(\mathbb{Z}/q\mathbb{Z})^{*}$

3 $(\alpha_{0},\ldots,\alpha_{n-1})\leftarrow(\alpha,\alpha^{D}\bmod
q,\ldots,\alpha^{D^{n-1}}\bmod q)$

4 $\lambda\leftarrow kT$

5 $\mu\leftarrow\lceil(\ell n\lg D)/\lg\lambda\rceil$

6 $L\leftarrow$ thread-safe list of integer triples

7 for  _$i=1,2,\ldots,\mu$ in parallel_ do

8 $p_{i}\leftarrow$ random prime in the range $[\lambda,2\lambda]$

9 $(D_{0},\ldots,D_{n-1})\leftarrow(1,D\bmod p_{i},\ldots,D^{n-1}\bmod p_{i})$

10 $f_{i}\leftarrow
f(\alpha_{0}z^{D_{0}},\ldots,\alpha_{n-1}z^{D_{n-1}})\bmod(z^{p_{i}}-1)\in(\mathbb{Z}/q\mathbb{Z})[z]$
using dense polynomial arithmetic

11 for  _$j=0,1,\ldots,p_{i}-1$_ do

12 $c_{ij}\leftarrow$ coefficient of $z^{j}$ in $f_{i}$

13 if  _$c_{ij}\neq 0$_ then Add $(c_{ij},j,p_{i})$ to $L$

14

15

16Sort $L$ in parallel by the coefficients $c_{ij}$ in each triple

17 $F\leftarrow$ thread-safe list of coefficient/exponent tuples

18 foreach  _Unique coefficient $c$ in $L$ in parallel_ do

19 $(d_{c1},p_{c1}),\ldots,(d_{cu},p_{cu})\leftarrow$ the $u\geq 1$ exponents
and primes appearing with coefficient $c$ in $L$

20 if  _$u\geq\mu/2$_ then

21 $u^{\prime}\leftarrow$ least integer s.t. $\prod_{1\leq i\leq
u^{\prime}}p_{ci}\geq D^{n}$

22 $E_{c}\leftarrow$ least integer s.t. $E_{c}\equiv d_{ci}\bmod p_{ci}$ for
$0\leq i\leq u^{\prime}$ via Chinese remaindering

23 $c^{*}\leftarrow c\alpha^{-E_{c}}\in\mathbb{Z}/q\mathbb{Z}$, stored as an
integer in the range $[-q/2,q/2]$

24 $(e_{1},\ldots,e_{n})\leftarrow$ $D$-adic expansion of $E_{c}$

25 Add $c^{*}$ and $(e_{1},\ldots,e_{n})$ to $F$

26

27

Sort $F$ in parallel by the exponents and return the resulting sparse
polynomial

Our algorithm, which is detailed in procedure 3, is based on the reduction
idea in Garg and Schost (2009), with the diversification method introduced by
Giesbrecht and Roche (2011) and the partial collisions handling of Arnold et
al. (2013). It depends crucially on the parameters $k$ and $\ell$, which in
theory grow as $O(n\log D)$ and $O(1)$, respectively, but according to our
heuristic method can both be treated as constants.

###  3.1 Algorithm overview

The main idea is first to choose a prime $q$ large enough to recover the
coefficients, and then to apply the Kronecker substitution so that we are
really interpolating a univariate polynomial $g(z)$ with coefficients in
$\mathbb{Z}/q\mathbb{Z}$:

$g=f(\alpha z,(\alpha z)^{D},\ldots,(\alpha
z)^{D^{n-1}})=\sum_{i=1}^{T}c_{i}z^{E_{i}}\in(\mathbb{Z}/q\mathbb{Z})[z].$

Each term $c^{*}x_{1}^{e_{1}}\cdots x_{n}^{e_{n}}$ of the original polynomial
$f$ maps uniquely to a term with exponent
$E_{i}=e_{1}+e_{2}D+\cdots+e_{n}D^{n-1}$ and coefficient
$c_{i}=c^{*}\alpha^{E_{i}}\bmod q$ in $g$.

The parameter $Q$ controls the size of the prime $q$, and so should always be
set larger than $2H$ in order to recover the full precision of the
coefficients. However, it may be necessary to set $Q$ even larger than this
when the height bound $H$ is very small. We comment that choosing a random
prime $q$ (rather than one with special properties, as in the “big primes”
algorithm) is important for the probabilistic analysis below.

The evaluation phase of the algorithm computes the polynomial $g$ modulo
$(z^{p}-1)$ using dense arithmetic, for many small primes $p$. The details of
how this evaluation is performed will depend on the particular application. In
our multivariate multiplication application below, the exponents of the
original multiplicands are reduced modulo $p$, followed by dense polynomial
arithmetic in the ring of polynomials modulo $z^{p}-1$. More generally, if the
unknown polynomial $f$ is given as a straight-line program or arithmetic
circuit, each operation in the circuit can be computed over that ring
$(\mathbb{Z}/q\mathbb{Z})[z]/\langle z^{p}-1\rangle$, as in Garg and Schost
(2009). In our complexity analysis below, we assume a cost of
$\widetilde{O}(p\log q)$ for each evaluation on this step, according to the
cost of dense degree-$p$ polynomial arithmetic over $\mathbb{Z}/q\mathbb{Z}$.

The next phase of the algorithm is to gather the images of each term, discard
those which appear infrequently (and thus were resulting from collisions in
the reduction modulo $z^{p}-1$), and use Chinese remaindering to recover the
exponents of each nonzero term in $f$.

Observe that the list $L$ can be implemented in any convenient way according
to the details of the parallel implementation; it serves only as an unordered
accumulator of coefficient-exponent-prime tuples. The output polynomial $f$
could also be considered as an unordered accumulator, followed by another
parallel sort before the final return statement.

###  3.2 Parallel complexity analysis

We state the parameterized running time of the algorithm as follows:

######  Theorem 1.

Given bounds $T,D$ on the sparsity and degree of an unknown polynomial
$f\in\mathbb{Z}[x_{1},\ldots,x_{n}]$, and parameters $k,\ell,Q$, Algorithm 3
has worst-case running time

$\displaystyle\widetilde{O}\bigg{(}$ $\displaystyle\log\ell+n\log D\log
Q+kT\log Q$ $\displaystyle+\left(\tfrac{1}{m}\ell n\log D\right)\left(n+\log
D+kT\log Q\right)\bigg{)}.$

###### Proof.

The computation of $\mu$ at the beginning incurs the $\widetilde{O}(\log\ell)$
cost in the complexity, which for our ultimate choices of $\ell$ will never
actually dominate the complexity.

The first loop executes $\lceil\mu/m\rceil\in\widetilde{O}(1+\tfrac{1}{m}\ell
n\log D)$ times in each thread. Computing the powers of $D$ modulo $p_{i}$ on
Step 3 requires $\widetilde{O}(\log D+n\log p_{i})$ bit operations. The
subsequent evaluations using dense arithmetic on Step 3 will usually dominate
the complexity of the entire algorithm, as each costs $\widetilde{O}(p_{i}\log
q)$, which is $\widetilde{O}(kT\log Q)$. (The addition of this term makes the
$\log p_{i}$ factor in the cost of Step 3 become $\widetilde{O}(1)$.)

Both parallel sorts are on lists of size at most $\mu T$, which means their
cost of $\widetilde{O}(\tfrac{1}{m}\ell nT\log D)$ does not dominate the
complexity.

The final for loop executes $O(1+\tfrac{1}{m}\mu T)$ times in each thread, but
the nested if statement can only be triggered $O(T)$ times overall. Within it,
the most expensive step is computing each $\alpha^{-E_{c}}\bmod q$, requiring
$\widetilde{O}(n\log D\log Q)$ bit operations. This contributes only
$\widetilde{O}(n\log D\log Q)$ to the overall complexity, as the term
$\frac{1}{m}Tn\log D\log Q$ is already dominated by the parallel cost of Step
3. ∎

The key feature of procedure 3 is that its potential parallel speedup, from
the previous theorem, is a factor of $O(\ell n\log D)$, depending on the
number of parallel processors $m$ that are available. This exceeds the $O(n)$
parallel speedup of previous approaches, and means that our algorithm should
scale better to a large number of processors when the number of variables
and/or degree are sufficiently large.

###  3.3 Correctness and probability analysis

We first use prior work to prove the bounds necessary to ensure correctness
with provably high probability.

######  Theorem 2.

If the parameters $k,\ell,Q$ and bounds $D,T,H$ satisfy

$\displaystyle Q$ $\displaystyle\geq\max\left(2H,\tfrac{1}{4}(\ell nT\lg
D)^{2}D\right),$ $\displaystyle k$ $\displaystyle\geq\max(21,\lceil 20n\ln
D\rceil),\text{and}$ $\displaystyle\ell$ $\displaystyle\geq 2,$

then with probability at least $1/2$, procedure 3 correctly computes the
coefficients and exponents of $f\in\mathbb{Z}[x_{1},\ldots,x_{n}]$.

###### Proof.

The condition $Q\geq 2H$ guarantees that positive or negative coefficients
whose absolute value is at most $H$ will all be distinct modulo $q$, for any
prime $q\geq Q$ as chosen in the algorithm.

Because $\lambda=kT$, Lemma 8 from Arnold et al. (2013) tells us that, for
each $p_{i}$, the probability that more than $T/3$ terms collide modulo
$z^{p_{i}}-1$ is at most $1/4$. Since the most number of terms that could
collide is $T$, this means the _expected_ number of collisions, for each
$p_{i}$, is at most $T/2$.

The total number of nonzero coefficients in all polynomials $f_{i}$ examined
on step 3 is at most $\mu T\leq\ell nT\lg D$. Many of these correspond to
single terms in $f$ itself, but some will be collisions of terms in $f$. Using
the reasoning of Lemma 4.1 in Arnold et al. (2014), and the proof of Theorem
3.1 from Giesbrecht and Roche (2011), every coefficient of a single term in
$f$, or a collision that appears in any of the $f_{i}$’s, is distinct modulo
$q$, due to the bound on $Q$ and because $q\geq Q\geq\tfrac{1}{2}(\mu
T)^{2}\lg D$.

We conclude that, with probability at least $1/2$, each term in $f$ appears
un-collided in at least $\mu/2$ of the polynomials $f_{i}$, and furthermore
that these coefficients are all distinct and are the only ones that repeat in
the list $L$. This guarantees that the correct coefficients and exponents are
recovered in the final loop, and the algorithm outputs the correct
interpolated polynomial. ∎

Applying the bounds in Theorem 2 to the analysis in Theorem 1 gives the
provable complexity bound for the algorithm as stated in (1).

As usual, the probability of success in either the provable or the heuristic
version of the algorithm (described below) can be increased arbitrarily high
by running the same algorithm repeatedly and choosing the most common
polynomial returned among all runs to be the most likely candidate for $f$.

###  3.4 Heuristic approach

The heuristic version of this procedure is simply to choose appropriate
constants for the parameters $k$ and $\ell$, with the intuition that there can
be some trade-off between the size of each prime (governed by $k$) and the
number of chosen primes (governed by $\ell$). Furthermore, the bound on $Q$
required in theory to obtain diversification between the coefficients in $f$
and in any collisions is unnecessarily high for most “typical” polynomials.
There do exist pathological counterexamples, but they require the degrees of
many terms to be equivalent modulo $q$. As the prime $q$ is also chosen
randomly in our approach, we have a good indication that this heuristic
approach will work with high probability for a randomly-chosen sparse
polynomial. We state this as a conjecture, which is also backed by the
experimental evidence reported in Section 5 below.

######  Conjecture 3.

For any sufficiently large height bound $H$, and using $Q=2H$, there exist
constants $k,\ell\geq 1$ such that, for a polynomial
$f\in\mathbb{Z}[x_{1},\ldots,x_{n}]$ chosen at random with at most $T$ nonzero
terms, the probability that algorithm 3 successfully interpolates $f$ is at
least 1/2.

The heuristic complexity under this conjecture is stated in (2).

###  3.5 Example

Consider an unknown bivariate polynomial $f(x,y)$ with 3 nonzero terms and
degree less than 10.

The primes that we are going to use are going to be small for the sake of this
example. The primes are 7, 13, and 17.

  1. 1.

Now we compute $f$ modulo $(z^{p}_{i}-1)$, we get:

$\displaystyle f\bmod(z^{17}-1)$ $\displaystyle=2z^{9}+7z^{8}+3z^{3}$
$\displaystyle f\bmod(z^{13}-1)$ $\displaystyle=10z^{7}+2z^{4}$ $\displaystyle
f\bmod(z^{7}-1)$ $\displaystyle=3z^{6}+7z^{3}+2z$

We notice here that a collision happened with the prime 13. This won’t affect
our calculation later because resulting coefficient $10$ doesn’t appear
anywhere else (For this example, we are going to assume that we have a good
coefficient if it appears twice or more).

Now we use these values to fill the list $L$. Every triple in $L$ consists of
the coefficient, the prime, and the exponent.

$\begin{array}[]{c|*{8}{c|}}\cline{2-9}\cr\text{coefficient}&2&7&3&10&2&3&7&2\\\
\cline{2-9}\cr\text{exponent}&9&8&3&7&4&6&3&1\\\
\cline{2-9}\cr\text{prime}&17&17&17&13&13&7&7&7\\\ \cline{2-9}\cr\end{array}$

We then sort $L$ based on coefficient values.

$\begin{array}[]{c|*{8}{c|}}\cline{2-9}\cr\text{coefficient}&10&7&7&3&3&2&2&2\\\
\cline{2-9}\cr\text{exponent}&7&8&3&3&6&9&4&1\\\
\cline{2-9}\cr\text{prime}&13&17&7&17&7&17&13&7\\\ \cline{2-9}\cr\end{array}$

We use the Chinese remainder theorem to get back the exponent that corresponds
to each coefficient.

$\left.\begin{array}[]{r@{\ }c@{\ }l}e_{1}\bmod 17&=&8\\\ e_{1}\bmod 7&=&3\\\
\end{array}\right\\}\Rightarrow e_{1}=59$

$\left.\begin{array}[]{r@{\ }c@{\ }l}e_{2}\bmod 17&=&3\\\ e_{2}\bmod 7&=&6\\\
\end{array}\right\\}\Rightarrow e_{2}=20$

$\left.\begin{array}[]{r@{\ }c@{\ }l}e_{3}\bmod 17&=&9\\\ e_{3}\bmod 13&=&4\\\
e_{3}\bmod 7&=&1\\\ \end{array}\right\\}\Rightarrow e_{3}=43$

This results in the univariate polynomial

$f(z)=3z^{20}+2z^{43}+7z^{59}.$

Finally, inverting the Kronecker map $f(z)=f(x,y^{10})$, we obtain

$f(x,y)=3y^{2}+2x^{3}y^{4}+7x^{9}y^{5}.$

##  4 Parallel implementation

We completed a low-level implementation of Procedure 3 written in the C
programming language. Our complete source code, as well as the exact source we
tested for the comparisons and benchmark problems listed later, is available
upon request by email.

We give a few details here on the choices of our implementation, in particular
the libraries that were utilized.

###  4.1 FLINT for sparse and dense polynomial arithmetic

The key advantage to the “small primes” approach which we employed is the
reliance on fast subroutines for dense polynomial arithmetic. The experiments
we ran always used a word-sized modulus $q$, and so the most expensive
computations involved computing with dense, low-degree polynomials with word-
sized modular coefficients.

FLINT (http://flintlib.org/) is a free, open-source C library for fast number
theoretic computations Hart et al. (2013). Our dense polynomial arithmetic,
which is the dominating cost both in theory and in practice in our
experiments, was performed using the `nmod_poly` data type.

In order to store the result of sparse interpolation and complete the
correctness testing in our experiments, we also added rudimentary support for
sparse integer polynomials on top of FLINT. We created a new data type,
`fmpz_sparse`, to represent sparse univariate polynomials in $\mathbb{Z}[x]$
for testing purposes, using FLINT’s multiple precision type `fmpz` as both the
coefficient and exponent storage for sparse polynomials.

Note that using multiple-precision integers for exponents is especially
important, as we have _not_ yet completed full multivariate polynomial support
within FLINT. Instead, for the purposes of our experiments, we always used a
standard Kronecker substitution to store $n$-variate, degree-$D$ multivariate
polynomials as univariate sparse polynomials with degree bounded by $D^{n}$.
Employing a multiple-precision data type allows for the largest possible
degree and number of variables, which is crucial as it is precisely this
_supersparse_ case in which our approach has the greatest potential parallel
speedup.

###  4.2 MPI for multi-processor parallelism

Message Passing Interface (MPI) allows us to parallelize our algorithm. As
stated earlier, since the slowest part of our algorithm involves calculating
the unknown polynomial $f$ modulo many polynomials $(x^{p_{i}}-1)$, we use MPI
to perform each of these evaluations in parallel.

The function `MPI_Init` is called at the beginning of the program to spawn an
arbitrary number of processes, as specified on the command line. Each of the
allocated processes will be executing separately with separate copies of all
variables in the original process. All processes will have unique id numbers.
The root process will have id number 0. By knowing processes id’s we can
separate what each process executes.

We used a master-slave model for our algorithm. The master process evenly
distributes how many primes each slave calculates. After getting the primes,
the slave process will compute the following for every prime:
$\prod_{i=1}^{m}f_{i}(x)\bmod(x^{p_{i}}-1)$, then traverse the resulting
univariate polynomial and save all nonzero terms, along with the prime
$p_{i}$, to an array. This array of triples is sent back to the master
process, which later sorts all the concatenated evaluation arrays and uses
Chinese remaindering to recover the full polynomial $f$, as described above.

While our experiments were performed on a multi-core machine, and hence using
a simpler threading library would have also worked, our goal in using MPI was
to demonstrate the full parallel potential of this approach by explicitly
detailing the inter-process communication. Furthermore, the MPI implementation
could also be used without modification on heterogeneous clusters or other
architectures besides multi-core.

##  5 Experimental results

We ran our tests on a machine using an Intel(R) Core(TM) i7-3930K CPU @
3.20GHz simulating 12 hyper-threaded cores on 6 physical cores, with 32 GB of
RAM. We used a Debian GNU system, running the “unstable” branch, with Linux
kernel version 3.16.0-4-amd64. This is a bleeding-edge system with the most
current versions of all software available within the Debian repositories.

###  5.1 Determining the parameters $k$ and $\ell$

The first task was to determine experimentally what kind of settings for the
parameters $k$ and $\ell$ would be appropriate for our heuristic interpolation
method. We found that $k=38$ and $\ell=2$ worked for a wide range of problem
sizes with almost no failures in the randomized algorithm.

The only exceptions we found here were that when the bound on the number of
terms $T$ in the output was very small, even setting $k=38$ the number of
primes $p$ in the range $[\lambda,2\lambda]$ was simply not sufficient to
ensure a high success probability. In these small-size extremes, the value of
$k$ was increased to accommodate; in particular, we settled on $k=50$ and
$\ell=2$ when $T<1000$, and when $T<100$ we had to select $k\geq 10000/T$.
Under these parameter settings, no failures were observed in any of our
experiments.

We comment that a smaller $k$ value and a larger $\ell$ value would be
preferable, because that would increase the number of primes $\mu$ and hence
the potential parallel speedup for the algorithm, while decreasing the size
$\lambda$ of each prime $p_{i}$. However, we found that modestly larger values
for $\ell$ did not allow for $k$ to be reduced and consistently result in
correct output. We consider the exploration of some better balance between the
parameters $k$ and $\ell$ as future work.

###  5.2 Single-threaded performance

$m$ variables terms max-degree $\mu$ $\lambda$ Mathemagix Ours (single thread)
Ours (12 threads) 1 20 3 40 14 50 000 0.078 0.035 0.029 2 20 9 80 18 5 000
0.155 0.151 0.048 3 20 27 120 18 6 900 0.305 0.329 0.116 4 20 81 160 18 4 900
0.598 0.323 0.085 5 20 243 200 17 9 900 2.156 0.785 0.175 6 20 729 240 15 39
900 5.053 3.084 0.814 7 20 2187 280 14 80 050 13.333 8.714 2.225 8 20 6561 320
13 321 300 41.070 43.911 10.605

Table 1: Sparse interpolation benchmark times (in seconds) Figure 1:
Comparison With Mathemagix Sparse Interpolation Program

The first experiment was to test our algorithm without parallelization against
the efficient “big primes” implementation in Mathemagix
(http://www.mathemagix.org/), as reported in van der Hoeven and Lecerf (2015).
We downloaded the source code for the Mathemagix sparse interpolation program,
then ran it with $m=1,2,3,4,5,6,7,8$ polynomials being multiplied. We kept
each polynomial at 20 variables, with degree 40, 3 terms and coefficients up
to $2^{30}$. Then we ran our algorithm with the same parameters. The results
are shown in Taible 1 and summarized in Figure 1. (We also attempted a
comparison against the sinterp function in Maple 2015, but it was more than an
order of magnitude slower in all experiments.)

Our non-parallelized algorithm seems to be on par with the Mathemagix
implementation for this range of problem sizes. As a comparison, and to
emphasize the main point of our paper, we also show the full parallel speedup
for the same problems in Figure 1.

###  5.3 Parallel speedup

Figure 2: Parallel speedup for our implementation

The second experiment was to test the parallel speedup for our implementation
by varying the degree size $D$ and leaving $k=38$, $l=2$, partial $f_{i}$, and
$m=6$ constant. $D$ was varied from $20^{25}$, $20^{50}$, $20^{100}$. This
test was performed 3 separate times and the resultant data was calculated from
median of the three trials. The results are seen in Figure 2.

Figure 2 shows a linear speedup increase as the number of parallel processes
used gets closer to the physical number of cores on the machine. Additionally,
we can see the most significant parallel speedup occurs for the highest degree
tested. Recall that on our machine, there are only 6 physical cores that are
hyper-threaded to 12 virtual cores. Furthermore, when running only six
threads, the “turbo mode” clock rate is increased to 3.8GHz. This may help to
explain the dip in performance seen around $m=6$ parallel processes.

Observe also that the parallel speedup is best, and continues the furthest, in
the most extremely sparse case with very high degree.

Our parallel speedup is demonstrated in Figure 2.

##  6 Conclusions and future work

We have shown that the “small primes” sparse interpolation algorithm is
competitive with the state of the art, even without parallelization,
especially for very sparse problem instances. Furthermore, there is greater
potential for parallelism in the small primes technique. These theoretical
results are borne out in practice in our experimental results compared to
other available software implementations.

There is significantly more work to be done, however, before we might suggest
widespread adoption of our heuristic sparse interpolation method. We would
like to understand the theory behind the heuristic approach in order to have a
less _ad hoc_ way of determining the parameters $k$ and $\ell$. On the other
hand, our implementation could be greatly enhanced with further
experimentation on a wider range of benchmark problems and incorporating true
multivariate sparse polynomial representations.

## Acknowledgments

This work was supported by the National Science Foundation under award number
#1319994. We thank Andrew Arnold and the anonymous PASCO 2015 reviewers for
their comments on an earlier draft of this paper.

## References

  * Arnold et al. [2013] Andrew Arnold, Mark Giesbrecht, and Daniel S. Roche.  Faster sparse interpolation of straight-line programs.  In Vladimir P. Gerdt, Wolfram Koepf, Ernst W. Mayr, and Evgenii V. Vorozhtsov, editors, _Proc. Computer Algebra in Scientific Computing (CASC 2013)_ , volume 8136 of _Lecture Notes in Computer Science_ , pages 61–74. Springer, September 2013.  doi:10.1007/978-3-319-02297-0_5.
  * Arnold et al. [2014] Andrew Arnold, Mark Giesbrecht, and Daniel S. Roche.  Sparse interpolation over finite fields via low-order roots of unity.  In _Proceedings of the 39th International Symposium on Symbolic and Algebraic Computation_ , ISSAC ’14, pages 27–34, New York, NY, USA, 2014. ACM.  doi:10.1145/2608628.2608671.
  * Arnold et al. [2015] Andrew Arnold, Mark Giesbrecht, and Daniel S. Roche.  Faster sparse multivariate polynomial interpolation of straight-line programs.  _CoRR_ , abs/1412.4088, 2015.  URL http://arxiv.org/abs/1412.4088.
  * Avendaño et al. [2006] Martín Avendaño, Teresa Krick, and Ariel Pacetti.  Newton-Hensel interpolation lifting.  _Found. Comput. Math._ , 6(1):81–120, 2006.  ISSN 1615-3375.  doi:10.1007/s10208-005-0172-3.
  * Ben-Or and Tiwari [1988] Michael Ben-Or and Prasoon Tiwari.  A deterministic algorithm for sparse multivariate polynomial interpolation.  In _Proceedings of the twentieth annual ACM symposium on Theory of computing_ , STOC ’88, pages 301–309, New York, NY, USA, 1988. ACM.  ISBN 0-89791-264-0.  doi:10.1145/62212.62241.
  * Blahut [1979] R.E. Blahut.  Transform techniques for error control codes.  _IBM Journal of Research and Development_ , 23(3):299–315, May 1979.  doi:10.1147/rd.233.0299.
  * Boyer et al. [2012] Brice Boyer, Matthew T. Comer, and Erich L. Kaltofen.  Sparse polynomial interpolation by variable shift in the presence of noise and outliers in the evaluations.  In _Electr. Proc. Tenth Asian Symposium on Computer Mathematics (ASCM 2012)_ , 2012.
  * Cuyt and Lee [2008] Annie Cuyt and Wen-shin Lee.  A new algorithm for sparse interpolation of multivariate polynomials.  _Theoretical Computer Science_ , 409(2):180–185, 2008.  ISSN 0304-3975.  doi:10.1016/j.tcs.2008.09.002.  Symbolic-Numerical Computations.
  * Garg and Schost [2009] Sanchit Garg and Éric Schost.  Interpolation of polynomials given by straight-line programs.  _Theoretical Computer Science_ , 410(27-29):2659–2662, 2009.  ISSN 0304-3975.  doi:10.1016/j.tcs.2009.03.030.
  * Giesbrecht and Roche [2011] Mark Giesbrecht and Daniel S. Roche.  Diversification improves interpolation.  In _Proceedings of the 36th international symposium on Symbolic and algebraic computation_ , ISSAC ’11, pages 123–130, New York, NY, USA, 2011\. ACM.  ISBN 978-1-4503-0675-1.  doi:10.1145/1993886.1993909.
  * Grenet et al. [2015] Bruno Grenet, Joris van der Hoeven, and Grégoire Lecerf.  Randomized root finding over finite FFT-fields using tangent Graeffe transforms.  In _Proc. 40th International Symposium on Symbolic and Algebraic Computation_ , ISSAC ’15, page to appear, 2015.
  * Hart et al. [2013] W. Hart, F. Johansson, and S. Pancratz.  FLINT: Fast Library for Number Theory, 2013.  Version 2.4.0, http://flintlib.org.
  * van der Hoeven and Lecerf [2015] Joris van der Hoeven and Grégoire Lecerf.  Sparse polynomial interpolation in practice.  _ACM Commun. Comput. Algebra_ , 48(3/4):187–191, February 2015.  doi:10.1145/2733693.2733721.
  * Javadi and Monagan [2010] Seyed Mohammad Mahdi Javadi and Michael Monagan.  Parallel sparse polynomial interpolation over finite fields.  In _Proceedings of the 4th International Workshop on Parallel and Symbolic Computation_ , PASCO ’10, pages 160–168, New York, NY, USA, 2010\. ACM.  ISBN 978-1-4503-0067-4.  doi:10.1145/1837210.1837233.
  * Kaltofen and Lee [2003] Erich Kaltofen and Wen-shin Lee.  Early termination in sparse interpolation algorithms.  _Journal of Symbolic Computation_ , 36(3-4):365–400, 2003.  ISSN 0747-7171.  doi:10.1016/S0747-7171(03)00088-9.  ISSAC 2002.
  * Kaltofen and Yagati [1989] Erich Kaltofen and Lakshman Yagati.  Improved sparse multivariate polynomial interpolation algorithms.  In P. Gianni, editor, _Symbolic and Algebraic Computation_ , volume 358 of _Lecture Notes in Computer Science_ , pages 467–474. Springer Berlin / Heidelberg, 1989.  doi:10.1007/3-540-51084-2_44.
  * Kaltofen et al. [1990] Erich Kaltofen, Y. N. Lakshman, and John-Michael Wiley.  Modular rational sparse multivariate polynomial interpolation.  In _Proceedings of the international symposium on Symbolic and algebraic computation_ , ISSAC ’90, pages 135–139, New York, NY, USA, 1990. ACM.  ISBN 0-201-54892-5.  doi:10.1145/96877.96912.
  * Kaltofen [2010] Erich L. Kaltofen.  Fifteen years after DSC and WLSS2: What parallel computations I do today [invited lecture at PASCO 2010].  In _Proceedings of the 4th International Workshop on Parallel and Symbolic Computation_ , PASCO ’10, pages 10–17, New York, NY, USA, 2010. ACM.  ISBN 978-1-4503-0067-4.  doi:10.1145/1837210.1837213.
  * Kaltofen et al. [2011] Erich L. Kaltofen, Wen-shin Lee, and Zhengfeng Yang.  Fast estimates of hankel matrix condition numbers and numeric sparse interpolation.  In _Proceedings of the 2011 International Workshop on Symbolic-Numeric Computation_ , SNC ’11, pages 130–136, New York, NY, USA, 2011\. ACM.  ISBN 978-1-4503-0515-0.  doi:10.1145/2331684.2331704.
  * Mansour [1995] Yishay Mansour.  Randomized interpolation and approximation of sparse polynomials.  _SIAM Journal on Computing_ , 24(2):357–368, 1995\.  doi:10.1137/S0097539792239291.
  * Zippel [1990] Richard Zippel.  Interpolating polynomials from their values.  _Journal of Symbolic Computation_ , 9(3):375–403, 1990.  ISSN 0747-7171.  doi:10.1016/S0747-7171(08)80018-1.  Computational algebraic complexity editorial.