simplicity we will combine these bits with the two uncontrolled evolutions,
giving a total of $(K+1)n_{t}+2$. Therefore we implement exponentials of $A$
controlled by the time registers and perform $\mathrm{PREP}_{B}$ and
$\mathrm{SEL}_{B}$ to block-encode $B$ with cost

$[(K+1)n_{t}+2]\ \textrm{controlled-}e^{isA}+2K\ \mathrm{PREP}_{B}+K\
\mathrm{SEL}_{B}.$ (308)

Note that this only represents a upper bound on the cost and can often be
improved by using specific properties of the target Hamiltonian; see Section
III for improvements for performing the first quantized quantum chemistry
simulation. This implements the interaction picture algorithm for time $\tau$.

The success amplitude of this procedure can be made larger than $1/2$. To make
the amplitude equal to $1/2$, we introduce an ancilla register with state

$\cos\theta\lvert 0\rangle+\sin\theta\lvert 1\rangle,$ (309)

so that the success amplitude is now multiplied by $\cos\theta$. Letting
$n_{\theta}$ be the number of bits used to represent $\theta$, we can then
prepare this state with $n_{\theta}-3$ Toffoli gates.

To boost the success amplitude close to unity, we need to perform a single
step of oblivious amplitude amplification, which triples the implementation of
the truncated Dyson series and introduces two additional reflections.
Specifically,

  1. 1.

We implement the pre-amplified truncated Dyson series and prepare the $\theta$
state with a complexity as described above.

  2. 2.

We perform a reflection on $n_{k}+K\left(n_{t}+n_{B}\right)+2$ qubits.

  3. 3.

We invert the truncated Dyson series and unprepare the $\theta$ state.

  4. 4.

We perform a reflection again on $n_{k}+K\left(n_{t}+n_{B}\right)+2$ qubits.

  5. 5.

We implement the truncated Dyson series without preparing the $\theta$ state.

The $+2$ accounts for the qubit rotated by $\theta$ here, as well as a qubit
rotated for the preparation of the equal superposition state for preparing
$k$. Each reflection can be implemented with Toffoli complexity

$n_{k}+K\left(n_{t}+n_{B}\right),$ (310)

where $n_{B}$ is the number of ancilla qubits for block-encoding operator $B$.

We now analyze the normalization factor introduced by the above steps. Such a
factor will come from three sources: (i) the preparation of the Dyson-
expansion coefficients; (ii) the preparation of the equal superposition state
on the $\lvert k\rangle$ register; and (iii) the implementation of the block
encoding of operator $B$. Same as in Section III.4, we have that the success
probability of Item (ii) is $p={\rm Ps}\left(\Sigma(0),b_{r}\right)$, whereas
Items (i) and (iii) together contribute a normalization factor of
$\beta=\sum_{k=0}^{K}c^{k}/(d^{k}k!)\leq\exp(c/d)$. To boost this close to
unity with a single step of oblivious amplitude amplification, we choose
coprime integers $c$ and $d$ so that

$\frac{{\rm Ps}\left(\Sigma(0),b_{r}\right)}{\exp(c/d)}\geq\frac{1}{2},$ (311)

which implies

$\mathcal{N}=\frac{\lambda_{B}td}{c}.$ (312)

The error of our circuit implementation is due to the truncation of the Dyson-
series expansion ($\epsilon_{K}$) and the finite binary representation of the
time steps ($\epsilon_{t}$), as well as the preparation of $\theta$ state
($\epsilon_{\theta}$). With $n_{\theta}$ digits to represent the angle
$\theta$, the error in the amplitude $\cos\theta$ is upper bounded by

$\epsilon_{\theta}=\frac{\pi}{2^{n_{\theta}}},$ (313)

which causes an overall error of at most $\epsilon_{K}+\epsilon_{\theta}$ in
the pre-amplified truncated Dyson series. After a single step of oblivious
amplitude amplification, the error is bounded as (Childs  _et al._ , 2018,
Lemma G.2)

$(\epsilon_{K}+\epsilon_{\theta})\frac{(\epsilon_{K}+\epsilon_{\theta})^{2}+3(\epsilon_{K}+\epsilon_{\theta})+4}{2}.$
(314)

To ensure that the entire simulation has accuracy $\epsilon$, we require that

$\mathcal{N}(\epsilon_{K}+\epsilon_{\theta})\frac{(\epsilon_{K}+\epsilon_{\theta})^{2}+3(\epsilon_{K}+\epsilon_{\theta})+4}{2}+\mathcal{N}\epsilon_{t}\leq\epsilon.$
(315)

We summarize the cost of the circuit implementation in the following theorem.

######  Theorem 12.

Using the interaction picture algorithm, the time evolution under the
Hamiltonian $H=A+B$ can be simulated for time $t$ that is a rational multiple
of $1/\lambda_{B}$ to within error $\epsilon$ by choosing positive integers
$K$, $n_{t}$, $n_{\theta}$, $b_{r}$, coprime $c$ and $d$, and using
$3[\mathcal{N}(K+1)n_{t}+2]\ \textrm{controlled-}e^{isA}$, $6\mathcal{N}K\
\mathrm{PREP}_{B}$, $3\mathcal{N}K\ \mathrm{SEL}_{B}$, and a number of Toffoli
gates

$\displaystyle
6\mathcal{N}(3n_{k}+2b_{r}-9)+3\mathcal{N}\left(-1+\sum_{k=0}^{K-1}\lceil\log(\Sigma(k))\rceil\right)+6\mathcal{N}Kn_{t}+12\mathcal{N}(n_{t}-1){\rm
St}\left(K\right)$
$\displaystyle+6\mathcal{N}(K-1)(n_{t}-1)+2\mathcal{N}(n_{\theta}-3)+2\mathcal{N}\left[n_{k}+K(n_{t}+n_{B})\right],$
(316)

where ${\rm St}\left(K\right)$ is the cost of sorting $K$ items,
$\tau=c/(d\lambda_{B})$, $\mathcal{N}=t/\tau$ is an integer, and

$\displaystyle\frac{{\rm
Ps}\left(\Sigma(0),b_{r}\right)}{\exp(c/d)}\geq\frac{1}{2},$ (317)
$\displaystyle\mathcal{N}(\epsilon_{K}+\epsilon_{\theta})\frac{(\epsilon_{K}+\epsilon_{\theta})^{2}+3(\epsilon_{K}+\epsilon_{\theta})+4}{2}+\mathcal{N}\epsilon_{t}\leq\epsilon,$
(318)
$\displaystyle\epsilon_{K}=\exp(c/d)-\sum_{k=0}^{K}\frac{c^{k}}{d^{k}k!},\qquad\qquad\epsilon_{\theta}=\frac{\pi}{2^{n_{\theta}}}$
, (319)
$\displaystyle\epsilon_{t}=\left(\frac{2\|A\|}{\lambda_{B}}+\frac{\|A\|^{2}}{\lambda_{B}^{2}}\right)\left[2^{n_{t}}\left(e^{\lambda_{B}\tau/2^{n_{t}}}-1\right)-\lambda_{B}\tau(1+\lambda_{B}\tau/2^{n_{t}+1})\right].$
(320)

###### Proof.

The complete interaction picture simulation is indicated in Figure 7, with the
interaction picture operator to be amplified shown in Figure 8 and the clock
preparation shown in Figure 9. The total costs for each of the parts of the
simulation are shown in Table 8. The sum of the costs in entries 1 to 7 gives
the cost in (12). The costs in terms of the evolution under $A$ and the block
encoding of $B$ are given in entries 8 and 9 of Table 8.

Procedure |  Cost
---|---
The preparation of the equal superposition over $\Sigma(0)$ numbers; see (161). |  $6\mathcal{N}(3n_{k}+2b_{r}-9)$
Inequality tests used for the preparation of the superposition over $k$; see (306). |  $3\mathcal{N}\left(-1+\sum_{k=0}^{K-1}\lceil\log(\Sigma(k))\rceil\right)$
Controlled Hadamards for preparing the superpositions over times; see (166). |  $6\mathcal{N}Kn_{t}$
The sorting of the time registers; see (167). |  $12\mathcal{N}n_{t}{\rm St}\left(K\right)$
The differences of the times; see (151). |  $6\mathcal{N}(K-1)(n_{t}-1)$
The preparation of $\theta$ state; see (309). |  $2\mathcal{N}(n_{\theta}-3)$
The reflection cost for oblivious amplitude amplification; see (310). |  $2\mathcal{N}\left[n_{k}+K(n_{t}+n_{B})\right]$
The exponential of $A$ needs to be performed $K+1$ times; see (308). |  $3[\mathcal{N}(K+1)n_{t}+2]\ \textrm{controlled-}e^{isA}$
Perform the controlled block encoding of operator $B$ $K$ times; see (308). |  $6\mathcal{N}K\ \mathrm{PREP}_{B}+3\mathcal{N}K\ \mathrm{SEL}_{B}$
Table 8: The costs of the interaction picture algorithm for time evolution of
a general Hamiltonian. The costs are Toffoli costs unless otherwise specified
(the last two entries).

Provided $\lambda_{B}t$ is rational, for any integer $\mathcal{N}$,
$\lambda_{B}t/\mathcal{N}$ must be rational, and we can choose $c$ and $d$ to
satisfy $c/d=\lambda_{B}\tau$ for $\tau=t/\mathcal{N}$. Then (317) comes from
(311), and constrains the value of $c/d$. Since $c/d$ is set by $\mathcal{N}$,
we need a way to select $\mathcal{N}$ so that $c/d$ satisfies (317). Choosing
some value of $b_{r}$ we can find a lower bound $P_{\min}$ on ${\rm
Ps}\left(\Sigma(0),b_{r}\right)$ (for example as given in Appendix J below).
Then we can set $\mathcal{N}=\lceil\lambda_{B}t/\ln(2P_{\min})\rceil$, and
$\lambda_{B}\tau=\lambda_{B}t/\mathcal{N}$ is strictly less than
$\ln(2P_{\min})$. Then $c/d=\lambda_{B}\tau<\ln(2P_{\min})$ implies that (317)
is satisfied.

The error bound in (318) is from (315). The quantities $\epsilon_{K}$,
$\epsilon_{\theta}$ and $\epsilon_{t}$ are bounds on the error due to
truncating the Dyson series, the rotation on the ancilla qubit for the
amplitude amplification, and the time discretisation, respectively. The norm
of each successive term in the Dyson series is bounded as
$(\tau\lambda_{B})^{k}/k!$, and $\lambda_{B}\tau=c/d$. Therefore the error in
truncating the Dyson series to order $K$ can be bounded by the error in
truncating the exponential of $c/d$, giving the expression for $\epsilon_{K}$
in (319). The expression for $\epsilon_{\theta}$ in (319) is from (313). The
bound in (320) is from (293), where we have used $\lambda_{B}$ as the upper
bound on $\|B\|$. ∎

The algorithm will be most efficient when $c/d$ is a close rational
approximation of $\ln 2$, so $\mathcal{N}$ is minimized. This is because $c/d$
is upper bounded by $\ln[2{\rm Ps}\left(\Sigma(0),b_{r}\right)]$ according to
(317), ${\rm Ps}\left(\Sigma(0),b_{r}\right)$ is close to 1, and
$\mathcal{N}=\lambda_{B}td/c$. It would be expected that for a real
application one can choose the exact time to be simulated, so that $c$ and $d$
are relatively small (for example $9/13$ is less than but close to $\ln 2$).
If evolution for a specific $t$ is needed where $\lambda_{B}t$ is not
rational, we prepare an ancilla register in the state

$\cos\theta_{B}\lvert 0\rangle+\sin\theta_{B}\lvert 1\rangle$ (321)

and use the amplitude $\cos\theta_{B}$ to artificially adjust the
normalization factor $\lambda_{B}$. This completes our costing of the
interaction picture algorithm for generic Hamiltonians.

|
---|---
(a) Interaction picture algorithm with 3 steps. | (b) Oblivious amplitude amplification.
Figure 7: Quantum circuits for the interaction picture algorithm for
simulating generic Hamiltonians.

Figure 8: Circuit implementation of the pre-amplified operator ($K=2$).

Figure 9: Quantum circuit for preparing the clock state ($K=2$). Note that the
last time difference $M-1-m_{2}^{\prime}$ can be computed with a copying
operation followed by flipping the qubits.

##  Appendix J Bound on the success probability for preparing equal
superposition

