Trace reconstruction from deletion traces
A complete and readable account of the quasipolynomial sample bound
========================================================================
1. THE PROBLEM
========================================================================
Imagine a long string of symbols: a piece of DNA, or a line of text
copied by a careless scribe.
Nature makes a copy by walking along the string and, independently at
each position, keeping the symbol with probability p and deleting it
with probability q = 1 - p. What is left is called a trace. A trace is
a subsequence of the original: the surviving symbols in their original
order, with no record of where the deletions happened, and with no
labels on the survivors saying which original positions they came from.
You are given many independent traces. You are also told p and the
original length n. You must recover the original string.
Question. How many traces do you need?
Two words to fix.
Worst case means the bound must hold for every possible original
string, not merely for typical ones. An adversary picks the hardest
string.
Sample complexity means the number of traces. It does not mean the
running time of whatever program processes them.
A tiny example. Take n = 4 and x = 0110, with p = 3/4. The retention
pattern (keep, keep, delete, keep) has probability (3/4)^3 (1/4) and
produces the trace 010. The pattern (delete, keep, keep, delete)
produces 11. The pattern (keep, keep, keep, keep) produces 0110. The
pattern (delete, delete, delete, delete) produces the empty trace.
If p is small, most traces are short, and each trace tells you very
little. The whole difficulty of the subject is that p may be as small
as 1/n, in which case a typical trace has length about 1.
Two features of the problem make it hard.
The traces do not tell you where the deletions were. If they did, one
trace would trivially suffice.
Two different strings can produce traces that are easy to confuse. A
string with many repeated symbols is harder than a random string.
========================================================================
2. THE ANSWER
========================================================================
Throughout, log means natural logarithm unless another base is stated,
and exp(t) means e to the power t.
THEOREM 1 (main bound). Let n >= 2 be an integer and let 0 < p <= 1.
For every binary string x of length n, there is a reconstruction rule
that recovers x from independent deletion traces with probability at
least 2/3, provided the number of traces is at least
exp( O( p^{-1} (log n)^3 (1 + log log(2n))^6 ) ).
The O(.) hides an absolute constant that does not depend on n, on p, or
on x. The same number of traces works for all x simultaneously. The
rule may depend on n and on p (both of which are known) but not on x.
THEOREM 2 (finer form). Put q = 1 - p and
ell_* = 1 + log(1/q)
X = 1 + log(n) / ell_*.
Then the logarithm of the number of traces needed is at most
O( p^{-1} (log n) X^2 (1 + log X)^6 ).
THEOREM 3 (low deletion). Suppose q <= n^{-epsilon} for some fixed
epsilon > 0. Then either q = 0, in which case one trace suffices, or
ell_* >= 1 + epsilon log n, so X <= 1 + 1/epsilon, while p >= 1/2 for
large n. Hence polynomially many traces suffice; the degree of the
polynomial may depend on epsilon.
THEOREM 4 (general alphabets). The same asymptotic bounds hold for
strings over an arbitrary alphabet, assuming the surviving symbols are
observed exactly.
Sanity check on sizes. Fix p = 1/2 and n = 10^9. Then log n is about
20.7, (log n)^3 is about 8900, log log(2n) is about 3.06,
(1 + log log(2n))^6 is about 4500, and the product is about 4.0 * 10^7.
With p^{-1} = 2 included, the theorem promises fewer than
exp(8 * 10^7) traces. That is a finite number, and the shape of the
exponent is the point. A trivial bound is 2^n, which here is
10^{3 * 10^8}: vastly larger.
Vocabulary. A function of n is polynomial if it is at most n^C for some
fixed C. It is quasipolynomial if it is at most exp((log n)^C) for some
fixed C. Quasipolynomial growth sits strictly between polynomial and
exponential: for every fixed C and every fixed c > 0, exp((log n)^C)
grows faster than n^c and slower than exp(n^c).
Three things the theorem does not say.
It is not a matching lower bound. The true number of traces needed
might be much smaller.
It is not an efficient algorithm. See section 17.
For fixed p in (0,1) the bound stays quasipolynomial. It does not
become polynomial in n. Getting polynomial for any fixed p is open.
========================================================================
4. NOTATION AND THE TOOLS WE BORROW
========================================================================
4.1 Notation.
log natural logarithm
exp(t) e to the power t
a^b a to the power b
|z| modulus of the complex number z
Re z, Im z real and imaginary parts
sum, prod sum and product
inf infinity
O(g), o(g) see 4.3
|x|_2 Euclidean length of a vector
|x|_inf largest coordinate in absolute value
||mu|| total variation of a complex measure
sup, inf supremum, infimum
1_E indicator of the event E
gcd(a,b) greatest common divisor
A polynomial here is always an ordinary polynomial in one complex
variable z. A quantity is analytic in a region if it is given there by
a convergent power series. A measure is a rule assigning a complex
number to each nice set, additively; its total variation is the largest
possible value of the sum of absolute values of the pieces, over all
finite partitions.
4.2 Seven standard facts.
FACT 1 (geometric series). If |z| < 1 then sum_{m>=0} z^m = 1/(1-z),
and sum_{m>=0} |z|^m = 1/(1-|z|).
FACT 2 (Cauchy-Schwarz). |sum_i x_i y_i|^2 <= (sum_i |x_i|^2)
(sum_i |y_i|^2).
FACT 3 (union bound). For any events A and B, P(A or B) <= P(A) +
P(B).
FACT 4 (Hoeffding). Let S be the average of N independent random
variables, each in [0,1]. Then for every t > 0,
P( |S - E S| >= t ) <= 2 exp(-2 N t^2).
So to resolve a gap of size t you need about log(1/delta)/t^2
samples, where delta is the error you are willing to tolerate.
FACT 5 (Poisson). If f is analytic and nonzero on a disk, then log|f|
is harmonic there, so its value at the centre is a weighted average
of its boundary values, with the Poisson weight
(1 - r^2) / (2 pi (1 + r^2 - 2 r cos(angle)))
at radius r. We only use one consequence: a large value at one
interior point forces a comparable value somewhere on any arc that
carries enough Poisson weight. Boundary zeros are harmless; apply the
ineuality to log(|f| + delta) and let delta tend to 0.
FACT 6 (Hadamard three-lines). If G is analytic and |G| <= 1 in the
vertical strip a <= Re s <= b, then log sup_{Re s = t} |G| is a
convex function of t in [a,b]. Consequence: a lower bound at one
interior line plus an upper bound at both edges gives a lower bound
at the other edge, and the price is the ratio of the distances.
FACT 7 (Fine and Wilf). If a word w has periods p and q and
|w| >= p + q - gcd(p,q), then w has period gcd(p,q).
A period of a word of length L is a positive integer delta such that
every pair of positions at distance delta carries the same symbol,
for as long as both positions lie inside the word. Every word
trivially has every period delta >= L.
4.3 Big-O and little-o.
"f(n) = O(g(n))" means there is a constant C, depending only on
quantities that are explicitly fixed and not on n, such that
|f(n)| <= C g(n) for all sufficiently large n.
"f(n) = o(g(n))" means f(n)/g(n) tends to 0 as n tends to infinity.
In this paper, every O(.) and o(.) statement is uniform over the
remaining parameters: the channel probabilities in the stated range,
the induction level, the choice of pair, and so on.
"For all sufficiently large n" means there is an absolute threshold
n_0 such that the conclusion holds for n >= n_0. These thresholds are
never computed, and they are always far smaller than the actual size
of n in any application.
========================================================================
6. DOUBLING, ANCHORS, PROBES, AND THE DIFFERENCE POLYNOMIAL
========================================================================
6.1 Small cases.
n = 1. About 1/p traces suffice, because each trace either shows the
single symbol or not.
2 <= n <= some fixed constant. Collect ceil(p^{-n} log 3) traces.
With probability at least 2/3 at least one trace has full length, and
then you can simply read off x.
p < 1/n with n >= 2. The same full-length strategy costs at most
exp(O(p^{-1} log n)). Indeed, with s = n p <= 1,
n log(1/p) = p^{-1} s (log n + log(1/s)) = O(p^{-1} log n).
q < 1/n with n >= 2. Then p^n >= (1 - 1/n)^n >= 1/4, so a constant
number of traces suffices by the same strategy.
So assume from now on that p >= 1/n and q >= 1/n.
6.2 Doubling.
Prepend a known block of zeros of length n to each candidate. A trace
of the extended input is simulated from a trace of the original input
by prepending an independent deletion trace of the zero block. From
here on the letter n denotes the doubled length. At the very end switch
back to the original length; the exponent changes only by an absolute
factor.
Fix two distinct extended binary candidates x and y. They share the
known zero prefix. Let d >= n/2 be their first mismatch, using 0-based
indices.
6.3 The constants.
All of these are absolute constants or simple expressions in n, p, q.
They look arbitrary. They are chosen once and never changed.
D = 10^14 a safety constant
epsilon = 10^{-5} the mask decay rate
ell_* = 1 + log(1/q)
b = ceil( log_2( 1 + log n / ell_* ) )
L_h = 2^{b-h} for 0 <= h <= b
M_h = 4^{b-h} for 0 <= h <= b
a_0 = exp(-10^16) / (D L_0)
t_h = a_{h+1} = D q L_h a_h^2
P_0 = ceil(1 / a_0)
R_h = D q (h+1)^2 M_h a_h^2
B_j = 10^100 ell_* 8^j (b+1)^3
f(u) = exp(u^2 / 20)
What each one is for.
ell_* measures how much a single trace carries. It is large exactly
when p is small. It appears everywhere as the cost of seeing one
symbol.
b is the number of induction stages. For fixed p and large n it
grows like log log n; more precisely, 2^b is about log(n)/ell_*.
a_h is the radius parameter at stage h. It shrinks super-
exponentially: setting v_h = D q L_h a_h gives v_0 = q exp(-10^16)
and v_{h+1} = v_h^2 / 2, so v_h = 2 (v_0/2)^{2^h}, and a_h decays
even faster. In particular a_b n < .0001.
t_h = a_{h+1} is the decay rate that level-h masks may use. Since a
is decreasing in the index, t_h >= a_j for every h < j.
P_0 is the length parameter of the period-avoiding word. It is
O(2^b), so log(2 P_0) = O(b).
R_h is the decay rate of the compact representation of a level-h
mask, and R_h / t_h = (h+1)^2 M_h / L_h >= 1.
B_j is the log-scale accuracy demanded at stage j. It multiplies by 8
at each stage, but 8^b is only about X^3.
Easy consequences, used repeatedly:
b >= 1;
ell_* (2^b - 1) >= log n, hence ell_* (b+1) -> infinity, hence
B_0 -> infinity uniformly;
a_{h+1}/a_h = v_h <= v_0 = q exp(-10^16), which is far below .001;
log(1/a_h) <= 10^18 (b + 2^h) ell_*;
P_0 = O(2^b), so log(2 P_0) = O(b), independent of ell_*.
6.4 Probes.
A tuple probe of order k >= 1 consists of
indices 0 <= i_1 < ... < i_k < n,
predicates chi_e : {0,1} -> {0,1}, for e = 1,...,k,
edges the pairs (i_e, i_{e+1}) for e = 1,...,k-1.
The value of the probe at that tuple is the product of the numbers
chi_e(x_{i_e}) over e.
Some edges are designated exact. An exact edge carries the extra
factor 1_{i_{e+1} - i_e = 1}; it forces the two indices to be
consecutive and contributes nothing otherwise. All other edges are
called loose.
The block count A is one plus the number of loose edges. Equivalently,
A is the number of maximal runs of consecutive forced indices. Note
that different blocks need not be separated by gaps of size greater
than one; they can touch.
6.5 Masks.
A probe at stage j also carries finitely many masks of levels h < j.
A mask of level h uses r-1 differences, with 2 <= r <= 4:
g_s = i_{b_s} - i_{c_s} for s = 1,...,r-1,
for some chosen indices 1 <= c_s <= b_s <= k. All of them are
nonnegative by the ordering, and one may be zero if the two endpoints
coincide.
A mask is a factor
K(g) * prod_{s=1}^{r-1} exp(i beta_s g_s),
where the beta_s are real with |beta_s| <= 10^4 a_h, and K may be
complex. K must satisfy, for all g_s >= 0,
(M1) |K(g)| <= 2 exp( - epsilon t_h (g_1 + ... + g_{r-1}) ),
and must also have a compact representation
(M2) K(g) = exp( - R_h (g_1 + ... + g_{r-1}) )
* int_{R^{r-1}} exp(i xi . g) mu(d xi),
with mu a complex measure supported on the cube |xi|_inf <= C R_h and
total variation
||mu|| <= exp( C (R_h / t_h) B_h ).
Here C is a large absolute constant, the same at every level. Its
existence follows from the construction, and one value serves for all
stages.
(M1) and (M2) in plain terms.
(M1) says that a mask is negligible unless all of its differences are
small, of order 1/t_h. This is what keeps the probe's sums bounded.
(M2) says that a mask is an average of pure phases with bounded
frequencies. In other words, a mask behaves like a finite combination
of honest phase factors, so that later steps can treat it as a tuple
with extra phases. The total variation is the price of that
approximation, and it is exponential in (R_h/t_h) B_h.
A difference g_s covers exactly the edges (i_e, i_{e+1}) inside its
interval, counting each edge once; a zero difference covers nothing.
The masks' difference intervals are how the induction remembers
things. At each stage, the loose edges of the new probe must be covered
by masks, either old ones or the one new mask created at that stage.
6.6 Latest-anchor arrays and the difference polynomial.
Let W_x(i_1,...,i_k) be the product of all bit, exact-gap, and mask
factors of a probe. Define
s_x(m) = sum over 0 <= i_1 < ... < i_k < n with i_k = d + m
of W_x(i_1,...,i_k),
and define s_y(m) the same way with y. Invalid latest positions
contribute zero. Now define the difference polynomial
H(z) = sum_m ( s_x(m) - s_y(m) ) z^m.
Key observation. For m < 0 the condition i_k = d + m forces every index
to be less than d. On that region x and y agree, so each term cancels
against its counterpart. Hence H has no negative powers. It is an
ordinary polynomial.
This is the whole point of the latest-anchor convention. The latest
index is either at the mismatch or to the right of it; anything to the
left has already been decided and cancels.
========================================================================
7. THE INDUCTION GOAL
========================================================================
At stage j, with 0 <= j <= b, we want a probe with the following
properties.
(G1) Order and blocks.
1 <= k <= 2 P_0 4^j and 1 <= A <= 4^j.
(G2) Masks. For every level h < j there are at most 4^{j-h-1} masks of
that level, and every loose edge is covered by at least one of
their differences.
(G3) A pivot. There is a real theta with |theta| <= 5 a_j such that
|H( exp(-a_j + i theta) )| >= exp( - B_j / f(theta/a_j) ).
Stop as soon as a_j n <= .001. This happens by j = b, because
a_b n < .0001. At the moment of stopping, the radius parameter a_j is at
most .001/n, i.e. very far inside the unit disk.
Interpretation of (G3). The point exp(-a_j + i theta) sits just inside
the unit circle, at distance about a_j from it. The factor
1/f(theta/a_j) is between exp(-25/20) and 1 when |theta| <= 5 a_j, so it
is a mild penalty for choosing a direction theta away from the positive
real axis.
Size of the sums. Every probe satisfying (G1) and (G2) has
sup_m |s_x(m)| <= Z_j, sup_m |s_y(m)| <= Z_j,
where
Z_j = 2^{4^j} * ( 2 / (epsilon a_j) )^{4^j}.
Why. Each loose edge contributes a decay factor at least
exp(-epsilon a_j Delta), where Delta is its gap, because every level
h < j has t_h >= a_j and the coverage condition forces at least one
applicable mask. Summing over Delta >= 1 gives at most
1/(1 - exp(-epsilon a_j)) <= 2/(epsilon a_j). There are at most
A <= 4^j loose edges. Bit predicates are at most 1. Each mask factor has
modulus at most 2 by (M1), and there are at most 4^j of them, giving
2^{4^j}. Exact edges are fixed at distance 1 and contribute 1.
Consequence, used many times below: the auxiliary magnitudes are tiny
compared to the accuracy budget. For 0 <= j < b,
log Z_j + log Z_{j+1} + log(1/a_j) + log(1/t_j)
<= 10^21 ell_* 4^j (b + 2^j)
<= 10^{-60} B_j.
Proof: each term is at most about 4^j log(1/a_j), and
log(1/a_h) <= 10^18 (b + 2^h) ell_*. So the left side is at most
10^18 * 4^j * (b + 2^j) * ell_* times a small constant. And
4^j (b + 2^j) <= 8^j (b + 1) <= 8^j (b + 1)^3.
Since B_j contains 10^100 * 8^j * (b + 1)^3 * ell_*, there is enormous
room.
========================================================================
8. THE BASE CASE: A PERIOD-AVOIDING WORD
========================================================================
We must produce a probe at stage 0. It will have all edges exact, no
masks, and order k = 2 P_0.
8.1 Finding the word.
By the padding, d >= 2 P_0 for all sufficiently large n. Let
w' = x_{d-2P_0+1} ... x_{d-1} = y_{d-2P_0+1} ... y_{d-1}
be the common last 2 P_0 - 1 symbols before the mismatch.
CLAIM. At least one of the two words w'0 and w'1 has no period strictly
smaller than P_0.
Proof. Suppose w'0 has period p_1 < P_0 and w'1 has period p_2 < P_0.
Then w', being a prefix of each, has both periods. Since
|w'| = 2 P_0 - 1 >= p_1 + p_2 - gcd(p_1, p_2)
(because p_1, p_2 <= P_0 - 1), Fine and Wilf gives that w' has period
g = gcd(p_1, p_2).
Look at positions 2 P_0 - p_1 and 2 P_0 - p_2 inside w'. They differ by
p_2 - p_1, which is a multiple of g, so they carry the same symbol.
But period p_1 of w'0 forces (w'0)_{2P_0 - p_1} to equal (w'0)_{2P_0},
which is 0. And period p_2 of w'1 forces (w'1)_{2P_0 - p_2} to equal
(w'1)_{2P_0}, which is 1. So the common symbol would be both 0 and 1.
Contradiction.
Choose w in {w'0, w'1} with no period below P_0.
8.2 Building the base probe.
Let k = 2 P_0. Take all edges exact. Take predicates chi_e that test
the e-th symbol of w. There are no masks, so A = 1.
The number of occurrences of w ending exactly at d differs between x
and y: by construction one of them has w there and the other does not.
Hence s_x(0) - s_y(0) is +1 or -1, so the two strings already disagree
at offset 0.
8.3 Occurrences are far apart.
Suppose w occurs in x ending at d + m_1 and again ending at d + m_2,
with 0 < m_2 - m_1 < P_0. Call delta = m_2 - m_1.
The two occurrences have length 2 P_0 and are shifted by delta. At every
position covered by both, the symbol is determined twice: once from the
first occurrence and once from the second. Reading off the two
determinations, we get w_j = w_{j+delta} whenever both indices lie in
[1, 2 P_0]. That is exactly the statement that w has period delta.
Since delta < P_0, this contradicts the choice of w. The same argument
applies in y. So occurrence endings in each individual string are
P_0-separated.
8.4 The shape of H.
Putting the two previous paragraphs together, up to an overall sign,
H(z) = 1 - eta z^u + sum_{m >= P_0} c_m z^m
with
1 <= u < P_0, eta in {0,1}, |c_m| <= 1,
and if eta = 0 then u is arbitrary, since the term -eta z^u is then
absent for every choice of u.
Explanation. s_x(0) - s_y(0) = +-1 gives the constant term. The other
string has its nearest occurrence ending at d + u for some 1 <= u < P_0,
or none at all, in which case eta = 0; this gives the -eta z^u term.
Every further occurrence, in either string, ends at least P_0 away from
d, giving the tail.
8.5 The pivot by the Poisson ineuality.
Let v = 10 log(P_0) / P_0. At the positive real point z = exp(-v):
The first two terms have modulus at least
1 - exp(-u v) >= 1 - exp(-v) >= v/2,
since v is tiny.
The tail is at most
sum_{m >= P_0} exp(-v m) = exp(-v P_0)/(1 - exp(-v))
<= 2 P_0^{-10} / v <= v/6,
because exp(-v P_0) = P_0^{-10} and 1 - exp(-v) >= v/2.
Therefore |H(exp(-v))| >= v/2 - v/6 = v/3.
On the circle |z| = exp(-a_0) the crude bound |H(z)| <= 2/a_0 holds.
There are at most about 1/a_0 contributing offsets, each coefficient at
most 2; since 1/a_0 is astronomically large, this is trivial.
The Poisson ineuality on the disk of radius exp(-a_0) centred at 0,
evaluated at the interior point exp(-v) which is much closer to 1, then
gives
max_{|theta| <= a_0} |H(exp(-a_0 + i theta))|
>= exp( - O( (v/a_0) log(1/a_0) ) ).
Indeed, the harmonic weight of that arc is at least a constant times
a_0/v (because v >= 2 a_0 and v < .01), and log|H| is bounded above by
log(2/a_0) elsewhere. Boundary zeros are removed by the usual
regularization log(|H| + delta), delta -> 0.
Size of the exponent. We have
log(1/a_0) = 10^16 + log(D L_0) <= 10^16 + 35 + b,
v/a_0 <= 10 log(2 P_0) <= 10 (10^16 + 36 + b),
so the exponent is at most about 10 (10^16 + 36 + b)^2, which is under
1.1 * 10^33 for any conceivable n. On the other hand B_0 is at least
10^100. So the exponent is at most 10^{-66} B_0, comfortably below the
requirement B_0 / f(theta/a_0) >= 0.94 B_0.
This gives (G1), (G2), and (G3) at j = 0: order 2 P_0 <= 2 P_0 4^0,
block count 1 <= 4^0, no masks, and |theta| <= a_0 <= 5 a_0.
========================================================================
9. THE STEP: A GAUSSIAN WINDOW AND AN ARC TEST
========================================================================
Suppose (G1), (G2), (G3) hold at stage j < b. Write
a = a_j, t = a_{j+1}, B = B_j, E = B / f(theta/a_j).
Recall t/a <= v_0, far below .01; a <= a_0, far below .001; and, during
a genuine step (one not at the stopping point), a > .001/n.
9.1 Step one: a window.
Define the window
sigma^2 = 2.4 B / t^2
ell = floor( 10^{-10} E )
v(m) = exp(-y^2/2) * sum_{k=0}^{ell} (y^2/2)^k / k!,
with y = m/sigma throughout. The function v is a polynomial in m^2 of
degree 2 ell; it satisfies 0 <= v <= 1; and it is a smooth window, equal
to 1 inside a window and 0 far outside.
CLAIM 9.1. For 0 <= m <= 3 E/a we have 1 - v(m) <= exp(-2E).
Proof. On this range,
y^2/2 <= 2 E (t/a)^2 <= 2 E exp(-2 * 10^16),
because t/a <= exp(-10^16). So w = y^2/2 is astronomically small, and
the Taylor remainder of the exponential series after ell terms,
multiplied by exp(-w), is at most (e w/(ell+1))^{ell+1}, which is far
below exp(-2E) because ell is a positive fraction of E.
Interpretation: the window radius 3E/a is astronomically smaller than
the Gaussian scale sigma, so within the window the correction is
negligible.
Project to real or imaginary part. Choose c(s) = Re(s) or c(s) = Im(s),
whichever keeps at least half of the pivot. Since H = Re H + i Im H,
one of the two has modulus at least half of |H|.
Define, for 0 <= m with d + m < n, and zero elsewhere,
f_x(m) = ( 2 + c( s_x(m)/Z_j ) ) v(m) / Y
f_y(m) = ( 2 + c( s_y(m)/Z_j ) ) v(m) / Y
Y = 3 sum_m v(m).
Two immediate facts. Since |s_x|, |s_y| <= Z_j, we have
0 <= 2 + c(s/Z_j) <= 3, so 0 <= f_x, f_y and sum_m f_x(m), sum_m f_y(m)
<= 1. Also Y >= 3, because the term m = 0 has v(0) = 1; so there is no
division by zero.
The normalizing factor is harmless: log Y <= 10^{-8} E. To see this,
split at H_* = ceil(1000 B / t). The sum of v within that range is at
most about 2 H_*, and outside it v(m) is at most
2^ell exp(-m^2/(4 sigma^2)), which decays like a Gaussian. Using
log(1/t) <= 10^{-60} B, B <= 4E, and log B = o(E), the claim follows.
9.2 Step two: a polynomial on the unit circle.
Define
J(z) = sum_m ( f_x(m) - f_y(m) ) z^m.
Two bounds.
On the whole unit circle, sup_{|z| = 1} |J(z)| <= 2, because
sum_m |f_x - f_y| <= sum_m f_x + sum_m f_y <= 2.
At the pivot z_0 = exp(-a + i theta),
|J(z_0)| >= exp( -1.0005 E ).
Reason: masking costs at most exp(-2E) on the window (Claim 9.1) and
is exponentially small beyond it; 1/Y costs at most exp(-10^{-8} E);
log Z_j costs at most 10^{-60} B. Adding these up gives the stated
loss.
The truncation of J to |m| <= H_* is denoted J_diamond. Because all
negative powers of z vanish in the difference, and the omitted tail has
l_1 norm below exp(-200E), we also get
|J_diamond(z_0)| >= exp( -1.0006 E ).
9.3 Step three: a discrete grid on the arc.
Put
N = ceil( (1000 B / t)^6 ),
h = 2 pi / (N a),
v_0^* = theta / a,
gamma = floor( 10^4 / h ),
F(l) = sum_m f_x(m) exp( 2 pi i l m / N ),
G(l) = sum_m f_y(m) exp( 2 pi i l m / N ).
The grid points are z_l = exp(-a) exp(2 pi i l / N) = exp(-a + i a l h),
for |l| <= gamma. Note that |F|, |G| <= 1.
CLAIM 9.2 (phase test). Put Lambda = exp(-6B). Then there exist
r in {2,3,4} and integers l_1,...,l_r with
|l_i| <= gamma and |l_1 + ... + l_r| <= 1,
such that
| prod_i F(l_i) - prod_i G(l_i) | > Lambda.
Proof. Introduce the deficits
m_l = max( |F(l)|, |G(l)| ) (so 0 <= m_l <= 1)
g_l = min( 50, E^{-1} log(1/m_l) )
q_l = min( 50, E^{-1} log(2 m_l / |F(l) - G(l)|) ),
with the conventions g_l = 50 when m_l = 0, and q_l = 50 when the
denominator vanishes. Both are nonnegative and symmetric, because F and
G are real on the grid by construction.
The grid weights are
k_l = h / ( pi ( 1 + (l h - v_0^*)^2 ) ) for |l| <= gamma,
and we define
G_cal = sum_l k_l g_l, Q = sum_l k_l q_l.
The Poisson ineuality, applied to the pivot bound on J_diamond, implies
the deficit budget
G_cal + Q <= 1.003.
Proof of the deficit budget, in two parts.
Continuous bound. Since J_diamond is a polynomial, log|J_diamond| is
subharmonic on the unit disk. Evaluating at the interior point z_0 (push
the disk radius just past exp(-a) and let it return; capping at 100E can
only lower an average) and using the pivot bound gives
Poisson average of log( 2 / |J_diamond| ) <= log 2 + 1.0006 E,
which is at most 1.003 E for E >= 300. Under theta = a x, the Poisson
kernel at z_0 has density at least (1 - 2a)/(pi (1 + (x - v_0^*)^2)),
so the continuous integral dominates the Riemann sum that produces
G_cal + Q.
Discretization. By layer cake, f(x) = int_0^{infinity} 1_{f(x) > c} dc
for f(x) = min(100, E^{-1} log(2/|J_diamond(exp(-a + i a x))|)). A
polynomial with H_* index terms meets a fixed positive modulus at most
O(H_*) distinct values of x. So each level set f > c is a union of
O(H_*) intervals. Replacing each interval's length by the grid sum of
mesh h adds O(H_* h) = o(1) uniformly in n.
Now the combinatorial core. For every integer j'shift with
1 <= j'shift <= S/h, every fractional vertex cover of the graph on the
integer points in [-gamma, gamma] with edges (l, l + j'shift) has total
k_l weight at least M. Here S and M are the numbers in the table below.
Proof. Restrict to one residue class mod j'shift: it is a chain. Extend
the weights to the full Cauchy kernel on the whole line; the tail costs
less than .0001. A fractional cover assigns a number in [0,1] to each
vertex, with endpoint sums at least 1. It can be rounded to a random
binary cover without increasing expected weight: pick one uniform
threshold U_0 in [0,1], use "vertex number >= U_0" on even vertices and
"vertex number >= 1 - U_0" on odd vertices. So it is enough to bound
binary covers.
In each chain, take the central vertex l nearest to v_0^*/h. The
weights decrease outward along both arms. A binary cover either omits
the centre, in which case both first neighbours are forced and the
disjoint edges (2,3), (4,5), ... force the odd weights on both arms; or
it includes the centre, in which case the disjoint edges (1,2), (3,4),
... force the even weights. So the minimum is the lesser of the two
parity masses.
The parity masses have closed forms. With d' = j'shift * h,
x = l h - v_0^*, and K_Cauchy(x) = 1/(pi (1 + x^2)),
Xi(x) = sum_{nu in Z} K_Cauchy(x + 2 nu d')
= (1/(2d')) * (1 + 2 sum_{m>=1} exp(-m pi/d') cos(m pi x/d')),
using the Fourier transform identity
int K_Cauchy(x) exp(i xi x) dx = exp(-|xi|) and summing over periods.
The even mass is h Xi(x). The odd mass is h Xi(x + d'), which is never
larger: write Xi in closed form with r_0 = exp(-pi/d'),
Xi(x) = (1 - r_0^2) / (2 d' (1 + r_0^2 - 2 r_0 cos(pi x/d'))),
and cos is maximal at x = 0.
Summing the odd masses over the j'shift central representatives, up to
O(h) error, gives the integral
int_{-d'/2}^{d'/2} Xi(x + d') dx
= 1/2 - (2/pi) arctan( exp(-pi/d') ).
This is smallest at the largest possible d', which is S. So the cover
weight is at least M.
Second, directly by Riemann sum,
sum_{|l| <= S/h} k_l > J_*.
Table. Round |v_0^*| up to w in {1,2,3,4,5}. Set
w 1 2 3 4 5
S 1.5 2.1 2.8 3.6 4.4
M .419 .357 .297 .246 .208
J_* .524 .453 .380 .335 .292
These are conservative, in the sense of the four ineualities
M < 1/2 - (2/pi) arctan( exp(-pi/S) ) - .0002
J_* < (1/pi)( arctan(S+w) + arctan(S-w) ) - .0002
T_lim(v) = 6 f(v) - .1
T_lim(w-1) > max{ 2.006/M,
1.005 / (J_* - (1 + 2 J_*/M)^2 M/16) }.
All of them can be checked with rational arithmetic, without any
sensitive numerical integration:
use 3.14159 < pi < 3.14160;
bound exp(-pi/S) from above by a truncated exponential series;
write arctan(S+w) = pi/2 - arctan(1/(S+w));
the remaining arctangent arguments have modulus at most .6, so
evaluating their alternating series to six terms costs at most
.6^{15}/15;
lower-bound exp(u) in f(w-1) by 1 + u + u^2/2.
The resulting worst-endpoint values are
w 1 2 3 4 5
max{...} 4.96 5.79 6.94 8.26 9.71
while
T_lim(w-1) = 6 exp((w-1)^2/20) - .1
equals about 5.90, 6.21, 7.23, 9.31, 13.25. Each exceeds the
corresponding maximum (4.96, 5.79, 6.94, 8.26, 9.71).
Now the contradiction. Put R = G_cal/M and T_0 = T_lim(v_0^*), which is
at least T_lim(w-1). Then 2R < T_0, and there exists an integer s with
|s| <= S/h, 2R + q_s < T_0, 2g_s + q_s < T_0.
Otherwise, pointwise on the whole interval,
q_s >= (T_0 - 2R) max( 0, 1 - 2 g_s/T_0 ),
so by the cover bound
Q >= (T_0 - 2 G_cal/M)(J_* - 2 G_cal/T_0),
and minimizing this quadratic in G_cal gives
G_cal + Q >= T_0 ( J_* - (1 + 2 J_*/M)^2 M/16 ) > 1.004,
contradicting G_cal + Q <= 1.003. Choose the witness with s >= 0 by
symmetry.
For each 0 <= j'shift <= s, choose a reference pair (e_j, e_j +
j'shift) whose end deficits sum to at most R. Such a pair exists: if
the minimum were R' > R, the clipped deficits min(1, g_l/R') would be a
fractional cover of cost G_cal/R' < M, contradicting the cover bound.
For j'shift = 0, minimize 2 g_l and use the total kernel mass > .999
together with M < .499.
Suppose the claim fails. At every pair (l, -l), the failure bounds the
modulus of F(l) - G(l) by Lambda/m_l. At all relevant indices
2 g_l < T_0 < 50, so the cap is inactive and m_l = exp(-g_l E). Write
U(l) = arg F(l) - arg G(l) modulo 2 pi.
By reality of F and G, U(-l) = -U(l). Set
d_j = U(e_j + j) - U(e_j), so d_0 = 0.
The failed quadruple product at frequencies
(e_j + j, -e_j, -e_{j-1} - j + 1, e_{j-1})
sums to 1, so it bounds d_j - d_{j-1} mod 2 pi by 100 Lambda exp(2 R E):
both products have modulus at least exp(-2RE)/16, and two complex
numbers differing by at most Lambda, both of modulus at least v, have
angular distance at most pi Lambda/v.
Summing over j gives
|d_s| <= 100 N Lambda exp(2 R E) modulo 2 pi.
Similarly, the failed triple product at (e_s + s, -e_s, -s) sums to 1
and gives
|U(s) - d_s| <= 100 Lambda exp( (R + g_s) E ).
But R + g_s + q_s <= max(2R + q_s, 2g_s + q_s) < T_0. Therefore
|F(s) - G(s)| / m_s
<= Lambda/m_s^2
+ 100 Lambda ( exp((R+g_s)E) + N exp(2RE) )
< exp(-q_s E),
using log N <= 10^{-8} E and R + g_s + q_s < T_0 < 50. That contradicts
the uncapped definition of q_s, which was supposed to stay below 50. The
claim holds.
Consequence. We now have r in {2,3,4} and angles theta_1,...,theta_r
with
|theta_s| <= 10^4 a and |theta_1 + ... + theta_r| <= .01 t,
and a product discrepancy larger than Lambda.
========================================================================
10. SHUFFLE: MERGING THE COPIES
========================================================================
The discrepancy is a difference of products. We expand it.
10.1 Components.
Each array copy contributing a factor (one F or one G) is a bounded
combination of three component types:
(C1) a constant "latest valid" probe with a single wildcard position
(value 1),
(C2) the stage-j probe,
(C3) the complex conjugate of the stage-j probe (which also
conjugates the masks).
Each copy also carries the Gaussian weight v(m_s), the phase
exp(i theta_s m_s), and a scalar coefficient of modulus at most 1 (from
the floor, from 1/Y, from 1/Z_j, and from the projection c).
10.2 Categories.
Group the expanded products into categories. A category records
for each copy, its type (C1), (C2), or (C3);
the total weak ordering of all selected source positions across all
copies, including all equalities, i.e. all ties.
The number P of nonempty categories satisfies log P = o(B), uniformly.
Why. Once the types are fixed and all exact adjacencies inside each copy
are honoured, each copy contributes at most 4^{j+1} contiguous blocks.
Specify the weak order of the block starts: at most m! 2^m possibilities
for m <= 4^{j+1} blocks. Then process blocks in start order. A new start
either lies strictly beyond the maximum end index seen so far, or its
offset from that maximum is an integer between 0 and k-1. This
determines the entire weak ordering. The total number of categories is
at most 81 (4k+1)^{4^{j+1}}, so
log P = O( 4^j log( 4k + 4^{j+1} + 2 ) ) = O(4^j b) = o(B),
including the case k = 1.
10.3 Merged tuples.
Fix a category. All of its exact adjacency constraints act across
adjacent distinct slots of the merged order: if an old exact edge
spanned an intervening distinct slot, the category would be empty.
Collect these as the merged exact edges, and declare the rest loose.
The merged order has k' <= 4k slots and at most as many blocks as the
sum of the old block counts (each old block is contained in one merged
block, and each merged block contains at least one old block). Tie
bit-predicates are multiplied. Retain every old mask unchanged; each
contributes at most 4 copies here. The wildcard component contributes
no mask.
Coverage. Every loose merged edge is covered either by an inherited
mask or by the new mask we are about to add. Consider a copy with a slot
at the left endpoint of a loose merged edge. If that slot is the copy's
latest, it is not the overall latest (since the copy is not the
component chosen last), so it contributes a new covering interval. If it
is not the copy's latest, then the copy has a consecutive pair
straddling the edge; that pair was loose in the copy, so an old mask
already covered it.
Latest component and differences. Choose the component whose latest
selected index is last in the category, with ties broken
deterministically, and relabel it as s = r. Let Y_0 be the overall
latest offset. Put
g_r = 0,
g_s = Y_0 - m_s >= 0 for s < r,
g_bar = (g_1 + ... + g_r)/r,
Q(g) = sum_s (g_s - g_bar)^2.
10.4 Phase bookkeeping.
The product of phases factorizes:
prod_{s=1}^r exp(i theta_s m_s)
= exp( i (theta_1 + ... + theta_r) Y_0 )
* prod_{s<r} exp( -i theta_s g_s ).
So the outer phase carries the tiny sum theta_1 + ... + theta_r, and the
inner phases live in the differences g_s. The inner phases become a mask
of level j. This is the only place where the relaxation of the angle
budget pays off.
========================================================================
11. CENTER FOURIER EXPANSION: AN ADMISSIBLE MASK
========================================================================
The Gaussian weight is not yet a mask in the sense of (M1) and (M2). We
convert it.
Momentarily vary Y_0 as a real variable and define
W_g(Y_0) = exp(t Y_0) * prod_{s=1}^r v(Y_0 - g_s).
CLAIM 11.1. There is a centred Gaussian probability nu(d omega) with
variance r / ((1 - zeta) sigma^2), where zeta = 10^{-5}, such that
W_g(Y_0) = exp(1.203 B) int_R exp(i omega Y_0) K(g, omega) nu(d omega),
where
K(g, omega) = exp(-1.203 B) * W_g_hat(omega) / (d nu / d omega),
and W_g_hat is the Fourier transform
W_g_hat(omega) = (2 pi)^{-1} int_R W_g(Y_0) exp(-i omega Y_0) dY_0.
Moreover, for all real omega and all g_s >= 0,
|K(g, omega)| <= exp( - epsilon t (g_1 + ... + g_r) ). (11.1)
Proof. Scale:
z = (Y_0 - g_bar) sqrt(r) / sigma
rho = t sigma / sqrt(r)
u' = omega sigma / sqrt(r)
d_s = (g_s - g_bar) / sigma
p_ell(x) = sum_{k=0}^{ell} (x^2/2)^k / k!.
Shifting the contour by -i u', which is legal because the integrand is a
polynomial times a Gaussian, gives
W_g_hat = (2 pi)^{-1/2}
* exp( (t - i omega) g_bar - Q(g)/(2 sigma^2)
+ rho^2/2 - (u')^2/2 - i rho u' )
* E_{X ~ N(rho,1)} prod_s p_ell( (X - i u')/sqrt(r) - d_s ).
The product of polynomials is at most zeta^{-r ell} times
exp( (zeta/2)( X^2 + (u')^2 + Q(g)/sigma^2 ) ). Integrating X gives
|W_g_hat(omega)| / (density of nu at omega)
<= zeta^{-r ell}/(1 - zeta)
* exp( t g_bar - (1 - zeta) Q(g)/(2 sigma^2)
+ 2.4 B / (2 r (1 - zeta)) ).
Now Q(g) >= r g_bar^2/(r-1), because g_r = 0. So the exponent above,
plus the epsilon t sum_s g_s term, which equals epsilon t r g_bar, is at
most
(1 + r epsilon)^2 (2.4 B) (r-1) / (2 (1 - zeta) r).
For r = 4 this is at most .9001 B. Add 2.4B/(2r(1-zeta)) <= .3001 B and
the negligible log term; the total is below 1.203 B. Hence (11.1).
Interpretation. K is a Gaussian bump in g. For now it decays like
exp(-epsilon t |g|_1), which is the (M1) behaviour, before we make it
compact.
========================================================================
12. RADIAL COMPACTIFICATION
========================================================================
12.1 Applying the expansion.
Each fixed-omega discrepancy contribution, after the expansion, is a
merged sum with new outer factor
exp( -t Y_0 + i (omega + theta_1 + ... + theta_r) Y_0 )
and inner phases exp(-i theta_s g_s), plus K(g, omega). Its modulus is
at most 2 Z_{j+1}/(1 - exp(-t)) <= exp(10^{-8} B), apart from the outer
factor exp(1.203 B), because every merged loose edge has coverage and
(11.1) applies. Latest offsets Y_0 < 0 cancel exactly, leaving Y_0 >= 0
and hence a geometric sum from the outer damping.
Where |omega + theta_1 + ... + theta_r| > 5t, the Gaussian factor is at
most
exp(1.2031 B) * 2 exp( -(4.99)^2 (1 - zeta) 2.4 B / 8 )
= o( exp(-6B) ),
uniformly over categories, after multiplying by P = exp(o(B)). So only
the range |omega + theta_1 + ... + theta_r| <= 5t matters.
12.2 Making the mask compact.
We must replace K(g, omega), for the retained omega, by an admissible
K^c(g, omega) with
| K(g) - K^c(g) | <= exp(-20B) exp( -epsilon t sum_s g_s ). (12.1)
Everything is uniform on |omega| <= 5.01 t and measurable in omega.
Fix omega. Put R = R_j; recall R/t = (j+1)^2 M_j/L_j >= 1. Choose a
large absolute constant H_0 and set
W = H_0 B / t
tau^2 = W / R
U(x) = (1/(tau sqrt(2 pi)))
* int_{-W}^{W} exp( -(x - v)^2/(2 tau^2) ) dv.
So U is a smoothed indicator of the interval [-W, W].
CLAIM 12.1.
| K(g) - K(g) prod_{s<r} U(g_s) |
<= (1/4) exp(-20B) exp( -epsilon t sum_s g_s ) for g_s >= 0.
Proof. If all g_s <= W/2, use
1 - U(g_s) <= 2 exp( -W^2/(8 tau^2) ) = 2 exp(-R W/8)
= 2 exp(-H_0 B/8),
and then (11.1). If some g_s > W/2, use Q(g) >= W^2/(4r), by Cauchy
together with g_r = 0. That makes the Gaussian factor in the derivation
of (11.1) at most exp(-c W^2/sigma^2) = exp(-c H_0^2 B/2.4); and
0 <= U <= 1.
Now form the entire function of d' = r-1 variables
J_cal(g) = exp( R sum_s g_s ) K(g) prod_{s<r} U(g_s),
and its inverse transform
J_cal_hat(lambda) = (2 pi)^{-d'}
int_{R^{d'}} J_cal(g) exp(-i lambda . g) dg.
CLAIM 12.2. There are absolute constants C_2, C_3, allowed to depend on
H_0, with
int |J_cal_hat(lambda)| dlambda <= exp( C_2 R W ), (12.2a)
int_{|lambda|_inf > C_3 R} |J_cal_hat(lambda)| dlambda
<= (1/4) exp(-20B). (12.2b)
Proof. Continue the integral formula for K to complex
g = x - i y, with x and y real vectors in R^{d'} and the last coordinate
set to zero. Using Re Q(x - i y) = Q(x) - Q(y), one gets
|K(x - i y)| <= exp( C_4 B + C_4 t sum_{s<r} (|x_s| + |y_s|)
+ C_4 |y|_2^2/sigma^2 ).
Also
|U(x_s - i y_s)| <= (2W/(tau sqrt(2 pi)))
* exp( ( y_s^2 - max(|x_s| - W, 0)^2 ) / (2 tau^2) ).
Shifting g = x - i y in J_cal_hat by y, which is legal because of
Gaussian tail decay in every real direction, and using
max(|x_s| - W, 0)^2 >= x_s^2/2 - W^2
to complete squares, gives
|J_cal_hat(lambda)| <= tau^{d'}
exp( C_5 R W + C_5 |y|_2^2 / tau^2 - lambda . y ).
Take y = tau^2 lambda/(2 C_5). This leaves a Gaussian envelope in
lambda. Integrating in lambda proves (12.2a) and (12.2b).
Finally define
K^c(g) = exp( -R sum_s g_s )
int_{|lambda|_inf <= C_3 R} exp( i lambda . g )
J_cal_hat(lambda) dlambda.
Fourier inversion, valid because both J_cal and its transform are
integrable, together with Claim 12.1 and (12.2b), gives (12.1).
Combining (11.1) and (12.1) shows K^c also obeys the quadrant decay
|K^c| <= (some absolute constant) exp(-epsilon t sum g_s), and (12.2b)
gives the complex-measure representation with |xi|_inf <= C R_j and
total variation at most exp(C (R_j/t_j) B_j), after enlarging C once.
This completes the construction of stage j+1 from stage j.
========================================================================
13. STOPPING: FROM MASKS TO PHASE PARAMETERS
========================================================================
We now have a_j n <= .001, j <= b, and the pivot point
z_0 = exp(-a_j + i theta).
13.1 One frequency choice.
Expand all retained masks by their compact representations (M2). The
total variation of the product measure is at most
exp( C sum_{h<j} 4^{j-h-1} (R_h/t_h) B_h )
<= exp( O( ell_* 8^b (b+1)^6 ) ). (13.1)
Check. We have R_h/t_h = (h+1)^2 2^{b-h} and
B_h = 10^100 ell_* 8^h (b+1)^3, so
4^{j-h-1} (R_h/t_h) B_h
<= 4^{j-1} 2^b (h+1)^2 10^100 ell_* (b+1)^3,
and summing over h <= j <= b gives 8^b (b+1)^3 times a polynomial in b.
At every set of frequencies in the product integral, the difference
contribution is still a finite sum. So by the pivot magnitude and
(13.1), some frequency choice has an unintegrated discrepancy of
magnitude at least
exp( -O( ell_* 8^b (b+1)^6 ) ).
Fix it.
13.2 Reading off a pure product probe.
A mask now contributes, for each of its differences g_s, a number
exp( -R_h g_s + i (beta_s + xi_s) g_s ),
with |xi_s| <= C R_h and |beta_s| <= 10^4 a_h.
Assign these factors edgewise. Each g_s is a sum of consecutive gap
distances in the merged tuple, so each factor becomes a product of
per-edge factors. Drop factors on exact edges; they are scalars of
modulus at most 1 and nonzero. What remains is a pure product probe:
bit predicates, exact constraints, and, on each loose edge e, a complex
number
z_e^* = exp( -eta_e + i phi_e ).
CLAIM 13.1. There is an absolute c_0 in (0,1] with
kappa(z_e^*) := Re ( (1 + z_e^*)/(1 - z_e^*) )
= (1 - |z_e^*|^2) / |1 - z_e^*|^2
>= c_0 q. (13.2)
Proof. Let the sum over covering differences be counted with
multiplicity. Because every loose edge is covered, eta_e, the sum of the
relevant R_h, is strictly positive. The frequency perturbations sum to
at most C eta_e in absolute value, and the beta's satisfy
(sum beta_s)^2 <= eta_e sum_h (10^4 a_h)^2/R_h
<= C_6 (eta_e/q) sum_h 1/(h+1)^2
= O( eta_e/q ),
using at most 3 * 4^{j-h-1} <= M_h coverings per level h. When
eta_e <= 1, we have
1 - |z_e|^2 = 1 - exp(-2 eta_e) >= (1 - e^{-2}) eta_e
and |1 - z_e|^2 <= C_7 eta_e/q, because
|1 - z|^2 = 1 - 2 exp(-eta) cos phi + exp(-2 eta)
~ eta^2 + phi^2 <= C (eta^2 + eta/q).
When eta_e > 1, use |1 - z_e| <= 2 and 1 - |z_e|^2 >= 1 - e^{-2}.
For the root itself, kappa(z_0) >= 4/p. Because |1 - z_0| <= 6 a_j and
1 - |z_0|^2 >= a_j, and a_j <= .001/n, we get
1/(36 a_j) >= 1000 n/36 >= 4 n >= 4/p,
since p >= 1/n.
========================================================================
14. INTERPOLATION AND THE TRACE STATISTIC
========================================================================
14.1 A one-parameter family of edge arguments.
Put
S = 4/(p c_0),
and let varsigma be a complex variable in the strip
1/2 <= Re varsigma <= S.
For each loose edge e, write
alpha_e + i delta_e = (1 + z_e^*)/(1 - z_e^*).
By (13.2), alpha_e = kappa(z_e^*) >= c_0 q. For each exact edge, set the
target z_e^* = 0, i.e. alpha_e = 1 and delta_e = 0. Define
z_e(varsigma) = (alpha_e varsigma + i delta_e - 1)
/ (alpha_e varsigma + i delta_e + 1).
Then z_e(1) = z_e^* and
kappa( z_e(varsigma) ) = alpha_e Re varsigma.
At Re varsigma = S this is at least c_0 q S = 4 q/p.
14.2 The interpolating function.
Let K_ex = k - A be the number of exact edges. Define G(varsigma) as the
difference sum at z_0 with the following substitutions:
each exact-edge selector on source distance Delta becomes
z_e(varsigma)^{Delta - 1}, so Delta = 1 gives 1, as before;
each loose edge at source distance Delta uses z_e(varsigma)^Delta;
the whole sum is normalized by n^{-A} (1 + varsigma)^{-2 K_ex};
all negative latest-offset terms are cancelled.
So at an exact edge the weight relaxes from a strict constraint to a
geometric sum over possible source distances, and the value 1 of that
sum at varsigma = 1 recovers the original constraint.
CLAIM 14.1.
(a) |G(varsigma)| <= 2 on the entire strip.
(b) At varsigma = 1,
log |G(1)| >= -A log n - 2 K_ex log 2
- O( ell_* 8^b (b+1)^6 )
= -O( ell_* 8^b (b+1)^6 ). (14.1)
Proof of (a). For an exact edge, the sum over Delta >= 1 of
|z_e|^{Delta - 1} equals 1/(1 - |z_e|). Since alpha = 1 and delta = 0,
1 - |z_e|^2 = 4 Re varsigma / |1 + varsigma|^2,
so this costs at most 1 against the normalization |1 + varsigma|^{-2},
per exact edge. For the other edges, use at most n choices per gap and
modulus at most 1.
Proof of (b). Use k <= 2 P_0 4^b = O(8^b), A <= 4^b, and
log n <= ell_* 2^b.
14.3 Hadamard three-lines.
Apply the three-lines theorem on the vertical strip
1/2 <= Re varsigma <= S, with bounds |G| <= 2 on both vertical edges.
The right-line weight of the point Re varsigma = 1 is
(1/2)/(S - 1/2). So there is a point varsigma of the strip with
Re varsigma = S and
|G(varsigma)| >= exp( -O( p^{-1} ell_* 8^b (b+1)^6 ) ). (14.2)
14.4 The disk of radius p centred at q.
The following algebraic fact is the bridge from analytic parameters to
actual traces.
FACT (product ratio). For u and v in the open unit disk with
kappa(u), kappa(v) > 0,
|1 - u v|^2 / (1 - |u v|^2)
<= |1 - u|^2/(1 - |u|^2) + |1 - v|^2/(1 - |v|^2). (14.3)
Proof: 1 - uv = (1 - u) + u(1 - v), then the weighted
Cauchy-Schwarz ineuality with weights (1 - |u|^2) and
|u|^2 (1 - |v|^2), whose sum is exactly 1 - |uv|^2.
Apply (14.3) with u = z_0, which has kappa >= 4/p, and v =
z_e(varsigma), which has kappa >= 4q/p. The right side is at most
p/4 + p/(4q) <= p/(2q).
So, for y_e = z_0 z_e(varsigma),
|1 - y_e|^2 / (1 - |y_e|^2) <= p/(2q).
Use the identity
p^2 - |y - q|^2 = p(1 - |y|^2) - q |1 - y|^2.
Combining with the previous line gives, for y = y_e,
p^2 - |y_e - q|^2 >= (p/2)(1 - |y_e|^2)
>= (p/2)(1 - |z_e|^2).
The same computation applies to z_0 itself.
So z_0 and y_1,...,y_{k-1} all lie in the closed disk of centre q and
radius p. Define
lambda_e = (y_e - q)/p, lambda_0 = (z_0 - q)/p.
All have modulus at most 1. Moreover
1 - |lambda_e| >= (1 - |z_e|^2)/(4p). (14.4)
Proof: p - |y_e - q| >= (p^2 - |y_e - q|^2)/(2p) >= (1 - |z_e|^2)/4,
then divide by p.
14.5 The trace statistic.
Given a trace U of the extended input, define
T(U) = n^{-A} (1 + varsigma)^{-2 K_ex}
* sum over 0 <= s_1 < ... < s_k < |U|
of lambda_0^{s_1}
* prod_{e=1}^{k-1} lambda_e^{ s_{e+1} - s_e - 1 }
* prod_{e=1}^{k} chi_e( U_{s_e} ).
CLAIM 14.2. |T(U)| <= 1 for every trace.
Proof. Expand the sum over the exact-edge gaps. For an exact edge the
geometric weight over possible source distances Delta is
sum_{Delta >= 1} |lambda_e|^{Delta - 1} = 1/(1 - |lambda_e|). By (14.4)
together with 1 - |z_e|^2 = 4 S/|1 + varsigma|^2,
1 - |lambda_e| >= (1 - |z_e|^2)/(4p) = S/(p |1 + varsigma|^2).
So the geometric sum is at most p |1 + varsigma|^2/S. After its share
|1 + varsigma|^{-2} of the normalization, each exact edge contributes at
most
p/S = p^2 c_0/4 <= 1.
The loose edges, of which there are A-1, contribute at most n each,
because there are at most n possible gap lengths. The first index s_1
ranges over at most n values. Together with the normalization n^{-A},
the choices of the A block starts cost at most 1.
The bit predicates are bounded by 1.
14.6 The expectation gap.
Condition on the survival of the source symbols at the tuple positions.
Each survives independently with probability p. If all k chosen positions
survive, the symbols observed there are exactly x's bits at those
positions. The counts of survivors before s_1 and between consecutive
s_e have independent binomial laws.
By linearity of expectation, and the identity
q + p lambda_0 = z_0,
q + p lambda_e = z_0 z_e for each loose edge e,
together with the fact that at an exact edge the exponent is zero (there
are no survivors strictly between two forced-consecutive positions, so
both sides contribute 1), we get
E_x T - E_y T = p^k z_0^{d - k + 1}
* ( prod over loose e of z_e(varsigma)^{-1} )
* G(varsigma).
The powers of z_0 and z_e reorganize the inner sum into exactly
G(varsigma) times an overall factor, via the elementary identity
z_0^{i_1} prod_{e=1}^{k-1} (z_0 z_e)^{i_{e+1}-i_e-1}
= z_0^{i_k-k+1} prod_{e=1}^{k-1} z_e^{i_{e+1}-i_e-1},
which shows where each power of z_0, each power of a z_e, and each
factor of G(varsigma) comes from.
Modulus of the prefactor. We have
|z_0|^{d - k + 1} >= exp(-a_j n) >= exp(-.001),
because a_j n <= .001 and d - k + 1 <= n;
|prod over loose e of z_e^{-1}| >= 1,
because each |z_e| < 1 at Re varsigma = S, and |z_e| > 0 by (14.2);
p^k >= exp(-k log(1/p)) >= exp(-O(8^b/p)),
because k = O(8^b) and log(1/p) <= 3/p.
So, by (14.2),
|E_x T - E_y T| >= exp( -O( p^{-1} ell_* 8^b (b+1)^6 ) ). (14.5)
========================================================================
15. COUNTING TRACES
========================================================================
We now have a bounded statistic T with an expectation gap (14.5). That
is an information statement. To turn it into a reconstruction rule, use
Hoeffding and a union bound.
15.1 One pair.
Fix x and y. Let delta = |E_x T - E_y T|. By Claim 14.2, T takes
values in [0,1] after the harmless shift T -> (T + 1)/2, which does not
change the gap.
Hoeffding's ineuality: with N independent traces, the empirical
average of T deviates from its mean by delta/2 with probability at most
2 exp(-N delta^2/2). So
N = O( delta^{-2} log(1/eps) ).
With eps = e^{-2n} this is
exp( O( p^{-1} ell_* 8^b (b+1)^6 ) + log(2n) ),
and since the first term grows like a fixed power of log n while
log(2n) does not, the whole exponent is O( p^{-1} ell_* 8^b (b+1)^6 ).
So this many traces distinguish x from y with error at most e^{-2n} per
pair, uniformly over pairs.
15.2 All pairs at once.
There are 2^n binary strings of the doubled length. By the union bound
over all 2^n challengers y for the true x,
total error <= 2^n * e^{-2n} = (2/e^2)^n < 1/3 for n >= 1.
So that many traces suffice to recover x from its traces with
probability at least 2/3, simultaneously for every x.
15.3 The final bound.
Substitute the constants:
N = exp( O( p^{-1} ell_* 8^b (b+1)^6 ) ).
Since 2^b <= 2X and ell_* X <= 1 + 2 log n, this is
exp( O( p^{-1} (log n) X^2 (1 + log X)^6 ) ),
which is Theorem 2, and in turn
exp( O( p^{-1} (log n)^3 (1 + log log(2n))^6 ) ),
which is Theorem 1.
For completeness, the steps in detail. We used the doubled length n
throughout. Switching back to the original length changes the exponent
by at most an absolute factor, since doubling n only changes log n by
an additive constant.
The exponent is p^{-1} ell_* 8^b (b+1)^6 times an absolute constant.
We have
2^b <= 2 X (by the definition of X),
ell_* X <= 1 + 2 log n (since ell_* = 1 + log(1/q) and q >= 1/n).
So
p^{-1} ell_* 8^b (b+1)^6
<= p^{-1} ell_* (2X)^3 (1 + log_2(1 + log n/ell_*) + 1)^6
<= O( p^{-1} ell_* X^3 (1 + log X)^6 )
<= O( p^{-1} (log n) X^2 (1 + log X)^6 ),
because ell_* X <= 3 log n for n >= 3. This is Theorem 2.
The cruder form, Theorem 1, follows because X <= 1 + log n, as
ell_* >= 1. Then X^2 (1 + log X)^6 <= (log n)^2 (1 + log log(2n))^6
times an absolute constant, and multiplying by log n gives
(log n)^3 (1 + log log(2n))^6.
The polynomial case, Theorem 3. Suppose q <= n^{-epsilon} for fixed
epsilon > 0. If q = 0, one trace suffices. Otherwise
ell_* >= 1 + epsilon log n, so X <= 1 + 1/epsilon while p >= 1/2 for
large n. Hence the exponent is at most O(log n), and polynomially many
traces suffice.
Two honest remarks about this counting.
The bound counts how many traces suffice in principle. The rule that
uses them is allowed to be infeasible to compute: it may need to
search over all candidates, evaluate the test on each pair, and so
on. That is consistent with the statement of the theorem, which only
speaks about information.
An equivalent, more elementary way to say the same thing: the
expectation gap of a bounded real test implies, by total variation, a
gap at least half as large for an optimal indicator test. The
optimal indicator test for a binary pair is simply the one that
counts how many times each trace appears as a subsequence of each
candidate. Each trace word's probability under x is its embedding
count as a subsequence of x, times a length-dependent factor that is
the same for x and y. So the sample complexity already contains the
subsequence-count test, and the bounded statistic is only an
analytical device for proving the gap.
========================================================================
16. GENERAL ALPHABETS
========================================================================
Let Sigma be the alphabet, with |Sigma| <= n after the small-case
reduction.
16.1 A preliminary pass to see all symbols.
Collect ceil(p^{-1} log(12n)) preliminary traces. By a union bound over
the n source positions, with probability at least 11/12 every source
symbol appears somewhere in these traces. Indeed a given position
survives none of them with probability
(1-p)^{ceil(p^{-1} log(12n))} <= exp(-log(12n)) = 1/(12n),
and there are n positions.
Let S be the set of observed symbols. Always |S| <= n.
16.2 Reducing to binary.
Condition on success.
On independent fresh traces, run the binary test for every ordered pair
of candidate strings in S^n. For a pair (a,b) whose first disagreement
is at position i, map symbols to bits by the indicator of the symbol
a_i, that is, send a_i to 0 and every other symbol to 1 (or do the
symmetric thing). Then a and b project to two distinct binary strings,
and the binary pair test distinguishes those projections, hence a from
b.
The Hoeffding sample size is multiplied by O(log n) so that the
per-pair error is at most e^{-3n log n}. Since there are at most
n^{2n} = e^{2n log n} ordered pairs, the union bound over all pairs gives
total error at most e^{-n log n} < 1/3. So a log n factor on top of the
same sample budget works.
Output a candidate that wins all its comparisons, or declare failure if
there is none.
The remaining cases are the bounded-size ones and the cases
min(p, 1-p) < 1/n; those are handled directly by the full-length-trace
strategy of section 6.1, and the increase does not affect the stated
asymptotic budget.
This extension concerns sample complexity only.
========================================================================
20. REFERENCES
========================================================================
Burudgunte, Valiant, Wang. "A polynomial-time algorithm for the
deletion channel." These are the BVW references in the text.
Chase. "Separating words and trace reconstruction." The
period-avoiding word base.
Chen, De, Lee, Servedio, Sinha. The low-deletion setting, cited as
CDLS.
Fine and Wilf. "Uniqueness theorems for periodic functions." The
period lemma.
Hoeffding. "Probability ineualities for sums of bounded random
variables." The concentration used in section 15.
Rubinstein. "Explicit separating words." The period and separation
ideas used in section 8.
For later developments, see the recent work on trace reconstruction
problem. arXiv:2207.11489.
========================================================================
END
========================================================================