1

Sub-string/Pattern Matching in Sub-linear Time
Using a Sparse Fourier Transform Approach
Nagaraj T. Janakiraman, Avinash Vem, Krishna R. Narayanan, Jean-Francois Chamberland
Department of Electrical and Computer Engineering

arXiv:1704.07852v1 [] 25 Apr 2017

Texas A&M University
{tjnagaraj,vemavinash,krn,chmbrlnd}@tamu.edu

Abstract
We consider the problem of querying a string (or, a database) of length N bits to determine all
the locations where a substring (query) of length M appears either exactly or is within a Hamming
distance of K from the query. We assume that sketches of the original signal can be computed off line
and stored. Using the sparse Fourier transform computation based approach introduced by Pawar and
Ramchandran, we show that all such matches can be determined with high probability in sub-linear time.
Specifically, if the query length M = O(N µ ) and the number of matches L = O(N λ ), we show that
for λ < 1 − µ all the matching positions can be determined with a probability that approaches 1 as
N → ∞ for K ≤ 61 M . More importantly our scheme has a worst-case computational complexity that

is only O max{N 1−µ log2 N, N µ+λ log N } , which means we can recover all the matching positions
in sub-linear time for λ < 1 − µ. This is a substantial improvement over the best known computational

complexity of O N 1−0.359µ for recovering one matching position by Andoni et al. [1]. Further, the
number of Fourier transform coefficients that need to be computed, stored and accessed, i.e., the sketching

complexity of this algorithm is only O N 1−µ log N . Several extensions of the main theme are also
discussed.

I. I NTRODUCTION AND P ROBLEM S TATEMENT
We consider the substring/pattern matching problem, which has been studied extensively in theoretical
computer science. In this problem, a signal ~x := (x[0], x[1], . . . , x[N − 1]) of length N symbols representing a string, library, or database is available. The objective is to answer queries regarding whether a
given string ~y := (y[0], y[1], . . . , y[M − 1]) of length M is a substring of ~x. We are especially interested
in the case where a sketch of ~x can be computed offline and stored, and where the one-time computational
complexity of creating the sketch can be amortized over repeated queries or ignored. Moreover, we focus

April 27, 2017

DRAFT

2

on the random setting in which the x[i]’s form a sequence of independent and identically distributed
(i.i.d.) random variables, each taking values in A ⊂ R. We begin our analysis by restricting our attention
to the binary case where A := {±1}.1 Within this context, we consider the following two settings:
Exact Pattern Matching: In the exact pattern matching problem, a substring of length M of ~x, namely
~y := ~x[τ : τ + M − 1], is obtained by taking M consecutive symbols from ~x and is presented as a query.

This pattern may appear within ~x in up to L different locations τ1 , τ2 , . . . , τL and the task is to determine
all the locations τi , ∀i ∈ [1 : L]. We consider the probabilistic version where our objective is to recover
the matching locations with a probability that approaches 1 as M, N → ∞.
Approximate Pattern Matching: In approximate pattern matching, ~y is a noisy version of a substring,
i.e., ~y = ~x[τ : τ + L − 1]

~b, where ~b is a noise sequence with b[i] ∈ {±1} such that dH (~y, ~x[τ :

τ + M − 1]) ≤ K . Function dH denotes the Hamming distance, and

represents component-wise

multiplication. The objective is to determine all locations τi such that dH (~y, ~x[τi : τi + M − 1]) ≤ K
with a probability that approaches 1 as K, M and N → ∞.
We evaluate our proposed algorithm according to two metrics - (i) the space required to store the sketch
of ~x, which we refer to as sketching complexity, and (ii) the computational complexity in answering the
query.
These problems are relevant to many applications including text matching, audio/image matching, DNA
matching in genomics, metabolomics, radio astronomy, searching for signatures of events within large
databases, etc. The proposed techniques are particularly relevant now due to the interest in applications
involving huge volumes of data. Our proposed approach is most useful in the following situations. (i)
The string ~x is available before querying and one time computations such as computing a sketch of ~x
can be performed offline and stored, and the complexity of computing the sketch of ~x can be ignored.
Then, when queries in the form of ~y appear, one would like to decrease the computational complexity
in searching for the string ~y. (ii) The string ~x is not centrally available, but parts of the string are sensed
by different data collecting nodes distributively and communicated to a central server. A search query is
presented to the server; and this server must decide whether the string appears in the data sensed by one
or more of the distributed nodes and, if present, it must also identify when the queried string appeared.
In this latter case, we wish to minimize the amount of data communicated by the nodes to the server
and the computational complexity in searching for the string. The proposed formulation is most useful
when the query ~y is not a pattern that can be predicted and, therefore, creating a lookup table to quickly
identify commonly-occurring patterns will not be effective.
1

Extensions to other alphabets A ⊂ R are straightforward. Extension to the non-i.i.d. case is also possible.

April 27, 2017

DRAFT

3

A naive approach to searching for substring ~y in ~x is to first compute the cross-correlation between ~x
and ~y, which we denoted by ~r = [r[0], r[1], . . . , r[N − 1]] with
r[m] = (~x ∗ ~y)[m] ,

M
X

x[m + i − 1]y[i],

(1)

i=1

and, subsequently, choosing the index m that maximizes r[m]. This strategy uses the N samples of ~x
and has a super-linear computational complexity of M N = O(N 1+µ ). A computationally more efficient
approach uses the fact that ~r can be computed by taking the inverse Fourier transform of the product
of the Fourier transforms of ~x and ~y∗ [−n], where ~y∗ [−n] is the conjugate of the time reversed version
of ~y. Such an approach still uses all the N samples of ~x, but reduces the computational complexity to
O(N log N ). Note that even though the Fourier transform of ~x can be precomputed, the N -point Fourier

transform of ~y still needs to be computed online resulting in the O(N log N ) computational complexity.
Both the exact pattern matching problem and the approximate pattern matching problem have been
extensively studied in computer science. Three recent articles [1], [2] and [3] provide a brief summary of
existing contributions. For the exact matching problem, the Rabin-Karp algorithm solves a more general
problem of finding the occurrence of a set of query strings. However, the algorithm has a computational
complexity that is at least linear in N . Boyer and Moore presented an algorithm in [4] for finding the first
occurrence of the match (only τ1 ) that has an expected complexity of O(N/M log N ) = O(N 1−µ log N ),
whereas the worst case complexity (depending on τ1 ’s) can be O(N log N ). For large M , the algorithm
indeed has an average complexity that is sub-linear in N . More recently, it has been shown that techniques
based on the Burrows-Wheeler transform can be used to solve the exact matching problem with sublinear time complexity [5] using a storage space of O(N ) bits. This problem is well studied under the
read alignment setting by the Bioinformatics community [6], [7]. For small |A|, it has the best known
complexity; however, the complexity increases with |A|. Further, extensions to approximate matching
setting [8] has a complexity that increases exponentially in K and, hence, appear to be infeasible for
K = O(M ). The Boyer and Moore algorithm has been generalized to the approximate pattern matching