Here we consider the probability of success for preparation of an equal
superposition state and place a lower bound on it for given $b_{r}$. The
probability for success was given in Sanders  _et al._ (2020), and is

${\rm Ps}\left(n,b_{r}\right)=\frac{n}{2^{\lceil\log
n\rceil}}\left[\left(1+\left(2-\frac{4n}{2^{\lceil\log
n\rceil}}\right)\sin^{2}\theta(n,b_{r})\right)^{2}+\sin^{2}(2\theta(n,b_{r}))\right],$
(322)

where $\theta(n,b_{r})$ is the angle of rotation on an ancilla qubit. It turns
out that there is a simple lower bound for the success probability, though the
expression for $\theta$ is adjusted slightly from the formula used before. The
result is as in the following theorem.

######  Theorem 13.

For all $n,b_{r}\in\mathbb{N}$, taking

$\displaystyle x$ $\displaystyle=\frac{n}{2^{\lceil\log n\rceil}},$ (323)
$\displaystyle\theta_{0}$
$\displaystyle=\arcsin\left(\frac{1}{2\sqrt{x}}\right),$ (324)
$\displaystyle\theta(n,b_{r})$ $\displaystyle=\frac{2\pi}{2^{b_{r}}}{\rm
round}\left[\frac{2\pi}{2^{b_{r}}}\left(\theta_{0}-\frac{(\pi/2^{b_{r}})^{2}(2x-1)}{\sqrt{4x-1}}\right)\right]$
(325)

gives the success probability for preparation of an equal superposition state
lower bounded as

${\rm Ps}\left(n,b_{r}\right)\geq
1-\left(\frac{3}{2}\right)^{2}\left(\frac{\pi}{2^{b_{r}}}\right)^{2}.$ (326)

###### Proof.

Now with $b_{r}=1$ or $b_{r}=2$ the result is trivial, because the right-hand
side is negative. We therefore aim to prove the result for $b_{r}\geq 3$.
Because the expression for ${\rm Ps}\left(n,b_{r}\right)$ depends only on the
ratio $x=n/2^{\lceil\log n\rceil}$, we can prove the lower bound for all
values of $x$ in the range $x\in(1/2,1)$, and that will give us the result for
all $n\in\mathbb{N}$. Then the probability of success can be written as

${\rm
Ps}\left(n,b_{r}\right)=x\left[\left(1+\left(2-4x\right)\sin^{2}\theta\right)^{2}+\sin^{2}(2\theta)\right].$
(327)

Let us take $\Delta\theta=\theta-\theta_{0}$. Using elementary trigonometric
identities, the probability of success can be rewritten as

$\displaystyle{\rm Ps}\left(n,b_{r}\right)$
$\displaystyle=1-4(1-x)(4x-1)\sin^{2}(\Delta\theta)$
$\displaystyle\quad+4(1-x)\left[2(4x-1-2x^{2})\sin^{2}(\Delta\theta)-(2x-1)\sqrt{4x-1}\sin(2\Delta\theta)\right]\sin^{2}(\Delta\theta).$
(328)

To prove lower bounds for $x\in(1/2,1)$, we use

$\displaystyle(2x-1)(3-x)$ $\displaystyle\geq 0$ $\displaystyle(2x-1)(3-x)+2$
$\displaystyle>0$ $\displaystyle 4(1-x)(2x-1)[(2x-1)(3-x)+2]$
$\displaystyle\geq 0$ $\displaystyle 4(1-x)(2x-1)[(2x-1)(3-x)+2]+1$
$\displaystyle>0$
$\displaystyle\left[2(4x-1-2x^{2})\right]^{2}-\left[(2x-1)\sqrt{4x-1}\right]^{2}$
$\displaystyle>0.$ (329)

Next, $4x-1-2x^{2}=(2x-1)(3/2-x)+1/2$ is positive for $x\in(1/2,1)$, and
$2x-1$ and $4x-1$ are non-negative. As a result, we have

$2(4x-1-2x^{2})>(2x-1)\sqrt{4x-1}.$ (330)

We also find that

$\sin^{2}\Delta\theta-\left[\sin(2\Delta\theta)-2\Delta\theta\right]=(\Delta\theta)^{2}+\mathcal{O}((\Delta\theta)^{3}),$
(331)

so it is non-negative for small $\Delta\theta$. Moreover this expression has
zeros of its derivative for $n\pi$ and $n\pi-\arctan(1/2)$. This means the
only zeros of its derivative for $\Delta\theta\in[-0.7,0.7]$ are at 0 and
$-\arctan(1/2)\approx-0.46$. This means that the function
$\sin^{2}\Delta\theta-\left[\sin(2\Delta\theta)-2\Delta\theta\right]$ can turn
around towards negative values for $\Delta\theta<-0.46$, but computing it at
$\Delta\theta=0.7$ gives about $0.000466$, which is still positive. This means
that for $\Delta\theta\in[-0.7,0.7]$ we have

$\sin^{2}\Delta\theta-\left[\sin(2\Delta\theta)-2\Delta\theta\right]\geq 0,$
(332)

because negative values would require an additional turning point. That
implies

$2(4x-1-2x^{2})\sin^{2}(\Delta\theta)\geq(2x-1)\sqrt{4x-1}\left[\sin(2\Delta\theta)-2\Delta\theta\right],$
(333)

which in turn implies

$2(4x-1-2x^{2})\sin^{2}(\Delta\theta)-(2x-1)\sqrt{4x-1}\sin(2\Delta\theta)\geq-2(2x-1)\sqrt{4x-1}\Delta\theta$
(334)

for $x\in(1/2,1)$ and $\Delta\theta\in[-0.7,0.7]$. This means that the
probability is lower bounded by

$\displaystyle{\rm Ps}\left(n,b_{r}\right)$ $\displaystyle\geq
1-4(1-x)(4x-1)\sin^{2}(\Delta\theta)-8(1-x)(2x-1)\sqrt{4x-1}\Delta\theta\,\sin^{2}(\Delta\theta)$
$\displaystyle=1-4(1-x)\left[(4x-1)+2(2x-1)\sqrt{4x-1}\Delta\theta\right]\sin^{2}(\Delta\theta).$
(335)

Next, we have for $x\in(1/2,1)$ and $\Delta\theta\in[-0.7,0.7]$

$\displaystyle 4(1-x)(3x-1)$ $\displaystyle\geq 0$ $\displaystyle
4x-1-3(2x-1)^{2}$ $\displaystyle\geq 0$ $\displaystyle\sqrt{4x-1}$
$\displaystyle\geq\sqrt{3}(2x-1)$ $\displaystyle 4x-1$
$\displaystyle\geq\sqrt{3}(2x-1)\sqrt{4x-1}$ $\displaystyle 4x-1$
$\displaystyle\geq-2(2x-1)\sqrt{4x-1}\Delta\theta.$ (336)

As a result, we can lower bound the success probability by

${\rm Ps}\left(n,b_{r}\right)\geq
1-4(1-x)\left[(4x-1)+2(2x-1)\sqrt{4x-1}\Delta\theta\right](\Delta\theta)^{2}.$
(337)

Note also that the derivative of this expression with respect to
$\Delta\theta$ is

$-8(1-x)\left[(4x-1)+3(2x-1)\sqrt{4x-1}\Delta\theta\right]\Delta\theta.$ (338)

Now $4x-1\geq\sqrt{3}(2x-1)\sqrt{4x-1}$ implies the expression in square
brackets must be non-negative for $\Delta\theta\in[-1/\sqrt{3},1/\sqrt{3}]$.
This implies that, provided we are considering
$\Delta\theta\in[-1/\sqrt{3},1/\sqrt{3}]$, we can lower bound the probability
of success in any region by the probabilities of success at the bounds of the
region.

Now, using the formula in the theorem for $\theta$, we have

$\Delta\theta+\frac{(\pi/2^{b_{r}})^{2}(2x-1)}{\sqrt{4x-1}}\in[-\pi/2^{b_{r}},\pi/2^{b_{r}}].$
(339)

Note that $\sqrt{4x-1}\geq\sqrt{3}(2x-1)$ implies

$\frac{2x-1}{\sqrt{4x-1}}\leq\frac{1}{\sqrt{3}}.$ (340)

Therefore the maximal value of $\Delta\theta$ with ${b_{r}}\geq 3$ is

$\frac{(\pi/8)^{2}}{\sqrt{3}}+\pi/2^{8}\approx 0.481734<\frac{1}{\sqrt{3}}.$
(341)

Because $\Delta\theta\in[-1/\sqrt{3},1/\sqrt{3}]$, we can lower bound the
probability by the lower bounds at the extremal points in the range of
$\Delta\theta$, which give

$\displaystyle{\rm Ps}\left(n,b_{r}\right)$ $\displaystyle\geq
1-4(1-x)(4x-1)\left(\frac{\pi}{2^{b_{r}}}\right)^{2}$
$\displaystyle\quad+(2x-1)^{2}(1-x)\left(\frac{\pi}{2^{b_{r}}}\right)^{4}\pm\frac{24(2x-1)^{3}(1-x)}{\sqrt{4x-1}}\left(\frac{\pi}{2^{b_{r}}}\right)^{5}+\frac{8(2x-1)^{4}(1-x)}{4x-1}\left(\frac{\pi}{2^{b_{r}}}\right)^{6}.$
(342)

Next, we have that $\sqrt{4x-1}\geq\sqrt{3}(2x-1)$ implies

$\displaystyle 20(2x-1)^{2}(1-x)\sqrt{4x-1}$ $\displaystyle\geq
20\sqrt{3}(2x-1)^{3}(1-x)$ $\displaystyle 20(2x-1)^{2}(1-x)$
$\displaystyle\geq\frac{24(2x-1)^{3}(1-x)}{\sqrt{4x-1}}$
$\displaystyle(2x-1)^{2}(1-x)\left(\frac{\pi}{2^{b_{r}}}\right)^{4}$
$\displaystyle\geq\frac{24(2x-1)^{3}(1-x)}{\sqrt{4x-1}}\left(\frac{\pi}{2^{b_{r}}}\right)^{5}$
(343)

because $\pi/2^{b_{r}}<1$. Hence we have

$(2x-1)^{2}(1-x)\left(\frac{\pi}{2^{b_{r}}}\right)^{4}\pm\frac{24(2x-1)^{3}(1-x)}{\sqrt{4x-1}}\left(\frac{\pi}{2^{b_{r}}}\right)^{5}\geq
0,$ (344)

and the sixth-order term in (J) is non-negative, so we obtain the lower bound

${\rm Ps}\left(n,b_{r}\right)\geq
1-4(1-x)(4x-1)\left(\frac{\pi}{2^{b_{r}}}\right)^{2}.$ (345)

Now $(1-x)(4x-1)$ has its maximal value of $9/16$ for $x=5/8$, so the lower
bound becomes

${\rm Ps}\left(n,b_{r}\right)\geq
1-\left(\frac{3}{2}\right)^{2}\left(\frac{\pi}{2^{b_{r}}}\right)^{2}$ (346)

as required. ∎

If we take $b_{r}=8$ for example, then the probability of success must be at
least $0.999661$. That is useful for the general interaction picture
simulations described in Appendix I, because then $\ln[2{\rm
Ps}\left(\Sigma(0),8\right)]>0.6928$. That is larger than $9/13$, so one can
use $c/d=9/13$ as a good approximation of $\ln 2$ and still satisfy (317).

##  Appendix K Simulating first quantized real space grid representations

Whether the goal is to perform non-Born-Oppenheimer dynamics or to solve the
electronic structure problem, in order to simulate molecular or material
systems on a quantum computer one must discretize the wavefunction and map it
to qubits. Perhaps the most obvious representation involves projecting the
wavefunction onto a real space grid and encoding the positions of each
particle explicitly in a quantum register. In this case, the computational
basis states that span the wavefunction of interest might be stored in the
quantum computer as $\lvert r_{1}r_{2}\cdots r_{\eta}R_{1}R_{2}\cdots
R_{L}\rangle$ where the $r_{i}$ and $R_{\ell}$ index the coordinates of the
$\eta$ electrons and $L$ nuclei, respectively, using binary registers with $N$
grid points (each requiring $\log N$ qubits). Thus, this encoding would
require $(\eta+L)\log N$ qubits.

