Efficiently constructing a quantum uniform superposition over bit strings near a binary linear code
Edward Farhi Stephen P. Jordan


unsrt

Performant near-term quantum combinatorial optimization
    Phillip C. Lotshaw
    April 29, 2024
=======================================================




    We demonstrate that a high fidelity approximation to |Ψ_b⟩, the quantum superposition over all bit strings within Hamming distance b of the codewords of a dimension-k linear code over ℤ_2^n, can be efficiently constructed by a quantum circuit for large values of n, b and k which we characterize. We do numerical experiments at n=1000 which back up our claims. The achievable radius b  is much larger than the distance out to which  known classical algorithms can efficiently find the nearest codeword. Hence, these states cannot be prepared by quantum constuctions that require uncomputing to find the codeword nearest a string. Unlike the analogous states for lattices in ℝ^n, |Ψ_b⟩ is not a useful resource for bounded distance decoding because the relevant overlap falls off too quickly with distance and known classical algorithms do better. Furthermore the overlap calculation can be dequantized. Perhaps these states could be used to solve other code problems. The technique used to construct these states is of interest and hopefully will have applications beyond codes.




§ INTRODUCTION



Given a probability distribution p over the bit strings of length n one can define a corresponding quantum state |p⟩ = ∑_𝐱∈{0,1}^n√(p(𝐱))|𝐱⟩. The problem of preparing these states when p can be efficiently sampled from classically was originally called qsampling in <cit.> but we call it Quantum Distributional State Construction, to emphasize that it is a problem of preparing coherent quantum states. Given |p⟩ one can sample from the classical probability distribution p by measuring |p⟩ in the computational basis. QDSC can also be used to do things that classical sampling cannot. In particular, given two probability distributions p and q on {0,1}^n, there is no way in general to determine whether p is close to q (e.g. in total variation distance) by drawing polynomially many samples from p and q. In contrast, using polynomially many copies of |p⟩ and |q⟩ one can answer this question using the Hadamard test.

Here, we consider the problem of QDSC for the uniform distribution over all bit strings within radius b of a binary linear code. Any binary linear code can be expressed in terms of a generator matrix B ∈ℤ_2^k × n as follows:

    C = {𝐱 B : 𝐱∈ℤ_2^k }.

The quantities n and k are known as the length and dimension of the code, respectively. There are 2^k codewords. (We are assuming that the rows of B are linearly independent.) For a given linear code C ⊂ℤ_2^n and a radius b ∈{0,1,…,n}  we will show how to efficiently produce the state

    |Ψ_b⟩ = 𝒩∑_𝐜∈ C∑_𝐳∈ℤ_2^nΘ(b-|𝐳|) |𝐳⊕𝐜⟩,