problem in [9] with an average case complexity of O(N K/M log N ), which provides a sub-linear time
algorithm only when K  M . In [1], Andoni, Hassanieh, Indyk and Katabi have given the first sub-linear
time algorithm with a complexity of O(N/M 0.359 ) for K = O(M ).
II. O UR M AIN R ESULTS AND R ELATION TO P RIOR W ORK
N
Assume that a sketch of ~x of size O( M
log N ) can be computed and stored. Then for the exact pattern

matching and approximate pattern matching (with K = ηM ) problems, with the number of matches L
scaling as O(N λ ), we show an algorithm that has

April 27, 2017

DRAFT

4

•


N
a sketching function for ~y that computes O( M
log N ) = O N 1−µ log N samples

a computational complexity of O max{N 1−µ log2 N, N µ+λ log N }

•

a decoder that recovers all the L matching positions with a failure probability that approaches zero

•

asymptotically in N

N
When L < O M
(i.e. λ < 1 − µ), which is typically the interesting case, our algorithm has a sub-linear
time and space complexity.
There are important differences between our paper and [10], [1], [4], [2]. First and foremost, the
algorithms used for pattern matching are entirely different. While their algorithms are combinatorial in
nature, our algorithm is algebraic and uses signal processing and coding theoretic primitives. Secondly,
the system model considered in our paper differs from the model in [10], [1], [4], [2] in that we allow
for preprocessing or creating a sketch of the data ~x. Our algorithm exploits this fact and results in a
computational complexity O(N/M ) which is better than that in [1] for the approximate pattern matching
problem. Finally, we also consider the problem of finding all matches of the pattern ~y instead of looking
for only one match.
Our paper is inspired by and builds on two recent works by Hassanieh et al. in [10] and Pawar and
Ramchandran [11]. In [10], Hassanieh et al., considered the correlation function computation problem
for a Global Positioning System (GPS) receiver and exploited the fact that the cross-correlation vector ~r
is a very sparse signal in the time domain and, hence, the Fourier transform of ~r need not be evaluated
at all the N points. In the GPS application, which was the focus of [10], the query ~y corresponds to the
received signal from the satellites and, hence, the length of the query was at least N . As a result, the
computational complexity is still O(N log N ) (still linear in N ) and only the constants were improved in
relation to the approach of computing the entire Fourier transform. In a later paper by Andoni et al., [1],
a sub-linear time algorithm for shift finding in GPS is presented; however, this algorithm is completely
combinatorial and eschews algebraic techniques such as FFT-based techniques.
In [11], Pawar and Ramchandran present an algorithm based on aliasing and the peeling decoder
for computing the Fourier transform of a signal with noisy observations for the case when the Fourier
transform is known to be sparse. This algorithm has a complexity of O(N log N ) and they do not
consider the pattern matching problem. Our algorithm is similar to that of Pawar and Ramchandran’s
algorithm with one important distinction. We modify their algorithm to exploit the fact that the peak
of the correlation function of interest is always positive. This modification allows us to compute the
Sparse Inverse Discrete Fourier Transform (SIDFT) with sub-linear time complexity of O(N 1−α log N ),
0 < α ≤ 1 . One of the main contributions of this paper is to show that signal processing and coding

theoretic primitives, i.e., Pawar and Ramchandran’s algorithm with some modifications can be used to

April 27, 2017

DRAFT

5

solve the pattern matching algorithm in sub-linear time.
III. N OTATIONS
The table below introduces the notations we adopt throughout this paper. We denote signals and vectors
Symbol

Meaning

N

Size of the string or database in symbols

M

Length of the query in symbols

L

Number of matches

µ

Smallest 0 < µ < 1 such that M = O(N µ )

λ

Smallest 0 < λ < 1 such that L = ON λ )

K

maxτ dH (~x[τ : τ + M − 1], ~y)

η

K
M

d

Number of stages in the FFAST algorithm

fi ≈ N

α

Length of smaller point IDFT at each stage-i

gi = N/fi

Sub-sampling factor in Fourier domain for stage-i

B

Number of shifts (or branches) in each stage
γ

Number of blocks (for parallel processing)

1−γ

Length of one block (for parallel processing)

G=N
Ñ = N

TABLE I: Parameters and various quantities involved in describing the algorithm

using the standard vector notation of arrow over the letter, time domain signals using lowercase letters
and the frequency domain signals using uppercase letters. For example ~x = (x[0], x[1], . . . , x[N − 1])
~ = FN {~x} denotes the
denotes a time domain signal with ith time component denoted by x[i], and X
N -point Fourier coefficients of ~x. We denote matrices using boldface upper case letters. We denote the

set {0, 1, 2 · · · , N − 1} by [N ].
IV. D ESCRIPTION O F T HE A LGORITHM
In this section, given the input string ~x and the query string ~y, we describe our algorithm that finds
the matching positions T := {τ1 , τ2 , · · · τL } with sample and time complexities that are sub-linear in
N . The main idea exploits the fact that the correlation vector ~r is sparse (upto some noise terms) with

dominant peaks at L matching positions denoted by T and noise components at N − L positions where
the strings do not match.
Consider the correlation signal ~r in the case of exact matching:

 M, if m ∈ T
r[m] =
 n , m ∈ [N ] − T
m

April 27, 2017

(2)

DRAFT

6

RSIDFT

Sketch of ~y computation

~y0 = ~y∗[−n]
~y0 ∈ RN

Folding

0
~yi,j
∈ Rfi

fi-FFT

Y~ 0[S]