Representing the potential part of the Hamiltonian is straightforward when
using a real space grid (the operators $U$, $V$ and $V_{\rm nuc}$ would take
the same diagonal forms as they do in (1)). However, representing the kinetic
operators will require some thought. One can see that if basis functions are
not overlapping in space (e.g., delta functions or finite elements with
disjoint support) the integrals from the Galerkin discretization would suggest
that the kinetic operator is also diagonal. To avoid this, one idea would be
to take a finite-difference approach and approximate the Laplacian using a
$k$-point stencil; this is exactly the approach taken in first quantization in
Kivlichan  _et al._ (2017) and described in second quantization in Appendix A
of Babbush  _et al._ (2018b). Another idea, similar to using what is referred
to as a “discrete value representation” (DVR) Light and Carrington Jr (2000),
would be to exactly define the kinetic operator in its eigenbasis and
approximate a transformation that relates the grid representation to the
eigenbasis of the kinetic operator. For example, if using a regularly spaced
grid then one might make the approximation that the quantum Fourier transform
(QFT) diagonalizes the kinetic operator. After all, the kinetic operator is
diagonal in the momentum basis (plane waves) and plane waves are obtained as
the continuous Fourier transform of delta functions. However, this is an
approximation when using a finite grid; the QFT is a discrete Fourier
transform (which is linear) and a linear combination of $N$ plane waves cannot
sum to a delta function. This sort of representation is employed by the first
paper to consider quantum simulations of molecular systems using a grid
representation by Kassal _et al._ Kassal  _et al._ (2008). Components of that
algorithm were assessed for fault-tolerance in the work of Jones _et al._ Cody
Jones  _et al._ (2012) but that work stopped short of estimating all of the
constant factors associated with the algorithm realization.

In addition to the necessity of the approximations already mentioned there are
other unfavorable characteristics of real space grid representations for use
in quantum simulation. For example, basis errors are no longer variationally
bounded when deviating from a Galerkin discretization and extrapolations to
the continuum limit might be less reliable. Furthermore, ways of softening the
potential and removing core electrons (akin to pseudopotentials for plane wave
methods) are less developed for grids. Thus, the majority of work on
simulating molecular systems focuses on Galerkin discretizations involving
basis functions rather than grids. Nonetheless, in this section we will
briefly describe how one can block encode a real space representation of the
Hamiltonian with the same asymptotic complexity as the main approaches pursued
in this paper. Since these grid representations are not periodic, it is
possible that they might converge faster than plane waves when simulating non-
periodic systems. Also, since the block encoding algorithms are entirely
different, it is possible that this approach might have lower constant
factors.

We can express the aforementioned first quantized real space Hamiltonian as

$\displaystyle H_{\rm BO}$
$\displaystyle=T+U+V+\frac{1}{2}\sum_{\ell\neq\kappa=1}^{L}\frac{\zeta_{\ell}\zeta_{\kappa}}{\left\|R_{\ell}-R_{\kappa}\right\|}\qquad\qquad
H_{\rm non-BO}=T+T_{\rm nuc}+U_{\rm non-BO}+V+V_{\rm nuc}$ (347)
$\displaystyle T$ $\displaystyle=\sum_{i=1}^{\eta}{\rm
QFT}_{i}\left(\sum_{p\in G}\frac{\left\|k_{p}\right\|^{2}}{2}\lvert
p\rangle\\!\langle p\rvert_{i}\right){\rm QFT}_{i}^{\dagger}\qquad T_{\rm
nuc}=\sum_{\ell=\eta+1}^{L+\eta}{\rm QFT}_{\ell}\left(\sum_{p\in
G}\frac{\left\|k_{p}\right\|^{2}}{2\,m_{\ell}}\lvert p\rangle\\!\langle
p\rvert_{\ell}\right){\rm QFT}_{\ell}^{\dagger}$ (348) $\displaystyle U$
$\displaystyle=-\sum_{i=1}^{\eta}\sum_{\ell=1}^{L}\sum_{p\in
G}\frac{\zeta_{\ell}}{\left\|R_{\ell}-r_{p}\right\|}\lvert p\rangle\\!\langle
p\rvert_{i}\qquad\quad U_{\rm non-
BO}=-\sum_{i=1}^{\eta}\sum_{\ell=\eta+1}^{L+\eta}\sum_{p,q\in
G}\frac{\zeta_{\ell}}{2\left\|r_{p}-r_{q}\right\|}\lvert p\rangle\\!\langle
p\rvert_{i}\lvert q\rangle\\!\langle q\rvert_{\ell}$ (349) $\displaystyle V$
$\displaystyle=\sum_{i\neq j=1}^{\eta}\sum_{p,q\in
G}\frac{1}{2\left\|r_{p}-r_{q}\right\|}\lvert p\rangle\\!\langle
p\rvert_{i}\lvert q\rangle\\!\langle q\rvert_{j}\qquad V_{\rm
nuc}=\sum_{\ell\neq\kappa=\eta+1}^{L+\eta}\sum_{p,q\in
G}\frac{\zeta_{\ell}\zeta_{\kappa}}{2\left\|r_{p}-r_{q}\right\|}\lvert
p\rangle\\!\langle p\rvert_{\ell}\lvert q\rangle\\!\langle q\rvert_{\kappa}$
(350)

where ${\rm QFT}_{i}$ is the usual quantum Fourier transform applied to
register $i$888Again, we emphasize that $T$ is only approximately given by the
expression involving the QFT. This relation is exact in the continuum limit.
For finite sized grids, it cannot be the case that the QFT completely
diagonalizes the momentum operator.. As in the main text, we have used atomic
units where $\hbar$ as well as the mass and charge of the electron are unity
and $\|\cdot\|$ denotes a function taking the 2-norm of its argument. In the
above expression $\ell$ and $\kappa$ index nuclear degrees of freedom; thus,
$R_{\ell}$ represents the positions of nuclei and $\zeta_{\ell}$ the atomic
numbers of nuclei. Furthermore, we have the following definition of grid
points and their frequencies in the dual space defined by the QFT:

$r_{p}=\frac{p\,\Omega^{1/3}}{N^{1/3}}\qquad\qquad k_{p}=\frac{2\pi
p}{\Omega^{1/3}}\qquad\qquad p\in G\qquad\qquad
G=\left[-\frac{N^{1/3}-1}{2},\frac{N^{1/3}-1}{2}\right]^{3}\subset\mathbb{Z}^{3}\,,$
(351)

where $\Omega$ is the volume of our simulation cell and $N$ is the number of
grid points in the cell. Again, this is the same representation used in the
first work on quantum simulating chemistry in first quantization, by Kassal
_et al._ Kassal  _et al._ (2008), well over a decade ago. We note that like
with the plane waves, it would be possible to use different grids for the
different particles but we avoid doing that here for simplicity.

We now show that one can simulate the Hamiltonian of (347) with the same
complexity as the approaches of this paper using either qubitization, or
interaction picture simulation frameworks. Here we will discuss only the Born-
Oppenheimer case for simplicity, but the techniques readily generalize to the
non-Born-Oppenheimer Hamiltonian. Focusing first on the qubitization approach,
the goal is to block encode the operators $T$, $U$ and $V$. The block encoding
of $T$ can be accomplished using almost the same techniques used for plane
waves. The only difference comes in the implementation of the $\mathrm{SEL}$
operator. Whenever the ancilla register flags that a term from $T$ should be
applied, one must apply the QFT, then phase the register by the appropriate
amount, then apply the inverse QFT. The QFT needs to be performed
independently on $\eta$ registers each of $\log N$ qubits, so has gate
complexity $\widetilde{\mathcal{O}}(\eta\log N)$, where here the
$\widetilde{\mathcal{O}}$ indicates terms double-logarithmic in $N$ are
omitted Cleve and Watrous (2000). That is complexity in terms of one and two-
qubit gates, and requires rotations up to $\mathcal{O}(\log N)$ bits of
precision. The complexity in terms of Toffoli gates is
$\widetilde{\mathcal{O}}(\eta\log^{2}N)$.

Here the block encoding of $U$ and $V$ will deviate from the other strategies
of this paper. Our approach will be to break up the terms as a sum of many
terms with equal magnitude coefficients that can be computed on-the-fly and
applied as phases, mirroring a strategy first introduced in Berry  _et al._
(2014b). First, we begin by re-expressing the $V$ operator as

$\displaystyle V$ $\displaystyle=\frac{v_{\rm max}}{M}\sum_{i\neq
j=1}^{\eta}\sum_{p,q\in G}\sum_{m=1}^{M}v\left(m,r_{p},r_{q}\right)\lvert
p\rangle\\!\langle p\rvert_{i}\lvert q\rangle\\!\langle q\rvert_{j}$ (352)
$\displaystyle v\left(m,r_{p},r_{q}\right)$
$\displaystyle=\begin{cases}1&\qquad\frac{m\,v_{\rm
max}}{M}<\frac{1}{2\left\|r_{p}-r_{q}\right\|}\\\ 0&\qquad{\rm
otherwise}\end{cases}$ (353) $\displaystyle v_{\rm max}$
$\displaystyle=\textrm{max}_{p\neq
q}\left(\frac{1}{2\left\|r_{p}-r_{q}\right\|}\right)=\frac{N^{1/3}}{2\,\Omega^{1/3}}\,,$
(354)

where in the inequality above we take advantage of the fact that the repulsive
Coulomb kernel is always greater than zero. From this definition we can see
that

$\left|\frac{1}{2\left\|r_{p}-r_{q}\right\|}-\frac{v_{\rm
max}}{M}\sum_{m=1}^{M}v\left(m,r_{p},r_{q}\right)\right|\leq\frac{v_{\rm
max}}{M}\,.$ (355)

Hence, we can make sure that the error in the Hamiltonian coefficients is less
than $\epsilon$ if $M>v_{\rm max}/\epsilon$. This error can be suppressed
exponentially since the cost of the comparison test is quadratic in the number
of qubits used to represent the number (because we need to compuare squares)
Sanders  _et al._ (2020), which is logarithmic in $M$, and thus logarithmic
$1/\epsilon$.

We can now block encode $V$ in the following way. We will start by applying
$\mathrm{PREP}$ as a series of Hadamard gates which realize the equal
superposition state over the ancilla registers $\lvert i,j\rangle\lvert
m\rangle$ (essentially all of the complexity of our approach will enter
through $\mathrm{SEL}$). Then, controlled on register $\lvert i\rangle$ we
will copy the momentum held in the $i^{\rm th}$ system register to an ancilla
register and translate it into $r_{p}$ using the relation in (351). We will
obtain $r_{q}$ in the same way, but this time controlled on the register
$\lvert j\rangle$. Then, we will compute $v(m,r_{p},r_{q})$ and record its
value in an ancilla register. Because this ancilla register is initialized to
$\lvert 0\rangle$ by $\mathrm{PREP}$, if it is changed to the $\lvert
1\rangle$ state by $\mathrm{SEL}$ then it will be removed from the linear
combination of unitaries. Indeed, we wish to remove the term from the linear
combination of unitaries whenever the inequality in (353) is false. Therefore,
the final step of our procedure should be to apply a bit flip operator to the
ancilla holding $v(m,r_{p},r_{q})$ and to uncompute the $r_{p}$ and $r_{q}$
ancillae. Thus, $\mathrm{SEL}$ will act with the following steps (with $p_{i}$
and $p_{j}$ equivalent to $p$ and $q$, respectively):