where | 𝐳 | denotes the Hamming weight of 𝐳, and Θ is the unit step function:

    Θ(y)  = {[         1 if  y ≥ 0;         0 if  y < 0 ]. .


For now, assume that b is less than half the distance of the code. (The distance of the code is the minimum Hamming weight of the difference between any two codewords.) In this case, the state |Ψ_b⟩ is a superposition of non-overlapping balls of radius b centered at the codewords of C. The normalization factor is then

    𝒩 = 1/√(2^k Vol(b))

where Vol(b) is the number of points contained within a ball of Hamming radius b:

    Vol(b) = ∑_j=0^b nj.

If the balls overlap the normalization is more complicated.

The uniform probability distribution over all bit strings within Hamming distance b of a binary linear code is easy to sample from classically. First, sample 𝐱 uniformly from bit strings of length k and sample 𝐳 uniformly from bit strings of length n and Hamming weight at most b. Then compute 𝐱 B ⊕𝐳. This suggests an approach to prepare the state |Ψ_b⟩. First, prepare

    ( 1/√(2^k)∑_𝐱∈ℤ_2^k |𝐱⟩) ( 1/√(Vol(b))∑_𝐳∈ℤ_2^n Θ(b-|𝐳|) |𝐳⟩) |0⟩^⊗ n.

Next, reversibly compute 𝐱 B ⊕𝐳 into the last register, yielding

    1/√(2^k)∑_𝐱∈ℤ_2^k |𝐱⟩ 1/√(Vol(b))∑_𝐳∈ℤ_2^n Θ(b-|𝐳|) |𝐳⟩ |𝐱 B ⊕𝐳⟩.

Last, we need to uncompute 𝐱 and 𝐳 from the first two registers to obtain |Ψ_b⟩. By assumption, b is less than half the distance of the code. Therefore, the balls do not overlap and 𝐱 and 𝐳 are uniquely determined by 𝐱 B ⊕𝐳. However, determining 𝐱 and 𝐳 from 𝐱 B ⊕𝐳 is in general a computationally hard problem. Specifically, it is a bounded distance decoding problem for the code C where the bound on distance is b.


    We are given a linear code C ⊂ℤ_2^n specified by its generator matrix, as in eq:generator. We are given a bit string 𝐯∈ℤ_2^n and are promised that 𝐯 has Hamming distance at most d from the nearest codeword: min_𝐜∈ C | 𝐯 - 𝐜 | ≤ d. The bounded distance decoding problem is to find 𝐜.


Bounded distance decoding for general binary linear codes can only be solved out to d logarithmic in n using known polynomial-time algorithms <cit.>. Consequently the method just described is only capable of efficiently preparing the state |Ψ_b⟩ for b = O(log n). Preparing |Ψ_b⟩ out to larger radius by this direct approach fails due to an uncomputation barrier. In this work, we show that by taking a more intrinsically quantum approach we can surpass the uncomputation barrier and efficiently prepare |Ψ_b⟩ out to b much larger than the radius to which we can solve bounded distance decoding.

In <cit.> it was shown that the ability to prepare superpositions over points within Euclidean distance b of the elements of a lattice in ℝ^n would yield efficient quantum algorithms to solve classically-intractible instances of the Bounded Distance Decoding (BDD) problem over these lattices. No efficient quantum circuits for this state preparation task are known. The states |Ψ_b⟩ that we prepare here are closely analogous to the states described in <cit.>, except over ℤ_2 instead of ℝ and using Hamming distance instead of Euclidean distance.

The state |Ψ_b⟩ can be used as a resource for BDD on binary linear codes using the same approach that was proposed in <cit.> for lattices over ℝ. For 𝐯∈ℤ_2^n let T_𝐯 be the corresponding unitary translation operator defined by T_𝐯|𝐱⟩ = |𝐱⊕𝐯⟩. If the overlap between |Ψ_b⟩ and its translation by T_𝐯 is bigger than 0 then 𝐯 is within distance 2b of the lattice. Given |Ψ_b⟩, one can use the Hadamard test to estimate the overlap ⟨Ψ_b| T_𝐯|Ψ_b⟩, as illustrated in figure <ref>. This overlap is determined by the distance from 𝐯 to the nearest codeword and decreases monotonically with this distance. (See sec:overlap.) By iteratively flipping a bit of 𝐯 and accepting the move if the distance is decreased one can find the nearest codeword.



The number of Hadamard tests required to detect that ⟨Ψ_b| T_𝐯|Ψ_b⟩ is nonzero and hence that 𝐯 is less than distance 2b from the nearest codeword is on the order of 1/|⟨Ψ_b| T_𝐯|Ψ_b⟩|^2. Using amplitude amplification, this cost can be reduced quadratically <cit.>. However, the overlap |⟨Ψ_b| T_𝐯|Ψ_b⟩| falls rapidly with |𝐯|. Consequently, as we show in sec:bddimp, using the states |Ψ_b⟩ in this way beats brute force search but fails to outperform the classical algorithm Information Set Decoding. Furthermore, in sec:dequant, we show that this approach can be dequantized.

It is interesting to see another example where the uncomputation barrier thwarts the naive approach to QDSC. As discussed in <cit.>, the graph isomorphism problem can be reduced to QDSC as follows. Given a graph G with n vertices, one wishes to prepare the uniform superposition over all of its permutations,

    |Φ_G⟩ = 1/√(N(G))∑_π∈ S_n|π(G)⟩

where N(G) is the appropriate normalization factor.
If G and G' are isomorphic graphs then ⟨Φ_G | Φ_G'⟩ = 1. Otherwise ⟨Φ_G | Φ_G'⟩ = 0. These cases are easily distinguished using the Hadamard test. Furthermore, the uniform probability distribution over permutations of a graph is easy to sample from classically by sampling a permutation at random and then applying it to G.

Now, consider the approach to quantum state preparation that is analogous to the above (failed) attempt to make |Ψ_b⟩. For simplicity, we restrict attention to the case that G has no automorphisms, so the number of distinct permutations is n!. The hoped for state construction process would be:

    1/√(n!)∑_π∈ S_n|π⟩|G⟩   →   1/√(n!)∑_π∈ S_n|π⟩|π(G)⟩
       →   1/√(n!)∑_π∈ S_n|identity⟩|π(G)⟩.

The second step of unentangling the first register from the second by reverting it to the identity is achievable if, with knowledge of G, π can be computed from π(G). In this case, one  computes π from π(G) reversibly in superposition and applies π^-1 to the first register. However, computing π from π(G) is exactly the graph isomorphism problem. So we have not made any progress.

In <cit.> reductions are also given from quadratic residuosity and discrete logarithm to QDSC. In the case of discrete logarithm, the uncomputation barrier arises just as in the above two examples. In the case of quadratic residuosity, the uncomputation task that arises is to uncompute a superposition rather than an individual bit string. These examples show that quantum approaches to several hard problems are thwarted by uncomputation barriers. Our results in this paper provide an example of how an uncomputation barrier can be surmounted.







§ PRELIMINARIES


We now give the background necessary for understanding our construction of |Ψ_b⟩. Our procedure begins by constructing a state over the dual code so let us review what the dual is.

 Given a code C ⊂ℤ_2^n, the dual code C^⊥ is defined as

    C^⊥ = {𝐝∈ℤ_2^n : 𝐝·𝐜 = 0   ∀𝐜∈ C},

where 𝐝·𝐜 = ⊕_i=1^n d_i c_i is the ℤ_2^n inner product.

The dual to a code of length n and dimension k has length n and dimension n-k. We call the n-k by n generator matrix of the dual code B_⊥. Given B we can efficiently construct B_⊥, as discussed in appendix <ref>. 

On the quantum side each codeword and each dual codeword can be viewed as a computational basis state in a 2^n-dimensional space. The normalized uniform superposition of the 2^k codewords is related to the normalized uniform superposition of the 2^(n-k) dual codewords by the Hadamard transform.  Given any function f:ℤ_2^n →ℂ, its Hadamard transform f̃:ℤ_2^n →ℂ is defined by

    f̃(𝐱) = 1/√(2^n)∑_𝐲∈ℤ_2^n (-1)^𝐱·𝐲 f(𝐲).

This transform can be applied to the amplitudes of a quantum state by applying the Hadamard gate H = 1/√(2)[ [  1  1;  1 -1 ]] to each qubit:


    H^⊗ n∑_𝐱∈ℤ_2^n f(𝐱) |𝐱⟩ = ∑_𝐱∈ℤ_2^nf̃(𝐱) |𝐱⟩.



Given a uniform superposition over a code C, its Hadamard transform is the uniform superposition over the dual code C^⊥. That is:

    If    f(𝐱) = 1/√(|C|)∑_𝐜∈ Cδ_𝐱,𝐜   then   f̃(𝐱) = 1/√(|C^⊥|)∑_𝐝∈ C^⊥δ_𝐱,𝐝.




Writing out the Hadamard transform yields

    f̃(𝐱) =1/√(2^n)1/√(|C|)∑_𝐜∈ C  (-1)^𝐱·𝐜.

For 𝐱∈ C^⊥

    1/√(|C|)∑_𝐜∈ C (-1)^𝐱·𝐜 = √(|C|).

Recall that  |C| = 2^k and accordingly |C^⊥| = 2^n-k. Thus f̃(𝐱) = 1/√(|C^⊥|) for any 𝐱∈ C^⊥. For all 𝐱∉ C^⊥ one must have f̃(𝐱) = 0 since any additional nonzero amplitudes would yield norm greater than one, which cannot occur because the Hadamard transform is unitary.


Our algorithm for preparing |Ψ_b⟩ begins by constructing a preliminary state |Ψ_b⟩ and then applying the quantum Hadamard transform H^⊗ n. The state |Ψ_b⟩ is the Hadamard transform of the target state  |Ψ_b⟩, 

    |Ψ_b⟩ = H^⊗ n|Ψ_b⟩,

so we see that acting on this state with the Hadamard transformation gives the target state.

By a slightly involved calculation which we defer to sec:hadamard,

    |Ψ_b⟩ = N ∑_𝐝∈ C^⊥ K_b^n-1(|𝐝| -1 )  |𝐝⟩,

where K_j^n(x) is the Krawtchouk polynomial defined by

    K_j^n(x) = ∑_r=0^j xrn-xj-r (-1)^r,

and N is the normalization factor. If the balls do not overlap then

    N = 1/√(2^n-kVol(b)),

which is also shown in section  sec:hadamard.

If we can construct |Ψ_b⟩ we can act on it with the Hadamard transform and get our desired state,


    |Ψ_b⟩ = H^⊗ n|Ψ_b⟩.

So our task now is to give an efficient construction of  |Ψ_b⟩.
Suppose we can make the state

    ||Ψ_b|⟩ = N ∑_𝐝∈ C^⊥| K_b^n-1(|𝐝| -1) |  |𝐝⟩.

For a given dual codeword 𝐝∈ C^⊥ it is easy to compute the sign of the amplitude K_b^n-1(|𝐝| -1) and apply it using phase kickback, as described in appendix <ref>, to make the state eq:psiperp from eq:unsigned. So the task now is to make this unsigned state. 

The unsigned state can be described in terms of a probability distribution as follows. Let B_⊥∈ℤ^k^⊥× n be the generator matrix for C^⊥. Again k^⊥ = n - k. Then,

    ||Ψ_b|⟩ = ∑_𝐮∈ℤ_2^k^⊥√(p(𝐮))|𝐮B_⊥⟩,

where p(𝐮) is the normalized probability distribution over ℤ_2^k^⊥ given by

    p(𝐮) = N^2  [ K_b^n-1(|𝐮 B_⊥| -1) ]^2

States of the form ∑_𝐱√(p(𝐱))|𝐱⟩ can in some cases be prepared efficiently using the method of conditional rotations, described in the next section.



§ CONDITIONAL ROTATIONS BY METROPOLIS MONTE CARLO



Given a probability distribution p(u_1, …, u_k^⊥), the method of conditional rotations, introduced in <cit.> and further developed in <cit.>, can be used to  efficiently prepare the state

    ∑_𝐮∈ℤ_2^k^⊥√(p(𝐮))|𝐮⟩

provided certain marginal conditional probabilities can be efficiently computed. Here, our goal is to prepare this state in the case that p(𝐮) is as given in eq:pb. Once we have done this, it is straightforward to perform the transformation ∑_𝐮√(p(𝐮))|𝐮⟩→∑_𝐮√(p(𝐮))|𝐮 B_⊥⟩, using classical reversible computation, to get the desired state eq:target. This is spelled out in sec:alg. Our method for computing the necessary conditional marginal probabilities is Metropolis Monte Carlo.

The method of conditional rotations works as follows. Start with k^⊥ qubits initialized to the zero state, |0⟩^⊗ k^⊥. Suppose there is an efficient way to calculate the marginals p(u_1=0) and p(u_1=1) which are the probabilities that the first bit is 0 or 1. Then perform a rotation on the first qubit, obtaining

    (√(p(u_1=0))|0⟩ + √(p(u_1=1))|1⟩) |0⟩^⊗ (k^⊥-1).

Assuming there is an efficient way to compute subsequent conditional marginals, perform a rotation on the second qubit conditioned on the value of the first qubit, obtaining 

    [              [ √(p(u_1=0))|0⟩(√(p(u_2=0|u_1=0))|0⟩+√(p(u_2=1|u_1=0))|1⟩) + .; . √(p(u_1=1))|1⟩(√(p(u_2=0|u_1=1))|0⟩ + √(p(u_2=1|u_1=1))|1⟩) ] |0⟩^⊗(k^⊥-2) ].

Similarly, rotate the third qubit conditioned on the first two according to p(u_3|u_2,u_1) and so on until all the bits are rotated. At the end of this process
one is left with the state 

    ∑_u_1 … u_k^⊥ a(u_1, …, u_k^⊥) |u_1 … u_k^⊥⟩

where

    a(u_1,…,u_n)     =    √(p(u_1))√(p(u_2|u_1))×…×√(p(u_k^⊥|u_k^⊥-1,…,u_1))
        =    √(p(u_1,…,u_k^⊥))

as desired.

For this process to be efficient, it must be possible to efficiently evaluate all the necessary conditional marginal probabilities in which the first m bits are fixed and the last k^⊥-m bits are summed over, keeping track of how often bit m+1 is 0 or 1 to get the marginal. Classically precomputing these sums would be infeasible because there are exponentially many of them. Instead, these evaluations are  done in superposition by reversible circuits. Thus, the key question is whether these conditional marginal probabilities can each be computed by an efficient classical algorithm.

The classical algorithm that we propose for estimating these conditional probabilities is the method of Markov chain Monte Carlo, using the Metropolis rule. We find that the necessary Markov chains converge remarkably rapidly, despite the exponentially large state space of codewords in C^⊥ for certain values of b and k.  We discuss the conditions that guarantee good convergence in sec:largeb.

Consider first the Markov chain used to estimate the rotation probability for the first bit. At each step in the Markov chain our random walker is located at some string 𝐮∈ℤ_2^k^⊥. We then choose one bit of 𝐮 unformly at random and flip it, yielding 𝐮'. Since we are going for eq:pb we then evaluate:

    f(𝐮') = [ K^n-1_b(|𝐮' B_⊥|-1)]^2.