~ i,j
R
~ i,j ]
= R[S

fi-IDFT

~zi,j

&

Decoder

{τ1, τ2 · · · , τL}

binning

X[S]
Sketch of ~x

Fig. 1: Schematic of the proposed scheme using sparse Fourier transform computation.
where nm is the noise component that is induced due to correlation of two i.i.d. sequence of random
variables each taking values from A := {+1, −1}. The sparse vector ~r can be computed indirectly using
Fourier transform approach as shown below:
−1
~r = FN
{FN {~x}
III

FN {~y0 }}

I

(3)

II

−1
{·} refer to N -point discrete Fourier transform and its inverse respectively,
where FN {·} and FN

is

the point-wise multiplication operation and y 0 [n] = y ∗ [−n]. Fig. 1 presents a notional schematic of our
Algorithm. As evident from Eq. (3), our algorithm for computing ~r consists of three stages:
•

~ = FN {~x} of ~x
Computing the sketch X

•

~ 0 = FN {~y0 } of ~y
Computing the sketch Y

•

~ =X
~
Computing the IDFT of R

~ 0 given X
~ and Y
~0
Y

A. Sparse Inverse Discrete Fourier Transform
In this section we present Robust Sparse Inverse Discrete Fourier Transform(RSIDFT) scheme that
exploits sparsity in the cross-correlation signal ~r and efficiently recovers its L dominant coefficients.
The architecture of RSIDFT is similar to that of the FFAST scheme proposed in [11], but the decoding
algorithm has some modifications to handle the noise model induced in this problem. We will see in
~ and Y
~ 0 are computed efficiently but for this section we will focus only
Sec. IV-B how the sketches X
~ and Y
~ 0.
on the recovery of the sparse coefficients in ~r given X
~ =X
~
Consider the RSIDFT framework shown in Fig. 2. Let R

~ 0 be the DFT of the cross-correlation
Y

signal of ~x and ~y. We begin by factoring N into d relatively prime factors {f1 , f2 , . . . , fd }, where d
is a parameter in the algorithm. The design scheme for choosing fi ’s for various values of µ such
that fi divides N and fi = N α + O(1)∀ i ∈ [d] are given in Sec. IV-A2. The RSIDFT algorithm
consists of d-stages with each stage corresponding to a sub-sampling factor of

N
fi .

In each stage, there

are B = O(log N ) branches with shifts from the set {s1 , s2 , · · · sB }, where s1 = 0 in the first branch
and the rest are chosen uniformly at random from [N ].

April 27, 2017

DRAFT

7

~ that is sparse in time
Fig. 2: RSIDFT Framework to compute inverse Fourier Transform of a signal R

domain.
~ , in branch j of ith stage of RSDIFT, referred to as branch (i, j) for simplicity,
Given the input R
~ at
RSIDFT sub-samples the signal R
Si,j := {sj , sj + gi , sj + 2gi , · · · sj + (fi − 1)gi },

where gi :=

N
fi

i ∈ [d], j ∈ [B]

(4)

~ i,j := R[S
~ i,j ]. The sub-sampling operation is followed by a fi -point IDFT in
to obtain R

branch (i,j) of stage-i to obtain ~ri,j . Notice that ~ri,j is an aliased version of ~r due to the property that
sub-sampling in Fourier domain is equivalent to aliasing in time domain.
Let ~zi,k ∈ RB , for k ∈ [fi ] + 1, be the k th binned observation vector of stage-i formed by stacking
~ri,j [k], j ∈ [B], together as a vector i.e.
h
iT
~zi,k = ri,1 [k], ri,2 [k], · · · , ri,B [k]

Note that this gives us a total of fi binned observation vectors in each stage-i. Using the properties of
Fourier transform, we can write the relationship between the observation vectors ~zi,k at bin (i, k) and

April 27, 2017

DRAFT

8

sparse vector to be estimated ~r as:
h
iT
~zi,k = Wi,k × r[k], r[k + fi ], · · · , r[k + (gi − 1)fi ]

(5)

where we refer to Wi,k as the sensing matrix at bin (i, k) and is defined as
h
i
Wi,k = w
~ k, w
~ k+fi , . . . , w
~ k+(gi −1)fi

(6)

h j2πks
iT
j2πksB
j2πks2
1
and w
~k = e N , e N , ··· , e N

 2πs k
i
wik := ej N

r[0]
~z1,1
r[1]
~z1,2
r[2]
~z1,3
r[3]
~z2,1
r[4]
~z2,2
r[5]


r1,1[1]
=
r [1]
 1,2
r1,1[2]
=
r [2]
 1,2
r1,1[3]
=
r [3]
 1,2
r2,1[1]
=
r [1]
 2,2
r2,1[2]
=
r2,2[2]

= r[0] + r[3]
= r[0]w20 + r[3]w23
= r[1] + r[4]
= r[1]w21 + r[4]w24
= r[2] + r[5]
= r[2]w22 + r[5]w25








= r[0] + r[2] + r[4]



= r[0]w20 + r[2]w22 + r[4]w24



= r[1] + r[3] + r[5]



= r[1]w21 + r[3]w23 + r[5]w25



Fig. 3: Example of a Tanner graph formed in a RSIDFT framework with system parameters being N = 6,
f1 = 2, f2 = 3 (i.e., d = 2) and B = 2. The variable nodes (gray circles) represent the cross-correlation

vector ~r and the bin nodes (white squares) represent the binned observation vector ~zi,k . The figure
illustrates the relationship between ~zi,k and ~r.

We represent the relation between the set of observation vectors {~zi,k , i ∈ [1 : d], k ∈ [fi ]} and ~r using
a Tanner graph, an example of which is shown in Fig. 3. The nodes on the left, which we refer to as
variable nodes, represent the N elements of vector ~r. Similarly the nodes on the right, which we refer to
P
as bin nodes, represent the i≤d fi sub-sensing signals. We will now describe the decoding algorithm
which takes the set of observation vectors {~zi,k , i ∈ [1 : d], k ∈ [fi ]}, each of length B , at dfi bins as
input and estimates the L-sparse ~r.

April 27, 2017

DRAFT

9

1) Decoder: Observe from the Tanner graph that the degree of each variable node is d and that of each
bin node at stage i is gi . A variable node is referred to as non-zero if it corresponds to a matching position
and as zero if it corresponds to a non-matching position. Note that even though the cross-correlation vector
value corresponding to a non-matching position is not exactly zero but some negligible noise value we
refer to them as zero variable nodes for simplicity. We refer to a bin node as zero-ton (or Hz ) if the
number of non-zero variable nodes connected to the bin node is zero. The singleton (Hs ), double-ton
(Hd ) and multi-ton (Hm ) bin nodes are defined similarly where the number of non-zero variable nodes
connected are one, two and greater than two, respectively. The peeling decoder has the following three
steps in the decoding process.
Bin Classification: In this step a bin node is classified either as a zero-ton or a singleton or a
multi-ton. At bin (i, j) the classification is done based on comparing the first observation zi,j [1], which
corresponds to zero shift, with a predefined threshold. For zi,j [1] = z , the classification hypothesis at bin
b i,j , can be written as follows:
(i, j), H

b i,j =
H



Hz






Hs



H

 d



H

m

z/M < γ1
γ1 < z/M < γ2

(7)

γ2 < z/M < γ3
z/M > γ3

3−4η 5−6η
where (γ1 , γ2 , γ3 ) = ( 1−2η
2 , 2 , 2 ). Note that for the case of exact matching η = K/M = 0.

Singleton decoding: If a bin node (i, j) is classified as a singleton in the bin classification step, we
need to identify the position of the non-zero variable node connected to it. This is done by correlating
the observation vector ~zi,j with each column of the sensing matrix Wi,j := [~
wj , w
~ j+fi , . . . , w
~ j+(gi −1)fi ]
and choosing the index that maximizes the correlation value.
k̂ = arg max ~z†i,j w
~k
k∈{j+lfi }

where † denotes the conjugate transpose. The value of the variable node connected to the singleton bin
is decoded as:
r̂[k̂] = M (1 − η).

Note that for the case of exact matching we know the value to be exactly equal to M . But in the
case of approximate matching, the actual value of r[k] ∈ {M (1 − 2η), . . . , M − 1, M } and our estimate
r̂[k̂] = M (1 − η) is only approximate. But this suffices for recovering the positions of matches i.e., the

indices of the sparse coefficients in ~r.

April 27, 2017

DRAFT

10