$\displaystyle\mathrm{SEL}\lvert i,j\rangle\lvert m\rangle\lvert
0\rangle^{\otimes{\cal O}(\log\left(N/\epsilon\right)}\underbrace{\lvert
p_{0}\cdots p_{i}\cdots p_{j}\cdots p_{\eta}\rangle}_{\lvert\psi\rangle}$
$\displaystyle\mapsto\lvert i,j\rangle\lvert m\rangle\lvert
r_{p_{i}},r_{p_{j}}\rangle\lvert 0\rangle\lvert p_{0}\cdots p_{i}\cdots
p_{j}\cdots p_{\eta}\rangle$ (356) $\displaystyle\mapsto\lvert
i,j\rangle\lvert m\rangle\lvert r_{p_{i}},r_{p_{j}}\rangle\lvert
v\left(m,r_{p_{i}},r_{p_{j}}\right)\rangle\lvert p_{0}\cdots p_{i}\cdots
p_{j}\cdots p_{\eta}\rangle$ $\displaystyle\mapsto\lvert i,j\rangle\lvert
m\rangle\lvert 0\rangle^{\otimes{\cal
O}(\log\left(N/\epsilon\right))}\lvert\neg
v\left(m,r_{p_{i}},r_{p_{j}}\right)\rangle\lvert p_{0}\cdots p_{i}\cdots
p_{j}\cdots p_{\eta}\rangle\,$

where $\neg$ implies the negation of the following bit.

One might be concerned that we would need to compute a reciprocal and a square
root as part of the 2-norm in the inequality for $v(m,r_{p},r_{q})$. However,
we can much more efficiently compute (353) by multiplying the norm through to
the other side and squaring it; in other words we can evaluate
$v(m,r_{p},r_{q})$ by testing

$4\,m^{2}v_{\rm max}^{2}\left\|r_{p}-r_{q}\right\|^{2}<M^{2}\,.$ (357)

The cost of evaluating this expression is also
$\mathcal{O}(\log^{2}M)=\mathcal{O}(\log^{2}(1/\epsilon))$ gates using the
result of Lemma 9, which again only contributes a polylogarithmic contribution
to the gate complexity. Apart from simplicity, this form of a comparison test
is much easier (will have lower constant factors) to implement because the
squared 2-norm is equivalent to just computing the sum of squares; whereas the
reciprocal square root can be a costly function to evaluate directly Soeken
_et al._ (2017); Häner  _et al._ (2018). We note that another approach to
avoiding the computation of reciprocal square roots in the context of quantum
simulating first quantized real space Hamiltonians is discussed in Poirier and
Jerke (2021), where it suggested that one take advantage of structure arising
from the inverse Laplace transformed version of the Coulomb kernel represented
in a Cartesian component-separated approach. While that approach may have some
advantages for product formula based methods, it is unclear whether it might
also be useful for further reducing the cost of block encoding the Coulomb
operator.

The block encoding of $U$ follows a similar pattern. In this case, we will re-
express the $U$ operator as

$\displaystyle U$ $\displaystyle=-\frac{u_{\rm
max}}{M}\sum_{i=1}^{\eta}\sum_{\ell=1}^{L}\zeta_{\ell}\sum_{p\in
G}\sum_{m=1}^{M}u\left(m,r_{p},R_{\ell}\right)\lvert p\rangle\\!\langle
p\rvert_{i}$ (358) $\displaystyle u\left(m,r_{p},R_{\ell}\right)$
$\displaystyle=\begin{cases}1&\qquad{\rm if}\qquad\frac{m\,u_{\rm
max}}{M}<\frac{1}{\left\|r_{p}-R_{\ell}\right\|}\\\ 0&\qquad{\rm
otherwise}\end{cases}$ (359) $\displaystyle u_{\rm max}$
$\displaystyle=\textrm{max}_{p,\ell}\left(\frac{1}{\left\|r_{p}-R_{\ell}\right\|}\right)={\cal
O}\left(\frac{N^{1/3}}{\Omega^{1/3}}\right)\,.$ (360)

where we will require that the $R_{\ell}$ do not overlap with any of the grid
points and we once again take advantage of the fact that the right-hand side
of the inequality in $u(m,r_{p},R_{\ell})$ is always greater than zero.
Similar to the strategy employed for $V$, here the $\mathrm{SEL}$ operation
will compute $u(m,r_{p},R_{\ell})$ into an ancilla that is then negated (thus,
excluding the cases where $u(m,r_{p},R_{\ell})=0$ from the linear combination
of unitaries). But there are two main differences: (1) the value $R_{\ell}$
will need to be loaded using a QROM and (2) we have not incorporated the
$\zeta_{\ell}$ into the definition of the inequality function that is computed
as part of $\mathrm{SEL}$.

Instead, we will include the $\zeta_{\ell}$ through the $\mathrm{PREP}$
operation. At the same time we will directly load the $\lvert R_{\ell}\rangle$
values into superposition (rather than computing it from $\ell$). Thus,
controlled on a flag ancilla specifying that the $U$ term is to be applied,
our $\mathrm{PREP}$ circuit will act as

$\mathrm{PREP}\lvert 0\rangle^{\otimes{\cal
O}\left(\log\left(N/\epsilon\right)\right)}\mapsto\sum_{\ell=1}^{L}\sum_{i=1}^{\eta}\sum_{m=1}^{M}\sqrt{\frac{\zeta_{\ell}}{\eta
M}}\lvert R_{\ell}\rangle\lvert i\rangle\lvert m\rangle\lvert 0\rangle\,.$
(361)

Using the coherent aliasing strategies developed in Babbush  _et al._ (2018a),
the weighted superposition over $R_{\ell}$ can be prepared with no more than
$\epsilon$ error in the amplitudes at Toffoli complexity $L+{\cal
O}(\log(1/\epsilon))$. That will be a negligible additive cost compared to the
rest of the $\mathrm{SEL}$ operation which will have Toffoli complexity ${\cal
O}(\eta)$ for charge neutral systems. Thus, up to exponentially small errors
in $\mathrm{PREP}$ we can see that

$\left|\sum_{\ell=1}^{L}\frac{\zeta_{\ell}}{\left\|r_{p}-R_{\ell}\right\|}-\sum_{\ell=1}^{L}\frac{u_{\rm
max}\zeta_{\ell}}{M}\sum_{m=1}^{M}u\left(m,r_{p},R_{\ell}\right)\right|\leq\frac{u_{\rm
max}\lambda_{\zeta}}{M}\,$ (362)

with $\lambda_{\zeta}=\sum_{\ell}\zeta_{\ell}$, and again we can make sure
that the error in the Hamiltonian coefficients is less than $\epsilon$ if
$M>u_{\rm max}\lambda_{\zeta}/\epsilon$. Once again, this $\epsilon$ will only
enter the gate complexity inside of logarithms.

Therefore, we see that we can implement the block encoding of
$(T+U+V)/\lambda$ with gate complexity $\widetilde{\cal O}(\eta)$. Given our
scheme, the $\lambda_{T}$ value is exactly the same as in (25), and

$\lambda_{V}=\eta^{2}v_{\rm max}={\cal
O}\left(\frac{\eta^{2}N^{1/3}}{\Omega^{1/3}}\right)\qquad\qquad\lambda_{U}=\eta\left(\sum_{\ell=1}^{L}\zeta_{\ell}\right)u_{\rm
max}={\cal O}\left(\frac{\eta^{2}N^{1/3}}{\Omega^{1/3}}\right)$ (363)

where these values of $\lambda_{V}$ and $\lambda_{U}$ involve slightly
different constant factors compared to the forms in (25), but the same
asymptotic scaling. Accordingly, we see that the overall asymptotic complexity
of the qubitization approach using this Hamiltonian would be

$\widetilde{\cal
O}\left(\frac{\left(\lambda_{T}+\lambda_{U}+\lambda_{V}\right)\eta}{\epsilon}\right)=\widetilde{\cal
O}\left(\frac{\eta^{2}N^{2/3}}{\epsilon\,\Omega^{2/3}}+\frac{\eta^{3}N^{1/3}}{\epsilon\,\Omega^{1/3}}\right)=\widetilde{\cal
O}\left(\frac{\eta^{4/3}N^{2/3}+\eta^{8/3}N^{1/3}}{\epsilon}\right)\,$ (364)

for $\eta\propto\Omega$, which is the same as in (26). It is likely that the
QFT will be the dominant cost in this method. The complexity for the plane
waves approach tends to be dominated by the cost of controlled selection of
the momentum registers, which has complexity $\mathcal{O}(\eta\log N)$ and is
the only part that is polynomial in $\eta$. The QFT has scaling as
$\widetilde{\mathcal{O}}(\eta\log^{2}N)$, which is larger and so should
dominate the complexity.

Finally, we discuss how the interaction picture scheme can be combined with
the real space first quantized representation in order to suppress the
complexity term proportional to $N^{2/3}$. To do this, we follow a similar
approach as in the main paper with the block encodings of $U$ and $V$ as
described in this appendix. The main difference is that before evolving under
the kinetic operator we must apply the QFT. The gate complexity of performing
the QFT in each time step is $\widetilde{\cal O}(\eta\log^{2}N)$ Toffoli
gates. A second difference is that it is no longer possible to keep the total
kinetic energy in an ancilla register and update it when we block-encode $U$
and $V$. At each step the kinetic energy will need to be recomputed, which has
complexity $\mathcal{O}(\eta\log^{2}N)$. That is almost as large as the
complexity scaling for the QFT. As a result there is the same asymptotic
scaling as the interaction picture algorithm in the main paper up to a
logarithmic factor,

$\widetilde{\cal
O}\left(\frac{\left(\lambda_{U}+\lambda_{V}\right)\eta}{\epsilon}\right)\widetilde{\cal
O}\left(\frac{\eta^{3}N^{1/3}}{\epsilon\,\Omega^{1/3}}\right)=\widetilde{\cal
O}\left(\frac{\eta^{8/3}N^{1/3}}{\epsilon}\right)\,.$ (365)

Again, the complexity will be increased over the plane wave approach primarily
due to the QFT.


$\text{\em positions}[i] = 1$

$i$ += 1

$i$ += $n - n\_0 - n\_1$

$i$ += 1

s $\coloneqq$ sequence

Algorithm to recover input sequence $s$.  We now give an overall algorithm (as shown in Algorithm <ref>) that recovers $s$ using the matching recovery algorithm and claims.
For the all 0 and all 1 input sequences, we can use $q^{(n+1)} = 0$, $q^{(n+2)} = 1$ to directly recover them (Algorithm  <ref>: line 2-5).
For the rest of the cases, we first recover the optimal isomorphic matching using the described algorithm (Algorithm  <ref>: line 6-15).
Let the recovered matching for $\mc{Q}$ be $m = x[1] \ldots x[n]$, ($x[i] \in \{0, 1\}$).
Denote the position of the leftmost 0 in $x$ to be $u$ ($1 \leq u \leq \ell$).
Then we know $s[1], \ldots, s[u-1] = 1$ by Proposition <ref>.
By using the sequence $q^{(n+2)} = 1$ to query $s$, we get the total number of 0's $(n_0)$ in $s$. Consider the substring $x[u, n]$, and delete the leading zeros in $x[u, n]$ until it has $n_0$ zeros. Suppose we obtain string $s_d$ after the deletion. We know that $s = 1^{u-1}s_d$  (Algorithm  <ref>: line 16-24).

If we are allowed to use $\mc{O}(n)$ extra characters in our queries, we have non-adaptive solutions with $\mc{O}(1)$ query complexity for DTW distance.
This assumption is stronger than the problem setting (where only $\mc{O}(1)$ extra characters are considered) throughout the paper.
For details of this complementary result, see Appendix <ref>.

§ RECOVERY WITH NON-ADAPTIVE DISTANCE ORACLE QUERIES

Consider two sequences $x$ and $y$ ($x \neq y$) defined on the binary alphabet $\bits$.
The query result from a distance oracle only gives very limited information, viz. 0 or 1 (which is more limited than the query from DTW oracle).
This 1-bit binary information restricts the power of sequence recovery with oracle.
Note that it is not possible to distinguish any sequences $x$ and $y$ under distance.
trivialTo see this and to see why the recovery problem is interesting for distance, we first define the concept of equivalent sequences under distance and revisit the problem from the perspective of equivalent sequences.

Given two sequences $x$ and $y$, we say $x$ and $y$ are equivalent if $y$ is obtained by taking any bit in $x$ and copying this bit contiguously any number of times.
For any pair of equivalent sequences, the distance between them is $0$.

A simple example of equivalent sequences under distance is two sequences, $1$ and $11$.
$11$ can be seen as copying the bit $1$ in the first sequence and the distance between $1$ and $11$ is $0$.
In addition, these two sequences cannot be distinguished by any query sequence.
This is because for the second sequence, the double $1$ characters can be matched to the same character in the query sequence as the single $1$ sequence.
This will not change the distance because the $l_\infty$ norm of the cost of matching edges is not changed.

From the perspective of equivalent sequences, for any two sequences $x$ and $y$, they are either in the same equivalence class (the distance is 0) or in different equivalence classes (the distance is 1).
Thus the distance between two sequences reflects whether or not they are equivalent.
Any equivalent sequences, therefore as suggested by its name, are not distinguishable, because all queries from the same equivalence class return 0 and all queries from different equivalence classes return 1.
Further, we can categorize all the equivalence classes under distance and then derive the lower bound of query complexity of recovering  non-equivalent sequences under distance, which is shown in the following theorem.

For a binary alphabet $\bits$, any algorithm to recover an arbitrary input sequence $s \in \bits^i$ revise_thm:frechet_lower_boundup to equivalence, where $0 \leq i \leq n$, by querying its distance to a non-adaptive set of sequences requires a query complexity of $\Omega(n)$.

classify_frechet We begin this proof of query complexity lower bound with a classification of all equivalence classes under the distance.
For each length $1 \leq i \leq n$, there exists two non-equivalent sequences under  distance, which are  $\underbrace{010101 \ldots}_{\text{of length} ~i}$ and $\underbrace{101010 \ldots}_{\text{of length} ~i}$, yielding $2n$ mutually non-equivalent sequences in total. As the distance oracle returns 0 when the input sequence and the query sequence are equivalent and 1 otherwise, we would need at least $2n-1$ queries to exactly recover the input sequence. If the number of queries is less than $2n-1$, we can always select 2 sequences from the $2n$ mutually non-equivalent sequences which are not covered by the queries, and these two sequences cannot be distinguished by the query sequences. This yields an $\Omega(n)$ lower bound on the query complexity.

In the analysis of non-adaptive strategies for DTW distance, we have shown that, with extra characters, we can obtain stronger results in recovering the exact sequence.
However, using queries from the extended alphabet (no matter how many extra characters are allowed) does not help increase the power of recovery under distance, proved in the following theorem.

Given two sequences $s$ and $s'$, if the distance $d_{F}(s, s') = 0$, then any query $q$ with extra characters cannot distinguish $s$ and $s'$.