If  f(𝐮') ≥ f(𝐮) this move is accepted and the walker moves from 𝐮 to 𝐮'. Otherwise, the move is accepted with probability f(𝐮')/f(𝐮). This is the Metropolis rule. Under mild conditions, this is guaranteed to converge to the limiting normalized distribution proportional to f.

One obtains a sequence of bit strings 𝐮, 𝐮', …∈ℤ_2^k_⊥ by running this Markov chain. The fraction of bit strings sampled in which the first bit is 1 is used as the marginal distribution for the first qubit rotation. 
For the mþ qubit rotation the same procedure is followed except that the first m bits of 𝐮 are frozen to specific values, and only flips of remaining k^⊥-m bits are proposed in the random walk, thereby sampling from the necessary conditional distribution. This is done reversibly in superposition over values that the first m bits are fixed to. Reversible implementation of the Markov chain is to be carried out using a reversible pseudorandom number generator.



§ QUANTUM STATE PREPARATION ALGORITHM FOR |Ψ_B⟩



We now present pseudocode for the construction of |Ψ_b⟩, which is an example of the Quantum Distributional State Construction problem.  It relies on sampling by Metropolis Monte Carlo which we described in the last section.
The key equation eq:psiperp is established later so if you take it on faith hopefully you have enough to follow the pseudocode.

Given a generator matrix for C, we construct |Ψ_b⟩ using the following steps.



    
  * Precompute K_b^n-1(h-1) for h=0, …, n so they need not be evaluated with the quantum device.  Note that they do not depend on the code.

    
  * Starting from the generator matrix for C, compute B_⊥, the k^⊥× n generator matrix for C^⊥. This can be done in polynomial time, see for example appendix <ref>.

    
  * Initialize a register of k^⊥ qubits and second register of n qubits to the all zero state,
    
    |0⟩^⊗ k^⊥|0⟩^⊗ n.


    
  * Use the method of conditional rotations, estimated using the Metropolis Monte Carlo method, to construct a weighted superposition over the coefficient vectors 𝐮,
    
    → N ∑_𝐮∈ℤ_2^k^⊥ | K_b^n-1 ( |𝐮 B_⊥|-1)  |  |𝐮⟩|0⟩^⊗ n.

    where N is the normalization constant.
    Note that our method automatically produces a normalized state so there is no need to know N in advance.

    
  * Reversibly compute the corresponding codewords of C^⊥ into the second register,
    
    → N ∑_𝐮∈ℤ_2^k^⊥ | K_b^n-1(|𝐮B_⊥|-1)  |  |𝐮⟩|𝐮 B_⊥⟩.


    
  * Uncompute the first register, returning it to |0⟩^⊗ k^⊥. We do this after having classically precomputed the right-inverse B_⊥^-1 for the (non-square) matrix B_⊥, that is, B_⊥ B_⊥^-1 = 𝕀.  Then in the first register add (𝐮 B_⊥) B_⊥^-1 to 𝐮. This yields
    
    → N ∑_𝐮∈ℤ_2^k^⊥ | K_b^n-1(|𝐮B_⊥|-1)  |  |0⟩|𝐮 B_⊥⟩.

    This can be rewritten as
    
    = N ∑_𝐝∈ C^⊥ | K_b^n-1(|𝐝|-1)  |  |𝐝⟩,

    where we have discarded the first register for clarity.

    
  * Use phase kickback (c.f. appendix <ref>) to put in the correct sign of K_b^n-1(|𝐝|-1),
    
    → N ∑_𝐝∈ C^⊥ K_b^n-1(|𝐝|-1)  |𝐝⟩.


    
  * Perform the quantum Hadamard transform H^⊗ n. As we show next in sec:hadamard, this yields the desired final state,
    
    →𝒩∑_𝐜∈ C∑_𝐳∈ℤ_2^nΘ(b-|𝐳|)  |𝐳⊕𝐜⟩ = |Ψ_b⟩.
 


The fidelity with which the final state is produced depends entirely on the precision to which the Metropolis Monte Carlo method can estimate the conditional probabilities defining the qubit rotations. This is step 4. All other steps of the algorithm are exact.



§ DERIVING THE FORM OF |Ψ_B⟩



In this section we collect and prove some facts regarding Hadamard transforms and then use them to derive eq:psiperp.


    
    
    H^⊗ n1/√(2^n)∑_𝐳∈ℤ_2^n∑_𝐲∈ℤ_2^n f(𝐲) g(𝐳⊕𝐲) |𝐳⟩ = ∑_𝐱∈ℤ_2^nf̃(𝐱) g̃(𝐱) |𝐱⟩

    where f̃ is the Hadamard transform of f
    
    f̃(𝐱) = 1/√(2^n)∑_𝐲∈ℤ_2^n (-1)^𝐱·𝐲 f(𝐲)

    and g̃ is the Hadamard transform of g.




    H^⊗ n1/√(2^n)∑_𝐳∈ℤ_2^n∑_𝐲∈ℤ_2^n f(𝐲) g(𝐳⊕𝐲) |𝐳⟩ = 1/2^n∑_𝐱,𝐲,𝐳∈ℤ_2^n (-1)^𝐱·𝐳 f(𝐲) g(𝐳⊕𝐲 ) |𝐱⟩

Let 𝐰 = 𝐳⊕𝐲. Then the previous expression can be rewritten as


    = 1/2^n∑_𝐱,𝐲,𝐰∈ℤ_2^n (-1)^𝐱· (𝐰⊕𝐲) f(𝐲) g(𝐰) |𝐱⟩
       = ∑_𝐱∈ℤ_2^n( 1/√(2^n)∑_𝐲∈ℤ_2^n (-1)^𝐱·𝐲 f(𝐲) ) ( 1/√(2^n)∑_𝐰∈ℤ_2^n (-1)^𝐱·𝐰 g(𝐰) ) |𝐱⟩
       = ∑_𝐱∈ℤ_2^nf̃(𝐱) g̃(𝐱) |𝐱⟩.




    
    Given a radius b ∈{0,1,2,…,n}, let g:ℤ_2^n →ℂ be
    
    g(𝐱) = Θ(b-|𝐱|),

    where Θ is the unit step function as defined in eq:unitstep. Then, the Hadamard transform of g is
    
    g̃(𝐱) = 1/√(2^n) K_b^n-1(|𝐱|-1),

    where K_b^n-1 is the Krawtchouk polynomial as defined in eq:Krawtchouk.



To compute the Hadamard transform of a ball, first rewrite
g(𝐲) = Θ(b-|𝐲|) as g(𝐲) = ∑_j = 0^b δ_j,|𝐲|. Then


    g̃(𝐱) = 1/√(2^n)∑_j=0^b ∑_𝐲∈ℤ_2^n (-1)^𝐱·𝐲δ_j,|𝐲|.

The summand (-1)^𝐱·𝐲δ_j,|𝐲| depends on 𝐱 and 𝐲 only through r, the number of bits which are equal to 1 in both 𝐱 and 𝐲, and s, the number of bits which are equal to 1 in 𝐲 but not 𝐱. Thus, 

    g̃(𝐱) = 1/√(2^n)∑_j=0^b ∑_r = 0^|𝐱|∑_s = 0^n-|𝐱||𝐱|rn - |𝐱|s (-1)^r δ_j,r+s


    = 1/√(2^n)∑_j=0^b ( ∑_r=0^j |𝐱|rn - |𝐱|j-r (-1)^r ).

The inner expression is recognizable as a Krawtchouk polynomial, which is defined by 

    K_j^n(x) = ∑_r=0^j xrn-xj-r (-1)^r

thus we can write 

    g̃(𝐱) = 1/√(2^n)∑_j=0^b K_j^n(|𝐱|).

Using the identity 

    ∑_j=0^b K_j^n(x) = K_b^n-1(x-1),

which we derive in appendix <ref>, we can further simplify our expression to 

    g̃(𝐱) = 1/√(2^n) K_b^n-1(|𝐱|-1).



With these facts in hand, we can derive eq:psiperp as follows. Examination of eq:psib shows that |Ψ_b⟩ is proportional to a convolution in the form of eq:conv. Namely,

    |Ψ_b⟩ = √(2^n/Vol(b))1/√(2^n)∑_𝐳∈ℤ_2^n∑_𝐲∈ℤ_2^n f(𝐲) g(𝐳⊕𝐲) |𝐳⟩,

where

    f(𝐱) = 1/√(2^k)∑_𝐜∈ Cδ_𝐜,𝐱

and

    g(𝐱) = Θ(b-|𝐱|).

The convolution theorem therefore tells us that

    H^⊗ n|Ψ_b⟩ = √(2^n/Vol(b))∑_𝐱∈ℤ_2^nf̃(𝐱) g̃(𝐱) |𝐱⟩.

Now the Hadamard transform of f(𝐱) is given by eq:cperp and  the Hadamard transform of g(𝐱) is given by  eq:gtilde so plugging in gives


    |Ψ_b⟩  = H^⊗ n|Ψ_b⟩= N ∑_𝐝∈ C^⊥ K_b^n-1(|𝐝| -1 ) |𝐝⟩

where N is given as

    N = 1/√(2^n-kVol(b)).

So we have established eq:psiperp.



§ IMPACT OF IMPERFECT CONVERGENCE ON FIDELITY



 In sec:convergence we show computer experiments to see how well a Metropolis Monte Carlo algorithm does at sampling from a Krawtchouk distribution at n=1000 with 10^10 samples.  In this section we quantify the effect of sampling with a finite number of samples versus the ideal limiting distribution. We calculate the fidelity between the ideal quantum state and one that arises from imperfect sampling assuming that the only source of error is the imperfect sampling. We do assume that the sampled distribution is uniform over codewords of a given Hamming weight, as it is in the ideal case. The result is eq:impover in terms of the ideal and sampled distributions and is evaluated in specific cases in the next section.

Let

    |h⟩ = 1/√(W(h))∑_𝐝∈ C^⊥
     |𝐝| = h|𝐝⟩,

where W(h) is the number of codewords in C^⊥ with Hamming weight h. Note that we are not assuming anything about W(h). Then

    ⟨ h | h' ⟩ = δ_h,h'.

By eq:unsigned,

    ||Ψ_b|⟩ = N ∑_h=0^n √(W(h))| K^n-1_b(h-1) |  |h⟩.

The probability distribution over Hamming weights defined by this state is

    p_ideal(h) = [ K^n-1_b(h-1) ]^2 W(h)/∑_h=0^n [ K^n-1_b(h-1) ]^2 W(h).

We can correspondingly rewrite eq:psiperph as

    ||Ψ_b|⟩ = ∑_h=0^n √(p_ideal(h))|h⟩.

Let p_sam(h) be the probability distribution over Hamming weights achieved from the product of conditional probabilities estimated by Markov chains, as in eq:markov1 and eq:markov2. Assuming no other sources of error, the unsigned state resulting from conditional rotations would be

    | | Ψ_b^sam |⟩ = ∑_h=0^n √(p_sam(h))|h⟩.

Here is where we assumed that the distribution over sampled dual codewords at fixed Hamming weight is uniform.
By eq:orthonormal, the inner product between 
| |Ψ_b^sam|⟩ and | |Ψ_b | ⟩ is

    ⟨ | Ψ_b^sam | | . |Ψ_b| ⟩ = ∑_h=0^n √(p_sam(h))√(p_ideal(h)).

The minus signs imposed by phase kickback are always exact, and so the signed state will have this same inner product

    ⟨Ψ_b^sam| . Ψ_b ⟩ = ∑_h=0^n √(p_sam(h))√(p_ideal(h)).

Let |Ψ_b^sam⟩ be the approximation to |Ψ_b⟩ obtained by taking the Hadamard transform of |Ψ_b^sam⟩,

    |Ψ_b^sam⟩ = H^⊗ n|Ψ_b^sam⟩.

Now | Ψ_b ⟩ is the Hadamard of |Ψ_b⟩ so the fidelity of the ideal state  with the sampled state is

    ⟨Ψ_b^sam | Ψ_b ⟩ = ∑_h=0^n √(p_sam(h))√(p_ideal(h)).

By doing numerical experiments where we classically sample to get p_sam(h) we can estimate the quantum fidelity between the less than perfect quantum state and the ideal.



§ EMPIRICAL CONVERGENCE OF METROPOLIS MONTE CARLO



To see if the Monte Carlo algorithm converges with not too many samples we consider the Markov chain used in the first step of the method of conditional rotations. That is, none of the bits of 𝐮 are fixed and the state space consists of all 2^k^⊥ bit strings 𝐮. Any subsequent Monte Carlo step with m bits of 𝐮 frozen is essentially a version of the same problem with m rows of the generator matrix for C^⊥ removed. As we will see in sec:largeb, failures of convergence can arise when one or both of the ball radius b and the Hamming distance f of the step size in the random walk are too large. The step size f is the Hamming weight of a row of the generator matrix B_⊥. As B_⊥ is in systematic form, these rows have Hamming weight approximately k/2. In the mþ random walk, we have removed m of these rows. This leaves b unchanged and f still has an expected size of k/2.  So we expect that if the first step converges, the subsequent ones do as well.

For our examples we  use random codes. Because this state space is exponentially large, it is not practical to verify the convergence of the Metropolis Markov chain by directly comparing the sampling frequencies to eq:pb. Instead, we measure the frequency at which elements of C^⊥ with given Hamming weight are sampled. If the algorithm is converged, these sampling frequencies should be proportional to

    W(h) [ K_b^n-1(h-1)]^2

where W(h) is the number of codewords in C^⊥ with Hamming weight h. It is well known that random linear codes have a weight distribution which is close to binomial

    W(h) ≃2^k^⊥/2^nnh

so we take

    p_ideal(h) ∝nh[ K_b^n-1(h-1)]^2.

We do not need the explicit form of the normalization for what follows.

We want to see if sampling can produce the distribution in eq:ph.  We show two numerical examples. The first is at n=1000, k=100, b=20 and is shown in figure <ref>. 
The results show close agreement between the measured sampling frequencies as a function of Hamming weight and the prediction of eq:ph. This close agreement suggests both that the Metropolis algorithm is well converged and that the binomial distribution is a good approximation to W(h).



Figure <ref> shows that  p_ideal(h) is very small outside of a limited range of Hamming weights. In the example shown p_ideal(h) is very small for Hamming weights less than roughly 340 or greater than roughly 660. This is significant, because finding codewords of low (or high) Hamming weight is an approximate version of the lightest codeword problem which is classically hard. Algorithms  that require  sampling from dual codewords of very low Hamming weight are unlikely to be successful.  In this example, for our task of constructing the state eq:psib, we succeed without going near low Hamming weight dual codewords.  

 For a linear code, we can estimate the minimum Hamming weight of any nonzero codeword. As shown in <cit.>, if a generator matrix is chosen uniformly at random then with high probability the resulting code will have minimum distance D saturating the Gilbert-Varshamov bound. That is, D is obtained by solving

    2^k = 2^n/Vol(D) = 2^n/∑_j=0^D nj.

In our example the dual is a random code, and we estimate that minimum Hamming weight of any nonzero dual codeword is approximately 14. We see that we are not sampling near the smallest dual codewords and we do not need to.

By choosing larger k and b we can construct examples in which the support of p_ideal(h) extends to dual codewords of  low and high Hamming weight which take exponential time to reach. In such cases, the Markov chain converges nicely to the correct relative probabilities for the dual codewords of Hamming weight in the middle range but the outer tails are cut off. We illustrate this in figure <ref> where n=1000, k=300, and b=60.  In figure <ref> we show the central region but we rescale the two distributions so they both sum to one in this limited range.  Note the excellent agreement between the sampled and the ideal in the middle range.  This results in high fidelity between the ideal quantum state and the sampled quantum state. A quantitative analysis of which values of n,k,b are expected to yield fully-converged Markov chains as in figure <ref> versus which are expected to yield cut-off distributions as in figure <ref> is given in sec:largeb.

Note that, even in this second example (n=1000, k=300, and b=60), the ideal  distribution over dual codewords has negligible probability at Hamming weight below approximately 263, whereas the shortest nonzero dual codeword for a random code of these parameters is, with high probability, approximately 65. So achieving full convergence would not require sampling  near the shortest nonzero codeword.

Another phenomenon that we frequently observe in the unconverged regime is that some codewords just at the boundaries of the regions of accessible Hamming weights are sampled with frequencies that appear to be anomalously high.  An example of this phenomenon can be seen in figure <ref> in which there are a few outlier red points showing sampling frequencies substantially higher than expected.

By the results of sec:imperfect, namely eq:impover, we can compute the fidelity between a state arising from the sampled probability distributions and the ideal state, under the assumption that the sampling of codewords of a given Hamming weight is uniform as in the ideal case. We are ignoring the error that occurs from successive conditional rotations.  For the well converged example at n=1000,k=100,b=20 from figure <ref> we obtain a fidelity of 0.999995. For the example at n=1000,k=300,b=60 shown in figure <ref> we obtain fidelity 0.67. These overlaps are high and show that in practice we should be able to make the desired states with good fidelity.









§ ANALYSIS OF CONVERGENCE OF METROPOLIS MONTE CARLO



In this section we derive a criterion to predict which values of n,k,b will yield good convergence, and which will not. We then compare these predictions against numerical experiments.

Based on the form of p_ideal(h) illustrated in figure <ref>, we identify two circumstances that cause concern regarding the convergence of our Markov chain to its limiting distribution. The first we will show is actually not problematic. The second can prevent full convergence for some values of n,k,b that we fully characterize.

First, one expects the random walk to fail if p_ideal(h) contains regions of high probability separated by a region of very low probability with width greater than the Hamming distance of the trial steps. We would expect failure in this case because the Metropolis walk would need to go from one region to the other only by stepping across a region of low probability. The Metropolis walk will not accept moves into very low probability regions so they will create “bottlenecks” for the Metropolis walk.

As shown in figure <ref>,  p_ideal(h) oscillates, and has many points of low probability. If the Metropolis walk were only allowed to make moves of Hamming distance one, these low probability codewords would create bottlenecks, which might cause the Metropolis walk to fail. However, the walk actually moves by taking steps of Hamming distance one on the coefficient vector 𝐮. A step 𝐮→𝐮' where |𝐮' - 𝐮| = 1 corresponds to a step from the dual codeword 𝐮B_⊥ to the dual codeword 𝐮'B_⊥. That is, the new dual codeword is obtained by adding one row of B_⊥ to the current dual codeword. The distance of the move will be the Hamming weight of the row. When B_⊥ is in systematic form, as defined in appendix <ref>, the typical Hamming weight of a row is is approximately (1/2)(n-k^⊥) = k/2. By taking steps of Hamming distance k/2, the Metropolis walk can easily step over regions of very low probability as is confirmed in our numerical examples.

It is known that the Krawtchouk polynomial K_b^n(x) has b real zeros which all lie in the range n/2-√(b (n-b))≤ x ≤ n/2+√(b (n-b)) <cit.>. So the average spacing between the  zeros is 2√((n-b)/b) and numerically we see that the zeros are roughly uniformly spaced. Many of the values of x where the exact zeros occur are not at integers but as we can see from figure <ref> at the integer values near the exact zeros the polynomial is tiny.
As discussed above, a bottleneck could occur if there were a “trough” in which p_ideal(h) is consistently low over some region whose width is larger than the steps taken by the random walker. So  we can see that unless 2√((n-b)/b)≫ k/2, the oscillations in p_ideal(h) are of short enough wavelength that any troughs will be much narrower than the step length for the walker and this type of bottleneck will not occur.  Accordingly in  our numerical experiments we do not see this failure mode.

The second circumstance in which we might expect the Markov chain to fail to converge is when the distribution p_ideal(h) assigns nonnegligible probability to codewords of Hamming weight near either zero or n, which become exponentially hard to reach by the sampling process. (Note that this is not the problem of finding short codewords which are difficult to discover.  It is a problem of failure to sample.)
If a random walker is sitting at a bit string of very low Hamming weight, a move that is proposed in the Metropolis algorithm which flips roughly k/2 bits will be much more likely to raise the Hamming weight than lower it. And the reverse holds at very high Hamming weight. If the edges of the nonnegligible region for p_ideal(h) extend beyond these entropic barriers, then we expect the random walk to converge incompletely. After any reasonable number of moves in the random walk, the number of times the walker reaches the regions beyond the entropic barrier will remain tiny, thereby cutting off the tails of the desired distribution p_ideal(h).  We see this in figure  <ref>.

Unlike bottlenecks, which we believe never occur, this second failure mode can occur when b or k is large. We can quantify this as follows. Consider a random walker at bit string 𝐝 that has Hamming weight h. At a given step of the Metropolis algorithm a move is proposed by adding a random row of B_⊥, which will have Hamming weight f ≃ k/2 since B_⊥ is in systematic form. We model this as flipping f bits of 𝐝 chosen uniformly at random.  Let j be the number of ones in 𝐝 that flip and accordingly the number of zeros that flip is f-j. The Hamming weight of the flipped string is (h-j)+(f-j) so to decrease the Hamming weight we need j > f/2. Counting the number of ways that this can happen we get that the probability for the proposed move to decrease the Hamming weight is

    P_down(h) = ∑_j=⌈ f/2 ⌉^f hjn-hf-j/nf.
 
If h becomes too small relative to n and f then this probability becomes exponentially small. The value of h below which P_down becomes insurmountably small relative to the number of steps in the Markov chain defines an entropic barrier which the random walker is unlikely to cross. One also has the analogous situation at h close to n for which  P_up becomes insurmountably small.

The existence of this entropic barrier only prevents convergence if the region to which p_ideal(h) assigns nonnegligible probability extends outside this entropic barrier. Using the results of <cit.> on the asymptotics of Krawtchouk polynomials, we find that the region on which p_ideal(h) is nonnegligible is

    n/2-√(b(n-b))≤ h ≤ n/2+√(b(n-b)).

Since this region is symmetric about n/2, the region of nonnegligible p_ideal(h) extends above the high-Hamming-weight entropic barrier  for the same values of n,k,b in which it extends below the low-Hamming-weight barrier. Hence, given n,k,b we can predict whether the Markov chain will converge by solving P_down(h) = ϵ for h using some appropriately small ϵ. If the resulting value h_barrier(ϵ) satisfies h_barrier(ϵ) < n/2-√(b(n-b)) then the Markov chain can be expected to converge, and otherwise not.

Let us now compare this prediction against computer experiments. The location of the entropic barrier is not highly sensitive to the choice of ϵ so we semi-arbitrarily set ϵ=10^-6. Using this criterion, for any fixed n we can map out the k-b plane, classifying each point into one of three categories. The first category is where the left entropic barrier h_barrier(ϵ) is below the Krawtchouk edge at n/2-√(b(n-b)). In this case we expect efficient convergence across the full range of nonnegligible p_ideal(h). The second region is where the entropic barrier is above the Krawtchouk edge. In this case we expect the sampled distribution p_sam(h) to look like p_ideal(h) with the outer tails cut off. The third region is where the ball radius exceeds the minimum distance of the code C, that is 2b > D. This is the regime of overlapping balls, which we do not consider in this paper. These regions are plotted for n=1000 in figure <ref>.

We next examine a slice through these regions along the line k=5 × b. We choose a set of seven evenly-spaced points along this line, and at each one generate ten random codes and carry out our Markov chain for each one. The analysis above predicts that the first four points, which lie in the fully-convergent region, should yield sampled distributions over Hamming weights with high fidelity to the ideal distribution. The last three points, which lie in the cut-off region, should have fidelity falling off due to an increasing fraction of the tails of the distribution being cut off. Two of these seven points, namely k=100,b=20 and k=300,b=60 are used in sec:convergence to illustrate in more detail the behavior in the convergent and cut-off regions, respectively.






The results are shown in figure <ref>, which show good qualitative agreement with these predictions.
For our trials at the three points in the cut-off region, the fidelity, as calculated with eq:impover, is greater than 0.5, as shown in figure <ref>. This is higher than one might naively expect since a look at figure  <ref>  suggests that the sampled and ideal distributions are not close.  However inside of the entropic barrier where p_sam is not tiny, the two distributions are in fact proportional, that is, one is scale factor times the other.  See figure <ref> where both have been normalized over the plotted region.  The quantum state overlap  eq:impover has support only inside the entropic region where the two distributions are proportional and the fidelity is then the square root of the scale factor which is not very small. 





§ OVERLAP OF |Ψ_B⟩ WITH T_𝐯|Ψ_B⟩



Our original motivation for finding a way to construct |Ψ_b⟩ is that the overlap of this state with its translation by 𝐯 gives information about how close 𝐯 is to the  lattice.  We do not claim that we have an  efficient quantum algorithm for BDD. Still we want to understand how this overlap depends on how far 𝐯 is from the lattice and in this section we derive an expression for the ideal state overlap ⟨Ψ_b| T_𝐯|Ψ_b⟩. This determines the runtime needed to solve BDD using Hadamard tests, which we discuss in sec:bddimp.

Let δ be the displacement of 𝐯 from the nearest codeword. That is

    δ = 𝐯⊕𝐜

where 𝐜 is the element of the code C that minimizes | 𝐯⊕𝐜|.
If each ball in |Ψ_b⟩ overlaps with only one ball from T_𝐯|Ψ_b⟩ then ⟨Ψ_b| T_𝐯|Ψ_b⟩ is equal to the overlap between a pair of single normalized  balls displaced by δ. This is because there are 2^k overlapping ball pairs and there is the same factor in the denominator from the overall normalization. We denote the overlap of a single ball with its translate by δ as 𝒜(δ).

Since the triangle inequality holds for Hamming distance, one can observe that if 2b + |δ| is less than the minimum distance of the code C then it is guaranteed that each ball in |Ψ_b⟩ will overlap with at most one ball from T_𝐯|Ψ_b⟩. Furthermore, if 4b is less than the minimum distance of C then this will hold for all 𝐯, as illustrated in figure <ref>. In these cases ⟨Ψ_b| T_𝐯|Ψ_b⟩ = 𝒜(δ).  
In our first example in sec:convergence we have n=1000, k=100 and  b=20. We estimate the length of the smallest codeword is 320 so we are safely in the regime of 4b < D. For our next example at n=1000, k=300 and b=60 we estimate that the smallest codeword has length 192 so for δ up to 72 we have 2b + |δ| < D and so each ball in |Ψ_b⟩ will overlap only one ball from T_𝐯|Ψ_b⟩. 



By definition

    𝒜(δ) = ⟨ B_1 | B_2 ⟩

where, for j=1,2

    |B_j⟩    =    1/√(|B_j|)∑_𝐱∈ B_j|𝐱⟩
    
        B_j     =    {𝐱∈ℤ_2^n : |𝐱⊕𝐳_j| ≤ b }

and 𝐳_1,𝐳_2 are any pair of bit strings satisfying

    𝐳_1 ⊕𝐳_2 = δ.

Now

    ⟨Ψ_b| T_𝐯|Ψ_b⟩ = ⟨ B_1 | B_2 ⟩ = |B_1 ∩ B_2|/√(|B_1| · |B_2|) = |B_1 ∩ B_2|/ |B_1| ,

where we used that |B_1| = |B_2|.
Next, observe that

    |B_1| = ∑_j=0^b nj.

To compute |B_1 ∪ B_2| we can assume without loss of generality that

    𝐳_1     =    00…000…0^n
    𝐳_2     =    11…1_δ00…0_n-δ

where δ = |δ|, which is equal to the distance from 𝐯 to the nearest codeword.

Starting from 𝐳_2 one generates all elements of B_2 by flipping s_1 bits from the first δ and s_2 bits from the last n-δ, where s_1 + s_2 ≤ b. The number of ways to do this is δs_1n-δs_2. The Hamming weight of the resulting string is δ-s_1+s_2. Therefore,

    | B_1 ∩ B_2 | = ∑_(s_1,s_2) ∈ Pδs_1n-δs_2

where

    P = { (s_1,s_2) : s_1 ≥ 0, s_2 ≥ 0, s_1+s_2 ≤ b , and δ-s_1+s_2 ≤ b }.

Putting it all together, we have

    𝒜(δ) = ∑_(s_1,s_2) ∈ Pδs_1n-δs_2/∑_j=0^b nj

with P as given in eq:P. In the case that each ball overlaps with at most one other,

    ⟨Ψ_b| T_𝐯|Ψ_b⟩ = ∑_(s_1,s_2) ∈ Pδs_1n-δs_2/∑_j=0^b nj,

where δ is the Hamming distance from 𝐯 to the nearest codeword.



§ DEQUANTIZING THE CALCULATION OF ⟨Ψ_B| T_𝐯|Ψ_B⟩



Here we show that there is a classical sampling procedure that one can run to estimate ⟨Ψ_b| T_𝐯|Ψ_b⟩ by sampling a ± 1 variable and averaging the result. The number of samples needed to distinguish this quantity from zero is of order 1/|⟨Ψ_b| T_𝐯|Ψ_b⟩|^2, which has the same scaling as the number of repetitions of the Hadamard test that are needed in the quantum setting. The quantity ⟨Ψ_b| T_𝐯|Ψ_b⟩ is exponentially small in the distance of 𝐯 from the lattice so both the quantum and dequantized algorithms require time exponential in this distance to distinguish ⟨Ψ_b| T_𝐯|Ψ_b⟩ from zero. The dequantization shown here is modeled on the dequantization results regarding lattices on ℝ^n shown in <cit.>.

To dequantize the estimation of ⟨Ψ_b| T_𝐯|Ψ_b⟩ we first write T_𝐯 in terms of Pauli X operators which flip bits where 𝐯 is 1,

    T_𝐯 = ∏_i   s.t.  v_i=1 X_i .

Next, we recall that |Ψ_b⟩ = H^⊗ n|Ψ_b⟩. Since H X_i H = Z_i we have

    ⟨Ψ_b| T_𝐯|Ψ_b⟩    =    ⟨Ψ_b|∏_i   s.t.  v_i=1 X_i |Ψ_b⟩
        =    ⟨Ψ_b| H^⊗ n∏_i   s.t.  v_i=1 X_i H^⊗ n|Ψ_b⟩
        =    ⟨Ψ_b|∏_i   s.t.  v_i=1 Z_i |Ψ_b⟩.

Recalling the expression eq:psiperp for |Ψ_b⟩ this yields

    ⟨Ψ_b| T_𝐯|Ψ_b⟩ = N^2 ∑_𝐝∈ C^⊥[ K_b^n-1(|𝐝|-1) ]^2 ⟨𝐝|∏_i   s.t.  v_i=1 Z_i |𝐝⟩.

The quantity N^2 [ K_b^n-1(|𝐝|-1) ]^2 is the probability distribution over dual codewords that we sample from in our Markov chain. See eq:pb. In a slight abuse of notation we define

    p_ideal(𝐝) = N^2 [ K_b^n-1(|𝐝|-1) ]^2,

where N is the normalization factor eq:Ndef. The previously defined probability p_ideal(h) over Hamming weights is related to this by

    p_ideal(h) = ∑_𝐝∈ C^⊥
     |𝐝|=h p_ideal(𝐝).

For any given codeword, the quantity ⟨𝐝|∏_i   s.t.  v_i=1 Z_i |𝐝⟩∈{+1,-1} is easy to compute. Namely

    ⟨𝐝|∏_i   s.t.  v_i=1 Z_i |𝐝⟩ = (-1)^𝐝·𝐯.

Consequently,

    ⟨Ψ_b| T_𝐯|Ψ_b⟩ = ∑_𝐝∈ C^⊥ p_ideal(𝐝) (-1)^𝐝·𝐯.


The formula eq:dequant yields the following dequantization of our procedure for estimating ⟨Ψ_b| T_𝐯|Ψ_b⟩. Use Metropolis Monte Carlo to sample dual codewords according to p_ideal. For each dual codeword sampled, compute (-1)^𝐝·𝐯. After S samples, the average of this list of ± 1 values will be an estimate of ⟨Ψ_b| T_𝐯|Ψ_b⟩ accurate to ± O(1/√(S)), provided n,k,b are in the regime where the Markov chain converges. This is the same regime for which our quantum state preparation works and is characterized in sec:largeb.



§ BOUNDED DISTANCE DECODING WITH HADAMARD TESTS



To distinguish an overlap of ϵ from an overlap of zero via Hadamard tests requires order 1/ϵ^2 repetitions, as dictated by standard sampling statistics. Thus, 1/|⟨Ψ_b| T_δ|Ψ_b⟩|^2 provides an estimate of runtime for this task as a function of δ = |δ|. Neglecting sampling error, |⟨Ψ_b| T_δ|Ψ_b⟩| is given by eq:ideal_overlap. Given an empirical distribution over Hamming weights of dual codewords observed in our Monte Carlo simulations, we can also calculate |⟨Ψ_b| T_δ|Ψ_b⟩| making only the assumption that dual codewords of given Hamming weight remain equiprobable. We omit this calculation, as it is much more involved than the calculation of the ideal overlap given in sec:overlap.

But we show the runtimes in the presence of empirically observed sampling errors in figures <ref> and <ref> alongside the ideal runtimes achievable without sampling errors.

In figures <ref> and <ref> we compare these runtimes against classical brute force search, which has runtime nδ. We also compare against Information Set Decoding, a more sophisticated classical randomized algorithm described in sec:ISD. We find that the quantum approach outperforms brute force search but falls short of Information Set Decoding. All three of these algorithms could be quadratically sped up using Grover search, leaving their relative performance unaffected. 









§ INFORMATION SET DECODING



In this section we discuss a classical algorithm for BDD, specifically the most basic version of Information Set Decoding. This provides a useful point of reference because it is a simple classical algorithm that provides nearly state of the art performance. More recently, other more complicated classical algorithms have been devised that achieve slightly improved (but still exponential) scaling relative to ISD, e.g. <cit.>. Some noteworthy features of ISD are that it uses only polynomial memory and it can be quadratically accelerated through Grover search <cit.>.

In a bounded distance decoding problem, we are given a generator matrix B ∈ℤ_2^k × n for a code C and a vector 𝐯 which is near an unknown codeword 𝐜, that is, 𝐯 = 𝐜⊕δ∈ℤ_2^n. Here δ is the error vector satisfying |δ| = δ with δ not too large. So 𝐮 B = 𝐜 for some 𝐮∈ℤ_2^k. The goal is to find 𝐜, or equivalently 𝐮.

The ISD algorithm works as follows.


    
  * Select k elements of {1,2,…,n} uniformly at random. Let 𝐯^(k) be the length-k string obtained by keeping only these bits from 𝐯 and let B^(k) be the k × k submatrix of B obtained by keeping only these columns of B.
    
  * Check whether B^(k) is an invertible matrix over the field ℤ_2. If not, go back to step 1 and try a different subset.
    
  * Solve the linear system 𝐮̃ B^(k) = 𝐯^(k) for 𝐮̃∈ℤ_2^k using Gaussian elimination.
    
  * Compute the Hamming distance | 𝐮̃ B ⊕𝐯 |. If this is equal to δ then 𝐮̃ B is a codeword of distance δ from 𝐯 and we are done. Otherwise, go back to step 1 and try again with a different subset.


Next, let us analyze the probability of success for each trial. The probability that a uniformly random k × k matrix over the field ℤ_2 is invertible, in the limit of large k, is approximately 0.29 <cit.>. For generic dense B this is also an accurate prediction of the approximate probability that the random submatrix B^(k) will be invertible in step 2.

Since 𝐯 = 𝐜⊕δ, the linear system we solve in step 3 can be rewritten as

    𝐮̃ B^(k) = 𝐜^(k)⊕δ^(k).

If we are lucky, the subset of bits we chose contains no errors. That is δ^(k) = 0. In this case,

    𝐮̃ B^(k) = 𝐜^(k).

Next, recall that since 𝐜 is a codeword there is a 𝐮 such that 𝐮 B = 𝐜. Given 𝐜 one could solve this system of linear equations to obtain 𝐮. In fact, this system of equations is overdetermined, having n constraints but only k unknowns. Throwing away all but k of these constraints will still provide a fully determined system of equations, provided that the remaining submatrix of B is invertible. Since 
 δ^(k)=0, the system we solve in step three is eq:goodsubset, which is precisely this same system of equations, which we can solve to obtain the unique solution 𝐮̃ = 𝐮. Now 𝐮̃ B is 𝐜, which we have checked in step four has Hamming distance δ from 𝐯 and so the algorithm terminates with success. How likely is this to happen?

The string δ has n bits, of which δ have errors and each trial randomly selects k of these n bits. The probability that the subset contains no errors is exactly

    n-δk/nk.

Hence, the probability of success in one trial is

    p_ISD≃ 0.29 n-δk/nk.

If δ is small compared to n and k then one can make the approximation

    p_ISD≃ 0.29(1-k/n)^δ.
 
The expected runtime of ISD is 1/p_ISD, which can be reduced to 1/√(p_ISD) if one replaces random sampling with Grover search <cit.>.



§ DISCUSSION


In this work, we show how to solve the Quantum Distributional State Construction problem for |Ψ_b⟩, the superposition over bit strings within Hamming distance b of a linear code of length n and dimension k. We map out the region of n,k,b where this procedure efficiently yields a high fidelity approximation to |Ψ_b⟩. We find that this region extends beyond the uncomputation barrier that thwarts a direct approach to constructing |Ψ_b⟩.

We restricted our investigations in several respects. We considered only balls whose radius is less than half the distance of the code, so that the balls do not overlap. We considered only random binary linear codes. And the only application we considered is the bounded distance decoding problem where we did not beat a known classical algorithm  and our procedure can be dequantized. Given that |Ψ_b⟩ can be constructed efficiently out to surprisingly large b, it may be interesting to consider potential algorithmic, cryptographic, or complexity-theoretic applications beyond what we considered here. Or, the state preparation method used here could perhaps be applied to overcome uncomputation barriers in other Quantum Distributional State Construction problems.



Acknowledgments:  We thank Madhu Sudan for discussions a long time ago. We thank Noah Shutty and Bill Huggins for more recent discussions. And we particularly thank Oded Regev for a key observation.





§ PHASE KICKBACK



Let f:{0,1}^n →{0,1} be a function for which we know an efficient classical circuit. A unitary U_f can always be implemented efficiently such that <cit.>:

    U_f |𝐱⟩|z⟩ = |𝐱⟩|z ⊕ f(𝐱)⟩

where the second register is a single qubit.
By initializing the second register to |0⟩ one obtains f(𝐱) written into the computational basis in the second register.  One can instead compute f into the phase by initializing the second register with 1/√(2)( |0⟩ - |1⟩) so

    U_f |𝐱⟩1/√(2)( |0⟩ - |1⟩) = (-1)^f(𝐱)|𝐱⟩1/√(2)( |0⟩ - |1⟩).

The second register remains unentangled with the first and can be discarded. The net effect is |𝐱⟩→ (-1)^f(𝐱)|𝐱⟩. By linearity, an arbitrary superposition undergoes the transformation

    ∑_𝐱 a(𝐱) |𝐱⟩→∑_𝐱 (-1)^f(𝐱) a(𝐱) |𝐱⟩.

This is known as phase kickback and  it is  one of the  steps in the state preparation algorithm.



§ SYSTEMATIC FORM FOR LINEAR CODES



The k × n generator matrix B for a code C is not unique. The basis vectors 𝐛_1,…,𝐛_k that form the rows of B can be replaced by any invertible ℤ_2-linear combination of 𝐛_1,…,𝐛_k without changing the code. One consequence of this is that B can be put into systematic form:

    B = [ 𝕀_k × k | R ],

where 𝕀_k × k is an identity matrix and R is a general k × (n-k) matrix. Thinking of ℤ_2 as a finite field, systematic form is recognizable as reduced row echelon form, which can be obtained in O(k^2 n) classical operations using Gaussian elimination. Let

    B_⊥ = [ R^T | 𝕀_n-k× n-k].

By construction, B B_⊥^T = R ⊕ R = 0. Hence, B_⊥ is a generator matrix for the dual code C^⊥. From this construction, we see that given a generator matrix for C, a generator matrix for C^⊥ can be obtained in polynomial time classically.   Also if we want to generate a random code and its dual at the same time we can just generate a random k × (n-k) matrix R and then get B and B_⊥ immediately.



§ KRAWTCHOUK IDENTITY



As stated in eq:ident, Krawtchouk polynomials obey the identity

    ∑_j=0^b K_j^n(x) = K_b^n-1(x-1).


To show this, we note that binomial coefficients have the property xr = x-1r + x-1r-1. Thus

    ∑_j=0^b K_j^n(x)     =    ∑_j=0^b ∑_r=0^j xrn-xj-r (-1)^r 
        =    ∑_j=0^b ∑_r=0^j [ x-1r+x-1r-1] n-xj-r (-1)^r 
        =    ∑_j=0^b ( [ ∑_r=0^j x-1rn-xj-r (-1)^r ] + [ ∑_r=0^j-1x-1rn-xj-r-1 (-1)^r+1] ) 
        =    ∑_j=0^b ( K_j^n-1 (x-1) - K_j-1^n-1(x-1) ) 
        =     K_b^n-1(x-1),

where the last line follows from K_-1^n-1(x-1) = 0. □