Peeling Process: The peeling based decoder we employ consists of finding a singleton bin, then
identifying the single non-zero variable node connected to the bin, decoding its value and removing
(peeling off) it’s contribution from all the bin nodes connected to that variable node. The main idea
behind this decoding scheme is that (for appropriately chosen parameters), at each iteration, peeling a
singleton node off will induce at least one more singleton bin and the process of peeling off can be
repeated iteratively. Although the main idea is similar for exact matching and the approximate matching
scenarios, there are some subtle differences in their implementation.
Exact Matching: In the case of exact matching, we remove the decoded variable node’s contribution from
all the bin nodes it is connected to.
Approximate Matching: In this case, similar to the approach in [12], we remove the decoded variable
node’s contribution only from bins that are originally a singleton or a double-ton. We do not alter the
bins which are classified to be multi-tons with degree more than two.
Note: The differences in peeling implementation for exact and approximate matching cases is because
unlike exact matching the approximate matching may result in non-zero error (e1 ) between the actual
correlation value r[k] and the estimate r̂[k̂], i.e. e1 = |r[k] − r̂[k̂]| with 0 ≤ e1 ≤ ηM . This may lead
to error propagation if we use the same decoding scheme as in exact matching. Hence to overcome this
problem we impose a constraint on the type of bin nodes participating in the peeling process.
We present the overall recovery algorithm, which comprises of bin classification, singleton decoding
and peeling process, in the Algorithm.1 pseudo-code. Note that N(k) denote the neighborhood for variable
node k i.e., the set of bins connected to k th variable node.
2) Choosing fi and α for various µ: For a given value of µ, we will describe how to choose the
Q
parameters d and fi . Find a factorization for signal length N = d−1
i=0 Pi such that the set of integers
{P0 , P1 , . . . , Pd−1 } are pairwise co-prime and all the Pi are approximately equal. More precisely, let
Pi = F + O(1) ∀i for some value F. We can add zeros at the end of the vector ~x and increase the length

of the vector until we are able to find a factorization that satisfies this property.
•

For µ < 0.5: Choose fi = N/Pi .
1
Exact Matching: Find d ∈ N\{1, 2} such that µ ∈ ( d1 , d−1
]
1
Approximate Matching: If µ ∈ ( 18 , 21 ), choose d = 8. Else find d ≥ 8 such that µ ∈ ( d1 , d−1
)

•

For µ > 0.5: Choose fi = Pi
Exact Matching: Find d ∈ N\{1, 2} such that µ ∈ (1 −

1
d−1 , 1

− d1 ]

1
Approximate Matching: If µ ∈ ( 12 , 18 ), choose d = 8. Else find d ≥ 8 such that µ ∈ (1 − d−1
, 1 − d1 )

April 27, 2017

DRAFT

11

Algorithm 1 Peeling based recovery algorithm
b i,j for i ∈ [d], j ∈ [fi ]. (See Eqn. (7))
Compute H
b i,j = Hs , do
while ∃ i, j : H
(k̂, r̂[k̂]) =Singleton-Decoder(~zi,j )

Assign r̂[k̂] to k̂ th variable node
for (i0 , j0 ) ∈ N(k̂) do
if Exact Matching then
~zi0 ,j0 ← ~zi0 ,j0 − r̂[k̂]~
wk̂

else
~zi0 ,j0 ← ~zi0 ,j0 − r̂[k̂]~
wk̂

only if

b i ,j = Hs or Hd
H
0 0

b i ,j
Re-do the bin classification for (i0 , j0 ) and compute H
0 0
Algorithm 2 Singleton-Decoder
Input: ~zi,j
Output: (k̂, r̂[k̂])
Estimate singleton index to be k̂ = arg max ~z†i,j w
~k
k∈{j+lfi }

Estimate the value to be:
r̂[k̂] =



M

Exact Matching case


M − K

Approximate Matching case

Thus, for both the exact and approximate matching cases, for any 0 < µ < 1, we choose the downsampling factors fi to be approximately equal to N α where α > 1 − µ.
3) Distributed processing framework: Given a database(or string) of length N , we divide the database
into G = N γ blocks each of length Ñ = N/G. Now each block can be processed independently
(in parallel) using the RSIDFT framework with the new database length reduced from N to Ñ . This
distributed framework has the following advantages
•

Firstly, this enables parallel computing and hence can be distributed across different workstations.

•

Improves the sample and computational complexity by a constant factor.

•

Sketch of the database needs to be computed only for a smaller block length and hence requires
computation of only a shorter Ñ point FFT.

•

Helps overcome implementation issues with memory and precision as the scale of the problem is

April 27, 2017

DRAFT

12

reduced.
~ and Y
~
B. Sketches of X
~
~
As we have already seen in Sec. IV-A the RSIDFT framework requires the values of R(=
X

~)
Y