not-triangleGiven sequences $s, s'$ (where $d_{F}(s, s') = 0$) and query $q$ with extra characters, our goal is to show $d_{F}(s, q) = d_{F}(s', q)$ for every possible $q$.
The technique of our proof is, for an optimal matching between $s$ and any query $q$, we can construct a matching between $s'$ and $q$ with the same cost, and vice versa.
In this way, we know that $d_{F}(s', q) \leq d_{F}(s, q)$ and $d_{F}(s, q) \leq d_{F}(s', q)$, so $d_{F}(s, q) = d_{F}(s', q)$ and $q$ cannot distinguish $s$ and $s'$.

Since $d_{F}(s, s') = 0$, $s$ and $s'$ have the same condensed expression. Suppose $s$ and $s'$ has $k$ runs.
In the optimal matching between $q$ and $s$, let $q_{s^{(i)}}$ denote the substring in $q$ which is matched to the $i$-th run of $s$ for every $i\in[k]$.
We can always match all $q_{s^{(i)}}$'s to the $i$-th run of $s'$ instead.
Note that the $i$-th runs of $s$ and $s'$ (denoted by $s^{(i)}$ and $s^{\prime(i)}$, resp.) are of the same character with maybe various length.
The distance between $q_{s^{(i)}}$ and $s^{(i)}$ only depends on the characters in $q_{s^{(i)}}$ and thus $d_F(q_{s^{(i)}}, s^{(i)}) = d_F(q_{s^{(i)}}, s^{\prime(i)})$.
Therefore we obtain a matching between $s'$ and $q$ with a cost of $d_{F}(s, q)$.
This matching between $q$ and $s'$ may be not optimal but is valid, and therefore we can conclude $d_{F}(s', q)\leq d_{F}(s, q)$.
Due to the symmetry of the statement, we can similarly obtain $d_{F}(s, q) \leq d_{F}(s', q)$.
This finishes the proof of this theorem.

Since extra characters are not helpful in recovering from distance queries, we conclude the analysis with a trivially interesting approach to recover sequences up to equivalence.
The approach uses up to $2n-1$ queries, which exactly matches our query complexity lower bound, as shown in the following theorem.

For a binary alphabet $\bits$ and two input sequences $s, s' \in \bits^i$ where $0 \leq i \leq n$ and $s$ and $s'$ are non-equivalent sequences under distance, there exists an algorithm to distinguish the input sequences $s$ and $s'$, given $2n-1 \in \mc{O}(n)$ query sequences $\revise{\mc{Q}}$ and the distance of $s$ and $s'$ to each query sequence $q \in \revise{\mc{Q}}$.

We first show that, for each length $0 \leq i \leq n$, there are only two non-equivalent sequences under distance, which are  $010101...$ and $101010...$ sequences, viz., we can identify two non-equivalent sequences by specifying the sequence length $i$ and the starting bit.
Therefore, for the maximum sequence length $n$, there are only $2n$ mutually non-equivalent sequences.

Given any two different sequences from this $2n$-sized collection of non-equivalent sequences under distance, we can use $\mc{O}(n)$ query sequences to distinguish them.
That is, we can utilize the exact set of $2n$ non-equivalent sequences as the query sequences.
If the query sequence $p$ is exactly the input sequence $q$, the distance between $p$ and $q$ is $d_F(p, q) = 0$.
If the query sequence $p$ is not equivalent to the input sequence $q$, then the distance between $p$ and $q$ is $d_F(p, q) = 1$ because it is impossible to skip over a bit without paying cost $1$.
Note that any one of the $2n$ queries can be skipped since we know the fact that there would be exactly one 0 among the $2n$ query results.
Therefore, $2n-1 \in \mc{O}(n)$ query sequences suffice to distinguish any two sequences from the non-equivalent sequence set and this finishes the proof.

This theorem shows that, if an input is in the collection of non-equivalent sequences under distance, we can use $\mc{O}(n)$ queries to exactly recover this sequence given the query results under the distance.

Remark: Extension to non-binary alphabets.
Our results are presented for input sequences from binary alphabet $\bits$.
These results can be extended to any non-binary alphabet $\Sigma$ by encoding the non-binary alphabet in a binary domain.
This will increase the query complexity by a constant factor from $|\Sigma|$ (one-hot encoding) to $\log(|\Sigma|)$ (binary encoding).
This extension works for the results for all distance metrics shown in this paper.
However, we note that this extension may not be optimal if one considers a large alphabet (e.g., larger than $n$).
In fact, calculating some of the distances themselves on a general alphabet is under SETH [Abboud et al., 2016, Bringmann and Künnemann, 2015], which is a much hard problem than on the binary case [Kuszmaul, 2021].
Obtaining optimal results on the extension of the non-decomposable distance recovery problem leaves room for future research.

§ RELATED WORK

A distance embedding [Cormode, 2003] embeds  sequences from the original distance metric space to other distance measures (usually $l_p$ norms), such that the distance measurements in the original space can be preserved up to a factor of $D$, namely the distortion rate.
The sequence distance embedding problem is related to our problem in the sense that, in our problem, we intend to recover the input sequence from a list of query results that are in the $l_p$ space, which can be regarded as finding a special distance embedding.
Existing works on the sequence distance embedding problem mainly focus on constructing such an embedding which can have a close approximation (viz., low distortion rate) and reduce the computational complexity (i.e., cost) on the new distance space.
[Andoni et al., 2003] shows a lower bound of $3/2$ on the distortion rate of embedding edit distance into $\ell_p$ norm spaces.
An improvement of $(\log n)^{\frac{1}{2} - o(1)}$ on this lower bound [Khot and Naor, 2005] has been further simplified and improved into $\Omega(\log n)$ by [Krauthgamer and Rabani, 2009].

Distance embeddings can be used to estimate the distance on the complex metric space because the evaluation and computations on the new (simpler metric) space can be significantly faster [Cormode, 2003].
Under the asymmetric query model (when estimating the edit distance between $x$ and $y$, the algorithm has unrestricted power accessing $x$ but limited power accessing $y$), [Andoni et al., 2010] proposes a $(\log n)^{\mc{O}(1/\epsilon)}$ approximation algorithm that runs in $n^{1 + \epsilon}$ time.
[Charikar et al., 2018] considers the alignment problem when estimating the edit distance (finding the sequence of edits between the estimated sequences) and presents an alignment with $(\log n)^{\mc{O}(1/\epsilon^2)}$ approximation in time $\Tilde{\mc{O}}(n^{1+\epsilon})$.
The sequence distance embedding problem has been investigated on other distance metrics as well, for example, the block edit distance [Cormode, 2003] and the Ulam distance [Charikar and Krauthgamer, 2006].
Existing work also shows embeddings from edit distance to the Hamming space [Belazzougui and Zhang, 2016, Chakraborty et al., 2016].
However, to the best of our knowledge, there is no prior work considering the embedding problem of the DTW distance and the exact recovery problem based on distance oracle query results.

§ OPEN PROBLEMS

We initiate an exact recovery problem of sequences using queries to a non-decomposable distance oracle.
We show recovery algorithms for edit distance, DTW distance, and distance, as well as a general adaptive algorithm for a wide class of distance oracles.
We envision the following directions for future work.

First, for the edit distance, there is still a quadratic gap between the non-adaptive query complexity upper and lower bounds without extra characters. Closing this gap requires a deeper understanding about the properties of edit distance.

Second, for the DTW distance, it remains unclear whether 1 extra character suffices for an $\mc{O}(n)$ non-adaptive upper bound, or we can have an $\Omega(n^2)$ non-adaptive lower bound with $1$ extra character (our proof uses $2$ extra characters).

Furthermore, as the initial work on non-decomposable distance recovery problem, we consider a simpler setting where input sequences are drawn from binary alphabet $\bits$.
While our results can be naturally extended to a non-binary alphabet, as stated in the paper, with a compensation of increasing the query complexity up to a constant factor, we notice that for some distances (e.g., DTW), the calculation on the general alphabet is much harder than on the binary case.
This spawns the open question for follow-up work to consider: Would there exist a strategy specifically designed for the non-binary alphabet with lower query complexity (than using encoding extensions to our results on the binary alphabet)?

Lastly, it would be interesting to consider the exact sequence recovery problem using the properties of specific distance metrics.
For example, the Edit distance with Real
Penalty (ERP) distance [Chen and Ng, 2004] which supports local time shifting in time series by the marriage of the $\ell_1$ norm and edit distance, would be of interest.
One can also consider other variants of our problem in terms of adaptive queries or the approximate recovery problem in the presence of noise.

§ ACKNOWLEDGEMENT

All authors thank the anonymous reviewers of ITCS 2023 and IEEE Transactions of Information Theory for their detailed comments which helped to improve the paper during the revision process.
David P. Woodruff would like to thank support from ONR grant N00014-18-1-2562 and a Simons Investigator Award.
Hongyang Zhang would like to thank support from NSERC Discovery Grant RGPIN-2022-03215, DGECR-2022-00357.

[Abboud et al., 2015]
Amir Abboud, Arturs Backurs, and Virginia Vassilevska Williams.
Tight hardness results for LCS and other sequence similarity
In Venkatesan Guruswami, editor, IEEE 56th Annual Symposium
on Foundations of Computer Science, FOCS 2015, Berkeley, CA, USA, 17-20
October, 2015, pages 59–78. IEEE Computer Society, 2015.
URL <https://doi.org/10.1109/FOCS.2015.14>.

[Abboud et al., 2016]
Amir Abboud, Thomas Dueholm Hansen, Virginia Vassilevska Williams, and Ryan
Simulating branching programs with edit distance and friends: or: a
polylog shaved is a lower bound made.
In Daniel Wichs and Yishay Mansour, editors, Proceedings of the
48th Annual ACM SIGACT Symposium on Theory of Computing, STOC 2016,
Cambridge, MA, USA, June 18-21, 2016, pages 375–388. ACM, 2016.
URL <https://doi.org/10.1145/2897518.2897653>.

[Afshani et al., 2019]
Peyman Afshani, Manindra Agrawal, Benjamin Doerr, Carola Doerr, Kasper Green
Larsen, and Kurt Mehlhorn.
The query complexity of a permutation-based variant of mastermind.
Discret. Appl. Math., 260:0 28–50, 2019.
URL <https://doi.org/10.1016/j.dam.2019.01.007>.

[Aldridge et al., 2019]
Matthew Aldridge, Oliver Johnson, and Jonathan Scarlett.
Group testing: An information theory perspective.
Found. Trends Commun. Inf. Theory, 150 (3-4):0
196–392, 2019.
URL <https://doi.org/10.1561/0100000099>.

[Amir et al., 2018]
Amihood Amir, Mika Amit, Gad M. Landau, and Dina Sokol.
Period recovery of strings over the hamming and edit distances.
Theor. Comput. Sci., 710:0 2–18, 2018.
URL <https://doi.org/10.1016/j.tcs.2017.10.026>.

[Andoni et al., 2003]
Alexandr Andoni, Michel Deza, Anupam Gupta, Piotr Indyk, and Sofya
Lower bounds for embedding edit distance into normed spaces.
In Proceedings of the Fourteenth Annual ACM-SIAM Symposium on
Discrete Algorithms, January 12-14, 2003, Baltimore, Maryland, USA, pages
523–526. ACM/SIAM, 2003.
URL <http://dl.acm.org/citation.cfm?id=644108.644196>.

[Andoni et al., 2010]
Alexandr Andoni, Robert Krauthgamer, and Krzysztof Onak.
Polylogarithmic approximation for edit distance and the asymmetric
query complexity.
In 51th Annual IEEE Symposium on Foundations of Computer
Science, FOCS 2010, October 23-26, 2010, Las Vegas, Nevada, USA, pages
377–386. IEEE Computer Society, 2010.
URL <https://doi.org/10.1109/FOCS.2010.43>.

[Aronov et al., 2006]
Boris Aronov, Sariel Har-Peled, Christian Knauer, Yusu Wang, and Carola Wenk.
Fréchet distance for curves, revisited.
In Yossi Azar and Thomas Erlebach, editors, Algorithms - ESA
2006, 14th Annual European Symposium, Zurich, Switzerland, September 11-13,
2006, Proceedings, volume 4168 of Lecture Notes in Computer Science,
pages 52–63. Springer, 2006.
URL <https://doi.org/10.1007/11841036_8>.

[Belazzougui and Zhang, 2016]
Djamal Belazzougui and Qin Zhang.
Edit distance: Sketching, streaming, and document exchange.
In Irit Dinur, editor, IEEE 57th Annual Symposium on
Foundations of Computer Science, FOCS 2016, 9-11 October 2016, Hyatt
Regency, New Brunswick, New Jersey, USA, pages 51–60. IEEE Computer
Society, 2016.
URL <https://doi.org/10.1109/FOCS.2016.15>.

[Braverman et al., 2019]
Vladimir Braverman, Moses Charikar, William Kuszmaul, David P. Woodruff, and
Lin F. Yang.
The one-way communication complexity of dynamic time warping
In Gill Barequet and Yusu Wang, editors, 35th International
Symposium on Computational Geometry, SoCG 2019, June 18-21, 2019, Portland,
Oregon, USA, volume 129 of LIPIcs, pages 16:1–16:15. Schloss
Dagstuhl - Leibniz-Zentrum für Informatik, 2019.
URL <https://doi.org/10.4230/LIPIcs.SoCG.2019.16>.

[Bressan et al., 2021]
Marco Bressan, Nicolò Cesa-Bianchi, Silvio Lattanzi, and Andrea
Exact recovery of clusters in finite metric spaces using oracle
In Mikhail Belkin and Samory Kpotufe, editors, Conference on
Learning Theory, COLT 2021, 15-19 August 2021, Boulder, Colorado, USA,
volume 134 of Proceedings of Machine Learning Research, pages
775–803. PMLR, 2021.
URL <http://proceedings.mlr.press/v134/bressan21a.html>.

[Bringmann and Künnemann, 2015]
Karl Bringmann and Marvin Künnemann.
Quadratic conditional lower bounds for string problems and dynamic
time warping.
In Venkatesan Guruswami, editor, IEEE 56th Annual Symposium
on Foundations of Computer Science, FOCS 2015, Berkeley, CA, USA, 17-20
October, 2015, pages 79–97. IEEE Computer Society, 2015.
URL <https://doi.org/10.1109/FOCS.2015.15>.

[Bshouty, 2009]
Nader H. Bshouty.
Optimal algorithms for the coin weighing problem with a spring scale.
In COLT 2009 - The 22nd Conference on Learning Theory,
Montreal, Quebec, Canada, June 18-21, 2009, 2009.

[Buchin et al., 2022]
Maike Buchin, Anne Driemel, Koen van Greevenbroek, Ioannis Psarros, and Dennis
Approximating length-restricted means under dynamic time warping.
In Parinya Chalermsook and Bundit Laekhanukit, editors,
Approximation and Online Algorithms - 20th International Workshop,
WAOA 2022, Potsdam, Germany, September 8-9, 2022, Proceedings, volume
13538 of Lecture Notes in Computer Science, pages 225–253. Springer,
URL <https://doi.org/10.1007/978-3-031-18367-6_12>.

[Cai et al., 2019]
Xingyu Cai, Tingyang Xu, Jinfeng Yi, Junzhou Huang, and Sanguthevar
DTWNet: a dynamic time warping network.
In Hanna M. Wallach, Hugo Larochelle, Alina Beygelzimer, Florence
d'Alché-Buc, Emily B. Fox, and Roman Garnett, editors, Advances
in Neural Information Processing Systems 32: Annual Conference on Neural
Information Processing Systems 2019, NeurIPS 2019, December 8-14, 2019,
Vancouver, BC, Canada, pages 11636–11646, 2019.

[Cantor and Mills, 1966]
David G. Cantor and W. H. Mills.
Determination of a subset from certain combinatorial properties.
Canadian Journal of Mathematics, 18:0 42–48, 1966.

[Chakraborty et al., 2016]
Diptarka Chakraborty, Elazar Goldenberg, and Michal Koucký.
Streaming algorithms for embedding and computing edit distance in the
low distance regime.
In Daniel Wichs and Yishay Mansour, editors, Proceedings of the
48th Annual ACM SIGACT Symposium on Theory of Computing, STOC 2016,
Cambridge, MA, USA, June 18-21, 2016, pages 712–725. ACM, 2016.
URL <https://doi.org/10.1145/2897518.2897577>.

[Charikar and Krauthgamer, 2006]
Moses Charikar and Robert Krauthgamer.
Embedding the ulam metric into l\({}_{\mbox{1}}\).
Theory Comput., 20 (11):0 207–224, 2006.
URL <https://doi.org/10.4086/toc.2006.v002a011>.

[Charikar et al., 2018]
Moses Charikar, Ofir Geri, Michael P. Kim, and William Kuszmaul.
On estimating edit distance: Alignment, dimension reduction, and
In Ioannis Chatzigiannakis, Christos Kaklamanis, Dániel Marx,
and Donald Sannella, editors, 45th International Colloquium on
Automata, Languages, and Programming, ICALP 2018, July 9-13, 2018, Prague,
Czech Republic, volume 107 of LIPIcs, pages 34:1–34:14. Schloss
Dagstuhl - Leibniz-Zentrum für Informatik, 2018.
URL <https://doi.org/10.4230/LIPIcs.ICALP.2018.34>.

[Chen and Ng, 2004]
Lei Chen and Raymond T. Ng.
On the marriage of $l_p$-norms and edit distance.
In Mario A. Nascimento, M. Tamer Özsu, Donald Kossmann,
Renée J. Miller, José A. Blakeley, and K. Bernhard Schiefer,
editors, (e)Proceedings of the Thirtieth International Conference on
Very Large Data Bases, VLDB 2004, Toronto, Canada, August 31 - September 3
2004, pages 792–803. Morgan Kaufmann, 2004.
URL <http://www.vldb.org/conf/2004/RS21P2.PDF>.

[Cohen et al., 2019]
Jeremy M. Cohen, Elan Rosenfeld, and J. Zico Kolter.
Certified adversarial robustness via randomized smoothing.
In Kamalika Chaudhuri and Ruslan Salakhutdinov, editors,
Proceedings of the 36th International Conference on Machine Learning,
ICML 2019, 9-15 June 2019, Long Beach, California, USA, volume 97 of
Proceedings of Machine Learning Research, pages 1310–1320. PMLR,
URL <http://proceedings.mlr.press/v97/cohen19c.html>.

[Coja-Oghlan et al., 2020]
Amin Coja-Oghlan, Oliver Gebhard, Max Hahn-Klimroth, and Philipp Loick.
Optimal group testing.
In Jacob D. Abernethy and Shivani Agarwal, editors, Conference
on Learning Theory, COLT 2020, 9-12 July 2020, Virtual Event [Graz,
Austria], volume 125 of Proceedings of Machine Learning Research,
pages 1374–1388. PMLR, 2020.
URL <http://proceedings.mlr.press/v125/coja-oghlan20a.html>.

[Cormode, 2003]
Graham Cormode.
Sequence distance embeddings.
PhD thesis, University of Warwick, Coventry, UK, 2003.
URL <http://wrap.warwick.ac.uk/61310/>.

[Dorfman, 1943]
Robert Dorfman.
The detection of defective members of large populations.
The Annals of Mathematical Statistics, 140
(4):0 436–440, 1943.

[Eiter and Mannila, 1994]
Thomas Eiter and Heikki Mannila.
Computing discrete fréchet distance.
Technical Report CD-TR 94/64, Christian Doppler Laboratory for Expert
Systems, TU Vienna, Austria, April 1994.

[Fernandez et al., 2019]
Manuel Fernandez, David P. Woodruff, and Taisuke Yasuda.
The query complexity of mastermind with l\({}_{\mbox{p}}\) distances.
In Dimitris Achlioptas and László A. Végh, editors,
Approximation, Randomization, and Combinatorial Optimization.
Algorithms and Techniques, APPROX/RANDOM 2019, September 20-22, 2019,
Massachusetts Institute of Technology, Cambridge, MA, USA, volume 145 of
LIPIcs, pages 1:1–1:11. Schloss Dagstuhl - Leibniz-Zentrum für
Informatik, 2019.
URL <https://doi.org/10.4230/LIPIcs.APPROX-RANDOM.2019.1>.

[Fredman and Willard, 1990]
Michael L. Fredman and Dan E. Willard.
BLASTING through the information theoretic barrier with FUSION
In Harriet Ortiz, editor, Proceedings of the 22nd Annual ACM
Symposium on Theory of Computing, May 13-17, 1990, Baltimore, Maryland,
USA, pages 1–7. ACM, 1990.
URL <https://doi.org/10.1145/100216.100217>.

[Hu et al., 2023]
Zhuangfei Hu, Xinda Li, David P. Woodruff, Hongyang Zhang, and Shufan Zhang.
Recovery from non-decomposable distance oracles.
IEEE Transactions on Information Theory, pages 1–1,

[Hu et al., 2023]
Zhuangfei Hu, Xinda Li, David P. Woodruff, Hongyang Zhang, and Shufan Zhang.
Recovery from Non-Decomposable Distance Oracles.
In Yael Tauman Kalai, editor, 14th Innovations in Theoretical
Computer Science Conference (ITCS 2023), volume 251 of Leibniz
International Proceedings in Informatics (LIPIcs), pages 73:1–73:22,
Dagstuhl, Germany, 2023b. Schloss Dagstuhl – Leibniz-Zentrum
für Informatik.
ISBN 978-3-95977-263-1.
URL <https://drops.dagstuhl.de/opus/volltexte/2023/17576>.

[Jiang and Polyanskii, 2019]
Zilin Jiang and Nikita Polyanskii.
On the metric dimension of cartesian powers of a graph.
J. Comb. Theory, Ser. A, 165:0 1–14, 2019.
URL <https://doi.org/10.1016/j.jcta.2019.01.002>.

[Khot and Naor, 2005]
Subhash Khot and Assaf Naor.
Nonembeddability theorems via fourier analysis.
In 46th Annual IEEE Symposium on Foundations of Computer
Science (FOCS 2005), 23-25 October 2005, Pittsburgh, PA, USA, Proceedings,
pages 101–112. IEEE Computer Society, 2005.
URL <https://doi.org/10.1109/SFCS.2005.54>.

[Knuth, 1976]
Donald E Knuth.
The computer as Master Mind.
Journal of Recreational Mathematics, 90 (1):0
1–6, 1976.

[Krauthgamer and Rabani, 2009]
Robert Krauthgamer and Yuval Rabani.
Improved lower bounds for embeddings
SIAM J. Comput., 380 (6):0 2487–2498, 2009.
URL <https://doi.org/10.1137/060660126>.

[Kremer et al., 1995]
Ilan Kremer, Noam Nisan, and Dana Ron.
On randomized one-round communication complexity.
In Frank Thomson Leighton and Allan Borodin, editors,
Proceedings of the Twenty-Seventh Annual ACM Symposium on Theory of
Computing, 29 May-1 June 1995, Las Vegas, Nevada, USA, pages 596–605.
ACM, 1995.
URL <https://doi.org/10.1145/225058.225277>.

[Kuszmaul, 2021]
William Kuszmaul.
Binary dynamic time warping in linear time.
CoRR, abs/2101.01108, 2021.
URL <https://arxiv.org/abs/2101.01108>.

[Lécuyer et al., 2019]
Mathias Lécuyer, Vaggelis Atlidakis, Roxana Geambasu, Daniel Hsu, and
Suman Jana.
Certified robustness to adversarial examples with differential
In 2019 IEEE Symposium on Security and Privacy, SP 2019,
San Francisco, CA, USA, May 19-23, 2019, pages 656–672. IEEE, 2019.
URL <https://doi.org/10.1109/SP.2019.00044>.

[Levenshtein, 1966]
Vladimir I Levenshtein.
Binary codes capable of correcting deletions, insertions, and
In Soviet physics doklady, volume 10, pages 707–710, 1966.

[Li and Vitányi, 1991]
Ming Li and Paul M. B. Vitányi.
Combinatorics and kolmogorov complexity.
In Proceedings of the Sixth Annual Structure in Complexity
Theory Conference, Chicago, Illinois, USA, June 30 - July 3, 1991, pages
154–163. IEEE Computer Society, 1991.
URL <https://doi.org/10.1109/SCT.1991.160256>.

[Rodríguez-Velázquez et al., 2014]
Juan Alberto Rodríguez-Velázquez, Ismael González Yero,
Dorota Kuziak, and Ortrud R. Oellermann.
On the strong metric dimension of Cartesian and direct products of
Discret. Math., 335:0 8–19, 2014.
URL <https://doi.org/10.1016/j.disc.2014.06.023>.

[Schaar et al., 2020]
Nathan Schaar, Vincent Froese, and Rolf Niedermeier.
Faster binary mean computation under dynamic time warping.
In Inge Li Gørtz and Oren Weimann, editors, 31st Annual
Symposium on Combinatorial Pattern Matching, CPM 2020, June 17-19, 2020,
Copenhagen, Denmark, volume 161 of LIPIcs, pages 28:1–28:13. Schloss
Dagstuhl - Leibniz-Zentrum für Informatik, 2020.
URL <https://doi.org/10.4230/LIPIcs.CPM.2020.28>.

[Selberg, 1949]
Atle Selberg.
An elementary proof of the prime-number theorem.
Annals of Mathematics, 500 (2):0 305–313,
ISSN 0003486X.
URL <http://www.jstor.org/stable/1969455>.

[Shapiro and Fine, 1960]
Harold S Shapiro and NJ Fine.
The American Mathematical Monthly, 670 (7):0
697–698, 1960.

[Sima and Bruck, 2021]
Jin Sima and Jehoshua Bruck.
Trace reconstruction with bounded edit distance.
In IEEE International Symposium on Information Theory, ISIT
2021, Melbourne, Australia, July 12-20, 2021, pages 2519–2524. IEEE,
URL <https://doi.org/10.1109/ISIT45174.2021.9518244>.

[Soderberg and Shapiro, 1963]
Staffan Soderberg and H. S. Shapiro.
A combinatory detection problem.
The American Mathematical Monthly, 700 (10):0
1066–1070, 1963.
ISSN 00029890, 19300972.
URL <http://www.jstor.org/stable/2312835>.

[Sunjaya and Sunjaya, 2020]
Angela Felicia Sunjaya and Anthony Paulo Sunjaya.
Pooled testing for expanding covid-19 mass surveillance.
Disaster Medicine and Public Health Preparedness, 140
(3):0 e42–e43, 2020.

[Vershynin, 2011]
Roman Vershynin.
Lectures in geometric functional analysis.
Unpublished manuscript. Available at http://www-personal.
umich. edu/romanv/papers/GFA-book/GFA-book. pdf, 30 (3):0
3–3, 2011.

[Wang et al., 2018]
Chao Wang, Qing Zhao, and Chen-Nee Chuah.
Optimal nested test plan for combinatorial quantitative group
IEEE Trans. Signal Process., 660 (4):0
992–1006, 2018.
URL <https://doi.org/10.1109/TSP.2017.2780053>.

[Yelin et al., 2020]
Idan Yelin, Noga Aharony, Einat Shaer Tamar, Amir Argoetti, Esther Messer, Dina
Berenbaum, Einat Shafran, Areen Kuzli, Nagham Gandali, Omer Shkedi, Tamar
Hashimshony, Yael Mandel-Gutfreund, Michael Halberthal, Yuval Geffen, Moran
Szwarcwort-Cohen, and Roy Kishony.
Evaluation of COVID-19 RT-qPCR Test in Multi sample Pools.
Clinical Infectious Diseases, 710 (16):0
2073–2078, 05 2020.
ISSN 1058-4838.
URL <https://doi.org/10.1093/cid/ciaa531>.

§ OTHER RELATED WORK

Recovery problems in metric spaces.
Our problem is related to the recovery or reconstruction problems over metric spaces.
[Amir et al., 2018] study the period recovery problem on strings, which is to find the primitive periods between two strings such that the periodic distance is below a threshold.
They present an $\mc{O}(n \log n)$-time algorithm for Hamming distance and an $\mc{O}(n^{4/3})$-time algorithm for edit distance.
[Sima and Bruck, 2021] investigate the approximate recovery problem over bounded edit distance spaces in the presence of noise and show $n^{\mc{O}(k)}$ noisy samples suffice for (approximate) reconstruction.
Interestingly, [Bressan et al., 2021] consider the exact recovery problem using oracle queries but the objective of their work is to exactly recover the clusters in Euclidean space, which is similar but orthogonal to our problem.

Learning problems: Coin-weighing and group-testing problems.
The related “decomposable” instance to our problem of querying a Hamming distance oracle is equivalent to the coin-weighing problem [Bshouty, 2009] and the quantitative group testing problem [Wang et al., 2018].
Both the coin-weighing problem and group-testing problems are well-studied learning problems in the literature and have many real-world applications [Soderberg and Shapiro, 1963, Cantor and Mills, 1966, Li and Vitányi, 1991, Sunjaya and Sunjaya, 2020, Yelin et al., 2020].
The coin-weighing problem is to determine the weight of each coin (of two distinct weights $w_1$ and $w_2$) by using a minimal number of weighings of a subset of $n$ total coins each time.
[Cantor and Mills, 1966] and [Bshouty, 2009] respectively present $2n/\log n$ weighing solutions which are optimal non-adaptive solutions to this problem.
Assuming the number of $w_1$ weight coins is known to be $d$, this $d$-coin weighing problem can be solved by an adaptive algorithm in time $2d\log \frac{n}{d}/\log d + \mc{O}(d/\log d + d(\log \log d) \log \frac{n}{d}/(\log d)^2)$ [Bshouty, 2009].
The major difference between our problem and these well-studied problems is that we consider distance metrics which cannot be aligned and represented as $\sum_i^n f(x_i - y_i)$ (i.e., the edit distance, DTW distance, and distances).

§ COORDINATE DESCENT ALGORITHM INSTANTIATION

Now we briefly discuss how we apply the Coordinate Descent algorithm to all three distances we consider in this paper by justifying the two conditions hold.

Edit distance.
For condition 2, we know that $\forall s, q, \dist(s, q) \leq n$ since the maximum length of $s$ or $q$ is $n$.
For condition 1, in each iteration, we consider a set $Q$ that contains all sequences that can be transformed from $q$ by inserting, deleting or substituting one character in $q$ (edit operations).
Note that $|Q|$ cannot exceed $(n+1) + n +n= 3n+1$.
We claim that there exists a $q'$ in $Q$ such that $\revise{\dist(s, q) > \revise{\dist}(s, q')}$.
Let $\revise{\dist(s, q) = d}$.
By the definition of edit distance, there exists a chain of edit operations of length $d$ that transforms $s$ to $q$, resulting in a list of intermediate sequences $q_1, ..., q_{d-1}$.
Note that $\revise{\dist(s, q) \geq \dist(q_1, q) + 1}$, otherwise we have $\revise{\dist(q_1, q) > d - 1}$.
However, the chain implies we can transform $q_1$ to $q$ in $d-1$ edit operations, which leads to a contradiction.
Since $q_1 \in Q$, we can find $q_1$ satisfying the condition in $3n+1$ searches.
Therefore, the algorithm is guaranteed to recover the input in $\mc{O}(n^2)$ steps.

DTW distance.
For DTW distance, condition 2 holds since $\forall s, q, \dist(s, q) \leq n$.
For condition 1, consider the $\#\textsc{runs}(x)$ in $s$ and $q$. If
$\#\textsc{runs}(x)$ of $q$ $<$ $s$, then either adding an (arbitrary length) run to the start or the end of $q$ will decrease the DTW distance from $s$. On the other hand, if
$\#\textsc{runs}(x)$ of $q$ $>$ $s$, then either deleting a run from the start or the end of $q$ will decrease the DTW distance from $s$.
If $\#\textsc{runs}(x)$ of $q$ $=$ $s$ and $\revise{\dist}(s, q) \neq 0$, we can still decrease the distance from $q$ by either adding/deleting a run to the start/end of the sequence.
Therefore, the algorithm is guaranteed to recover the input in $\mc{O}(n^2)$ steps.

Condition 2 holds since $\forall s, q, \dist(s, q) \leq 1$.
For condition 1, enumerating $2n$ non-equivalent sequences, (i.e., 010101... and 101010...) guarantees to find $q'$ such that $\revise{\dist(s, q) > \revise{\dist}(s, q') = 0}$.
Therefore, the algorithm terminates in $\mc{O}(n)$ steps.