~ and Y
~ only at the indices in set S of cardinality dBfi .
only at indices S or in other words we need X
~
We assume that the sketch of ~x, X[S]
= {X[i], i ∈ S} is pre-computed and stored in a database.
~ 0 [Si,j ], i ∈ [d], j ∈ [B]} needs to be
Computing the sketch of ~y: : For every new query ~y, only {Y

computed where the subsets Si,j , defined in Eq. (4), are
Si,j := { sj + kgi : k ∈ [fi ]},

i ∈ [d], j ∈ [B]

(8)

~ 0 and the required
of cardinality fi . Naively, the FFT algorithm can be used to compute N -pt DFT of Y

subset of coefficients can be taken but this is inefficient and would be of O(N log N ) complexity. Instead,
~ 0 [Si,j ] is Y
~ 0 shifted by sj and sub-sampled by a factor of gi . Thus for a given (i, j)
we observe that Y

this corresponds to, in time domain, a point-wise multiplication by [1, wsj , ws2j , . . . , wsNj −1 ] followed by
folding the signal into gi (=

N
fi )

signals each of length fi and adding them up resulting in a single

0 . Formally the folding operation can be described as follows:
length-fi signal denoted by ~yi,j
0
~yi,j
=

gX
i −1

~y0 [mfi : (m + 1)fi − 1]

i
i +1
i −1
[wsmf
, wsmf
, . . . , ws(m+1)f
],
j
j
j

m=0

where, wsj = e−

j2πsj
N

0 produces Y
~ 0 [Si,j ] i.e.
. Taking fi -point DFT of ~yi,j
0
~ 0 [Si,j ] = Ff {~yi,j
Y
}
i

which is what we need in branch (i, j). To obtain all the samples in S required for the RSDIFT framework,
the folding technique followed by a DFT needs to be carried out for each (i, j), for i ∈ [d], j ∈ [B], a
total of dB times N α -point DFT.
V. P ERFORMANCE A NALYSIS
In this section, we will analyze the overall probability of error involved in finding the correct matching
positions. This can be done by analyzing the following three error events independently and then using
a union bound to bound the total probability of error.
•

E1 -Bin Classification: Event that a bin is wrongly classified.

•

E2 -Position Identification: Given a bin is correctly identified as a singleton, event that the position

of singleton is identified incorrectly.

April 27, 2017

DRAFT

13

•

E3 -Peeling Process: Given the classification of all the bins and the position identification of singletons

in each iteration is accurate, event that the peeling process fails to recover the L significant correlation
coefficients.
A. Bin Classification
Lemma 1. The probability of bin classification error at any bin (i, j) can be upper bounded by
P[E1 ] ≤ 6e−

N µ+α−1 (1−6η)2
16

Proof.
b i,j = Hz ] +
P[E1 ] = ≤ P[E1 |H
≤ e−

N µ−α (1−2η)2
8

≤ 6e−

+ 2e−

b i,j = Hs ]
P[E1 |H
N µ−α (1−4η)2
16

+ 2e−

+

b i,j = Hd ∪ Hm ]
P[E1 |H

N µ+α−1 (1−6η)2
16

+ e−

N µ+α−1 (1−6η)2
16

N µ+α−1 (1−6η)2
16

where the inequalities in the third line are due to Lemmas 8, 9, 10 and 11 respectively provided in
Appendix A.
B. Position Identification
Lemma 2. Given that a bin (i, j) is correctly classified as a singleton, the probability of error in
identifying the position of the non-zero variable node can be upper bounded by




N µ+α−1 (1 − 2η)2 (c21 − 1)
N µ+α−1 (c1 (1 − 2η) − 1)2
P[E2 ] ≤ exp −
+ exp −
8(c21 + 1)
8(1 + c21 )

Proof. The detailed proof is provided in Appendix A.
C. Peeling Process
To analyze the peeling process alone independently, we refer to a oracle based peeling decoder which
has the accurate classification of all the bins and can accurately identify the position of the singleton
given a singleton bin at any iteration. In other words, oracle based peeling decoder is the peeling part of
our decoding scheme but with the assumption that the bin classification and position identification are
carried out without any error.

April 27, 2017

DRAFT

14

d

2

3

4

5

6

7

8

δ

1.000

0.4073

0.3237

0.2850

0.2616

0.2456

0.2336

dδ

2.000

1.2219

1.2948

1.4250

1.5696

1.7192

1.8688

TABLE II: Constants for various error floor values
Lemma 3 (Exact Matching). For the exact matching case, choose F d−1 = δN α where δ is chosen as
given in Table. II. Then the oracle based peeling decoder:
•

successfully uncovers all the L matching positions if L = Ω(N α ) and L ≤ N α , with probability at
1

least 1 − O(1/N d )
•

successfully uncovers all the L matching positions, if L = o(N α ), with probability at least 1 −
2

e−βε1 N

α/(4l+1)

for some constants β, ε1 > 0 and l > 0.

Proof. We borrow this result from Pawar and Ramchandran’s [11]. Although our RSDIFT framework
and their robust-FFAST scheme have three main differences:
•

We are computing smaller IDFT’s to recover a sparse bigger IDFT whereas in [11] the same is true
for DFT instead of IDFT.

•

Our problem model is such that the sparse components of the signal space has only positive amplitude
and thus our bin processing part (bin classification and position identification) is different when
compared to [11].

•

The sparsity of the signal L to be recovered is exactly known in the case of [11] whereas we have
no information about L not even the order with which the quantity scales in N .

Irrespective of these differences, the Tanner graph representation of the framework and the peeling part
of the decoder are identical to that of the robust-FFAST scheme. And thus the limit of the oracle based
peeling decoder for our scheme is identical to that in the robust-FFAST scheme [11]. With respect to
the third difference, in robust-FFAST scheme the authors choose F d−1 = δk where k is the sparsity of
the signal (which is assumed to be known) and show the first assertion of the lemma. They also showed
that upto a constant fraction (1 − ε) of k -variables node can be recovered with probability of failure that
decays exponentially in N . In our case, since L = o(N α ), this result translates to recovering all the L
non-zero variable nodes with an exponentially decaying failure probability.
In any iteration, given a singleton bin, the peeling process, in the case of approximate matching, runs
the Singleton-Decoder algorithm on the bin only if it was either originally a singleton or originally
a double-ton with one of the variable nodes being peeled off already. This is in contrast to the exact
matching case where the peeling decoder runs the Singleton-Decoder on the bin irrespective of it’s original

April 27, 2017

DRAFT

15

degree. Hence we need to analyze the oracle based peeling decoder for the approximate matching case
separately compared to the exact matching case.
Lemma 4 (Approximate Matching). For the approximate matching case, choose parameter d ≥ 8 as
described in Sec. IV-A2 and F d−1 = 0.7663N α . Then the oracle based peeling decoder:
•

successfully uncovers all but a small fraction ε = 10−3 of the L matching positions, if L = Ω(N α )
and L ≤ N α with a failure probability that decays exponentially in N

•

successfully uncovers all the L matching positions, if L = o(N α ), with probability at least 1 −
2

e−βε1 N

α/(4l+1)

for some constants β, ε1 > 0 and l > 0.

Proof. As mentioned earlier, the key difference in the approximate matching case is peeling off variable
nodes from only singleton and double-tons. An identical peeling decoder is used and analyzed in
the problem of group testing [12] by Lee, Pedarsani and Ramchandran which the authors refer to as
SAFFRON scheme. In SAFFRON, the authors claim that for a graph ensemble which has a regular
degree of d on the variable nodes and a Poisson degree distribution on the bins, this peeling decoder with
a left degree of d = 8 and a total number of bins at least equal to 6.13k recovers at least (1 − ε) fraction
of the k non-zero variable nodes with exponentially decaying probability. Note that

6.13
8

≈ 0.7663 is

approximately the number of bins per stage for d = 8. We also leverage the result from [11] that the
Tanner graph representation of the robust-FFAST (or equivaently RSDIFT framework) has a Poisson
degree distribution on the bins. Combining these two results gives us the required results.
Theorem 5 (Overall Probability of Error). The RSIDFT framework succeeds with high probability
asymptotically in N if the number of samples in each branch fi = N α + O(1) satisfies the condition
α > 1 − µ.

Proof. The overall probability of error P[Etotal ] can be bounded using an union bound on the three error
events E1 , E2 and E3 given by
P[Etotal ] ≤ P[E1 ] + P[E2 ] + P[E3 ]
Using the expressions for error probabilities from Lemmas 1, 2, 3 and 4, we can see that all the terms
vanish to zero as N → ∞ if µ + α − 1 ≥ 0, i.e. α ≥ 1 − µ.
VI. S AMPLE AND C OMPUTATIONAL C OMPLEXITY
In this section, we will analyze the sketching complexity which is the number of samples we access
from the sketch of the signal ~x stored in the database and the computational complexity as a function of
the system parameters.

April 27, 2017

DRAFT

16

A. Sample Complexity
In each branch of the RSDIFT framework we down-sample the N samples by a factor of

N
fi

to get

fi ≈ N α samples. We repeat this for a random shift in each branch for B = O(log N ) branches in each

stage thus resulting in a total of O(N α log N ) samples per block per stage. We repeat this for d =

1
1−α

such stages resulting in a total of dN α log N samples per block. So, the total number of samples is given
by
S = O (dN α log N ) = O(N 1−µ log N )

B. Computational Complexity
As described in Eq. 3, the computation of ~r involves three steps:
1) Operation - I: Since we assume that the sketch of database ~x, FN {~x}, is pre-computed, we do not
include this in computational complexity.
2) Operation - II: As described in Sec. IV-B, in each branch (i, j), we use a folding based technique
to compute the sketch of ~y, FN {~y0 } at points in the set Si,j . The folding technique involves two
steps: folding and adding (aliasing) which has a complexity of O(M ) computations , and computing
fi -point DFTs that takes O(N α log N α ) computations. So, for a total of dB branches the number

of computations in this step is given by

CII = dB ( |{z}
N µ + N α log N α )
|
{z
}
Folding

Shorter FFTs

= O(max(N 1−µ log2 N, N µ log N )).

Note: Folding and adding, for each shift, involves adding N 1−α vectors of length N α . We know
that the length of the query is M = N µ , i.e., the number of non-zero elements in ~y (zero-padded
version of the query) is M and hence we only need to compute M additions instead of length of
the vector N .
−1
3) Operation - III: Computing FN
{~x} involves two parts:

RSIFT framework and the decoder. The RSIFT framework involves computing smaller fi point
IDFTs, which takes approximately O(N α log N α ) computations in each branch. For a total of
dB branches, we get a complexity of O(dBN α log N α ). In the decoding process, the dominant

computation is from position identification. Each position identification process involves correlating
the observation vector of length B with

N
fi

≈ N 1−α column vectors, which amounts to BN 1−α

computations. There will be a maximum of dL such position identifications, which gives a complexity

April 27, 2017

DRAFT

17

ofO(dLBN 1−α ). Now, plugging in α = 1 − µ (condition for vanishing probability of error) the total
number of computations involved in this step, CI , is given by
1−α
CIII = dB (O(N α log N α ) + |L N
{z })
|
{z
}
Shorter IFFTs /block/stage

Correlations

= O(max(N 1−µ log2 N, N µ+λ log N )).

Thus, the total number of computations, C = max{CII , CIII }, is given by
C = O(max{N 1−µ log2 N, N µ+λ log N })

10−1
10−2
10−3
10−4
10−5
100 200 300 400 500 600 700 800
Sample Gain
(a) M = 105 (µ = 0.41), Ñ = 107 , G = 105

Prob of Missing a Match

Prob of Missing a Match

VII. S IMULATION R ESULTS

100
10−1
10−2
10−3
10−4
10−5

2

4

6

8 10 12
Sample Gain

14

16

(b) M = 103 (µ = 0.25), Ñ = 106 , G = 106

Fig. 4: Plots of probability of missing a match vs. sample gain for exact matching of a query of length
M from a equiprobable binary {+1,-1} sequence of length N = 1012 , divided into G blocks each of

length Ñ . The substring was simulated to repeat in L = 106 (λ = 0.5) locations uniformly at random.

Simulations were carried out to test the performance of RSIDFT framework for exact matching scenario
on a database of length N = 1012 for two different query lengths M = 105 (µ = 0.41) and M = 103
(µ = 0.25). The database was generated as a equiprobable {+1, −1} sequence of length N . A substring
of length M from the generated database is presented as a query. Also the chosen query was repeated
at L = 106 randomly chosen locations in the database.
The sample gain, defined as the ratio of N to the number of samples used from the sketch of database,
was varied and the probability of RSIDFT framework to miss a match (Pe ), as defined below, was
measured.
Pe =

April 27, 2017

# of correctly identified locations
L
DRAFT

18

18

The plots of Pe vs. sample gain, is presented in Fig. 4 for two different query lengths: M = 105 (µ = 0.41)
in Fig 4(a) and M = 103 (µ = 0.25) in Fig 4(b). As can be inferred from the plots we achieve a sample
gain of 200-300 (depending on the tolerable error probability) for the query length corresponding to
µ = 0.41 and a sample gain of 2-4 for µ = 0.25. This sample gain results from an average number

of samples per branch fi ≈ 9.25 × 107 (α = 0.66) for µ = 0.41, and fi ≈ 6.94 × 109 (α = 0.82) for
µ = 0.25. The trend in the results almost matches with the theoretical findings of α = 1 − µ. We also

notice a sharp threshold in the sample gain, below which the RSIDFT framework succeeds with very
high probability.
A PPENDIX
Lemma 6 (Hoeffding tail bound). Let X1 , X2 , . . . , Xn be a sequence of independent random variables
such that Xi has mean µi and sub-Gaussian parameter σi . Then for any δ > 0:


hX
i
δ2
Upper Tail : P
(Xi − µi ) ≥ δ ≤ exp − P 2
2 σi


hX
i
δ2
Lower Tail : P
(Xi − µi ) ≤ −δ ≤ exp − P 2
2 σi
Note that for bounded random variables Xi ∈ [a, b] the sub-Gaussian parameter is σi =
the upper tail Hoeffding bound can be simplified to
" n
#

X
P
(Xi − µi ) ≥ δ ≤ exp −
i=1

2δ 2
n(b − a)2

b−a
2

whereupon


.

(9)

Similarly the lower tail bound can be simplified.
Lemma 7 (Tail bounds for noise terms). Let us consider r[θ0 ], r[θ1 ], . . . , r[θgi −1 ] where θj = θ0 + jfi
and θj ∈
/ {τ1 , . . . , τL } is not one of the matching positions for any j . Then for any δ > 0:
Upper Tail:






2
X X
1
M
δ
x[θj + k]y[k] ≥ δ  ≤ exp −
P 
M
2gi
j∈[gi ] k∈[M ]

Lower Tail:






2
X X
1
M
δ
P 
x[θj + k]y[k] ≥ δ  ≤ exp −
M
2gi
j∈[gi ] k∈[M ]

Recall that [gi ] is used to denote the set {0, 1, . . . , gi − 1}.

April 27, 2017

DRAFT

19

Proof. Since θj is not one of the matching positions for any j , x[θj + k] 6= y[k] and more importantly
x[θj +k] ⊥
⊥ y[k] ∀j, k . This implies that x[θj +k]y[k] = ±1 with equal probability and E[x[θi +k]y[k]] = 0.

Let the set of random variables corresponding to a position θj be Sj := {x[θj + k]y[k], k ∈ [M ]}. It is
clear that the random variables in the set Sj are all independent with respect to each other due to our
i.i.d assumption on the database ~x and θj being a non-matching position. For the case of µ < α we have
M < fi for large enough N thus resulting in non-overlapping parts of ~x participating in the correlation

coefficients r[θi ] and r[θj ]. Hence it can be shown that Si ⊥
⊥ Sj ∀i, j and we can apply the bounds from
Eqn. (9) achieve the required result.
For the case of µ ≥ α, M > fi for large enough N which results in a coefficient x[j] participating in
multiple correlation coefficients r[θj ]. Therefore we pursue an alternate method of proof by defining
X
pj,l := x[θj + l]
y[l + kfi ] for l ∈ [fi ],
k∈[ M
]
f
i

where θj = θ0 + jfi . W.L.O.G we assume that fi divides M evenly although the proof can be extended
on similar lines for the case where fi does not divide M evenly. Now we can show that the required
sum
X X

x[θj + l]y[l] =

j∈[gi ] l∈[M ]

X X

pj,l .

j∈[gi ] l∈[fi ]

From the above equivalent representation of the required sum, we need the following to achieve the
required result:
•

•