§ PROOFS OF CLAIM <REF> IN LEMMA <REF>

Proof of Claim <ref>. We prove this by contradiction. Suppose $\exists(i,j)$ where $i\in[n]$ and $j\in[\ell]$ such that  $\revise{\deg}(q[i])>1$ and $\revise{\deg}(s[j])>1$.
Let $\mb{X} = \{x \in [n] \mid \revise{\deg}(q[x]) > 1\}$, $\mb{Y} = \{y \in [\ell] \mid \revise{\deg}(s[y]) > 1 \}$ and $\mb{Z} = \{z \in [\ell] \mid \exists x \in \mb{X} \text{ such that edge }(q[x], s[z])\in M\}$.
According to Claim <ref>, we know that $\mb{Y}\bigcap \mb{Z} = \emptyset$. Let $d = \min_{y\in \mb{Y}, z\in \mb{Z}}|y-z|$. We would have $d>0$. Suppose we have $x_0\in \mb{X}, y_0\in \mb{Y}, z_0\in \mb{Z}$ such that edge $(q[x_0], s[z_0])\in M$ and $|y_0-z_0| = d$.
There are two cases to discuss. As we have already solved the case $y_0<z_0$ in the proofs of Claim <ref> in Lemma <ref>, here we only discuss the case $y_0 > z_0$.