q
pj,l is sub-Gaussian with parameter σi = M
fi since, from Eqn. (9),
)
(
−δ 2
P [pj,l ≤ δ] ≤ exp
2M
fi
The random variables {pj,l , ∀j, l} are independent of each other due to the i.i.d assumption on the
database

Now we can apply the tail bounds for sum of gi fi sub-Gaussian random variables each with sub-Gaussian
q
parameter M
fi from Lem. 6 and arrive at the required result.
We employ classification rules based only on the first element of the measurement vector at bin (i, j)
which can be given by
gi −1 M −1
P P


nl,k



`=0 k=0


gP

−1
i −2 M
P
M
+
nl,k
1
Z[1] =

`=0 k=0



gP
−1
i −3 M

P


nl,k
M1 + M2 +

if

H = Hz

if

H = Hs

if

H = Hd

(10)

`=0 k=0

April 27, 2017

DRAFT

20

where nl,k = x[θ` +k]y[k] and θ` ∈
/ {τ1 , τ2 , . . . , τL }. Also for the case of exact matching M1 = M2 = M
whereas in the case of approximate matching the values of M1 , M2 ∈ [M (1 − 2η) : M ].
Lemma 8 (zero-ton). Given that the bin (i, j) is a zero-ton, the classification error can be bounded by
P[E1 |Hz ] ≤ e−

N µ+α−1 (1−2η)2
8

Proof. The above expression can be derived by observing that a bin is not classified as zero-ton if
Z[1]
M

≥

1−2η
2 .

Let us denote the probability of this event as pz1 which can be bounded as:


Z[1]
1 − 2η
≥
pz1 =P
M
2
−

≤e

≈ e−

M gi (1−2η)2
8g 2
i
N µ+α−1 (1−2η)2
8

where the second bound is due to Eqn. (10) and Lemma. 7. The approximation in the third line is from
our design that all the gi are chosen such that gi ≈ N 1−α and M = N µ .
Lemma 9 (singleton). Given that the bin (i, j) is a singleton, the classification error can be bounded by
P[E1 |Hs ] ≤ 2e−

N µ+α−1 (1−4η)2
16

Proof. We observe that a bin is not classified as singleton if

Z[1]
M

≤

1−2η
2

or

Z[1]
M

≥

3−4η
2 .

Let us denote

the probability of the two events as ps1 and ps2 respectively which can be bounded as:
" gi −2 M −1
#
1 X X
1 − 2η M1
ps1 = P
nl,k ≤
−
M
2
M
`=0 k=0
#
" gi −2 M −1
1 − 2η
1 X X
nl,k ≤ −
≤P
M
2
`=0 k=0

−

≤e

≈ e−

M gi (1−2η)2
16g 2
i
N µ+α−1 (1−2η)2
16

where we used lower tail of Lemma 6 and gi ≈ N 1−α and the lower bound on

M1
M

≥ (1 − 2η). Similarly

ps2 can be upper bounded by:
gi −2 M −1
1 X X
3 − 4η M1
=P
nl,k ≥
−
M
2
M
`=0 k=0
"
#
gi −2 M −1
1 X X
1 − 4η
≤P
nl,k ≥ −
M gi
2gi

"

ps2

#

`=0 k=0

−N

≈e

April 27, 2017

µ+α−1 (1−4η)2
8

DRAFT

21

Thus the overall probability of error for classifying a singleton can be obtained by combining ps1 and
ps2 .

Lemma 10 (double-ton). Given that the bin (i, j) is a double-ton, the classification error can be bounded
by
P[E1 |Hd ] ≤ 2e−

N µ+α−1 (1−6η)2
16

Proof. We observe that a bin is not classified as double-ton if

Z[1]
M

≤

3−4η
2

or

Z[1]
M

≥

5−6η
2 .

Let us denote

the probability of these two events as pd1 and pd2 respectively which can be bounded similar to Lemma
9.
gi −3 M −1
1 X X
3 − 4η M1 + M2
=P
−
nl,k ≤
M
2
M
`=0 k=0
#
" gi −3 M −1
1 − 4η
1 X X
nl,k ≤ −
≤P
M
2

"

pd1

#

`=0 k=0

≤e

−

≈ e−

M (gi −2)(1−4η)2
16(gi −2)2
N µ+α−1 (1−4η)2
16

.

Similarly pd2 can be bounded as
gi −3 M −1
5 − 6η M1 + M2
1 X X
nl,k ≥
−
=P
M
2
M
`=0 k=0
" gi −3 M −1
#
1 X X
1 − 6η
≤P
nl,k ≥
M
2

"

pd2

#

`=0 k=0

−

≤e

≈ e−

M (gi −2)(1−6η)2
8(gi −2)2
N µ+α−1 (1−6η)2
8

where we use the lower bounds M1 , M2 ≤ M .
Lemma 11 (multi-ton). Given that the bin (i, j) is a multi-ton, the classification error can be bounded
by
P[E1 |Hm ] ≤ e−

April 27, 2017

N µ+α−1 (1−6η)2
16

DRAFT

22

Proof. We observe that a bin is not classified as multi-ton if

Z[1]
M

≤

5−6η
2 .

Let us denote the probability

of this event as pm1 which can be bounded as:
" gi −3 M −1
#
1 X X
5 − 6η Mm
pm1 = P
nl,k ≤
−
M
2
M
`=0 k=0
" gi −m M −1
#
1 X X
1 − 6η
≤P
nl,k ≤ −
M
2
`=0 k=0

−

≤e

M (gi −m)(1−4η)2
16(gi −m)2

− M (1−6η)
16n

≤e

2

i

≈ e−

N µ+α−1 (1−6η)2
16

.

We will analyze the singleton identification in two separate cases:
•

E21 : Event where the position is identified incorrectly when the bin is classified correctly a singleton

•

E22 : In the case of approximate matching, event where the position is identified incorrectly when

the bin is originally a double-ton and one of the non-zero variable nodes has already been peeled
off
Definition 12 (Mutual Incoherence). The mutual incoherence µmax (W) of a matrix W = [~
w1 w
~2 ···w
~i···w
~ N]
is defined as

µmax (W) , max
∀i6=j

|~
wi† w
~ j|
||~
wi ||.||~
wj ||

Lemma 13 (Mutual Incoherence Bound for sub-sampled IDFT matrix [[11],Proposition A.1). The mutual
incoherence µmax (Wi,k ) of the sensing matrix Wi,k (defined in Eq. 6), with B shifts, is upper bounded
by
r
µmax < 2

log(5N )
B

Proof. The proof follows similar lines as the proof for Lemma V.3. in [11].
Lemma 14. For some constant c1 ∈ R and the choice of B = 4c21 log 5N , the probability of error in
identifying the position of a singleton at any bin (i, j) can be upper bounded by


N µ+α−1 (1 − 2η)2 (c21 − 1)
P[E21 ] ≤ exp −
8(c21 + 1)

April 27, 2017

DRAFT

23

Proof. Let jp be the variable node participating in the singleton (i, j). Then the observation vector ~zi,j
is given by


h
z i,j = w
~ j2 , · · ·
~ j1 , w

= r[jp ] w
~ jp +

X

w
~ jp , · · ·



 n1 
 . 
 .. 




r[jp ]


i 
.. 

w
~ j gi ×  . 



 n 
 j 
 . 
 . 
 . 


ngi

nk w
~ jk

k6=p

~ jk = w
~ jk as defined in Eq. and
where for convenience we use a simpler notation jk = j + (k − 1) N
fi , w
M
−1
P
nl =
x[θ` + k]y[k] as defined in Eq. (10).
k=0

The estimated position p̂ is given by
p̂ = arg max
l

w
~ j†l z i,j
B

(11)

where † denotes the conjugate transpose of the vector. Also note that ||~
wjk || = B for any j and k . From
Eq. (11) we observe that the position is wrongly identified when ∃p0 such that
r[jp ] †
1 X
1 X
~ jk ≤
nk w
~ j†p w
w
~ jp0 w
~ jp + n p 0 +
nk w
~ j†p0 w
~ jk
r[jp ] +
B
B
B
k6=p
k6=p,p0
!
X
w
~ j†p0 w
~ jp
↔
αk nk + βnp0 ≥ r[jp ] 1 −
≥ M (1 − 2η)(1 − µmax )
B
0
k6=p,p

where αk and β are constants and can be shown to be in the range αk ∈ [−2µmax , 2µmax ] and β ∈
[1 − µmax , 1 + µmax ]. Now using the bound given Chernoff Lemma in Lem. 7 we obtain


2M (1 − 2η)2 (1 − µmax )2
P[E21 ] ≤ exp −
16gi µ2max + 4(1 + µmax )2


2M (1 − 2η)2 (1 − µmax )2
≤ exp −
16(gi µ2max + 1)


2M (1 − 2η)2 (c1 − 1)2
≤ exp −
16(gi + c21 )


N µ+α−1 (1 − 2η)2 (c21 − 1)
≈ exp −
8(c21 + 1)

April 27, 2017

DRAFT

24

The second inequality follows by the definition of µmax ≤ 1. We choose B = 4c21 log 5N , and substituting
q
µmax ≤ 2 logB5N = 1/c1 (Lemma 13) we get the third inequality.

Lemma 15. For some constant c1 ∈ R and the choice of B = 4c21 log 5N , the probability of error in
identifying the position of second non-zero variable node at a double-ton at any bin (i, j), given that the
first position identification is correct, can be upper bounded by


N µ+α−1 (c1 (1 − 2η) − 1)2
P[E22 ] ≤ exp −
8(1 + c21 )
Proof. E22 :
Let jp and jp̃ be the two variable nodes participating in the doubleton (i, j). Then the observation
vector ~zi,j is given by


h

z i,j = w
~ j1 , w
~ j2 , · · ·

X













r[jp ]


 . 
 .. 

i 



w
~ j gi ×  n j 

 . 
 .. 






r[jp̃ ]


 .. 
 . 


ngi

w
~ jp , · · ·

= r[jp ] w
~ jp + r[jp̃ ] w
~ jp̃ +

n1
..
.

nk w
~ jk

k6=p

Let the contribution from jp̃ be peeled off from the doubleton at some iteration, then we get
~zi,j = r[jp ] w
~ jp +

X
e1
w
~ jp̃ +
nk w
~ jk
B
k6=p

where e1 ∈ [−ηM, ηM ] is an extra error term induced due to peeling off.
Now the estimated second position p̂ is calculated using Eq. (11). We can observe that the position is
wrongly identified when ∃p0 such that
w
~ j†p z i,j
B

April 27, 2017

≤

w
~ j†p0 z i,j
B

DRAFT

25

=⇒ r[jp ] +

1 X
e1 †
nk w
~ j†p w
~ jk + w
~ w
~j
B
B jp p̃
k6=p,p̃

≤

↔

r[jp ] †
1
w
~ jp0 w
~ jp + n p 0 +
B
B

X

X

nk w
~ j†p0 w
~ jk +

k6=p,p0 ,p̃

αk nk + βnp0 ≥ r[jp ] 1 −

w
~ j†p0 w
~ jp

k6=p,p0 ,p̃

!

B

−

e1 †
w
~ w
~j
B jp0 p̃

2ηM †
w
~ jp 0 w
~ jp̃
B

≥ M (1 − 2η)(1 − µmax ) − 2ηM µmax = M (1 − 2η − µmax )

where αk and β are constants and can be shown to be in the range αk ∈ [−2µmax , 2µmax ] and β ∈
[1 − µmax , 1 + µmax ]. Now using the bound given by Chernoff Lemma in Lem. 7 we obtain


2M (1 − 2η − µmax )2
P[E22 ] ≤ exp −
16gi µ2max + 4(1 + µmax )2


2M (1 − 2η − µmax )2
≤ exp −
16(gi µ2max + 1)


M (c1 (1 − 2η) − 1)2
≤ exp −
8(gi + c21 )


N µ+α−1 (c1 (1 − 2η) − 1)2
≤ exp −
8(1 + c21 )

q
where for the choice of B = 4c21 log 5N , µmax ≤ 2 logB5N = 1/c1 .

R EFERENCES
[1] A. Andoni, H. Hassanieh, P. Indyk, and D. Katabi, “Shift finding in sub-linear time,” in Proceedings of the Twenty-Fourth
Annual ACM-SIAM Symposium on Discrete Algorithms, pp. 457–465, Society for Industrial and Applied Mathematics,
2013.
[2] A. Amir, M. Lewenstein, and E. Porat, “Faster algorithms for string matching with k mismatches,” Journal of Algorithms,
vol. 50, no. 2, pp. 257–275, 2004.
[3] G. Navarro, “A guided tour to approximate string matching,” ACM computing surveys (CSUR), vol. 33, no. 1, pp. 31–88,
2001.
[4] R. S. Boyer and J. S. Moore, “A fast string searching algorithm,” Communications of the ACM, vol. 20, no. 10, pp. 762–772,
1977.
[5] P. Ferragina and G. Manzini, “Indexing compressed text,” Journal of the ACM (JACM), vol. 52, no. 4, pp. 552–581, 2005.
[6] H. Li and R. Durbin, “Fast and accurate short read alignment with burrows–wheeler transform,” Bioinformatics, vol. 25,
no. 14, pp. 1754–1760, 2009.
[7] H. Li and R. Durbin, “Fast and accurate long-read alignment with burrows–wheeler transform,” Bioinformatics, vol. 26,
no. 5, pp. 589–595, 2010.

April 27, 2017

DRAFT

26

[8] N. Zhang, A. Mukherjee, D. Adjeroh, and T. Bell, “Approximate pattern matching using the burrows-wheeler transform,”
in Proceedings of the Conference on Data Compression, p. 458, IEEE Computer Society, 2003.
[9] W. I. Chang and T. G. Marr, “Approximate string matching and local similarity,” in Annual Symposium on Combinatorial
Pattern Matching, pp. 259–273, Springer, 1994.
[10] H. Hassanieh, F. Adib, D. Katabi, and P. Indyk, “Faster GPS via the sparse Fourier transform,” in mobicom, pp. 353–364,
ACM, 2012.
[11] S. Pawar and K. Ramchandran, “A robust R-FFAST framework for computing a k-sparse n-length DFT in O(k log n)
sample complexity using sparse-graph codes,” in Proc. IEEE Int. Symp. Information Theory, pp. 1852–1856, 2014.
[12] K. Lee, R. Pedarsani, and K. Ramchandran, “Saffron: A fast, efficient, and robust framework for group testing based on
sparse-graph codes,” arXiv preprint arXiv:1508.04485, 2015.

April 27, 2017

DRAFT