In this case, we can assume that $s[y_0]$ is matched to

$$\{q[w], q[w+1], ..., q[w+\revise{\deg}(s[y_0])-1]\}$$
and $q[x_0]$ is matched to
$$\{s[z_0-\revise{\deg}(q[x_0])+1], s[z_0-\revise{\deg}(q[x_0])+2],..., s[z_0]\}.$$

We remove $d+1$ edges $E = \{(q[x_0], s[z_0]), (q[x_0+1], s[z_0+1]), ..., (q[w], s[y_0])\}$and add $d$ new edges

$$E'=\{(q[x_0+1], s[z_0]), (q[x_0+2], s[z_0+1]), ..., (q[w], s[y_0-1])\}$$
to construct a new matching $M'$. Since $\revise{\deg}(s[y_0])>1$ and $\revise{\deg}(q[x_0])>1$ in $M$, $M'$ would still be a valid matching. Computing the sum of two sets of edges $E$ and $E'$, respectively, would yield
Equation. <ref> (see the cross-column equations).

In this case, we can assume that $s[y_0]$ is matched to
$\{q[w], q[w+1], ..., q[w+\revise{\deg}(s[y_0])-1]\}$ and $q[x_0]$ is matched to $\{s[z_0-\revise{\deg}(q[x_0])+1], s[z_0-\revise{\deg}(q[x_0])+2],..., s[z_0]\}.$
We remove $d+1$ edges $E = \{(q[x_0], s[z_0]), (q[x_0+1], s[z_0+1]), ..., (q[w], s[y_0])\}$and add $d$ new edges

$$E'=\{(q[x_0+1], s[z_0]), (q[x_0+2], s[z_0+1]), ..., (q[w], s[y_0-1])\}$$
to construct a new matching $M'$. Since $\revise{\deg}(s[y_0])>1$ and $\revise{\deg}(q[x_0])>1$ in $M$, $M'$ would still be a valid matching. Computing the sum of two sets of edges $E$ and $E'$, respectively, would yield the following.

\begin{align}
\Cost(E) &  = |s[y_0]-q[w]|+\sum_{i = 1}^{d}|s[z_0+i-1]-q[x_0+i-1]| \label{eqn_dbl_2} \\
& > |q[w]-q[x_0]|+ \sum_{i = 1}^{d}|s[z_0+i-1]-q[x_0+i-1]|  \hfill \tag{Equation.~\ref{equ:single_direction_o1}} \\
& = \sum_{i = 1}^{d}|q[x_0+i-1]-q[x_0+i]|+ \sum_{i = 1}^{d}|s[z_0+i-1]-q[x_0+i-1]|  \hfill \tag{Monotonicity of $q$}  \\
& =  \sum_{i = 1}^{d}(|q[x_0+i-1]-q[x_0+i]|+|s[z_0+i-1]-q[x_0+i-1]|) \nonumber \\
& \geq \sum_{i = 1}^{d}|q[x_0+i]-s[z_0+i-1]| \hfill \tag{Triangle Inequality}\\
& = \revise{\Cost}(E'). \nonumber
\end{align}

Hence, $M'$ would be a better matching than $M$, a contradiction. Combining 1) and 2) completes the proof of Claim <ref>.

\begin{align}
\Cost(E) &  = |s[y_0]-q[w]|+\sum_{i = 1}^{d}|s[z_0+i-1]-q[x_0+i-1]| \label{eqn_dbl_2} \\
& > |q[w]-q[x_0]|+ \sum_{i = 1}^{d}|s[z_0+i-1]-q[x_0+i-1]|  \hfill \tag{Equation.~\ref{equ:single_direction_o1}} \\
& = \sum_{i = 1}^{d}|q[x_0+i-1]-q[x_0+i]|+ \sum_{i = 1}^{d}|s[z_0+i-1]-q[x_0+i-1]|  \hfill \tag{Monotonicity of $q$}  \\
& =  \sum_{i = 1}^{d}(|q[x_0+i-1]-q[x_0+i]|+|s[z_0+i-1]-q[x_0+i-1]|) \nonumber \\
& \geq \sum_{i = 1}^{d}|q[x_0+i]-s[z_0+i-1]| \hfill \tag{Triangle Inequality}\\
& = \revise{\Cost}(E'). \nonumber
\end{align}

§ RECOVERY USING NON-ADAPTIVE DTW DISTANCE ORACLE WITH ON EXTRA CHARACTERS

Define a sequence of $n$ elements, each of which has $\mc{O}(\log n)$ bit complexity, as a query sequence.
For a binary alphabet $\bits$ and an input sequence $s := \bits^{\ell}$ where $0 \leq \ell \leq n$, there exists an algorithm to recover the input sequence $s$, given $4 \in \mc{O}(1)$ query sequences $\revise{\mc{Q}}$ and the $d_{\revise{\DTW}}(s, q)$ to each query sequence $q \in \revise{\mc{Q}}$.

We note that, if we remove the constraint of $\mc{O}(\log n)$ bit complexity, we can give a straightforward solution by leveraging the query string $q = \{ \frac{1}{(n+1)}, \frac{1}{(n+1)^2}, \dots, \frac{1}{(n+1)^n} \}$ to encode much more information in a single query.
With the word RAM bit complexity requirement [Fredman and Willard, 1990] on the queries though, namely that each entry fits into a single $O(\log n)$-bit word, such solutions are not allowed.

Proof of Theorem <ref>.
Note that we can still use query sequences $0$ and $1$ to recover input sequences consisting of only $0$s or $1$s. For simplicity, we assume in the rest of the proof that the input sequence $s$ contains both $0$ and $1$ and let $s = s[1]s[2]...s[\ell]$.

We give our proof by constructing $2$ query sequences $q$ and $q'$ and presenting an algorithm to recover an input sequence $s$ from its DTW distance to these $2$ query sequences.

Query Sequences Construction.
Let $P_{prime} = \{p_1, p_2, ..., p_n\}$ be the first $n$ primes not including $2$. By the prime number theorem [Selberg, 1949], we have that $p_n = \mc{O}(n\log n)$.  Note that for any prime number $p_i > 2$, $\exists 1\leq x_i < p_i$ such that $\frac{1}{4} < \frac{x_i}{p_i} < \frac{1}{2}$. We obtain $q$ by selecting such a $\frac{x_i}{p_i}$ for each $p_i\in P_{prime}$ and rearranging them in increasing order. Then we construct $q$ as $q = q[1]q[2]...q[n]$ where $\frac{1}{4} < q[1] < q[2] < ... < q[n] < \frac{1}{2}$. Let $q'[i] = 1-q[i], 1\leq i \leq n$, and let $q' = q'[1]q'[2]...q'[n]$. We would have $\frac{3}{4} > q'[1] > q'[2] > ... > q'[n] > \frac{1}{2}$. Since $p_n = \mc{O}(n\log n)$, it is easy to verify that each $q[i]$ and $q'[i]$ does have bit complexity  $\mc{O}(\log n)$.

According to Lemma <ref>, each element in $q$ would be involved exactly once in $d_{\revise{\DTW}}(q, s)$, and a similar argument would hold for $q'$. We hereby present an algorithm to determine the value of the matched element for each element in $q$, and the same algorithm can also be applied to $q'$.

Algorithm to determine matched elements for a query sequence. Suppose $q[i] = \frac{x_{t_i}}{p_{t_i}}$ where $\{t_j\}$ is a permutation of $[n]$. Letting $m_i$ be the value matched to $q[i]$ in the optimal DTW matching for $(q, s)$ (different $m_i$s could correspond to the same element in $s$), $m_i \in\{0, 1\}$ , we would have

\begin{align*}
d_{\revise{\DTW}}(q, s) = \sum_{i=1}^{n}|m_i-q[i]| =
\sum_{i=1}^{n}\left|\frac{m_{i}p_{t_i}-x_{t_i}}{p_{t_i}}\right| \\
=\frac{\sum_{i=1}^{n} \left( |m_{i}p_{t_i}-x_{t_i}| \cdot \Pi_{j\neq t_i}p_j \right) }{{\Pi_{i=1}^{n}p_{i}}}.
\end{align*}

Let $d_{\revise{\DTW}}(q, s) = \frac{u}{v}$, where $u$ and $v$ are co-primes.
We have $u = \sum_{i=1}^{n}\left(|m_{i}p_{t_i}-x_{t_i}|\cdot\Pi_{j\neq t_i}p_j \right)$ and $v = \Pi_{i=1}^{n}p_{i}$.
Consider $u \bmod p_{t_k}$ for a specific $k$.
As each term in the summation has a factor $p_{t_k}$ except $|m_{k}p_{t_k}-x_{t_k}|\cdot\Pi_{j\neq t_k}p_j$, we have
$a \equiv |m_{k}p_{t_k}-x_{t_k}|\cdot\Pi_{j\neq t_k}p_j \mod{p_{t_k}}.$
Note that $p_{t_k}-x_{t_k}\not\equiv x_{t_k}\mod{p_{t_k}}$, so $(p_{t_k}-x_{t_k})\cdot \Pi_{j\neq t_k}p_j\not\equiv x_{t_k}\cdot \Pi_{j\neq t_k}p_j\bmod p_{t_k}$. Thus, we can determine $m_k$ by checking whether $(p_{t_k}-x_{t_k})\cdot \Pi_{j\neq t_k}p_j\equiv u\mod p_{t_k}$ or $x_{t_k}\cdot \Pi_{j\neq t_k}p_j\equiv u\mod p_{t_k}$.

Furthermore, we have the following claim for the optimal DTW matching between $q$, $q'$ and $s$.

For any given input sequence $s$ and optimal DTW matching $M$ and $M'$ for $(q,s)$ and $(q', s)$ respectively, we have $\revise{\deg}(s_[i]) = 1$ in $M$ if $s[i]=1$ and $\revise{\deg}(s[i]) = 1$ in $M'$ if $s[i]=0$ .

Proof of claim. We give a proof by contradiction. Given an optimal DTW matching $M$ for $(q,s)$, suppose $\exists 1\leq i\leq \ell$ such that $s[i] = 1$ and $\revise{\deg}(s[i]) > 1$. Suppose $s[i]$ is matched to $q[j], q[j+1], \ldots, q[j+\revise{\deg}(s[i])-1]$.

First, we show that we can “swap” $s[i]$ with its neighboring element while maintaining optimality of the matching. If one of the neighboring elements of $s[i]$ is $1$, without loss of generality, suppose $s[i+1] = 1$, then we can construct an alternate optimal matching $M^*$ where $\revise{\deg}(s[i]) = 1$ and $\revise{\deg}(s[i+1])>1$. According to Lemma <ref>, $s[i+1]$ cannot be matched with any of $q[j], q[j+1], \ldots, q[j+\revise{\deg}(s[i])-1]$ in $M$, otherwise there would exist $j+1\leq k\leq j+\revise{\deg}(s[i])-1$ such that $\revise{\deg}(q[k]) = 2$. Thus by matching $q[j+1], \ldots, q[j+\revise{\deg}(s[i])-1]$ to $s[i+1]$ instead of $s[i]$, we would obtain a new optimal matching $M^*$ where $\revise{\deg}(s[i]) = 1$ and $\revise{\deg}(s[i+1])>1$.

As there exists at least one $0$ in $s$, we know that there exists an optimal DTW matching $M_0^{*}$ for $(q, s)$ where $\exists s[i]$ such that $s[i] = 1$, $\revise{\deg}(s[i]) > 1$ and one of the neighboring elements of $s[i]$ is $0$. Without loss of generality, suppose $s[i+1] = 0$. Similarly, according to Lemma <ref>, $s[i+1]$ cannot be matched with any of $q[j], q[j+1], \ldots, q[j+\revise{\deg}(s[i])-1]$ in $M_0$. Here we construct a new matching $M_0'$ by matching $q[j+1], \ldots , q[j+\revise{\deg}(s[i])-1]$ to $s[i+1]$ instead of $s[i]$. Fig <ref> illustrates an example of such a construction. Considering the total cost of differing edges in both matchings, we have
> \sum_{k=j+1}^{j+\revise{\deg}(s[i])-1}\frac{1}{2}
> \sum_{k=j+1}^{j+\revise{\deg}(s[i])-1}|s[i+1]-q[k]|.$
Thus $M_0^*$ would be a better matching than $M_0$, causing a contradiction and thus finishing the proof.
A similar proof can be derived for query sequence $q'$ and the case $s[i]=0$.

Algorithm to recover $s$.  We now give an overall algorithm that recovers $s$ using the above algorithm and claim.
Applying the above algorithm gives the matched elements of $q$ and $q'$.
Let the matching result for $q$ and $q'$ be $m = m_1...m_n$ and $m' = m_1'...m_n'$, ($m_i, m_i' \in \{0, 1\}$) respectively.
We break $m$ and $m'$ into blocks such that each block is the longest substring that contains either $0$ or $1$. By also breaking $s$ into such blocks, we know that $m$ has the same number of blocks as $s$ according to Lemma <ref>. Similarly $m'$ has the same number of blocks as $s$. Let $l$ be the number of blocks that $m$ and $m'$ have. Then we can represent $m$ and $m'$ as $m=A_{1}...A_{l}$ and $m'=B_{1}...B_{l}$. Note that if $A_i$ contains only $1$, then $s$ must have the same number of $1$'s in the i-th-2$i$-th block, otherwise there will be some $s[k]$ for which  $\revise{\deg}(s[k])>1$, which contradicts Claim <ref>. Similarly, if $B_j$ contains only $0$, then $s$ has the same number of $0$'s in the i-th-3$j$-th block.
Then we can fully recover $s$ as $h(A_1)...h(A_l)$ where $h(X_i) = A_i$ if $X_i$ contains $1$ and $h(X_i) = B_i$ if $X_i$ contains $0$.